Image analysis using the direct exponential curve resolution algorithm

Abstract

A method of enhancing a series of images comprising the steps of: providing x digital image files, each image file including a matrix of n rows.times.m columns of pixels; producing a matrix of the x digital image files of x rows of m.times.n pixels per row, where each row is an unfolded image of an x digital image; splitting the produced matrix into two data sets A and B, each data set with an equal amount of images less than x; applying singular value decomposition (SVD) to data set A and limiting the SVD space to a number of expected components; calculating the projections A.sup.p and B.sup.p of data set A and data set B into the SVD space of A, with the eigenvalues of the SVD representing the projection of A; calculating the eigenvalues and eigenvectors of the product of A.sup.p and the inverse of B.sup.p ; and transforming the results of the latter calculation back to the original data set and calculating the resolved images and their contributions into the original images.

Description

FIELD OF THE INVENTION

This invention relates to a method designed to enhance a series of images, such as clinical Magnetic Resonance Images (MRI's) possessing an exponential dependency by extracting underlying images with pure exponential changes using a fast algorithm. This method compresses the information of the original images into fewer images between which the differences are maximized, while preserving all of the information in the original data set.

BACKGROUND OF THE INVENTION

Magnetic Resonance Imaging (MRI)

MRI is an imaging modality which is used predominantly in a clinical setting for obtaining spatial information of human anatomy. It is a powerful technique for visualizing tissue pathology. MRI is used to a lesser extent for materials imaging. The non-invasive nature of MRI makes it particularly useful to probe diffusion processes, material components, and materials properties without disturbing the sample. MRI is an extension of the spectroscopic technique generally referred to as Nuclear Magnetic Resonance (NMR) which provides structural detail of molecules in solution.

In clinical MRI, image contrast of various tissue types is accomplished due to the fact that the hydrogen nuclei (.sup.1 H) within these tissues posses intrinsic properties which result in unique signal response. These properties are described by two constants of the .sup.1 H nuclei, T.sub.1 and T.sub.2, which define the rate about which the signal decays back to an equilibrium value after excitation. The relationship between the constant and the signal is exponential. Another tissue property involves the spatial density of the .sup.1 H nuclei, .rho.. A higher value for .rho. within a given region of space will give a proportionately (linear) higher signal intensity in the image. Water self-diffusion (D.sub.w) within a tissue type is a fourth property which also gives rise to image contrast and is also related to the signal intensity with an exponential. Several parameters are varied in a typical MRI experiment to obtain an optimal contrast for a specific anatomy. For materials imaging, image contrast may occur for different material types or for different material properties, such as localized stress or component diffusion.

FIG. 1 represents a schematic of the main components of an MRI system 10, which includes magnet and detector assembly 12, system control 14, and computer 16. MRI utilizes several devices to create an image. A magnet 18 is used to create a large, static, homogeneous field which polarizes the .sup.1 H nuclei within the subject 20 to be imaged. A shim coil assembly 22 is generally used to trim the magnetic field for the optimum homogeneity. A radio frequency (RF) antenna 24 is used to irradiate the subject 20, thereby exciting the .sup.1 H nuclei, and collecting an RF signal as the nuclei relax to an equilibrium state. A series of three orthogonal magnetic field gradient coils 26 provide a means to obtain spatial information. These devices are contained within the magnet and detector assembly 12 and are powered with their respective power supplies 23 and 25. A system control section 14 contains several devices which control the electronic properties of the RF and gradient pulses. A pulse controller 28 connects to three gradient power supplies and amplifiers 30, 32, 34, which control the output of all gradient coils. The pulse controller 28 also controls the output of the RF transceiver 36 with power supply 37. The RF amplifier section 38 is controlled for both excitation power and receiver gain. A computer system 16 is used to interface with the pulse controller 28 using acquisition parameters, and reconstruct, display, and store the resultant images. System 16 includes computer/image processor 40 and main console/image display 42. Full instrument control is accessible through the main console 42.

The MRI operator selects from a variety of MRI experiments. These experiments, called pulse sequences, control the hardware to provide precisely controlled RF and gradient pulses which are used to condition and excite the .sup.1 H nuclei within the object An image is created by processing the collected, digitized RF signals.

The MRI signal acquired through use of a pulse sequence called pulsed gradient spin echo (PGSE) is related to the four constants, T.sub.1, T.sub.2, .rho., and D.sub.w, by the following equation:

S.varies..rho.e-TE/T.sub.2 (1-e.sup.-TR/T.sbsp.1)e.sup.-D.sbsp.w.sup.(g.gamma..delta.).spsp.2.sup..DELTA. ( 1) where .gamma. is a physical constant of the .sup.1 H nucleus. The parameters TE, TR, g, .delta., and .DELTA. are timing values within the pulse sequence which can be programmed to give the desired signal intensity variation and contrast within an image. These parameters may also be varied in a systematic way to probe the value of one or more of the constants, T.sub.1, T.sub.2, .rho., and D.sub.w. For example, by choosing a short TE (TE<<T.sub.2), a long TR (TR >>T.sub.1), and g=0, one can obtain an image which is weighted according to T.sub.2. If one obtains a series of images where the value of TE is varied, an image containing a true spatial distribution of T.sub.2 can be, in principle, created.

Anomalous (pathological, etc.) tissue usually gives different contrast compared to healthy tissue for a variety of reasons. A single image produced with an appropriate set of acquisition parameters can provide qualitative information of components within a region of space. A series of multivariate images, processed appropriately, can provide quantitative information of components within a region of space by utilizing differences in the constants T.sub.1, T.sub.2, .rho., or D.sub.w for the various tissue types. For example, in a clinical setting, quantifying tissue regions will lead to the classification of specific tissue. This process is useful in diagnosing the type and spatial extent of a specific pathology.

FIG. 2 shows a flow diagram describing the MRI process. As shown, the process includes the following steps:

1) place object into the MRI instrument (box 50);

