Method of cancelling ghosts from NMR images

A method of cancelling ghosts from NMR images. The method involves estimating a phase difference function .DELTA. (n.sub.1, n.sub.2) 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.sub.1,n.sub.2) and x(n.sub.1,n.sub.2 +N/2), respectively, where the image has dimensions N.times.N.sub.s. Experimental values of .DELTA. (n.sub.1, n.sub.2) for a variety of objects indicate that its variation along n.sub.1 is considerably larger than along n.sub.2. Thus, for each column n.sub.1, the phase difference function .DELTA. (n.sub.1, n.sub.2) can be modelled as a one-dimensional function of n.sub.2 with two parameters .alpha. (n.sub.1) and .beta. (n.sub.1), which are estimated from the pixels in the 2-D FFT processed reconstructed image Y(n.sub.1,n.sub.2). These parameters are then used to estimate .DELTA. (n.sub.1, n.sub.2), which is ultimately used to de-ghost the image.

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.sub.x -k.sub.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.sub.x -k.sub.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 .phi. (n.sub.1,n.sub.2) and .theta. 
(n.sub.1,n.sub.2). More specifically, if x(n.sub.1,n.sub.2) denotes the 
true density distribution of the object under consideration, and 
Y(n.sub.1,n.sub.2) 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.sub.even (n.sub.1,n.sub.2) and 
Y.sub.odd (n.sub.1,n.sub.2) can be modeled as: 
##EQU1## 
where the dimensions of the reconstructed image are N.times.N.sub.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, .phi. 
(n.sub.1,n.sub.2).ident..theta. (n.sub.1,n.sub.2)). 
Having modeled the ghosted image, the objective can be stated as estimation 
of the true object density distribution x(n.sub.1,n.sub.2) from the 
observed ghosted image Y(n.sub.1,n.sub.2). 
From equations (1) and (2), it is clear that if the phase functions .phi. 
(n.sub.1,n.sub.2) and .theta.(n.sub.1,n.sub.2) are known for all values of 
n.sub.1 and n.sub.2, then x(n.sub.1,n.sub.2) and x(n.sub.1,n.sub.2 +N/2) 
can be determined from Y.sub.even (n.sub.1,n.sub.2) Y.sub.odd 
(n.sub.1,n.sub.2) by solving a 2.times.2 linear system of equations. In 
practice, the difference between .phi. (n.sub.1,n.sub.2) and .theta. 
(n.sub.1,n.sub.2) 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: 
EQU x(n.sub.1,n.sub.2 +N/2)=0 
Substituting this into equations (1) and (2): 
##EQU2## 
The phase difference between Y.sub.even (n.sub.1,n.sub.2) Y.sub.odd 
(n.sub.1,n.sub.2) can be used to obtain 
##EQU3## 
Similarly, by placing the object in the lower half of the FOV: 
EQU x(n.sub.1,n.sub.2)=0 
##EQU4## 
The phase difference between Y.sub.even (n.sub.1,n.sub.2) and Y.sub.odd 
(n.sub.1,n.sub.2) can be used to obtain 
##EQU5## 
Thus, experimental values of .DELTA. (n.sub.1,n.sub.2) and .DELTA. 
(n.sub.1,n.sub.2 +N/2) can be used to determine A(n.sub.1,n.sub.2) and 
B(n.sub.1,n.sub.2) by solving the above linear system of equations. Once A 
and B are determined, their magnitudes can be used to find 
x(n.sub.1,n.sub.2) and (n.sub.1,n.sub.2 +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 .DELTA. (n.sub.1,n.sub.2) is somewhat object 
dependent. More specifically, although the general shape of .DELTA. 
(n.sub.1,n.sub.2) 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 .DELTA. (n.sub.1,n.sub.2) of a 
test object for all values of n.sub.1 and n.sub.2 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.sub.1,n.sub.2); 
(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 .alpha. (n.sub.1) and .beta. (n.sub.1) for each remaining 
column of data; i.e. n.sub.1 =0, . . . N.sub.s -1, by: 
(i) finding the phase difference function .DELTA. (n.sub.1,n.sub.2) for all 
ghosting pixels of the column; and 
(ii) solving the following simultaneous equations to find linear least 
square estimates of .alpha. (n.sub.1) and .beta. (n.sub.1): 
##EQU7## 
(e) using .alpha. (n.sub.1) and .beta. (n.sub.1) in the above equation to 
find the phase difference .DELTA. (n.sub.1,n.sub.2) for 0.ltoreq.n.sub.2 
&lt;N; and 
(f) using .DELTA. (n.sub.1,n.sub.2) in equation (11) to find 
A(n.sub.1,n.sub.2) and B(n.sub.1,n.sub.2) for 0.ltoreq.n.sub.2 &lt;N. 
From equations (9) and (10), the true object density distribution at 
(n.sub.1,n.sub.2) and (n.sub.1,n.sub.2 +N/2) are found by taking magnitude 
of A(n.sub.1,n.sub.2) and B(n.sub.1,n.sub.2), respectively.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT 
Referring first to FIG. 1, a simplified block diagram of the method of 
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 .DELTA. 
(n.sub.1,n.sub.2). As shown therein, the phase difference function is 
appropriately symmetrical along the center of the n.sub.2 axis. The phase 
difference function is defined by the offset .alpha. and slope .beta. 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 .DELTA. (n.sub.1,n.sub.2) 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 .DELTA. (n.sub.1,n.sub.2) is 
considerably smaller along n.sub.2 than n.sub.1. In fact, for fixed 
n.sub.1, variation along n.sub.2 is symmetric about n.sub.2 =N/2 and 
.DELTA. (n.sub.1,n.sub.2) 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 .alpha. (n.sub.1) and .beta. (n.sub.1) for 
each column of the data, i.e. n.sub.1 =0, . . . ,N.sub.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.sub.1,n.sub.2). 
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 .alpha. (n.sub.1) and .beta. 
(n.sub.1) for all n.sub.1 .epsilon. S. The estimation procedure will be 
discussed at length later. 
5. Use .alpha. (n.sub.1) and .beta. (n.sub.1) in equation (12) to find 
A(n.sub.1,n.sub.2) and B(n.sub.1,n.sub.2) for 0.ltoreq.n.sub.2 &lt;N. From 
equations (9) and (10), the true object density distribution at 
(n.sub.1,n.sub.2) and (n.sub.1,n.sub.2 +N/2) are found by taking magnitude 
of A(n.sub.1,n.sub.2) and B(n.sub.1,n.sub.2), 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.sup.2 tolerance: discard pixels whose llsq.sup.3 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: 
EQU 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.sub.1 th column is defined to be: 
##EQU10## 
where Y(n.sub.1,n.sub.2) denotes the value of the ghosted image at 
location (n.sub.1,n.sub.2). 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.times.10.sup.8 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 .DELTA. 
(n.sub.1,n.sub.2) 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.sub.1,n.sub.2) with 
##EQU11## 
appears at (n.sub.1,n.sub.2 +N/2), and vice versa. Thus, the location of 
the ghost of a pixel at (n.sub.1,n.sub.2) is (n.sub.1,(n.sub.2 +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 .DELTA. 
(n.sub.1,n.sub.2) can be exploited to find the parameters .alpha. and 
.beta. 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.sub.1,n.sub.2) is a ghosting one, then by definition: 
1. It corresponds to a high energy point in the object. Therefore 
EQU .vertline.x(n.sub.1,n.sub.2).vertline.&gt;&gt;0 
2. The pixel at location (n.sub.1,(n.sub.2 +N/2) mod N) corresponds to a 
low energy part of an object or empty space in the FOV. That is 
EQU .vertline.x(n.sub.1,(n.sub.2 +N/2) mod N).vertline..apprxeq.0 
Step 4 of the automatic ghost correction algorithm derives the parameters 
associated with the n.sub.1 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 .alpha. (n.sub.1) and 
.beta. (n.sub.1). 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.sub.1,n.sub.2) are identical. Specifically, if Y(n.sub.1,n.sub.2) is a 
ghosting pixel, from equation (11): 
EQU .vertline.Y.sub.even (n.sub.1,n.sub.2).vertline.=.vertline.Y.sub.odd 
(n.sub.1,n.sub.2).vertline.=.vertline.A(n.sub.1,n.sub.2).vertline.=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.sub.1,n.sub.2) are equal or somewhat close 
to each other, then Y(n.sub.1,n.sub.2) 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.sub.1,n.sub.2) and Y(n.sub.1,n.sub.2 +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 .alpha. (n.sub.1) and .beta. (n.sub.1) 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 
.alpha. (n.sub.1) and .beta. (n.sub.1). 
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 .theta. (n.sub.1,n.sub.2)=.phi. (n.sub.1,n.sub.2) 
or equivalently 
.DELTA. (n.sub.1,n.sub.2)=0 
the observed image, Y(n.sub.1,n.sub.2), becomes ghost free and is identical 
to the true object density function x(n.sub.1,n.sub.2). Recall that if the 
pixel at location (n.sub.1,n.sub.2) is a ghosting one, then by definition, 
the magnitude of x(n.sub.1,n.sub.2) must be large (i.e. not at the noise 
level) and the magnitude of x(n.sub.1,n.sub.2 +N/2) must be very small 
(i.e. not at the noise level). From equations (15) and (16), if the 
difference between .theta. (n.sub.1,n.sub.2) and .phi. (n.sub.1,n.sub.2) 
is small (or equivalently .DELTA. (n.sub.1,n.sub.2) is small), then a 
ghosting pixel at (n.sub.1,n.sub.2) results in a large value of the 
following ratio: 
##EQU13## 
Specifically, as .DELTA. (n.sub.1,n.sub.2) changes from 0 to .pi., the 
above ratio changes from .infin. to 0. It has been found experimentally 
that .DELTA. (n.sub.1,n.sub.2) is smaller for columns closer to the center 
of the magnet (i.e. around n.sub.1 =N.sub.s /2). This implies that the 
ratio of equation (17) becomes larger as the index of the column under 
investigation becomes closer to N.sub.s /2. Thus, the second criterion for 
detecting ghosting pixels of the n.sub.1 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 columns become closer to N.sub.s /2. In the program listed 
in Appendix (A), the threshold for the n.sub.1 th column is: 
the internal variable "threshold" set by "-t" option, if n.sub.1 =N.sub.s 
/2. 
decreases linearly with .vertline.n.sub.1 -N.sub.s /2.vertline.for 
.vertline.n.sub.1 -N.sub.s /2 .vertline.&lt;15. 
is equal to 1 for .vertline.n.sub.1 -N.sub.s /2.vertline.&gt;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.sub.1,n.sub.2). If this ratio is close to one, then either 
Y(n.sub.1,n.sub.2) or Y(n.sub.1,n.sub.2 +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.sub.1,n.sub.2). 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 .DELTA. (n.sub.1,n.sub.2) 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 .DELTA. (n.sub.1,n.sub.2), i.e. .alpha. and 
.beta., 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 .beta. in equation 
(12) becomes rather large. This error manifests itself as large magnitude 
for .beta., resulting in unrealistic values of the phase difference 
function. Since experimental results indicate that the magnitude of .beta. 
is small for most columns, the present invention sets to .beta.=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, .beta. 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.sub.1,i.sub.2, . . . , i.sub.upper &lt;N/2and 
j.sub.1,j.sub.2, . . . , j.sub.lower .gtoreq.N/2 denote the indices of the 
ghosting pixels of the n.sub.1 th column, taking into account equation 
(12), the linear least-squares estimation of .alpha. and .beta. consists 
of solving the following overdetermined linear system of equations: 
##EQU14## 
If .alpha. (n.sub.1) and .beta. (n.sub.1) denote the solution of the above 
linear least-squares problem, the residue of the ghosting pixels at 
(n.sub.1,i.sub.1) and (n.sub.1,j.sub.1) are defined to be: 
EQU res(n.sub.1,i.sub.1).ident..DELTA.(n.sub.1,i.sub.1)-(.alpha.(n.sub.1)+.beta 
.(n.sub.1)i.sub.1) 
EQU res(n.sub.1,j.sub.1).ident..DELTA.(n.sub.1,j.sub.1)-(.alpha.(n.sub.1)+.beta 
.(n.sub.1)N-.beta.(n.sub.1)j.sub.1) (19) 
and the means square error is defined to be 
##EQU15## 
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.sub.1,n.sub.2). If this ration is close to 1 (or more 
specifically is in the range [1/eoratio, eoratio]) then either 
Y(n.sub.1,n.sub.2) or Y(n.sub.1,(n.sub.2 +N/2) mod N) are classified as 
ghosting pixels. If the phase difference function is less than .pi. 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 .pi., 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 .pi.) 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.sub.1,n.sub.2) is 
the ghosting pixel and Y(n.sub.1,n.sub.2 +N/2) is the ghosted pixel. On 
the other hand if the ratio of equation (17) is smaller than 1/eoratio, 
then Y(n.sub.1,n.sub.2 +N/2) is the ghosting and Y(n.sub.1,n.sub.2) 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 .pi., which in our situation corresponds to the central 
columns of the image. For columns away from the center, the phase 
difference function might become larger than .pi., 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 .pi.. 
Fortunately, this is a simple visual task since in the areas of image with 
large phase difference function (i.e. larger than .pi.) 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.sub.1 th column is different from those of the past four 
columns by more than an amount specified by "smooth", then .alpha. 
(n.sub.1) and .beta. (n.sub.1) are discarded and the n.sub.1,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 .alpha. and .beta. 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 .alpha. and .beta.. 
The "-b" option sets the internal double precision parameter "betasmth" 
which smoothes the values of .beta. 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: 
EQU 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: 
EQU 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: 
EQU 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: 
EQU 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: 
EQU 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: 
EQU correct -g -t1. -r1.5 -s5. -m3. -f -d.6 -b.06 -a1000 e4 lvrnum lvrproc 
A seventh example 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: 
EQU 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: 
EQU 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##