USRE36679E

USRE36679E - Method of cancelling ghosts from NMR images

Info

Publication number: USRE36679E
Application number: US08/198,866
Authority: US
Inventors: Avideh Zakhor; Richard R. Rzedzian
Original assignee: Aurora Imaging Technology Inc
Current assignee: Aurora Imaging Technology Inc
Priority date: 1989-08-07
Filing date: 1994-02-18
Publication date: 2000-05-02
Anticipated expiration: 2009-08-07
Also published as: US5089778A

Abstract

A method of cancelling ghosts from NMR images. The method involves estimating a phase difference function Δ (n₁, n₂) and using that function to solve a linear system of equations to find the magnitudes of the true object densities at the true image and ghost locations x(n₂,n₂) and x(n₁,n₂ +N/2), respectively, where the image has dimensions N×N_s. Experimental values of Δ (n₁, n₂) for a variety of objects indicate that its variation along n₁ is considerably larger than along n₂. Thus, for each column n₁, the phase difference function Δ (n₁, n₂) can be modelled as a one-dimensional function of n₂ with two parameters α (n₁) and β (n₁), which are estimated from the pixels in the 2-D FFT processed reconstructed image Y(n₁,n₂). These parameters are then used to estimate Δ (n₁, n₂), which is ultimately used to de-ghost the image.

Description

BACKGROUND OF THE INVENTION

The present invention relates to a method for improving reconstructed NMR images and, more specifically, to a method for cancelling ghosts from NMR images.

Under ideal conditions, reconstructed NMR images must be positive and real. This is because, in NMR experiments, the observed signal is Fourier transform of density distribution of the object under consideration, which by definition is a real positive quantity. In practice however, even the most straightforward ways of scanning the k_x -k_y space, (e.g., row by row or column by column scanning without date reversal), result in complex images. This is partly due to a shift in the data from the notional origin in k space. If readout gradient is constant and uniform time samples of data are inverse Fourier transformed to generate the image, delay in the time data translates into a linear phase shift in the reconstructed image. These phase shifts can be easily determined and eliminated, since they do not affect the magnitude of the reconstructed images.

When NMR data is obtained by scanning the k_x -k_y space, with data reversal on alternate y lines (where y is the horizontal axis across which readout occurs), time delays between the start of data acquisition and the start of the readout pulse are different for even and odd lines. The effect of this on the image manifests itself as a ghost separated by half the image size. Under these conditions, if readout gradient is constant and data is sampled uniformly in time, then the ghost image can be entirely removed by a first-order phase difference between odd and even lines. However, when the readout gradient is sinusoidal and the image is reconstructed by inverse Fourier transforming non-uniform samples (see, e.g., the method disclosed in U.S. Pat. No. 4,740,748) the difference between even and odd line delays degrades introduces ghosts and thus the quality of the resolution.

This, together with asymmetry of the sinusoidal readout gradient for even and odd lines can be modeled by multiplying even and odd parts of the NMR image by two separate phase functions φ (n₁,n₂) and θ (n₁,n₂). More specifically, if x(n₁,n₂) denotes the true density distribution of the object under consideration, and Y(n₁,n₂) denotes the 2-D inverse discrete Fourier transform of the time data (i.e., the reconstructed ghosted image), then the even and odd parts of the observed image, Y_even (n₁,n₂) and Y_odd (n₁,n₂)can be modeled as: ##EQU1## where the dimensions of the reconstructed image are N×N_s, and the number of echoes is N. (Note that if the sinusoidal readout y gradient was identical for even and odd lines, then the even and odd phase functions would have been identical. That is, φ (n₁,n₂).tbd.θ (n₁,n₂)).

Having modeled the ghosted image, the objective can be stated as estimation of the true object density distribution x(n₁,n₂) from the observed ghosted image Y(n₁,n₂).

From equations (1) and (2), it is clear that if the phase functions φ (n₁,n₂) and θ(n₁,n₂) are known for all values of n₁ and n₂, then x(n₁,n₂) and x(n₁,n₂ +N/2) can be determined from Y_even (n₁,n₂) Y_odd (n₁,n₂) by solving a 2×2 linear system of equations. In practice, the difference between φ (n₁,n₂) and θ (n₁,n₂) can be determined experimentally by placing a test object in the upper and lower half of the field of view (FOV) and measuring the difference between even and odd parts of the resulting images. More specifically, when the object is in the upper half of the FOV, by definition:

x(n.sub.1,n.sub.2 +N/2)=0

Substituting this into equation (1) and (2): ##EQU2##

The phase difference between Y_even (n₁,n₂) Y_odd (n₁,n₂) can be used to obtain ##EQU3##

Similarly, by placing the object in the lower half of the FOV: ##EQU4##

The phase difference between Y_even (n₁,n₂) and Y_odd (n₁,n₂) can be used to obtain ##EQU5##

Thus, experimental values of Δ (n₁,n₂) and Δ (n₁,n₂ +N/2) can be used to determine A(n₁,n₂) and B(n₁,n₂) by solving the above linear system of equations. Once A and B are determined, their magnitudes can be used to find x(n₁,n₂) and (n₁,n₂ +N/2) respectively.