2) scan the object using an appropriate MRI experiment to obtain the image plane of interest (box 52);

3) acquire a series of images using an appropriate MRI experiment which will provide image contrast based upon either T.sub.2 or T.sub.1 relaxation or that based upon diffusion properties (box 54);

4) store the series of images in a directory on a local computer (box 56); and

5) visualize and process the images (box 58).

Prior art methods have involved iterative fitting procedures or principal component analysis (PCA) to obtain the information and generally are time consuming, subjective, and approximate.

NMR spectral analysis

In order to illustrate the data processing approach, the analysis of NMR spectra is briefly described in this section and subsequently followed with a description of the MR image analysis. The exponential character of NMR signals has been utilized to visualize related variables (ppm values) in a series of spectra, e.g. the DOSY technique (K. F. Morris, C. S. Johnson Jr., J. Am. Chem. Soc. 115, 1993, 4291-4299). The next step in data analysis is to resolve the spectra into the pure components and their respective exponential contributions in the original spectra, which is done by the CORE (COmponent REsolved) technique (P. Stilbs, K. Paulsen, P. C. Griffiths, J. Phys., Chem., 100, 1996, 8180-8189.). This technique fits (linear combinations of) exponential functions. Once the proper exponential contributions in each of the mixture spectra are known, the spectra can be resolved.

A recent technique developed by Antalek and Windig (B. Antalek, W. Windig, J.Am.Chem, Soc,, 118, 1996, 10331-10332) enables resolution of MR spectral data in seconds. Furthermore, the method obtains a direct solution for the exponentials; no iterative methods are used to fit the data. The principles are outlined below.

Kubista showed that mixture data can be resolved when two data sets are available with the following relation (K. S. Booksh and B. R. Kowalski, J. Chemom. 8 (1994) 287-292; I. Scarminio and M. Kubista, Anal. Chem. 65 (1993) 409-416):

A=CP.sup.T ( 2) B=C.alpha.P.sup.T ( 3)

A and B are data matrices of size s*v; s is the number of spectra and v is the number of variables (e.g., wave numbers). C (size s*n) and P (size n*v) are the concentration and pure spectra matrices of n components. The matrix a is a diagonal matrix (size n*n). In these two data sets, the pure spectra and the concentrations have a correlation coefficient of one. They only differ by a certain scaling factor as defined by the diagonal matrix .alpha..

Booksh and Kowalski showed that Kubista's method can be expressed in terms of the Generalized Rank Annihilation Method (GRAM), which eliminates some problems and restrictions with Kubista's method (K. S. Booksh and B. R. Kowalski, J. Chemom. 8 (1994) 287-292). For details about the GRAM algorithm, see (E. Sanchez and B. R. Kowalski, Anal. Chem. 58 (1988) 496-499; B. Wilson, E. Sanchez and B. R. Kowalski, J. Chemom. 3 (1989) 493-398).

