Geometric distortion correction in magnetic resonance imaging

[0002] This patent application claims the benefit of the Apr. 24, 1998 filing date of U.S. Provisional Application No. 60/083,010, which is herein incorporated by reference.

[0003] The present invention relates generally to the field of magnetic resonance imaging (MRI). More particularly, the present invention relates to the field of geometric distortion correction in MRI.

[0004] Neurosurgical stereotaxis refers to a collection of techniques used to superimpose radiological volume images intra-operatively to the patient's skull, brain, and cerebrovasculature. Typical stereotactic procedures begin with the collection of a three dimensional (3D) computer tomography (CT) image containing fiducial markers. The resulting 3D image space may then be used to determine an optimal craniotomy or skull opening site and to model virtually the potential intra-operative orientation of a surgical tool, such as a biopsy needle for example. Using such images helps reduce the size of the craniotomy and intracranial traverse of surgical tools.

[0005] For a typical procedure, a stereotactic frame or frameless fiducial imaging markers are mounted on the patient's head, often by screwing either into the skull. These fiducial markers are then used by a three dimensional (3D) digitizer to register the image-based surgical plan to the patient. Using the fiducials, the region of interest can be mapped to the coordinate system in the pre-operative image. The digitizer can be used to guide the craniotomy and initial insertion of surgical tools virtually from the pre-operative image.

[0006] Typical stereotactic surgical procedures emphasize accurate location of anatomical structures within pre-operative images. The accuracy of a neurosurgical navigational system depends on, for example, the mechanical accuracy of the 3D digitizing probe, the image registration algorithm, geometrical distortion of the pre-operative images, movement of the patient during scanning, movement of the patient with respect to the patient space determined with the 3D digitizing system, and movement of the brain between scanning and surgery.

[0007] Using CT for stereotactic pre-operative imaging helps ensure relatively high spatial fidelity is obtained between the patient's images and the patient's anatomy in the operative stereotactic guidance device. CT techniques use a beam of high-energy X-ray photons to probe the patient and display an image based on the internal influences on the passage of photons through the patient. CT may be used to generate a series of images containing a distribution of the X-ray attenuation co-efficient, which relates to the electron density of the atoms, of a 3D volume. Such image reconstruction is based on the principles of line-of-sight ray optics. CT therefore provides images of relatively high positional accuracy.

[0008] CT is useful in imaging bony structures yet has limited ability to differentiate between components of inhomogeneous soft tissue structures, such as the brain for example. CT is capable of differentiating tissues only insofar as the region of interest's X-ray attenuation coefficient differs from its surroundings. Diverse adjacent soft tissue structures, such as muscles, mucosa, nerves, and blood vessels for example, typically appear as a solid mass when imaged. As CT also involves X-ray exposure, its frequency of use must be limited to avoid tissue burning and genetic damage or mutation, for example, as a result of ionizing radiation overexposure.

[0009] Magnetic resonance (MR) is based on the presence of magnetic moments in some nuclei of atoms of the patient or other sample being imaged. The protons in many nuclei behave as tiny bar magnets or magnetic dipoles. In the absence of any external influences, these magnetic dipoles are randomly oriented with zero net magnetization. When placed in an external magnetic field, the nucleus precesses around that field at a gravitational field. The precession frequency, ω, depends on the strength of the external magnetic field B₀as follows:

[0010] ω=γB₀  Equation 1

[0011] where ω is the resonance or Larmor frequency of the nuclei of interest, γ is the gyromagnetic ratio of the nuclei of interest, and B₀is the applied magnetic field strength.

[0012] In magnetic resonance imaging (MRI), the sample is exposed to a radio frequency (RF) signal of the same precession frequency. The sample absorbs energy and changes its alignment relative to the applied field. Upon removal of the excitation signal, the magnetic dipoles relax and fall back to their original entropic alignment. As this occurs, each material in the sample emits a characteristic resonant frequency which can be detected and measured. The relaxation time required for each tissue sample depends on its proton environment as protons that may pass on the energy to their neighboring atoms with relative ease can return to their original state relatively faster.

[0013] MRI maps the proton density distribution weighted by T1 and T2 relaxation effects in the sample. Hydrogen, an abundant element in live tissue, has a nucleus comprising a single proton. Because of the freedom of hydrogen-containing-molecules to precess, MR images of living samples highlight different time dependent development of recovery of magnetization of the heterogeneous soft tissues in such samples.

