Magnetic Resonance Imaging (MRI) has become a widely accepted and commercially available technique for obtaining digitized visual images representing the internal structure of objects (such as the human body) having substantial populations of atomic nuclei that are susceptible to nuclear magnetic resonance (NMR) phenomena. In MRI nuclei in a body to be imaged are polarized by imposing a strong main magnetic field Ho on the nuclei. Selected nuclei are excited by imposing a radio frequency (RF) signal at a particular NMR frequency. By spatially distributing the localized magnetic fields, and then suitably analyzing the resulting RF responses from the nuclei, a map or image of relative NMR responses as a function of the location of the nuclei can be determined. Following a Fourier analysis, data representing the NMR responses in space can be displayed on a CRT.
As shown in FIG. 1, the NMR imaging system typically includes a magnet 10 to impose the static magnetic field, gradient coils 12 for imposing spatially distributed magnetic fields in the three orthogonal coordinates, and RF coils 14 to transmit and receive RF signals to and from the selected nuclei. The NMR signal received by the coil 14 is transmitted to a computer 16 which processes the data into an image displayed on the display 18. The magnetic resonance image is composed of picture elements called "pixels." The intensity of a pixel is proportional to the NMR signal intensity of the contents of a corresponding volume element or "voxel" of the object being imaged. The computer 16 also controls the operation of the RF coils 14 and gradient coils 12 through the RF amplifier/detector 20 and gradient amplifiers 22, respectively.
Each voxel of an image of the human body contains one or more tissues. The tissues of the human body are comprised primarily of fat and water. Fat and water have many hydrogen atoms (the human body is approximately 63% hydrogen atoms). Since hydrogen nuclei have a readily discernible NMR signal, magnetic resonance imaging of the human body primarily images the NMR signal from the hydrogen nuclei.
In NMR, a strong static magnetic field is employed to align atoms whose nuclei have an odd number of protons and/or neutrons, such that the nuclei have a spin angular momentum and a magnetic dipole movement. A second RF magnetic field, applied as a single pulse transverse to the first field, pumps energy into the nuclei, which causes the angular orientation of the nuclei to flip by, for example, 90.degree. or 180.degree.. After this excitation, the nuclei precess and gradually relax into alignment with the static field. As they relax and precess, the nuclei emit energy in the form of weak but detectable free induction decay (FID). These FID signals and/or RF or magnetic gradient refocused "echoes" thereof (collectively referred to herein as MR signals) sensed by an NMR imaging system are used by a computer to produce images of the body.
The excitation frequency and the FID frequency are related by the Larmor relationship. This relationship states that the angular frequency, .omega..sub.0, of the precession of the nuclei is the product of the magnetic field, B.sub.0, and the so-called magnetogyric ratio, .gamma., a fundamental physical constant for each nuclear species: EQU .omega..sub.0 =B.sub.0.multidot..gamma.
By superimposing a linear gradient magnetic field, B.sub.Z =Z.multidot.G.sub.Z on the static uniform field, B.sub.0, which is typically defined as the Z axis, for example, nuclei in a selected X-Y plane may be excited by proper choice of the frequency of the transverse excitation field applied along the X or Y axis. Similarly, a gradient magnetic field can be applied in the X-Y plane during detection of the MR signals to spatially localize emitted MR signals from the selected X-Y plane according to their frequency (and/or phase).
An example of an operation whereby the various coils produce the MR signal in a 3-D MRI is shown in FIG. 2A, which is a graphical representation of an example MRI acquisition sequence called a "field-echo" sequence. First, a gradient field, G.sub.slice, is superimposed along the main field to sensitize a "slab" of nuclei in the body to be imaged to a particular RF resonance frequency. An RF "nutation" pulse is then applied at the particular frequency to force some of the nuclei within the slab to precess in a perpendicular direction with respect to the main field. Thereafter, pulsed magnetic gradient fields of changing magnitudes, G.sub.pe and G.sub.slice, are used to phase encode the nuclei by inducing frequency differences, and hence phase differences, between nuclei in different locations along a direction within the slab. A gradient field applied perpendicular to the pe direction, in a "ro" (or "read out") direction, G.sub.ro, first de-phases and then rephases the nuclei to form "field-echo" MR signals. During the field-echo, the applied gradient field, G.sub.ro, frequency encodes the selected slab of nuclei in the readout direction. The resultant MR signal (also called "raw data" or "k-space data") is then read and analyzed by Fourier analysis. The frequency (domain) plot of that analysis is then scaled to render information about the nuclei population in Fourier space, which corresponds to an X-Y-Z position.
FIG. 2B shows another example of a 2-D MRI surface. In this case, a slice of nuclei are selected by a 90.degree. followed by a 180.degree. RF to generate a spin-echo. When the time between the 180.degree.-center and the echo-center is the same as that between the 90.degree.-center and the 180.degree.-center, a symmetric spin-echo is generated. However, the echo-center can be shifted by adjusting the field gradients so that an asymmetric spin-echo is produced instead.
In addition to using the frequency information content of an MR signal to generate images, the phase of an MR signal in frequency domain can be utilized to provide information indicative of some physical quantity. For example, depending on the type of pulse sequence used, the MR phase can represent a main B.sub.0 field inhomogeneity or can be proportional to the velocity of the moving spins.
As mentioned above, the main magnetic field can be altered by gradient magnetic fields created in the X, Y, and Z directions of the imaging volume. For the purpose of simplifying the descriptive mathematics, a rotating reference coordinate system X'-Y'-Z', that rotates at the nominal Larmor frequency about the Z' axis, is often used to describe nuclear phenomenon in NMR. Selected nuclei, which precess in alignment with the B.sub.0 field, are influenced (nutated) by the perpendicular magnetic field of an RF pulse at the Lamor frequency, causing a population of such nuclei to tip from the direction of the magnetic field B.sub.0. Thus, certain nuclei start aligned with the "Z'" axis by the static B.sub.0 field and then are rotated to the X'-Y' plane as a result of the RF signal being imposed on them. The nuclei then precess in the X'-Y' plane.
The RF spin-nutating signal will, of course, tip more than one species of the target isotope in a particular area. For the purpose of simplifying the description in analyzing an MRI sequence, each RF pulse can be characterized by its center where nuclei can be considered all in-phase and whereafter each species of nuclei precess at their own characteristic speed. The phase of the processing nuclei species will gradually differ (de-phase) as a result of such parameters as the physical or chemical environment that the nuclei are located in. Nuclei in fat, for example, precess at a different rate than do nuclei in water due to their difference in chemical shift. Field gradients and inhomogenieties also contribute to this de-phasing.
MRI sequences are almost always designed to have the field gradients fully refocused at the echo-center, either by reversing the gradients to create field-echoes or by applying RF refocusing pulses to create spin-echoes. Inhomogeneities in the B.sub.0 field and chemical shift are refocused in symmetric spin-echo, but not in field echoes and asymmetric spin-echoes.
Once the spins are disturbed from their equilibrium, processes known as relaxation cause the Z'-component to recover to its equilibrium magnitude, M.sub.0, and the phase-coherent component of magnetic moments in the X'-Y' plane to decay. These processes are usually characterized by exponentials whose time constants are called the spin-lattice relaxation time, T.sub.1, and the spin-spin relaxation time, T.sub.2, respectively. When magnetic resonance signals are observed through magnetic flux oscillation in a plane coexistent with the X'-Y' plane, both of these processes cause the signal strength to decrease as a function of time. An apparent relaxation time, T.sub.2 *, is used to characterize the transverse signal decay due to both the spin-spin relaxation and B.sub.0 inhomogeneities.
A. Water and Fat Pixels
Although MR images of both water and fat may contain the same or different diagnostic information, they often interfere with each other's interpretation when overlapped in an MRI image, and thus make it difficult to properly interpret the composite MR image. The difference in one tissue versus another can be used as information to separate MRI images of fat component and tissue from fluids or water-based tissue (for these purposes, "water-based tissue" and "fluids" are used interchangeably).
When water and fat signals are separated quantitatively in each image pixel, the resultant images are called true water-only and fat-only images. A less vigorous approach is to separate the image into pixels which contain more water than fat (herein referred to as water pixels) and pixels which have more fat than water (herein referred to as fat pixels). Such separation may already be sufficient for certain kinds of diagnostic evaluations.
B. Three-Point Dixon Methods
At high-field strengths, the separation of water and fat images or suppression of one of these two components can be achieved using selective excitation or non-excitation approaches. However, at mid- or low-field strengths, approaches based on chemical shift selectivity become impractical, if not impossible. At all-field strengths, the difficulties of water/fat image separation are further exacerbated when there are large magnetic field inhomogeneities.
One promising group of techniques, known as "Three-Point Dixon" methods, have attractive features for mid- or low-field strength applications. These methods require three images to obtain enough information for water/fat separation with correction of the effect of B.sub.0 field inhomogeneities. The Images can be acquired using spin-echo as well as field-echo sequences, in three separate scans as described in "Three-Point Dixon Technique for True-Water-Fat Decompositions with B.sub.0 Inhomogeneity Corrected," by Glover et al., Magnetic Resonance in Medicine 18, 371-383 (1991); or in two scans, "Separation of True Fat and Water Images by Correcting Magnetic Field Inhomogeneity In Situ", by Yenng et al., Radiology 159, 783-786 (1986), or in a single scan, "True Water and Fat MR Imaging with Use of Multiple-Echo Acquisition", by William et al., Radiology 173, 249-253 (1989) and "Separation of Water and Fat MR Images in a Single Scan at 0.35 T Using `Sandwich` Echoes," by Zhang et al, JMRI 6, 909-917 (1996), all the above of which are incorporated herein by reference.
In accordance with the above Three-Point Dixon method, acquisitions of the three images are controlled so that the difference in phase between the water image and the fat image changes by .pi. consecutively between three images. The image data may be represented as follows: EQU S.sub.--.pi. =(W-Fe.sup.-i.psi.) e.sub.0.sup.-i(.phi.-.phi.) EQU S.sub.0 =(W+Fe.sup.-i.psi.) e.sub.0.sup.-i.phi. EQU S.sub..pi. =(W-Fe.sup.i.psi.) e.sub.0.sup.-i(.phi.+.phi.)
where S.sub.--.pi., S.sub.0 and S.sub..pi. represent the three MRI images; W and F are the water and fat components, respectively; .psi. is the difference in phase between the water and fat signals in S.sub.O ; .phi..sub.o is the phase in S.sub.O due to B.sub.0 -offset induced and other system sources; .phi. is the B.sub.0 -offset induced phase change between S.sub.--.pi. and S.sub.O and between S.sub.O and S.sub..pi..
Images are then processed to obtain information to remove the effects of B.sub.0 inhomogeneities and finally to separate the water and fat images: B.sub.0 inhomogeneities are determined from the three images by a process of "phase unwrapping" as will be described further herein in greater detail: EQU .phi.=1/2 unwrap {arg(S.sub..pi..multidot.S.sub.-.pi.)}
where arg produces the phase angle of a complex number.
Water images and fat images can then be reconstructed in accordance with the following relationships: EQU W=S.sub.O +0.5 S.sub..pi. e.sup.+i.phi. +0.5 S.sub.--.pi. e.sup.-i.phi. EQU F=S.sub.O -0.5 S.sub..pi. e.sup.i.phi. -0.5 S.sub.--.pi. e.sup.-i.phi.
C. One-Point Dixon Methods
One-point Dixon methods, which form the subject matter of the present invention, acquire only one image in which the water and fat signals are out-of-phase, represented by: EQU S=(W-F)e.sup.i.phi.
Moreover only one image is available, attempts may still be made to determine the B.sub.0 inhomogeneities by phase unwrapping before separating the water and fat pixels: EQU .phi.=1/2 unwrap {arg(S.multidot.S)}
Water pixels (I.sub.water-pixel) and fat pixels (I.sub.fat-pixels) are then separated according to: EQU I.sub.water-pixel =.vertline.s.vertline.+se.sup.i.phi. EQU I.sub.fat-pixel =.vertline.s.vertline.-se.sup.i.phi.
D. Phase Unwrapping
Basically, phase unwrapping is a process of determining the absolute phase of a complex signal given the measurement of its principal phase value.
Since the phase angle of a complex number is unambiguous only between -.pi. and .pi., the phase of an MRI signal cannot be unambiguously determined from its argument. Any phase values beyond -.pi. and .pi. will be wrapped back around into values between -.pi. and .pi..
Preferred algorithms for phase unwrapping as implemented in the present invention involve a combination of modeling the static magnetic field using polynomial functions and a guided phase unwrapping by region-growing.
i. Polynomial Fitting
The magnetic field is modeled using a polynomial function: ##EQU1##
To find the coefficients a.sub.n and b.sub.n, partial spatial derivatives of the phase value .phi. are calculated and fit to polynomial functions as follows: ##EQU2##
Fitting is performed using a weighted least-square with the weighting factors determined according to: ##EQU3##
where S(x,y) is the pixel value in the in-phase image and S.sub.max is the maximum of that image.
Using p.sub.n and q.sub.n, a.sub.n and b.sub.n are calculated according to the following equations: ##EQU4##
ii. Binary Phase Unwrapping
If it can be assumed that the magnetic field fitting is relatively accurate within a small error range, for example, .+-.0.2.pi., then unwrapping can be performed by simply comparing the measured phase (with the predicted phase .phi..sub.p : EQU .DELTA..phi.=.phi..sub.p -.phi.
If .vertline..DELTA..phi..vertline.&gt;.pi., then .phi..sub..function. fused for water and fat image reconstruction is determined by ##EQU5##
where integer promotes the resulting quotient to whole number.
iii. Unwrapping by Region Growing
However, the field fitting may contain large errors which will cause cumulative errors in phase-unwrapping and consequently result in water/fat mutual contamination in the final images. To unwrap in a more fool-proof way, a "region growing" algorithm is implemented as follows:
(a) A pixel in the image is chosen as the subseed for unwrapping and the measured phase value is assigned to the final phase value used for water and fat image reconstruction. EQU .phi..sub.f (x.sub.o,y.sub.o)=.phi.(x.sub.o,y.sub.o) PA1 (b) The subseed is selected so that all pixels in a 6.times.6 region centered at the subseed have sufficient signal strength. The four immediately neighboring pixels of the subseed are first unwrapped by comparing the phase values to the subseed value. If the difference is larger than EQU .phi..sub..function. =.phi.+sign(.DELTA..phi.).times.2.pi. PA1 a predetermined threshold, a 2T unwrapping is executed: PA1 (c) A 3.times.3 pixel region centered at the subseed is then built based on the five already-determined pixels. A 2.pi. phase unwrapping is executed whenever the phase difference between a pixel and its already-unwrapped, immediate neighboring pixel is larger than the threshold. PA1 (d) The 3.times.3 region is then expanded to 4.times.4 region with a three-pixel prediction: EQU .phi..sub.p =1/3(.phi..sub.f.sup.-3 +.phi.f.sup.-2 +.phi..sub.f.sup.-1)+1/3(3.DELTA..phi..sup.-1 +2.DELTA..phi..sup.-2 +.DELTA..phi..sup.-3) PA1 where .phi..sub.p is the predicted phase value of a pixel; .phi..sub.f.sup.-i (i=1, 2, 3) are the phase values of the already unwrapped first (i=1), second (i=2), and third (i=3) neighboring pixels; .DELTA..phi..sup.-i are the phase spatial derivatives of the already-unwrapped, neighboring pixel along the direction of the prediction. PA1 (e) Continuing from the 4.times.4 seed region, a four column by four rows "cross" region is built using a four-pixel prediction: EQU .phi..sub.p =1/4(.phi..sub.f.sup.-4 +.phi..sub.f.sup.-3 +.phi..sub.f.sup.-2 +.phi..sup.-1)+1/4(.DELTA..phi..sub.f.sup.-3 +3.DELTA..phi..sub.f.sup.-2 +2.DELTA..phi..sub.f.sup.-3 +.DELTA..phi..sub.f.sup.-4) PA1 (f) Using the cross, the four quadrants of the image are unwrapped using the same 4-pixel prediction approach, but in two directions. Unwrapping is executed when both directions show the same execution for unwrapping. In other situations, the average of the predicted values is used. When the pixel value is below the intensity threshold, the phase value is again set to the predicted average value.
where: .DELTA..phi.=.phi.-.phi.(x.sub.0,y.sub.0)
Unwrapping is executed if .DELTA..phi.=.phi.-.phi..sub.p is larger than the threshold.