GRAM is mainly used as a calibration technique for techniques where a single sample results in a data matrix. The analysis of two mixture samples then results into the matrixes A and B with the relation described in Eqs. (2) and (3) for certain analytical techniques. Applications of GRAM and related techniques have been described for the excitation-emission matrices resulting from fluorescence spectroscopy (C. -N. Ho, G. D. Christian and E. R. Davidson, Anal. Chem. 50 (1978) 1108-1113; C. -N. Ho, G. D. Christian and E. R. Davidson, Anal. Chem. 52 (1980) 1071-1079; C. -N. Ho, Q. D. Christian and E. R. Davidson, Anal. Chem. 53 (1981) 92-96; M. M. C. Ferreira, M. L. Brandes, I. M. C. Ferreira, K. S. Booksh, W. C. Dolowy, M. Gouterman and B. R. Kowalski, Appl. Spectrosc. 49 (1995) 1317-1325); hyphenated techniques (L. S. Ramos, E. Sanchez and B. R. Kowalski, J. Chromatogr. 385 (1987) 165-180; J. Ohman, P. Geladi and S. Wold, Residual Bilinearization, J. Chemom. 4 (1990) 135-146; M. Josefson, L. Tekensbergs-Hjelte, Multicomponent analysis: a case report, J. Pharm. Biomed. Anal. 9 (1991) 643-654; M. J. P. Gerritsen, H. Tanis, B. G. M. Vandeginste, G. Kateman, Anal. Chem. 64 (1992) 2042-2056; M. J. P. Gerritsen, H. Peters, B. G. M. Vandeginste, G. Kateman, Chemometrics Intell. Lab. Syst. 18 (1993) 205-119) and time resolved spectral information of reactions (J. M. Henshaw, L. W. Burgess, K. S. Booksh and B. R. Kowalski, Anal. Chem. 66 (1994) 3328-3336). A major problem for practical applications is that the latter experiments need to be synchronized almost exactly and this type of calibration cannot be used as a routine method. For the application of GRAM to resolve mixture data, the same strict requirements need to be fulfilled. Furthermore, in an industrial environment only one sample of a certain class is often available, and to request another sample to fulfill the requirements of Eqs. (2) and (3) is not realistic.

As Antalek and Windig recently showed, problems such as synchronization of analyses and the availability of two samples can be avoided totally when data is available, where the signal of every component decays exponentially (B. Antalek, W. Windig, J. Am. Chem, Soc., 118, 1996, 10331-10332).

A simulated example of the exponential decays of two components with a different self-diffusion coefficient is given in following Table 1.

TABLE 1______________________________________D A B______________________________________27 8 27 8 9 49 4 9 4 3 23 2 3 2 1 11 1______________________________________ D represents a data set of which the exponentially decaying concentrations profiles of the two components are given in the first and second column under D. The first column has a faster decay than the second column. The data set A is created using the first three `spectra` of data set D, and data set B is created by using the last three spectra of data set D. Because of the exponential character of the profiles, the first column of A is a constant (three) times the first column of column B, and the second column of A is a constant (two) times the second column of B. As a consequence, the data sets A and B fulfill the requirements as expressed in Eqs. (2) and (3). In this case the matrix .alpha. in Eq. (3) contains 0.33 and 0.50 as diagonal elements. This means that in the case of exponentially decaying functions, one can use two different parts of the data set to create the two data sets necessary to apply Kubista's method or GRAM. MR image analysis

The DOSY technique mentioned above is not suited for image analysis, since the visual interpretation relies heavily on the spectral character of the data. The CORE technique is in principle suited for image analysis, but the computation time is prohibitive. It is in the range of hours for full spectra, for images the processing time would be at least an order of magnitude higher.

Sandor et. al. (T. Sandor, A. R. Bleier, P. W. Ruenzel, D. F. Adamns, F. A. Jolesz, Magnetic Resonance Imaging, 6, 1988, 27-40) described a maximum likelihood method to separate linear combinations of T.sub.2 signals. Since resolving single exponential curves takes 40 to 60 seconds this method is too slow for the analysis of complete images. This method requires operator-interaction to identify linear segments in graphic displays of the measured data.

