US RE36679 E Abstract A method of cancelling ghosts from NMR images. The method involves estimating a phase difference function Δ (n
_{1}, n_{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_{2},n_{2}) and x(n_{1},n_{2} +N/2), respectively, where the image has dimensions N×N_{s}. Experimental values of Δ (n_{1}, n_{2}) for a variety of objects indicate that its variation along n_{1} is considerably larger than along n_{2}. Thus, for each column n_{1}, the phase difference function Δ (n_{1}, n_{2}) can be modelled as a one-dimensional function of n_{2} with two parameters α (n_{1}) and β (n_{1}), which are estimated from the pixels in the 2-D FFT processed reconstructed image Y(n_{1},n_{2}). These parameters are then used to estimate Δ (n_{1}, n_{2}), which is ultimately used to de-ghost the image.Claims(12) 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 _{1},n_{2});(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 _{1}) and β (n_{1}) for.]. .Iadd.identifying ghosting pixels in .Iaddend.each remaining column of data; i.e. n_{1} =0, . . . N_{s} -1, .[.by:.]..[.(i).]. .Iadd.(e) .Iaddend.finding the .Iadd.actual .Iaddend.phase difference .[ƒ]. Δ.Iadd. _{actual} .Iaddend.(n_{1},n_{2}) 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 _{1}) and β (n_{1}): ##EQU18## .[.(e).]. .Iadd.(g) .Iaddend.using .Iadd.the linear least square estimates of .Iaddend.α (n_{1}) and β (n_{1}) .[.in the above equation.]. to find .[.the.]. .Iadd.an estimated .Iaddend.phase difference .Iadd.function .Iaddend.Δ (n_{1},n_{2}) for 0≦n_{2} <N.[.; and.]. .Iadd.using .Iaddend. ##EQU19## .[.(f).]. .Iadd.(h) .Iaddend.using .Iadd.the estimated phase difference function .Iaddend.Δ (n_{1},n_{2}) in ##EQU20## to find A(n_{1},n_{2}) and B(n_{1},n_{2}) for 0≦n_{2} <N, where the dimensions of the .[.reconstructed.]. image .Iadd.data x(n_{1},n_{2}) 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
_{1},n_{2}) 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
_{1}, n_{2}) 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_{1},n_{2}) is calculated using:Δ _{actual} (n_{1},n_{2})=phase (Y_{even} (n
_{1},n_{2}))-phase (Y_{odd} (n_{1},n_{2}))..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 _{1},n_{2}) 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_{1},n_{2}) between Y_{even} (n_{1},n_{2}) and Y_{odd} (n_{1},n_{2})..Iaddend..Iadd.20. The method of claim 19 wherein the actual phase difference Δ_{actual} (n_{1},n_{2}) is calculated using: Δ_{actual} (n_{1},n_{2})=phase (Y_{even} (n_{1},n_{2}))-phase (Y_{odd} (n_{1},n_{2}))..Iaddend..Iadd.21. The method of claim 19 wherein the step of processing the NMR data comprises using the transform data Y(n_{1},n_{2}) and the actual phase difference Δ_{actual} (n_{1},n_{2}) to calculate the image data, wherein the image data are denoted by x(n_{1},n_{2})..Iaddend..Iadd.22. The method of claim 21 wherein the step of processing the NMR data comprises estimating a phase difference function Δ(n_{1},n_{2}) from the actual phase difference Δ_{actual} (n_{1},n_{2})..Iaddend..Iadd.23. The method of claim 22 wherein the estimated phase difference function Δ(n_{1},n_{2}) is determined by(i) calculating a plurality of values of the actual phase difference Δ _{actual} (n_{1},n_{2});(ii) using the plurality of values of the actual phase difference Δ _{actual} (n_{1},n_{2}) to find linear least square estimates of α(n_{1}) and β(n_{1}) using ##EQU23## (ii) using the linear least square estimates of α(n_{1}) and β(n_{1}) to find the estimated phase difference function Δ(n_{1},n_{2}) for 0≦n_{2} <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
_{1},n_{2}) and the estimated phase difference function Δ(n_{1},n_{2}) to calculate the image data x(n_{1},n_{2})..Iaddend..Iadd.25. The method of claim 24 wherein the image data x(n_{1},n_{2}) are calculated by: (i) using the estimated phase difference function Δ(n_{1},n_{2}), Y_{even} (n_{1},n_{2}), and Y_{odd} (n_{1},n_{2}) to find A(n_{1},n_{2}) and B(n_{1},n_{2}) using ##EQU25## where ##EQU26## (ii) using the respective magnitude of A(n_{1},n_{2}) and B(n_{1},n_{2}) to find x(n_{1},n_{2}) and x(n_{1},n_{2} +N/2)..Iaddend..Iadd.26. 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 _{1},n_{2});(b) identifying ghosting pixels in columns of Y(n _{1},n_{2}); i.e. n_{1} =0 . . . N_{s} -1,(c) finding the actual phase difference Δ _{actual} (n_{1},n_{2}) for the ghosting pixels(d) solving the following simultaneous equations to find linear least square estimates of α(n _{1}) and β(n_{1}); ##EQU27## (e) using the linear least square estimates of α(n_{1}) and β(n_{1}) to find an estimated phase difference function Δ(n_{1},n_{2}) for 0≦n_{2} <N using ##EQU28## (f) using the estimated phase difference function Δ(n_{1},n_{2}) in ##EQU29## to find A(n_{1},n_{2}) and B(n_{1},n_{2}) for 0≦n_{2} <N, where the dimensions of the image data x(n_{1},n_{2}) 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
_{1},n_{2});(b) calculating Y _{even} (n_{1},n_{2}) and Y_{odd} (n_{1},n_{2}) using ##EQU31## (c) estimating a phase difference function Δ(n_{1},n_{2}) between Y_{even} (n_{1},n_{2}) and Y_{odd} (n_{1},n_{2}); and(d) using the estimated phase difference function Δ(n _{1},n_{2}) in ##EQU32## to find A(n_{1},n_{2}) and B(n_{1},n_{2}) for 0≦n_{2} <N, where the image data corresponding to the NMR image are denoted by x(n_{1},n_{2}), 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
_{1},n_{2})..Iaddend.Description 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 When NMR data is obtained by scanning the k 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 Having modeled the ghosted image, the objective can be stated as estimation of the true object density distribution x(n From equations (1) and (2), it is clear that if the phase functions φ (n
x(n Substituting this into equation (1) and (2): ##EQU2## The phase difference between Y Similarly, by placing the object in the lower half of the FOV: ##EQU4## The phase difference between Y Thus, experimental values of Δ (n 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 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. 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 (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 (i) finding the phase difference function Δ (n (ii) solving the following simultaneous equations to find linear least square estimates of α (n (e) using α (n (f) using Δ (n From equations (9) and (10), the true object density distribution at (n 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. 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 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 The algorithm of the present invention takes advantage of the above approximation by estimating α (n 1. If raw sampled data is input, take 2-D inverse Fourier transform of the time data to obtain the ghosted image Y(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 5. Use α (n 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 -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 Step 4. of the algorithm is now described in more detail. As set forth in the Background section above, the phase difference function Δ (n 1. It corresponds to a high energy point in the object. Therefore
|x(n 2. The pixel at location (n
|x(n Step 4 of the automatic ghost correction algorithm derives the parameters associated with the n 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
|Y Thus, if the ratio between the magnitude of even and odd parts of the pixel (the "eoratio") at location (n 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 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 or equivalently Δ (n the observed image, Y(n Specifically, as Δ (n 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 the internal variable "threshold" set by "-t" option, if n decreases linearly with |n is equal to 1 for |n 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 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 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 If α (n 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 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 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## Patent Citations
Referenced by
Classifications
Legal Events
Rotate |