Experimentally, there are two major drawbacks with the above approach. The first drawback has to do with the fact that the phase difference function is a function of the parameters for the NMR experiments. Some of these parameters are the strength of the x, y and z gradients, and the static magnetic field or the RF. Therefore, to be able to apply this method successfully, a different look up table is needed for different experimental set ups. The second drawback has to do with the fact that the phase difference function Δ (n₁,n₂) is somewhat object dependent. More specifically, although the general shape of Δ (n₁,n₂) does not vary drastically from one object to the next, the change is large enough to introduce considerable amount of ghosts. The third drawback of the above approach has to do with the fact that obtaining the phase difference function Δ (n₁,n₂) of a test object for all values of n₁ and n₂ is a non-trivial task from an experimental point of view. This has to do with factors such as susceptibility effects.

In short, it has been found that the performance of the above scheme is inadequate for most ghosted images. Accordingly, it would be highly desirable to process NMR signals using a method of ghost correction in the form of an algorithm which is automatic and does not require a look up table.

SUMMARY OF THE INVENTION

The object of the present invention is therefore to provide a method for automatically eliminating ghosts in NMR signals resulting from the difference between even and odd line delays in a traversal of k-space using a sinusoidal readout gradient without using a look-up table.

The present invention achieves the foregoing objective by a method involving the steps of:

(a) taking a two dimensional inverse Fourier transform of the time data to obtain the ghosted image Y(n₁,n₂);

(b) computing the signal energy for each column using ##EQU6##

(c) discarding the columns whose signal energy level are below a predetermined threshold;

(d) estimating α (n₁) and β (n₁) for each remaining column of data; i.e. n₁ =0, . . . N_s -1, by:

(i) finding the phase difference function Δ (n₁,n₂) for all ghosting pixels of the column; and

(ii) solving the following simultaneous equations to find linear least square estimates of α (n₁) and β (n₁): ##EQU7##

(e) using α (n₁) and β (n₁) in the above equation to find the phase difference Δ (n₁,n₂) for 0≦n₂ <N; and

(f) using Δ (n₁,n₂) in equation (11) to find A(n₁,n₂) and B(n₁,n₂) for 0≦n₂ <N.

From equations (9) and (10), the true object density distribution at (n₁,n₂) and (n₁,n₂ +N/2) are found by taking magnitude of A(n₁,n₂) and B(n₁,n₂), respectively.

BRIEF DESCRIPTION OF THE DRAWINGS

These and other features and advantages of the present invention will become apparent when the following text is read in conjunction with the accompanying drawings, in which:

FIG. 1 shows a flow diagram of the ghost correction method of the present invention;

FIG. 2 depicts a reconstructed image with ghosts.

FIG. 3 graphically depicts the phase difference function; and

FIGS. 4a-11b show before/after examples of NMR images processed with the algorithm of the present invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

Referring first to FIG. 1, a simplified block diagram of the method of the invention is shown. First, the raw NMR signal is converted to an image by performing a two dimensional inverse Fast Fourier Transform (FFT) 2. The resultant image has ghosts, which are eliminated by processing the image with the ghost elimination algorithm 4 of the present invention.

FIG. 2 depicts an image with ghosts. As shown therein, a bright line 6 represents the true image, while a less bright line 8 is a ghost of true image 6. Assuming that the traversal through k-space involves the sampling of 128 lines, the ghost will be separated by the true image by 64 lines, i.e. half the image size.

FIG. 3 graphically depicts the phase difference function Δ (n₁,n₂). As shown therein, the phase difference function is appropriately symmetrical along the center of the n₂ axis. The phase difference function is defined by the offset α and slope β of a straight line 10 obtained by taking the least square estimation of a number of points;

certain points

12,14 are eliminated (ignored) as being out-of-range since Δ (n₁,n₂) is assumed to vary smoothly.

The ghost correction algorithm of the present invention takes advantage of the approximate shape of the phase difference function, which as described above, can be obtained experimentally for test objects. Experiments indicate that variation of the function Δ (n₁,n₂) is considerably smaller along n₂ than n₁. In fact, for fixed n₁, variation along n₂ is symmetric about n₂ =N/2 and Δ (n₁,n₂) can be closely approximated by a piecewise linear function of the form: ##EQU8##

The algorithm of the present invention takes advantage of the above approximation by estimating α (n₁) and β (n₁) for each column of the data, i.e. n₁ =0, . . . ,N_s -1. It consists of the following steps:

1. If raw sampled data is input, take 2-D inverse Fourier transform of the time data to obtain the ghosted image Y(n₁,n₂).

2. Determine signal energy for each column of the data by computing ##EQU9##

3. Discard columns whose signal energy level is below a fixed threshold. This threshold is an input to the program of Appendix (A), and is denoted by the double precision variable "snr".

4. Let S denote the set of indices of the columns whose signal energy level is larger than the threshold snr. Estimate α (n₁) and β (n₁) for all n₁ εS. The estimation procedure will be discussed at length later.