[0014] Typical MRI techniques apply three orthogonal linearly varying magnetic field gradients to the nuclear magnetic dipoles of the sample in the presence of a homogeneous static magnetic field to define the location of the proton density. the sample is divided into consecutive planes or slices. The protons in an image slice can be selectively excited by a combination of a frequency selective pulse and a field gradient perpendicular to the image plane, termed the slice select gradient. The bandwidth of the RF pulse and the amplitude of the slice select gradient determine the thickness of the slice as follows: $\begin{array}{cc}\Delta z=\frac{{\mathrm{BW}}_{\mathrm{rf}}}{{\mathrm{\gamma G}}_z}& \mathrm{Equation}2\end{array}$

[0015] where Δz is the slice thickness, BW_rfis the bandwidth of the RF pulse, and G_z, is the slice select gradient. The slice thickness can be reduced by increasing the gradient strength or decreasing the RF bandwidth.

[0016] After the slice selection, a second magnetic field gradient can frequency encode the position of the proton density signal in one direction within the plane. The frequency encoding gradient, also known as the read-out gradient, causes the magnetic field to vary linearly from one end of the image plane to the other in the applied direction. As an extension of Equation 1, the precession frequency of the net magnetization vector varies from point to point along this direction as follows:

[0017] ω(y)=γ(B₀+yG_y)  Equation 3

[0018] where y is the position along the frequency encoded direction and G_yis the frequency encoding gradient.

[0019] The second dimension in the image plane is encoded by phase encoding. Phase of the magnetization vector can be defined as the state of precession of the vector at a given point in time. The application of a third orthogonal gradient, also known as the phase encoding gradient, causes the magnetization vectors in that direction to precess at different frequencies. When the phase encoding gradient is removed after a short time, the vectors precess again at the same frequency although the phase of each vector has been changed in proportion to the applied gradient.

[0020] Within a slice plane, the phase encoding gradient and frequency encoding gradient are applied in orthogonal directions resulting in a data matrix of amplitude points. The phase and frequency at point (x,y) are related to their spatial position and gradient strength as follows:

[0021] Φ_x=γ_xG_xΔt_xx  Equation 4

[0022] where Φ_xis the phase at point (x,y), G_xis the phase encoding gradient, Δt_xis the phase encoding period, and x is the position along the phase encoded direction, and

[0023] ω_y=γ_yG_yγ  Equation 5

[0024] where ω_yis the frequency at point (x,y), G_yis the frequency encoding gradient, and y is the position along the frequency encoded direction.

[0025] For one typical MRI technique, a human subject is placed relative to a magnet that orients the proton magnetic dipoles. A field gradient is imposed along the z-axis in the direction of the main magnetic field such that a narrow plane of protons resonate within a band of frequencies. A phase encoding gradient along the x-axis is activated for a short time during which the dipoles acquire a different phase. A frequency encoding gradient along the y-axis is then activated to frequency encode the positions of the dipoles while a receiver coil is activated to record the signal. Typically, 128, 256, or 512 data points are recorded along the frequency encode axis. As each recorded data point corresponds to a respective pixel in the image to be generated, the number of recorded data points determines the resolution of that image.

[0026] A Fourier Transformation (FT) algorithm is used to decode frequency information contained in a proton signal at each location in the imaged plane to corresponding intensity levels. The resulting image is then displayed as shades of gray in a matrix arrangement of pixels.

[0027] For 3D FT imaging, slice selection is obviated by simultaneously acquiring data from the entire imaging volume. The slice selection excitation pulse of the 2D imaging technique is replaced by a weaker gradient pulse which is simultaneously recorded from the entire image volume. Partitioning into individual slices is then achieved by stepping an additional phase encoding gradient, analogous to the phase encoding for 2D imaging. This additional gradient is applied orthogonal to the plane defined by the conventional phase and frequency encoding gradients.

[0028] As MR radio frequency electromagnetic radiation can penetrate the skull relatively well, MRI may be used for neurosurgery. MRI may be used to image the interior of the cranial vault, for example, without significant attenuation.

[0029] In contrast to CT, MRI involves no ionizing radiation and helps provide better definition to inhomogeneous soft tissue structures due to the dependence of the MR signal on nuclear proton density, T1 and T2 relaxation times, and fluid flow through tissue. Through instrument design and operator control, the influence of these parameters on the MR image may be emphasized to better discriminate between tissues having similar density. As one example, gray and white matter may be discriminated on the basis of their significantly different T1 values. As another example, tumor, edema, cerebrospinal fluid, gray matter, and white matter may be contrasted in T2 weighted images.