For image analysis several approaches have been described. A technique developed by Geladi utilizes PCA (H. Grahn, N. M. Szeverenyi, M. W. Roggenbuck, F. Delaglio, P. Geladi, Chemometrics and intelligent Laboratory Systems, 5, 1989, 311-322; H. Grahn, N. M. Szeverenyi, M. W. Roggenbuck, P. Geladi, Chemometrics and Intelligent Laboratory Systems, 7, 1989, 87-93; H. F. Grahn J. Saaf, Chemometrics and Intelligent Laboratory Systems, 14, 1992, 391-396). PCA is a technique that enables one to determine an efficient base for data. For example, when we have a data set of 15 mixture spectra of 3 components, where each spectrum has a 1000 wavenumbers, principal component analysis creates a base of only three vectors which described the whole data set. This is possible since the basic dimensionality of the data set is determined by the number of components, which is 3. Similarly, in a series of, for example, 15 images, the basic dimensionality may only be 2. In such case, the 15 images can be completely described by only two images. In this lower dimensional space, it is possible to determine clusters of pixels, i.e., pixels with are similar. These clusters of pixels can be defined interactively and these pixels can than be enhanced in images.

Hornak et al described a technique related to Geladi's. However, instead of using principal components as a base for the data set, Hornak uses the spin-lattice relaxation time T.sub.1, the spin-spin relaxation time T.sub.2 and the spin density p as a basis for the data set to determine clusters of pixels (L. M. Fletcher, J. B. Barsotti, J. P. Hornak, Mag. Res. Med., 29, 1993, 623-630). The techniques of Geladi and Hornak require user interaction to determine the clusters. Since the cluster are generally not well defined, these methods are subjective and time consuming.

There is thus a need for an image analysis method which obviates the disadvantages of the above-described techniques.

SUMMARY OF THE INVENTION