5. Use α (n₁) and β (n₁) in equation (12) to find A(n₁,n₂) and B(n₁,n₂) for 0≦n₂ <N. From equations (9) and (10), the true object density distribution at (n₁,n₂) and (n₁,n₂ +N/2) are found by taking magnitude of A(n₁,n₂) and B(n₁,n₂), respectively.

The algorithm of the present invention has been implemented in "C" programming language, and its listing is included in Appendix (A). The name of the "C" program is "correct.c". The usage of the program is described by simply invoking its executable without any parameters:

usage: correct

-g (if set, will do inverse two-dimensional FFT)

-s (snr: avoid processing columns with small signal component)

-t (ratio between--Y(n1,n2)--and--Y(n1,n2+N/2)--: must be larger for central columns)

-r (even to odd ratio: must be close to 1 for ghosting pixels)

-m (mse² tolerance: discard pixels whose llsq³ mse is large)

-h (ghost larger than image: 1 left, 2 right, 3 both)

-d (parameter controlling smoothness of alpha corrections)

-f (parameter controlling smoothness of beta corrections)

-a (scale to see the ghosts clearly)

-e (expand by an integer factor in each direction)

-c (convert to hexadecimal for halftone printout)

inputfile

outputfile

As it is seen, the program takes an input file representing the sampled time data or its two-dimensional (2-D) inverse Fourier transform, and generates an output file containing either the binary or hexadecimal version of the ghost-free image. Clearly, the parameters used by the program are set by the operator in such a way that the performance of the algorithm is optimized. Some of these parameters such as "-g", "-c", and "-f" are only used as switches to set flags, while others are used to set internal variables (either integer or floating point) to specific values. An example of the usage of the program is as follows:

correct -g -t100. -r1.5 -s5. -m2. -f -d.7 -b.04 -a750 -e4 inputfile outputfile

In the above example, the internal variable associated with "-t" is 100, the one associated with "-r" is 1.5, etc.

The function of various parameters used in step 1. through step 3. of the algorithm is as follows:

The "-g" flag determines whether or not the input ASCII file is the raw time data or its 2-D inverse Fourier transform. Specifically, if the "-g" option is set, then the program assumes the input file to contain raw time data, and computes the 2-D inverse Fourer transform of the data by successive application of the "ifft05" subroutine. This subroutine takes inverse FFT of one-dimensional sequences, and its listing is included in Appendix (B).

The "-s" option sets the internal double precision variable "snr" used in steps 2. and 3. above. A large number of the columns in the ghosted input image correspond to empty space in the magnet, and therefore have no signal component. To prevent the ghost correction algorithm from processing these columns, the signal energy for each column of the ghosted image is computed and compare it to a fixed threshold. This fixed threshold is the "snr" variable which is set by the "-a" option. The signal energy for the n₁ th column is defined to be: ##EQU10## where Y(n₁,n₂) denotes the value of the ghosted image at location (n₁,n₂). Thus, the AGC algorithm only processes columns whose energies exceed the threshold set by the variable "snr". An appropriate value of "snr" is 5 for input files containing raw time data, and 2×10⁸ for the ones containing the 2-D inverse Fourier transform of the raw time data.

Step 4. of the algorithm is now described in more detail. As set forth in the Background section above, the phase difference function Δ (n₁,n₂) can be found experimentally by placing a test object in the lower or upper half of the FOV. Specifically, the ghost of a test object at location (n₁,n₂) with ##EQU11## appears at (n₁,n₂ +N/2), and vice versa. Thus, the location of the ghost of a pixel at (n₁,n₂) is (n₁,(n₂ +N/2) mod N). In general, if an object only fills out half of the FOV, its ghost does not overlap with itself in the reconstructed image. Under these conditions, the phase difference function associated with the object and the particular experimental set up can be determined empirically for regions of the reconstructed image which correspond to the object rather than its ghost. (Two special cases of this were discussed in the Background section. These cases correspond to the object being in the lower or upper half of the FOV). On the other hand, when an object fills out more than half of the FOV, its phase difference function can only be determined for pixels whose ghosts are not superimposed on pixels corresponding to other parts of the object. This information about Δ (n₁,n₂) can be exploited to find the parameters α and β as defined by equation (12). Pixels which correspond to bright (high intensity) parts whose locations correspond to either empty space in the FOV or parts of the object with very little or no energy, will be referred to as "ghosting pixels". Specifically, if the pixel at location (n₁,n₂) is a ghosting one, then by definition:

1. It corresponds to a high energy point in the object. Therefore

|x(n.sub.1,n.sub.2)|>>0

2. The pixel at location (n₁,(n₂ +N/2) mod N) corresponds to a low energy part of an object or empty space in the FOV. That is

|x(n.sub.1,(n.sub.2 +N/2) mod N)|≈0

Step 4 of the automatic ghost correction algorithm derives the parameters associated with the n₁ th column in two steps. Specifically, it first finds the value of the phase difference function for all the "ghosting pixels" of the column, and then solves an overdetermined linear system of equations to find linear least square estimates of α (n₁) and β (n₁). At this point, the key question which remains to be answered is the way ghosting pixels are detected. The present invention uses two criteria for classifying pixels as ghosting ones.