[0030] MRI may also be used to reconstruct images for any standard orientation, such as transverse, coronal, or sagittal slices for example, or for any oblique orientation. CT techniques typically collect images in the axial plane, although spiral-scanned volumes may be reformatted to generate images in other orientations without repositioning the patient within the CT scanner. Some image quality may be lost during such reformatting to avoid more than the minimum X-ray exposure.

[0031] Typical MRI techniques generate images having a 256×256 pixel resolution. Higher image resolutions may lead to a lower signal-to-noise ratio (SNR). Typical CT techniques generate images having a 512×512 pixel resolution with relatively good SNR.

[0032] Scan time in typical MRI techniques is relatively longer than for an equivalent 3D CT spiral scan. Scan time T_scandepends on the repetition time T_R, the number of phase encodings N_y, and the number of signal acquisitions NSA as follows.

[0033] T_scan=T_RN_yNSA(*N_zin 3D)  Equation 6

[0034] Depending on the chosen acquisition parameters, such as T1, T2, and proton density for example, the scan time T_scancan range from seconds to minutes and is typically significantly greater than for an equivalent 3D CT spiral scan. Fast MRI techniques, such as Echo Planar Imaging for example, have been proposed with most having significantly lower SNR. Slower MR scan times increase the likelihood of blood flow and patient motion artefacts.

[0035] CT techniques may produce within-plane beam artefacts from dense metal objects while MRI is not suitable for imaging subjects that contain ferrous objects such as dental, cardiac pacemaker, or other common surgical implants.

[0036] MRI depends on relatively highly uniform magnetic fields. Any non-uniformity in these magnetic fields manifests as distortion in the reconstructed image. Because an accurate representation of location and size of organ sub-structures in images is crucial for stereotactic applications, such as multimodality-image-superposition and image-guided-surgery for example, geometric distortions inherent to MRI are to be considered a for such spatially quantitative applications. Geometric distortions in MR images result from non-uniformity introduced into the magnetic field from different sources, such as gradient field non-linearities, chemical shift resonance frequency offsets, and magnetic field inhomogeneity resonance frequency offsets for example.

[0037] Typical MRI systems generate gradient magnetic fields with currents passing through three orthogonal coils. Despite the relatively high precision fabrication of the coils, the resulting gradient fields have some residual non-linearity which may manifest differently depending on the imaging method. Gradient field non-linearities manifest in three ways in 2D MRI, as illustrated in FIGS. 2A, 2B, and 2C.

[0038] As illustrated in FIG. 2A, a transverse distortion in the slice plane results in a barrel aberration. As illustrated in FIG. 2B, gradient field non-linearities can cause gradual increase in slice thickness as one moves away from the gradient iso-center, resulting in a bow tie effect. As illustrated in FIG. 2C, a gradient field non-linearity in the slice selection direction results in a potato chip effect where the reconstructed image planes are wrapped into potato chip like shapes and may be interpreted as a volume.

[0039] Barrel aberration and bow tie distortions may be corrected with suitable mathematical techniques. Correction of the potato chip effect distortion, however, is more complex as the slices are not contiguous in most scans to avoid RF saturation effects.

[0040] Distortion in MRI also results due to chemical shift which occurs when materials, perhaps different portions of a single scanned object, appear to have moved relative to each other in an image. Water and fat are two basic types of molecules whose hydrogen atoms give rise to MR signals. Due to the different chemical environments, hydrogen atoms attached to different molecules resonate at slightly different frequencies while in the same applied magnetic field B₀. Inhomogeneous shifts between fat and water-rich soft tissue regions, for example, may occur in the image. Selective saturation or selective excitation methods may be used to account for chemical shift. Selective suppression of either fat or water-rich soft tissue regions yields a water-only or fat-only image, respectively. Selective excitation of either fat or water-rich soft tissue regions yields a fat-only or water-only image, respectively. Chemical shift generally does not pose a significant problem for stereotactic neurosurgical localization procedures.

[0041] Magnetic field inhomogeneity distortion is more complex as it depends on both the apparatus and the imaged object. Scanner based inhomogeneity can be caused by residual imperfections which cause non-linearities in the generated magnetic field. Shimming is a process whereby appropriate corrective currents are applied in shim coils to restore the homogeneity of the magnetic field. Improvements in magnet technology and shimming techniques have significantly reduced this form of inhomogeneity.