According to the present invention, there is provided a solution to the needs described above. Objects of the present invention are to provide an improved method for processing Magnetic Resonance Images (MRI's) to obtain separate images of different tissue/material type; to provide such a method based on exponential profiles; to provide quantitative information about the derived images; to compress the information of the original MRI's into a minimal number of images between which the differences are maximized, preserving all the information in the original data set; to provide a fast image processing method where the analysis of a series of 15 256.times.256 images should take less than one minute; to provide an image processing method that resolves the exponentials directly, rather than through iterative linear or non linear methods; and to provide an image processing method that does not require user interaction.

According to a feature of the present invention, there is provided a method of enhancing a series of images comprising the steps of: providing x digital image files, each image file including a matrix of n rows x m columns of pixels; producing a matrix of said x digital image files of x rows of m x n pixels per row, where each row is an unfolded image of an x digital image; splitting the produced matrix into two data sets A and B, each data set with an equal amount of images less than x; applying singular value decomposition (SVD) to data set A and limiting the SVD space to a number of expected components; calculating the projections A.sup.p and B.sup.p of data set A and data set B into the SVD space of A, with the eigenvalues of the SVD representing the projection of A; calculating the eigenvalues and eigenvectors of the product of A.sup.p and the inverse of B.sup.p ; and transforming the results of the latter calculation back to the original data set and calculating the resolved images and their contributions into the original images.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a block diagram of an MRI system which can incorporate the present invention.

FIG. 2 is a flow diagram of the MRI experiment.

FIG. 3 is a flow diagram of the method of the present invention.

FIG. 4 is a diagrammatic view of the test phantom described in the illustrative example.

FIG. 5 are graphical views of signal profiles derived from the original images.

FIG. 6 are graphical views of resultant signal profiles showing extracted component A.sup.----- ; Extracted component B.sup.- -.sup.- -; Reconstructed signal profile for image #1--; Signal profile from image #1----.

DETAILED DESCRIPTION OF THE INVENTION

Referring to FIG. 2, the entire process of performing the image acquisition and obtaining results is illustrated. The following steps are taken:

1. place object into the MRI instrument (box 50); 2. scan the object using a routine localizing pulse sequence to locate the image plane of interest (box 52); 3. obtain the image series a. choose an appropriate pulse sequence to perform the imaging experiment b. select an appropriate range of sequence parameters which will cause the signal variance within the set of images to depend upon either of the properties, T.sub.1, T.sub.2, or D.sub.w. c. collect a series of images based upon the sequence parameters (box 54); 4. store the images in a directory on a local computer (box 56); and 5. use the method of the invention illustrated in FIG. 3 and described to process the data (box 58).

In order to apply the method of the present invention for resolving MR spectra, it is necessary first to represent each image by an array rather than by a matrix. This can be done by so called unfolding. For the commonly used 256.times.256 matrix this is done by taking the first row of 256 elements, followed by appending the second, third etc. column, so that an array of 256.times.256 elements is obtained. This array is then stored as a row in the data matrix. After the unfolding procedure the data analysis developed for spectral data can be applied directly to the images.

Referring now to FIG. 3, there is shown a flow diagram of the method of the present invention.

1. Read the file with unfolded images; each unfolded image is stored as a separate row in the matrix. As a consequence, the size of the matrix for 15 images of each 256.times.256 pixels is 15.times.65536 (rows.times.columns) (box 60).

2. Split the file into two parts with an equal amount of images. For a data set as described above a good way to split is in a matrix A with images 1 to 14 and a matrix B with the images 2:15 (box 62).

3. Apply Singular Value Decomposition (SVD) to the first data set A in order to create a space into which both data files can be projected. Limit the dimensionality of the SVD space to the number of components expected (box 64).

4. Project both data sets into the common space. The eigenvalues of the SVD represent the projection of A. The notation of the projections is A.sup.p and B.sup.p (box 66).

5. Calculate the eigenvalues and eigenvectors of the product of A.sup.p and the inverse of B.sup.p (box 68).

6. Transform the results from the SVD based results back to the original space into the resolved images and their contributions into the original images (box 70).

The latter method in terms of equations can be described as follows.

Eqs. (2) and (3) can be rewritten as:

C=A(P.sup.T).sup.+ (4) C.alpha.=B(P.sup.T).sup.+ (5) where (P.sup.T).sup.+ represents the pseudoinverse of the matrix P.sup.T.

Post-multiplying the left and right side of Eq. (4) by .alpha. results in:

C.alpha.=A(P.sup.T).sup.+ .alpha. (6)

Combining Eqs. (5) and (6) results in:

A(P.sup.T).sup.+ .alpha.B(P.sup.T).sup.+ (7)

This can be rewritten as:

AZ.alpha.=BZ (8) where Z=(P.sup.T).sup.+ (9)

The expression in Eq. (8) is known as the generalized eigenvector problem, where Z contains the eigenvectors and .alpha. contains the eigenvalues

In order to solve Eq. (8), the matrices A and B need to be square. This can be achieved by using a common space in which both data files are projected. Several common spaces have been proposed in literature, depending on the properties of the data sets of A and B:

(2) The latent variables of A or B.

(3) The latent variables of A+B.

(4) The latent variables of the joint spaces of (A.vertline.B) and ##EQU1## Approach (2) is the appropriate one for this study. A singular value decomposition was applied to A, after which only the significant results were kept: U, S and V so that USV.sup.T reproduces the original data set within the experimental noise range.

Substituting USV.sup.T for A in Eq. (8) results in:

USV.sup.T Z.alpha.=BZ (10)

Substitute Z=VS.sup.-1 Z.sup.*, where Z.sup.* =SV.sup.T Z. (This substitution uses the relation Z=VV.sup.T Z, which can be proven by using the fact that P is a linear combination of V and using properties of the generalized inverse.)

USV.sup.T VS.sup.-1 Z.sup.* .alpha.=BVS.sup.-1 Z.sup.* (11) UZ.sup.* .alpha.=BVS.sup.-1 Z.sup.* (12)

Left multiplying by U.sup.T results in

Z.sup.* .alpha.=(U.sup.T BVS.sup.-1)Z.sup.* (13)

Since U.sup.T BVS.sup.-1 is a square matrix, this is the standard eigenvalue problem, where Z.sup.* are the eigenvectors and .alpha. are the eigenvalues.

From the eigenvectors Z.sup.* the spectra can be resolved as follows:

Z.sup.* =SV.sup.T Z=SV.sup.T (P.sup.T).sup.+ (14) V.sup.T (P.sup.T).sup.+ =S.sup.-1 Z.sup.* (15) (P.sup.T).sup.+ =VS.sup.-1 Z.sup.* (16) P.sup.T =(VS.sup.-1 Z.sup.*).sup.+ (17)

By use of the following definition:

A=CP.sup.T =USV.sup.T (18)

The concentrations can be derived as follows:

C=USV.sup.T (P.sup.T).sup.+ (19) C=USV.sup.T Z (20) C=USV.sup.T VS.sup.- Z.sup.* (21) C=UZ.sup.* (22)

Another representation of the method of the invention is given in the Appendix in the form of the MATLAB code. Some of the figures have been generated using this code.

It can be shown that the decay rate of the exponentials can be derived directly from the eigenvalue a. This is often done in order to determine parameters in eq. 1.

Data Transformations

The method of the invention requires profiles with the behavior expressed in eqs. (2) and (3). As shown in Table 1 above, this is the case for decaying exponentials. In MRI, other exponentials are used which do not behave according to eqs. (2) and (3). In these cases transformations can be applied to transform the data into a decaying exponential, so that the method of the invention can be applied. An example is the results of T1 experiments, which results in data with the following general equation:

.alpha..sub.1 (1-e.sup.-b.sbsp.1.sup.x) (23)

This signal does not have the properties described by eqs. (2) and (3). It is possible to transform this function into a sum of two exponentials by multiplying it with

.alpha..sub.2 e.sup.-b.sbsp.2.sup.x (24)

This results in:

.alpha..sub.1 .alpha..sub.2 e.sup.-b.sbsp.2.sup.x -.alpha..sub.1 .alpha..sub.2 e.sup.(-b.sbsp.1.sup.-b.sbsp.2.sup.)x (25)

A data set transformed this way will be rank deficient since we have basically only one equation with two unknown exponentials. In order to make it possible to resolve the data set an extra variable with the exponential character described by eq. (23) needs to be added to the data set. This transformation is also valid for data sets with more than one component The mathematical technique described above will resolve two exponentials; one with a decay rate of b.sub.1, and one with a decay rate of b.sub.1 +b.sub.2. The resolved images obtained with the transformed data set can be used to calculate the T1 profiles of the original (untransformed) data set. In practice this is done by transforming the data set by a decaying exponential, and then add one column to the data set with this exponential. This data set can now be solved. For the final results the introduced resolved image needs to be ignored, and the remaining images can be used to calculate the contributions in the original data set.

EXAMPLE

Description MRI

A test phantom was constructed to test the method of the invention and is shown in FIG. 4. This consisted of a PVC plastic box 80 with the dimensions, 18 cm.times.10 cm.times.1 cm, which contained two compartments 82 and 84 separated by a thin plastic sheet 86 of thickness, 0.18 cm. The compartments 82,84, therefore, formed two wedge-shaped spaces. The angle of the wedge is 6.2.degree.. These spaces were filled with water which contained a measured amount of paramagnetic salt. One space contained 8 millimolar NiCl.sub.2 and the other 0.5 millimolar MnCl.sub.2. The separation sheet 86 was glued so as to completely seal one compartment 82 from the other 84. Large bubbles were removed and a top sealed with glue was affixed. The presence of the paramagnetic salt influences the T.sub.1 and T.sub.2 relaxation behavior exhibited by the water .sup.1 H nuclei. The extent of this influence depends upon the amount of salt present.

Using a standard spin echo pulse sequence (i.e., g=0) 15, 256.times.256 pixel, images were acquired of a plane passing through long axes of the object. The parameter, TE, in the spin echo pulse sequence was varied starting at 15 milliseconds and incremented by 15 milliseconds for each image. TR was a constant of 2000 milliseconds and the field of view was 20 centimeters. Varying TE resulted in a dependence in signal intensity for each image based upon T.sub.2 relaxation. The phantom was constructed so that one may obtain an image with three general regions, two regions where the signal represents the pure decay behavior of each component (close to the ends) and a middle region which represents a weighted sum of the two.

Data analysis

The images are unfolded and each image is stored in the matrix data as a row. Since the size of the original 15 images was 256.times.256, the matrix data has a size of 15 rows of each 65536. The pseudocode to analyze this data file is given in the appendix.

Signal profiles are obtained from the original images by averaging 30 columns within each image comprising signal from within the phantom. These profiles are displayed in FIG. 5. Signal decays faster for the compartment containing MnCl.sub.2 (i.e. the left hand portion of the curves in FIG. 5). The resolved images are obtained and their respective signal profiles plotted in FIG. 6 along with the signal profile reconstructed from the extracted components and that from image #1.

The analysis of this data set with a Matlab function on a 90 MHz Pentium computer with 64 Mb of RAM memory takes about 40 seconds.

The method of the present invention can be stored in a computer program product containing a computer readable storage medium having a computer program stored thereon of the method. The computer readable storage medium can include, for example: magnetic storage media (magnetic floppy disc, magnetic hard drive, magnetic tape); optical storage media (optical disc, optical tape, machine readable visual pattern, i.e., bar code); solid state electronic storage device (RAM, ROM); or any other physical device or medium which might be employed to store a computer program.

Although the image processing method of the invention has been described with respect to MRI images, it will be understood by those skilled in the art that the method is also applicable to other applications where a series of digital images in general, or medical images in particular, are to be processed.

The invention has been described in detail with particular reference to certain preferred embodiments thereof, but it will be understood that variations and modifications can be effected within the spirit and scope of the invention.

______________________________________PARTS LIST______________________________________10 MRI system12 magnet and detector assembly14 system control16 computer18 magnet20 subject22 shim coil assembly23 power supply24 radio frequency antenna25 power supply26 gradient coils28 pulse controller30,32,34 power supplies and amplifiers36 transceiver37 power supply38 amplifier section40 computer/image processor42 main console/image display80 PVC plastic box82,84 compartments86 separation sheet______________________________________ __________________________________________________________________________APPENDIX__________________________________________________________________________%Pseudo code to calculate images.%The information that needs to be available is:%data: a matrix in which each row contains an unfolded image% when the original image of 256*256 pixels is in a matrix M, the unfolded version of M% in the first row of data is: M(1,1) M(2,1) M(3,1). . . M(256,1) M(1,2) M(2,2). . .% M(256,256)%ncom: the number of components (pure images) to be resolved%npixx: the number of pixels in the x direction (needed to fold resultsinto images again to plot)%In the notations used, A(1,:) means the first row of the matrix A.%If the matrix range contains the values 1:4 (1,2,3,4), the matrixA(:,range) represents the first%four columns of AINITIALIZE�nrows,ncols!=size(data);range1=1:nrows-1;range2=2:nrows;.%SINGULAR VALUE DECOMPOSITION OF FIRST DATA SET%LIMIT THE OUTPUT TO NCOM%u CONTAINS SCORES, v CONTAINS EIGENVALUES, v CONTAINS EIGENVECTORS�u,s,v!=svd(data(range1,:));s=s(1:ncom,1:ncom);u=u(:,1:ncom);v=v(:,1:ncom);%CALCULATE EIGENVECTORS ev AND EIGENVALUES aa�ev,aa!=eig(u'*data(range2,:)*v*inv(s);%CALCULATE PURE IMAGESpurimage=pseudoinverse(v*inv(s)*ev);%PLOT RESULTS IN FORM OF IMAGEfor i=1:ncom; plotimage(purimage(i,:));end;__________________________________________________________________________

Claims

1. A method of enhancing a series of images comprising the steps of:

providing x digital image files, each image file including a matrix of n rows.times.m columns of pixels;

producing a matrix of said x digital image files of x rows of m.times.n pixels per row, where each row is an unfolded image of an x digital image;

splitting the produced matrix into two data sets A and B, each data set with an equal amount of images less than x;

applying singular value decomposition (SVD) to data set A and limiting the SVD space to a number of expected components;

calculating the projections A.sup.p and B.sup.p of data set A and data set B into the SVD space of A, with the eigenvalues of the SVD representing the projection of A;

calculating the eigenvalues and eigenvectors of the product of A.sup.p and the inverse of B.sup.p ; and

transforming the results of the latter calculation back to the original data set and calculating the resolved images and their contributions into the original images.

2. The method of claim 1 wherein said providing step provides x medical digital images.

3. The method of claim 1 wherein said providing step provides x MRI digital images.

4. A computer program product, comprising:

a computer readable storage medium having a computer program stored thereon for performing the steps of:

providing x digital image files, each image file including a matrix of n rows.times.m columns of pixels;

producing a matrix of said x digital image files of x rows of m.times.n pixels per row, where each row is an unfolded image of an x digital image;

splitting the produced matrix into two data sets A and B, each data set with an equal amount of images less than x;

applying singular value decomposition (SVD) to data set A and limiting the SVD space to a number of expected components;

calculating the projections A.sup.p and B.sup.p of data set A and data set B into the SVD space of A, with the eigenvalues of the SVD representing the projection of A;

calculating the eigenvalues and eigenvectors of A.sup.p and the inverse of B.sup.p ; and

transforming the results of the latter calculation back to the original data set and calculating the resolved images and their contributions into the original images.