The first criterion is a direct consequence of equation (11) and the second part of the definition of ghosting pixels. It takes advantage of the fact that the magnitude of the even and parts of a ghosting pixel at location (n₁,n₂) are identical. Specifically, if Y(n₁,n₂) is a ghosting pixel, from equation (11):

|Y.sub.even (n.sub.1,n.sub.2)|=|Y.sub.odd (n.sub.1,n.sub.2)|=|A(n.sub.1,n.sub.2)|=x(n.sub.1,n.sub.2)                                                (14)

Thus, if the ratio between the magnitude of even and odd parts of the pixel (the "eoratio") at location (n₁,n₂) are equal or somewhat close to each other, then Y(n₁,n₂) can be classified as a ghosting pixel. In the program listing of the automatic ghost correction algorithm in Appendix (A), if the ratio between the magnitudes of even and odd parts of a pixel are between [1/eoratio, eoratio], then the pixel is classified as a ghosting one. As seen, "eoratio" denotes a double precision variable which is input to the program by the operator via parameter -r. It is important to note that under this criterion, the conditions for Y(n₁,n₂) and Y(n₁,n₂ +N/2) to qualify as ghosting pixels are identical Thus, if equation (14) is exactly or approximately satisfied, then there is an ambiguity as to which pixel is the ghosting one. As discussed below, the second criteria for ghosting pixel detection helps to resolve this ambiguity.

The appropriate values of "eoratio" is anywhere between 1 and 2. In most of the examples of the next section, the "-r" option (the internal variable "eoratio") is set to 1.5. If it is set to small values (i.e. too close to 1), then the number of ghosting pixels will be too small, and therefore estimation of α (n₁) and β (n₁) will not be robust. On the other hand, if it is set to larger values such as 2 or even 3, then the chosen pixels might not necessarily be the ghosting ones. This will increase the error in the observations for the linear least-squares estimation, and therefore will result in less accurate estimation of α (n₁) and β (n₁).

The second criteria takes advantage of the definition of ghosting pixels. To describe this condition, equations (1) and (2) are rewritten in the following way: ##EQU12## As expected, if θ (n₁,n₂)=φ (n₁,n₂)

or equivalently

Δ (n₁,n₂)=0

the observed image, Y(n₁,n₂), becomes ghost free and is identical to the true object density function x(n₁,n₂). Recall that if the pixel at location (n₁,n₂) is a ghosting one, then by definition, the magnitude of x(n₁,n₂) must be large (i.e. not at the noise level) and the magnitude of x(n₁,n₂ +N/2) must be very small (i.e. not at the noise level). From equations (15) and (16), if the difference between θ (n₁,n₂) and φ (n₁,n₂) is small (or equivalently Δ (n₁,n₂) is small), then a ghosting pixel at (n₁,n₂) results in a large value of the following ratio: ##EQU13##

Specifically, as Δ (n₁,n₂) changes from 0 to π, the above ratio changes from ∞ to 0. It has been found experimentally that Δ (n₁,n₂) is smaller for columns closer to the center of the magnet (i.e. around n₁ =N_s /2). This implies that the ratio of equation (17) becomes larger as the index of the column under investigation becomes closer to N_s /2. Thus, the second criterion for detecting ghosting pixels of the n₁ th column consists of

1. Computing the quantity shown in equation (17) for each pixel.

2. Comparing this ratio to a fixed threshold associated with that column.

Clearly, this threshold is column dependent and becomes larger as the indices of the column become closer to N_s /2. In the program listed in Appendix (A), the threshold for the n₁ th column is:

the internal variable "threshold" set by "-t" option, if n₁ =N_s /2.

decreases linearly with |n₁ -N_s /2|for |n₁ -N_s /2|<15.

is equal to 1 for |n₁ -N_s /2|>15.

In most of the following examples, the "-t" option (the internal double precision "threshold") is set to 100. In general, the appropriate value for the "-t" option depends on the amount of ghosting in the central part of the image, and lies between 1 and 10000. If the reconstructed image suffers from considerable amount of ghosting in the central columns, then "threshold" must be set at a small value e.g. 1. Otherwise, it should be set at a larger value, say 100, so that the center 30 columns of the image remain more or less unchanged by the algorithm. If the central columns of an image are ghost-free, setting the "-t" option at small values might result in unnecessary distortions in these columns.

To summarize, the first ghost detection criterion checks the ratio between the magnitudes of even and odd parts of the pixel at location (n₁,n₂). If this ratio is close to one, then either Y(n₁,n₂) or Y(n₁,n₂ +N/2) are classified as ghosting pixels. To resolve this ambiguity and to improve the detection procedure, a second criterion is used which computes the ratio shown in equation (17). For columns close to the center of the magnet, large values of this ratio imply a ghosting pixel at location (n₁,n₂). However, for columns further away from the center, the second criterion becomes more or less inconclusive, and other ways must be found to overcome the ambiguity problem of the first criterion. The present invention uses the apriori knowledge about the approximate shape of the phase difference function in order to resolve this ambiguity. Detailed experimental procedures for obtaining the phase difference function was described in the Background section above. Unlike that "look up" table approach, the present invention does not require detailed and exact values of the phase difference function. In fact, the algorithm of the present invention only needs to know as much about Δ (n₁,n₂) as to make binary decisions.