[0042] Magnetic field inhomogeneities due to magnetic susceptibility differences in the imaged object, however, present a complex source of geometric distortion error in MR images. Magnetic susceptibility errors are also referred to as object-induced magnetic field inhomogeneities. The materials in a scanned neurosurgical patient's head have different susceptibility to applied magnetic fields and structurally have different proportions and distribution of constituents with differing magnetic susceptibility, especially fat, water, and air/tissue interfaces. The shape of the head and the proportions of these constituents remodel with age and health status and vary from individual to individual. As a result, each case has a different susceptibility error and therefore different geometric distortion due to inhomogeneity in the scanned object's susceptibility to magnetism. This error is likely difficult to correct using only prescan calibration data. Differential susceptible to the magnetic field of materials possibly in the environment around the imaged object also present a source of geometric distortion error in MR images.

[0043] Overall, magnetic field inhomogeneities are dependent on the material present in the image volume, the shape of the structure being imaged, and its pose within the B₀field. Like other resonance offsets, such as chemical shift for example, these distortions occur only along the read-out gradient direction. Distortion of this nature occurs in gradient recalled echo and spin echo image sequences. The distortion may be significantly large when a small read-out gradient and a large magnetic field are used as in many stereotactic neurosurgical imaging protocols. Also, these distortions may be significant at interfaces of materials with largely differing magnetic susceptibility, such as the air/scalp, scalp/skull, skull/dura, and gray/white matter interfaces in the head for example.

[0044] In accordance with the present invention, a method corrects a three dimensional (3D) magnetic resonance (MR) image in a 3D magnetic resonance imaging (MRI) system.

[0045] For the method, a 3D MRI pulse sequence is performed to acquire a first 3D MR image of an examination region containing at least a portion of a subject. The 3D MRI pulse sequence is performed a second time to acquire a second 3D MR image of the examination region.

[0046] In performing the 3D MRI pulse sequence, the examination region may be spatially encoded with a first phase encoding gradient along a first axis in the examination region, a second phase encoding gradient along a second axis in the examination region, and a frequency encoding gradient along a third axis in the examination region. A Fast Low Angle SHot (FLASH) 3D MR-scan protocol may be performed.

[0047] A voxel error map is computed based on a phase difference between the first and second 3D MR images. Voxel positions are corrected in one of the 3D MR images in accordance with the voxel error map. A forward warping technique may be performed in correcting the one 3D MR image.

[0048] Frameless fiducial markers may be applied to at least a portion of the subject in the examination region for registering the corrected 3D MR image with at least a portion of the subject in the examination region.

[0049] Also in accordance with the present invention, a 3D MRI system comprises suitable means for performing the method. The 3D MRI system may comprise a static magnet, pulse sequence apparatus, pulse sequence control apparatus, and a computer system.

[0050] The static magnet produces a magnetic field along a predetermined axis relative to at least a portion of a subject in an examination region. The pulse sequence apparatus creates magnetic field gradients in the examination region, applies radio frequency (RF) pulsed magnetic fields in the examination region, and receives RF magnetic resonance (MR) signals from the examination region. The pulse sequence control apparatus controls the pulse sequence apparatus to apply a 3D MRI pulse sequence in the examination region a first time and to apply the 3D MRI pulse sequence in the examination region a second time. The pulse sequence control apparatus may control the pulse sequence apparatus to encode spatially the examination region with a first phase encoding gradient along the predetermined axis in the examination region, a second phase encoding gradient along a second axis in the examination region, and a frequency encoding gradient along a third axis in the examination region. The pulse sequence control apparatus may control the pulse sequence apparatus to perform a Fast Low Angle SHot (FLASH) 3D MR-scan protocol.

[0051] The computer system acquires a first 3D MR image from received RF MR signals resulting from the first time application of the 3D MRI pulse sequence in the examination region, acquires a second 3D MR image from received RF MR signals resulting from the second time application of the 3D MRI pulse sequence, computes a voxel error map based on a phase difference between the first and second 3D MR images, and corrects voxel positions in one of the 3D MR images in accordance with the voxel error map. The computer system may perform a forward warping technique in correcting the one 3D MR image.

[0052] The MRI system may also comprise frameless fiducial markers applied to at least a portion of the subject in the examination region for registering the corrected 3D MR image with at least a portion of the subject in the examination region.

[0053] The present invention is illustrated by way of example and not limitation in the figures of the accompanying drawings, in which like references indicate similar elements and in which:

[0054] FIG. 1 illustrates directions for image encoding in magnetic resonance imaging (MRI);