From the description of the automatic ghost correction algorithm, it is clear that the ghosting pixel detection part is rather hueristic. To decrease the sensitivity of the algorithm to this part, and to improve the estimation of the phase difference function, the following measures are preferably employed:

From classical results in estimation theory, it is clear that the error in estimating the parameters of Δ (n₁,n₂), i.e. α and β, is reduced as the number of observations is increased. In this case, the observations are the empirical value of the phase difference function for the ghosting pixels. Experimental results seem to indicate that for columns whose number of observations (or equivalently the number of ghosting pixels) is small, the error in estimating β in equation (12) becomes rather large. This error manifests itself as large magnitude for β, resulting in unrealistic values of the phase difference function. Since experimental results indicate that the magnitude of β is small for most columns, the present invention sets to β=0 for columns whose number of ghosting pixels is less than a fixed integer. In the program listing of Appendix (A), this integer has been chosen to be 8. To make the estimation part of the algorithm even more robust, β is set to zero when the magnitude of its estimated value exceeds a certain threshold. The optimal value for this threshold was found empirically from the approximate shape of the phase difference function for various test objects. For the program listed in Appendix (A), this threshold was set to 0.08.

A second measure taken to improve the robustness of the algorithm is to discard ghosting pixels whose least-square residue is too large. Specifically, if i₁,i₂, . . . , i_upper <N/2and j₁,j₂, . . . , j_lower ≧N/2 denote the indices of the ghosting pixels of the n₁ th column, taking into account equation (12), the linear least-squares estimation of α and β consists of solving the following overdetermined linear system of equations: ##EQU14##

If α (n₁) and β (n₁) denote the solution of the above linear least-squares problem, the residue of the ghosting pixels at (n₁,i₁) and (n₁,j₁) are defined to be: ##EQU15## and the means square error is defined to be ##EQU16##

To reduce the likelihood of false alarm (i.e. declaring non-ghosting pixels as ghosting ones), pixels for which the ratio between their residue and the mean squared error exceeds a predetermined. threshold "mse" are discarded.

The "-m" option sets the threshold "mse". Appropriate values for the "-m" option can range anywhere between 1 and 3. Clearly, if "mse" is too small, then most of the already detected ghosting pixels will be discarded for the second estimation process. On the other hand, if "mse" is too large, picels which have been misclassified as ghosting ones are not discarded and therefore result in large amounts of error in observations used for the linear least-squares estimation process. In most of the examples of the next section, the "-m" option is set to 2.

Another measure to improve the robustness of the algorithm is the -h option. Recall that the ghost detection part of the algorithm first checks the ratio between the magnitudes of even and odd parts of the pixel at location (n₁,n₂). If this ration is close to 1 (or more specifically is in the range [1/eoratio, eoratio]) then either Y(n₁,n₂) or Y(n₁,(n₂ +N/2) mod N) are classified as ghosting pixels. If the phase difference function is less than π then the magnitude of the ghost (i.e. the ghosted pixel) is smaller than that of the object (i.e. the ghosting pixel). On the other hand, if the phase difference function is larger than π, the magnitude of the ghost becomes larger than that of the object causing it. Since we expect the phase difference function to be small (less than π) in the center of the magnet, the differentiation between ghosting and ghosted pixels is straightforward for the central columns of the image, and therefore there is not any ambiguity in computing the phase difference function. In fact, as we mentioned earlier, the "-t" option takes advantage of this in order to detect/differentiate ghosting and ghosted pixels for the central 30 columns. Specifically, if the ratio shown in equation (17) is larger than the "eoratio" parameter set by the "-t" option, then Y(n₁,n₂) is the ghosting pixel and Y(n₁,n₂ +N/2) is the ghosted pixel. On the other hand if the ratio of equation (17) is smaller than 1/eoratio, then Y(n₁,n₂ +N/2) is the ghosting and Y(n₁,n₂) is the ghosted pixel. However, this clear distinction between ghosting and ghosted pixels is possible only if the phase difference function is known to be less than π, which in out situation corresponds to the central columns of the image. For columns away from the center, the phase difference function might become larger than π, thus creating an ambiguity about the ghosting and ghosted pixels. To resolve this ambiguity, the operator has to provide the algorithm with a clue as to whether or not the phase difference function is larger than π. Fortunately, this is a simple visual task since in the areas of image with large phase difference function (i.e. larger than π) the ghosts are brighter than the actual object causing the ghost. Therefore, if the ghost has larger intensity than the object in the areas to the left of central columns, "-h" option must be set to 1. Similarly, if the ghost has larger intensity than the object in the areas to the right of central columns, "-h" option must be set to 3. As shown below in most imaging situations, the ghost has lower intensity than the actual object itself, and thus there is no need to set the "-h" option to any value. However, in few of the heart images shown in FIGS. 4-11, the ghosts become larger than the objects and the "-h" option is needed to guide the program to correct them.

The "-d" option sets the internal double precision parameter "smooth" which smoothes the values of for neighboring columns. Specifically, if the value of for the n₁ th column is different from those of the past four columns by more than an amount specified by "smooth", then α (n₁) and β (n₁) are discarded and the n₁,th column appears unchanged in the output image Appropriate values for "smooth" can range between 0.3 and 2. Clearly, if "-d" option is set to a small value, e.g. 0.3, then α and β for most columns will be discarded and the final output will resemble the ghosted input image to a great extent. On the other hand, if "-d" option is set to a large value like 2, then the amount of smoothing is minimized, and the final image might have some columns which stand out amont their neighboring columns because of their drastically, different values of α and β.

The "-b" option sets the internal double precision parameter "betasmth" which smoothes the values of β for neighboring columns. The functional description of this variable is similar to that of the "-d" option. Appropriate values for "betasmth" however, can range between 0.01 and 1.

The "-f" flag determines whether or not the processed ghost-free image needs to be rearranged so that the center of the magnet coincides with the center of the image. This flag is set for all the examples shown in the next section.

The "-a" option sets the internal double variable "scale" which scales the processed image before it is written in the output file. This parameter is normally set anywhere between 500 and 1000.

The "-e" option sets the internal variable "expand" which expands the final output image by an integer in each direction. The expansion is basically done by repeating each pixel a fixed number of times in x and y directions.

The "-c" flag determines whether the processed image is written in binary or hexadecimal format. The hexadecimal format is necessary for generating halftone images with PostScript commands and the Apple LaserWriter.

"Inputfile" denotes the name of the input file which contains the ASCII characters representing the time raw data or its two-dimensional inverse Fourier transform. If such an input file does not exist, the program exists while printing "input file does not exist".

"Outputfile" denotes the name of the output file which contains numbers between 0 and 255 representing the processed image. The format of these numbers is either in binary or hexadecimal depending on whether or not the "-c" flag is set.

Examples of images processed by the ghost elimination algorithm are shown in FIGS. 4-10. For each example, the ghosted image, the processed ghost-free image, and the parameters used with the algorithm will be shown.

FIG. 4a shows a ghosted heart image which was processed by the ghost elimination algorithm with the following parameters:

correct -g -t100. -r1.5 -s5. -m3. -h1 -f -d2. -b1. -a1000 -e4 htnum htproc

The processed ghost-free version is shown in FIG. 4b. Clearly, the AGC algorithm has done an excellent job of removing the ghosts. As it is seen, the magnitude of the ghost in the left side of FIG. 4a is larger than that of the object causing it. Therefore, the parameter "-h" had to be set at 1 in order to remove the ambiguity in the phase difference function for the columns to the left of the image.

The second example of the ghost elimination algorithm is shown in FIG. 5. The original ghosted heart image is shown in FIG. 5a, and its processed ghost-free version is shown in FIG. 5b. The parameters of the algorithm were:

correct -g -t100. -r1.5 -s5. -m2. -f -d1. -b.05 -a1000 -e4 heart5num heart5proc

A third example of the algorithm is shown in FIG. 6a, and its processed ghost-free version is shown in FIG. 6b. The parameters of the algorithm were:

correct -g -t100. -r1.5 -s5. -m2. -f -d2. -b.1 -h1 -a750 -e4 heart4num heart4proc

Similar to FIG. 4a, the magnitude of the ghost in the left side of FIG. 6a is larger than that of the object causing it. Thus, "-h" option had to be set to 1, to remove the ambiguity of the phase difference function for the columns to the left of the ghosted image.

A fourth example of the algorithm is shown in FIG. 7. The original ghosted heart image is shown in FIG. 7a, and its processed ghost-free version is shown in FIG. 7b. The parameters of the algorithm were:

correct -g -t100. -r1.5 -s5. -m2. -f -d.7 -b.03 -h3 -a500 -e4 heart3num heart3proc

As it is seen in FIG. 7a, the magnitude of the ghost in the right and left part of the ghosted image is slightly larger than that of the object causing it. Therefore, the "-h" option is set at 3, to remove the ambiguity associated with the phase difference function.

A fifth example of the algorithm is shown in FIG. 8. The original ghosted heart image is shown in FIG. 8a, and its processed ghost-free version is shown in FIG. 8b. The parameters of the algorithm were:

correct -t10. -r1.3 -s200000000. -m2. -f -d.5 -b1. -a1000 -e4 newhtnum newhtproc

Note that in the above example, the input file contained the 2-D inverse Fourier transform of the raw time date. Therefore, the "-s" option had to be set at a different value from the previous examples.

A sixth example of the algorithm is shown in FIG. 9. The original ghosted heart image is shown in FIG. 9a, and its processed ghost-free version is shown in FIG. 9b. The parameters of the algorithm were:

correct -g -t1. -r1.5 -s5. -m3. -f -d.6 -b.06 -a1000 e4 lvrnum lvrproc

A seventh example of the algorithm of the algorithm is shown in FIG. 10. The original ghosted heart image is shown in FIG. 10a, and its processed ghost-free version is shown in FIG. 10b. The parameters of the algorithm were:

correct -g -t100. -r1.5 -s5. -m2. -f -d1. -b.1 -a1000 -e4 body2num body2proc

An eighth example of the algorithm is shown in FIG. 11. The original ghosted heart image is shown in FIG. 11a, and its processed ghost-free version is shown in FIG. 11b. The parameters of the algorithm were:

correct -g -t100. -r1.5 -s5. -m3. -f -d.7 -b.06 -a1000 -e4 body3num body3proc

Although the present invention has been described in connection with a preferred embodiment thereof, many other variations and modifications will now become apparent to those skilled in the art without departing form the scope of the invention. It is preferred, therefore, that the present invention be limited not by the specific disclosure herein, but only by the appended claims ##SPC1##

Claims

What is claimed is:

1. A method of cancelling .[.ghosts.]. .Iadd.a ghost .Iaddend.from .Iadd.a .Iaddend.NMR .[.images.]. .Iadd.image.Iaddend., comprising the steps of:

(a) taking a two dimensional inverse Fourier transform of a raw NMR signal to obtain a ghosted image Y(n₁,n₂);

(b) computing the signal energy for each column of said Fourier transformed signal using ##EQU17## (c) discarding the columns whose signal energy level are below a predetermined threshold;

(d) .[.estimating α (n₁) and β (n₁) for.]. .Iadd.identifying ghosting pixels in .Iaddend.each remaining column of data; i.e. n₁ =0, . . . N_s -1, .[.by:.].

.[.(i).]. .Iadd.(e) .Iaddend.finding the .Iadd.actual .Iaddend.phase difference .[ƒ]. Δ.Iadd._actual .Iaddend.(n₁,n₂) for .[.all.]. .Iadd.the .Iaddend.ghosting pixels .[.of the column; and.].

.[.(ii).]. .Iadd.(f) .Iaddend.solving the following simultaneous equations to find linear least square estimates of α (n₁) and β (n₁): ##EQU18## .[.(e).]. .Iadd.(g) .Iaddend.using .Iadd.the linear least square estimates of .Iaddend.α (n₁) and β (n₁) .[.in the above equation.]. to find .[.the.]. .Iadd.an estimated .Iaddend.phase difference .Iadd.function .Iaddend.Δ (n₁,n₂) for 0≦n₂ <N.[.; and.]. .Iadd.using .Iaddend. ##EQU19## .[.(f).]. .Iadd.(h) .Iaddend.using .Iadd.the estimated phase difference function .Iaddend.Δ (n₁,n₂) in ##EQU20## to find A(n₁,n₂) and B(n₁,n₂) for 0≦n₂ <N, where the dimensions of the .[.reconstructed.]. image .Iadd.data x(n₁,n₂) corresponding to the NMR image .Iaddend.are N×N_s .Iadd.and where .Iaddend. ##EQU21##

2. The method of claim 1, wherein the ratio of the magnitude of even and odd parts of a pixel at location (n₁,n₂) is used to determine whether a pixel is a true image or a ghost.

3. The method of claim 2, wherein the pixel is classified as a ghost if said ratio is approximately equal to one.

4. The method of claim 2, wherein the following ratio is computed, indicating that a pixel near the center of the magnetic cord is a ghost when the ratio extends to a predetermined value.

5. The method of claim 1, wherein columns with small signal components are not processed.

6. The method of claim 1, wherein pixels at location (n₁, n₂) whose linear base squares mean square error is larger than a predetermined amount are discarded when determining the phase difference function. .Iadd.7. The method of claim 1 wherein the actual phase difference Δ_actual (n₁,n₂) is calculated using:

Δ_actual (n₁,n₂)=phase (Y_even

(n₁,n₂))-phase (Y_odd (n₁,n₂))..Iaddend..Iadd.8. A method of producing a magnetic resonance image, comprising the steps of:

applying a bipolar readout gradient field;

sampling NMR data by scanning N lines of the k_x -k_y space in connection with the application of the bipolar readout gradient field, with reversal of the direction of data sampling on odd and even lines;