[0055] FIG. 2A illustrates barrel aberration in an MR image;

[0056] FIG. 2B illustrates bow tie effect distortion in an MR image;

[0057] FIG. 2C illustrates potato chip effect distortion in an MR image;

[0058] FIG. 3 illustrates in block diagram form an exemplary three-dimensional (3D) MRI system for performing object-induced geometric distortion correction in accordance with the present invention; and

[0059] FIG. 4 illustrates, for one embodiment, a flow diagram for correcting object-induced geometric distortion in 3D MRI.

[0060] A three-dimensional (3D) magnetic resonance imaging (MRI) system performs object-induced geometric distortion correction.

[0061] FIG. 3 illustrates an exemplary three-dimensional (3D) magnetic resonance imaging (MRI) system 100 for performing object-induced geometric distortion correction in accordance with the present invention.

[0062] The operation of 3D MRI system 100 is controlled by a computer system 110. A console 120 comprising a control panel 122 and a display 124 communicates with computer system 110 to enable an operator to control the production and display of 3D MR images on display 124.

[0063] To produce images with 3D system 100 at least a portion of a subject 102 of interest is placed within an examination region 104. A static magnet 132 produces a substantially uniform, temporally constant magnetic field B₀along a desired z-axis such that selected nuclear magnetic dipoles of subject 102 within examination region 104 preferentially align with the magnetic field. Although illustrated as a human subject, subject 102 may be an animal subject or any other suitable sample.

[0064] Computer system 110 communicates with a pulse program generator 142 to control a set of G_x, G_y, G_zgradient amplifiers and coils 134, a radio frequency (RF) transmitter 152, and an RF receiver 154 so as to carry out a desired 3D MRI scan sequence.

[0065] RF transmitter 152 transmits RF pulses into examination region 104 using RF coils 156 to cause magnetic resonance of the preferentially aligned dipoles of subject 102 within examination region 104. RF receiver 154 receives RF magnetic resonance (MR) signals detected by RF coils 156 from the resonating dipoles of examination region 104. Pulse program generator 142 also controls a transmit/receive (T/R) switch 158 selectively connecting RF transmitter and RF receiver 154 to RF coils 156. Separate transmit and receive RF coils may also be used, obviating any need for T/R switch 158. Computer system 110 comprises an analog-to-digital converter 144 to receive the RF MR signals from RF receiver 154 in digital form and processes the digitized RF MR signals to reconstruct an image representation for display on display 124. Gradient amplifiers and coils 134 impose magnetic field gradients with the static magnetic field along mutually orthogonal x, y, z-axes to spatially encode the received RF MR signals.

[0066] As 3D MRI system 100 scans a 3D image volume in examination region 104 in accordance with a 3D MRI pulse sequence, 3D image data of the entire image volume is simultaneously acquired. For one embodiment, a relatively weak gradient pulse is imposed along the z-axis in the direction of magnetic field B₀. A phase encoding gradient is then activated along the z-axis to encode contiguous parallel image planes p, q, . . . , z in the 3D image volume. An additional phase encoding gradient is activated along the x-axis for a short time during which the dipoles in the 3D image volume acquire a different phase. A frequency encoding gradient is activated along the y-axis to frequency encode the positions of the dipoles. Each recorded data point of the scanned 3D image volume corresponds to a respective voxel in the 3D image to be generated.

[0067] Computer system 110 performs a suitable 3D Fourier Transformation (FT) algorithm to decode frequency information contained in a dipole signal at each recorded data point to corresponding intensity levels. Computer system 110 may then display any desired 2D image in any desired orientation from the reconstructed 3D image data as shades of gray in a matrix arrangement of pixels on display 124.

[0068] Computer system 110 loads software or program code defining different 3D MRI pulse sequences into writable control storage areas of pulse program generator 142. Pulse program generator 142 executes program code corresponding to a given 3D MRI pulse sequence to provide suitable signals that control the operation of RF transmitter 152, RF receiver 154, T/R switch 158, and gradient amplifiers and coils 134 and thereby effectuate the given 3D MRI pulse sequence. Computer system 110 can specify and effectuate any suitable 3D MRI pulse sequence for 3D MRI system 110 as desired.

[0069] Although described in connection with 3D MRI system 100 as illustrated in FIG. 3, any suitable 3D MRI system may be used. One suitable MRI system for 3D scanning is a 1.5 T Siemens Magnetom Vision scanner.

[0070] Ideally, 3D MRI system 100 should maintain an homogeneous magnetic field B₀throughout examination region 104. Because of the magnetic susceptibility differences of subject 102in examination region 104, however, the magnetic field B₀becomes perturbed. Such perturbations lead to off-resonant frequencies at each voxel of the resulting image, thereby distorting the image. In accordance with the present invention, computer system 110 performs an object-induced geometric distortion correction algorithm to correct geometric distortion due to magnetic field inhomogeneities induced by subject 102 in imaging at least a portion of subject 102 with 3D MRI system 100.

[0071] Water and lipid are two major constituents in human and animal subjects, for example, that generate MR signals. An MR spectrum showing a distribution of resonant frequencies of various compounds in an object show that hydrogen atoms attached to water and fat molecules resonate at two distinct frequencies separated by about 205 Hertz (Hz) at B₀=1.5 Tesla (T). At a first approximation, then, the MR spectrum may be considered as two impulses in the frequency domain. The MR spectrum ρ(y,ω) at location y can be characterized as follows:

[0072] ρ(y‘,ρ)=ρ_ω(y)δ(ω)m+ρ_f(y)δ(ω−ω_f  Equation 7

[0073] where ρ_ω(y) and ρ_f(y) are the distributions of water and lipid, respectively, and ω_fis the chemical shift separation of lipid with respect to water.

[0074] Considering the one dimensional case where a line is scanned along the read-out direction y using a Gradient Recalled Echo (GRE) pulse, the shift Ω(y) in resonant frequency at location y is as follows:

[0075] Ω(y)=γB₀χ(y)  Equation 8

[0076] where χ(y) is the magnetic susceptibility of the material at location y. The shift Δ(y) that the voxel at location y undergoes is as follows. $\begin{array}{cc}\Delta \left(y\right)=\frac{\Omega \left(y\right)}{{\mathrm{\gamma G}}_y}& \mathrm{Equation}9\end{array}$

[0077] Extracting this misregistration profile Δ(y) is therefore helpful for correcting any geometric distortion due to object-induced magnetic field inhomogeneities.

[0078] When the pulse sequence is applied, the signal demodulated by water resonant frequency W_ω at echo time T is given by: $\begin{array}{cc}s\left(t,T\right)=\int \int \rho \left(y,\omega \right)e^{j{\mathrm{\gamma G}}_y\mathrm{yt}}{e}^{j\left[\omega +\Omega \left(y\right)\right]\left(t+T\right)}{e}^{-\frac{t+T}{T_2\left(y\right)}}dyd\omega & \mathrm{Equation}10\end{array}$

[0079] where T₂(y) is the spin-spin relaxation time. Substituting Equation 7 in Equation 10, assuming t<<T, yields: $\begin{array}{cc}s\left(t,T\right)=\int p_{\omega }\left(y\right)e^{j\Omega \left(y\right)\left(t+T\right)}{e}^{-\frac{T}{T_2\left(y\right)}}{e}^{j{\mathrm{\gamma G}}_y\mathrm{yt}}dy+\int p_f\left(y\right)e^{j\left[\left({\omega }_w-{\omega }_f\right)\left(t+T\right)\right]}{e}^{-\frac{T}{T_2\left(y\right)}}{e}^{{\mathrm{j\gamma G}}_y\mathrm{yt}}dy& \mathrm{Equation}11\end{array}$

[0080] or: $\begin{array}{cc}s\left(t,T\right)=\int p_{\omega }\left(y\right)e^{\left[\mathrm{j\Omega }\left(y\right)-\frac{1}{T_2\left(y\right)}\right]T}{e}^{{\mathrm{j\gamma G}}_{y}t\left[y+\frac{\Omega \left(y\right)}{{\mathrm{yG}}_y}\right]}dy+e^{{\mathrm{j\omega }}_{f}T}\int p_f\left(y\right)e^{j\left[\Omega \left(y\right)-\frac{1}{T_2\left(y\right)}\right]T}{e}^{{\mathrm{j\gamma G}}_{y}t\left[y+\frac{\left({\omega }_{\omega }-{\omega }_j\right)+\Omega \left(y\right)}{{\mathrm{yG}}_y}\right]}dy& \mathrm{Equation}12\end{array}$

[0081] Where G_y=3.13 mT/m, the chemical shift ω_f/γG_y=1.5 mm Using Equation 8 and an upper bound for χ(y) as 9.05 ppm, the magnetic susceptibility distortion Ω(y)/γG_yis less than 4 mm for B₀=1.5 T. As the voxel location y generally satisfies y>>4 mm, ω_f/γG_y<<y and Ω(y)/γG_y<<y. The chemical shift and magnetic susceptibility distortions are therefore negligibly small compared to the voxel location. Therefore: $\begin{array}{cc}s\left(t,T\right)=\int \left[p_{\omega }\left(y\right)+p_f\left(y+\frac{\left({\omega }_{\omega }-{\omega }_f\right)}{{\mathrm{\gamma G}}_y}\right)e^{j\left({\omega }_{\omega }-{\omega }_f\right)T}\right]e^{j\left[\Omega \left(y\right)-\frac{1}{T_2\left(y\right)}\right]T}{e}^{{\mathrm{j\gamma G}}_y\mathrm{yt}}dy& \mathrm{Equation}13\end{array}$

[0082] The reconstructed signal in the space-domain is then: $\begin{array}{cc}p\left(y,T\right)=\frac{{\mathrm{\gamma G}}_y}{2\pi }\int s\left(t,T\right)e^{-{\mathrm{j\gamma G}}_y\mathrm{yt}}dt& \mathrm{Equation}14\end{array}$

[0083] or: $\begin{array}{cc}p\left(y,T\right)=\frac{{\mathrm{\gamma G}}_y}{2\pi }\left[p_{\omega }\left(y\right)+p_f\left(y\right)e^{-j\left({\omega }_{\omega }-{\omega }_f\right)T}\right]e^{\left[\mathrm{j\Omega }\left(y\right)-\frac{1}{T_2\left(y\right)}\right]T}& \mathrm{Equation}15\end{array}$

[0084] Using the same equations for an echo time (T+ΔT) yields: $\begin{array}{cc}p\left(y,T+\Delta T\right)=\frac{{\mathrm{\gamma G}}_y}{2\pi }\left[{\rho }_{\omega }\left(y\right)+{\rho }_f\left(y\right)e^{{\mathrm{j\omega }}_f\left(T+\Delta T\right)}\right]e^{\left[\mathrm{j\Omega }\left(y\right)-\frac{1}{T_2\left(y\right)}\right]\left(T+\Delta T\right)}& \mathrm{Equation}16\end{array}$

[0085] Choosing ΔT such that ω_fΔT=2π yields: $\begin{array}{cc}p\left(y,T+\Delta T\right)=\frac{{\mathrm{\gamma G}}_y}{2\pi }\left[{\rho }_{\omega }\left(y\right)+{\rho }_f\left(y\right)e^{{\mathrm{j\omega }}_{f}T}\right]e^{\left[\mathrm{j\Omega }\left(y\right)-\frac{1}{T_2\left(y\right)}\right]\left(T+\Delta T\right)}& \mathrm{Equation}17\end{array}$

[0086] This choice of ΔT brings the fat signal for the echo time T=ΔT in phase with the fat signal for the echo time T. Dividing Equation 17 by Equation 15 yields: $\begin{array}{cc}\frac{p\left(y,T+\Delta T\right)}{p\left(y,T\right)}=e^{\mathrm{j\Delta }\mathrm{T\Omega }\left(y\right)}{e}^{-\frac{\Delta T}{T_2\left(y\right)}}& \mathrm{Equation}18\end{array}$

[0087] Considering the amplitude and phase of both sides: $\begin{array}{cc}{T}_2\left(y\right)=\frac{\Delta T}{\ln \Vert \frac{p\left(y,T\right)}{p\left(y,T+\Delta T\right)}\Vert }& \mathrm{Equation}19\\ \Omega \left(y\right)=\frac{\mathrm{phase}\{\frac{p\left(y,T+\Delta T\right)}{p\left(y,T\right)}\}}{\Delta T}& \mathrm{Equation}20\end{array}$

[0088] Hence, from Equation 9: $\begin{array}{cc}\Delta \left(y\right)=\frac{\mathrm{phase}\left\{\frac{p\left(y,T+\Delta T\right)}{p\left(y,T\right)}\right\}}{{\mathrm{\gamma G}}_y\Delta T}& \mathrm{Equation}21\end{array}$

[0089] Equation 21 suggests that a position error profile for the voxels can be obtained by acquiring two scans of the 3D image volume and calculating the phase difference. Knowing this error profile, the voxels in the resulting 3D image data can be shifted to their original locations.

[0090] Computer system 110 for one embodiment performs object-induced geometric distortion correction in accordance with flow diagram 200 of FIG. 4.

[0091] For step 202 of FIG. 4, 3D MRI system 100 acquires a pair of 3D MR images with similar scan parameters yet with different echo recall times T, where ΔT is the shift in echo recall time T between the two 3D image acquisitions. The off-resonant frequencies at each voxel accrue additional phase in the time interval ΔT. 3D MRI system 100 acquires a pair of 3D anatomical images M₁and M₂and a corresponding pair of respective phase volumes φ₁and φ₂at echo recall times T₁and T₂.

[0092] 3D MRI system 100 may perform any suitable imaging protocol to acquire the pair of 3D MR images. One exemplary imaging protocol to acquire the pair of 3D MR images using a 1.5 T Siemens Magnetom Vision scanner is a Fast Low Angle SHot (FLASH) 3D MR-scan protocol to obtain contiguous axial slices of 1 mm thickness with an in-plane resolution of 256×256 pixels, a Field of View (FOV) of 225 mm leading to voxels of 0.89×0.89×1.0 mm, G_y=6.52 mT/m, TR=10 ms, tip angle=20°, echo recall times T of 4 and 8.4 ms for the first and second 3D images, respectively, and thus ΔT=4.4 ms. The FOV may be increased to 256 mm to obtain contiguous axial slices of 1.2 mm thickness, sagittal orientation, for another exemplary imaging protocol.

[0093] In practice the phase difference between the two scans may exceed 2π in which case aliasing occurs leading to phase wrapping. Computer system 110 therefore computes the phase difference at each voxel location between the two 3D images for step 204. If the phase difference exceeds 2π as determined for step 206, computer system 110 for step 208 performs a suitable phase unwrapping technique to determine a suitable phase. As one example, computer system 110 performs a phase unwrapping technique by adding or subtracting multiples of 2π by tracking the trend of the variation of the phase, as suggested in Schneider, E., Glover, G., “Rapid in-vivo shimming”,Magn. Reson. Med.,Vol. 18, pp. 335-347 (1991).

[0094] Because 3D MRI system 100 scans the 3D image volume in accordance with a 3D MRI pulse sequence, the effects of gradient field non-linearities, such as the potato chip and bow tie effect distortions for example, are significantly reduced as a weak slab select gradient is used to excite the entire image volume simultaneously. Such effects may be further reduced by rapidly acquiring relatively high resolution images, that is with small voxel dimensions, using a relatively small field of view (FOV). Computer system 110 for one embodiment seeks to correct only the geometric distortion resulting from magnetic susceptibility differences and therefore presumes any voxel shifts occur only in the frequency encoding direction as a result of magnetic susceptibility differences. Computer system 110 for step 210 computes a voxel position error map in accordance with Equation 21. Computer system 110 for one embodiment computes error only for voxels above a predetermined voxel intensity threshold to eliminate contributions from background noise, such as random noise and static field inhomogeneity for example.

[0095] Computer system 110 then effectuates correction of one of the acquired 3D anatomical images for step 212 by shifting the voxels in that image to their correct locations in accordance with the voxel error map followed by an interpolation and resampling of gray scale values. Computer system 110 for one embodiment computes a voxel error map for correction by forward interpolation through the 3D imaged volume, that is by forward warping, accounting for sub-voxel size shifts in voxel position. Fractional contributions of gray levels from all neighboring voxels of the original 3D image are determined for each voxel position in the corrected 3D image data. Gray scale values across the 3D imaged volume are interpolated and resampled.

[0096] For step 214, computer system 110 may then display on display 124 any desired 2D image in any desired orientation from the corrected 3D image data.

[0097] 3D MRI system 100 may be used for frameless stereotactic neurosurgical procedures. Although significant object-induced geometric distortions are typically encountered at the skin/air and scalp/skull interfaces, where frameless fiducial markers are commonly applied, the voxel-level correction of such object-induced geometric distortions by 3D MRI system 100 allows for better registration of the corrected 3D MR image with at least a portion of the subject in examination region 104.

[0098] Although described as being determined by computer system 110, any suitable computer system, such as a personal computer for example, may be used to perform object-induced geometric distortion correction. As one example, a separate computer system may be connected to computer system 110 and programmed to perform object-induced geometric distortion correction based on image information transmitted from computer system 110.

[0099] In the foregoing description, the invention has been described with reference to specific exemplary embodiments thereof. It will, however, be evident that various modifications and changes may be made thereto without departing from the broader spirit or scope of the present invention as defined in the appended claims. The specification and drawings are, accordingly, to be regarded in an illustrative rather than a restrictive sense.