processing the NMR data to generate image data representative of a magnetic resonance image such that there is substantially no phase difference between even and odd parts of the image data..Iaddend..Iadd.9. The method of claim 1 wherein the phase difference between even and odd parts of the transform data is a function that varies in at least one of the x and y directions..Iaddend..Iadd.10. The method of claim 8 wherein the step of sampling NMR data comprises introducing a time delay in sampling the data, the time delay being selected to provide substantially no phase difference between even and odd parts of the image data..Iaddend..Iadd.11. The method of claim 10 wherein the step of sampling NMR data comprises introducing a time delay in sampling the data, the time delay being selected to provide substantially no first order phase difference between even and odd parts of the image data..Iaddend..Iadd.12. The method of claim 8 wherein the step of processing the NMR data comprises transforming the NMR data into transform data..Iaddend..Iadd.13. The method of claim 12 wherein the step of processing the NMR data comprises determining a function representative of the phase difference between even and odd parts of the transform data..Iaddend..Iadd.14. The method of claim 12 wherein the step of transforming the NMR data comprises taking an inverse Fourier transform of the NMR data..Iaddend..Iadd.15. The method of claim 14 wherein the step of transforming the NMR data comprises taking a two-dimensional inverse Fourier transform of the NMR data..Iaddend..Iadd.16. The method of claim 12 wherein the transform data are denoted by Y(n₁,n₂) and where: ##EQU22##.Iadd.17. The method of claim 16 further comprising the step of computing the signal energy for each column of the transform data using .Iadd.18. The method of claim 17 further comprising the step of discarding columns

with a signal energy below a predetermined threshold..Iaddend..Iadd.19. The method of claim 16 wherein there is an actual phase difference Δ_actual (n₁,n₂) between Y_even (n₁,n₂) and Y_odd (n₁,n₂)..Iaddend..Iadd.20. The method of claim 19 wherein the actual phase difference Δ_actual (n₁,n₂) is calculated using: Δ_actual (n₁,n₂)=phase (Y_even (n₁,n₂))-phase (Y_odd (n₁,n₂))..Iaddend..Iadd.21. The method of claim 19 wherein the step of processing the NMR data comprises using the transform data Y(n₁,n₂) and the actual phase difference Δ_actual (n₁,n₂) to calculate the image data, wherein the image data are denoted by x(n₁,n₂)..Iaddend..Iadd.22. The method of claim 21 wherein the step of processing the NMR data comprises estimating a phase difference function Δ(n₁,n₂) from the actual phase difference Δ_actual (n₁,n₂)..Iaddend..Iadd.23. The method of claim 22 wherein the estimated phase difference function Δ(n₁,n₂) is determined by

(i) calculating a plurality of values of the actual phase difference Δ_actual (n₁,n₂);

(ii) using the plurality of values of the actual phase difference Δ_actual (n₁,n₂) to find linear least square estimates of α(n₁) and β(n₁) using ##EQU23## (ii) using the linear least square estimates of α(n₁) and β(n₁) to find the estimated phase difference function Δ(n₁,n₂) for 0≦n₂ <N using ##EQU24##

.Iadd.24. The method of claim 22 wherein the step of processing the NMR data comprises using the transform data Y(n₁,n₂) and the estimated phase difference function Δ(n₁,n₂) to calculate the image data x(n₁,n₂)..Iaddend..Iadd.25. The method of claim 24 wherein the image data x(n₁,n₂) are calculated by: (i) using the estimated phase difference function Δ(n₁,n₂), Y_even (n₁,n₂), and Y_odd (n₁,n₂) to find A(n₁,n₂) and B(n₁,n₂) using ##EQU25## where ##EQU26## (ii) using the respective magnitude of A(n₁,n₂) and B(n₁,n₂) to find x(n₁,n₂) and x(n₁,n₂

+N/2)..Iaddend..Iadd.26. A method of reducing a ghost from a NMR image, comprising the steps of:

(b) identifying ghosting pixels in columns of Y(n₁,n₂); i.e. n₁ =0 . . . N_s -1,

(c) finding the actual phase difference Δ_actual (n₁,n₂) for the ghosting pixels

(d) solving the following simultaneous equations to find linear least square estimates of α(n₁) and β(n₁); ##EQU27## (e) using the linear least square estimates of α(n₁) and β(n₁) to find an estimated phase difference function Δ(n₁,n₂) for 0≦n₂ <N using ##EQU28## (f) using the estimated phase difference function Δ(n₁,n₂) in ##EQU29## to find A(n₁,n₂) and B(n₁,n₂) for 0≦n₂ <N, where the dimensions of the image data x(n₁,n₂) corresponding to the NMR image are N×N_s and where ##EQU30##

.Iadd.27. A method of reducing a ghost from a NMR image, comprising the steps of: (a) taking a two dimensional inverse Fourier transform of a raw NMR signal to obtain a ghosted image Y(n₁,n₂);

(b) calculating Y_even (n₁,n₂) and Y_odd (n₁,n₂) using ##EQU31## (c) estimating a phase difference function Δ(n₁,n₂) between Y_even (n₁,n₂) and Y_odd (n₁,n₂); and

(d) using the estimated phase difference function Δ(n₁,n₂) in ##EQU32## to find A(n₁,n₂) and B(n₁,n₂) for 0≦n₂ <N, where the image data corresponding to the NMR image are denoted by x(n₁,n₂), and where ##EQU33##

.Iadd.28. The method of claim 27 wherein the estimated phase difference function is predetermined..Iaddend..Iadd.29. The method of claim 28 wherein the estimated phase difference function is estimated from at least a portion of a Fourier transform of the NMR signal..Iaddend..Iadd.30. The method of claim 29 wherein the estimated phase difference function is estimated from the ghosted image Y(n₁,n₂)..Iaddend.