Correction of Artifacts Caused by Maxwell Terms in Slice Offset Echo Planar Imaging

The field of the invention is nuclear magnetic resonance imaging methods and systems. More particularly, the invention relates to the correction of image artifacts caused by "Maxwell terms" produced by imaging gradients in MRI systems.

When a substance such as human tissue is subjected to a uniform magnetic field (polarizing field B.sub.0), the individual magnetic moments of the spins in the tissue attempt to align with this polarizing field, but precess about it at their characteristic Larmor frequency. If the substance, or tissue, is subjected to a magnetic field (excitation field B.sub.1) which is in the x-y plane and which is near the Larmor frequency, the net aligned moment, M.sub.z, may be rotated, or "tipped", into the x-y plane to produce a net transverse magnetic moment M.sub.t. A signal is emitted by the excited spins, and after the excitation signal B.sub.1 is terminated, this signal may be received and processed to form an image.

When utilizing these signals to produce images, magnetic field gradients (G.sub.x, G.sub.y, and G.sub.z) are employed. Typically, the region to be imaged is scanned by a sequence of measurement cycles in which these gradients vary according to the particular localization method being used. The resulting set of received NMR signals are digitized and processed to reconstruct the image using one of many well known reconstruction techniques.

It is well known that imperfections in the linear magnetic field gradients (G.sub.x, G.sub.y, and G.sub.z) produce artifacts in the reconstructed images. It is a well known problem, for example, that eddy currents produced by gradient pulses will distort the magnetic field and produce image artifacts. Methods for compensating for such eddy current errors are also well known as disclosed, for example, in U.S. Pat. Nos. 4,698,591; 4,950,994; and 5,226,418.

It is also well known that the gradients may not be perfectly uniform over the entire imaging volume, which may lead to image distortion. Methods for compensating this non-uniformity are well known, and for example, are described in U.S. Pat. No. 4,591,789.

Other than uncompensated eddy current errors and gradient non-uniformity errors that escape correction, it can be assumed that the magnetic field gradients (G.sub.x, G.sub.y, and G.sub.z) produce linear magnetic fields exactly as programmed, thus spatially encoding the NMR data accurately. With these gradients, the overall magnetic field at location (x,y,z) is conventionally given as B.sub.0 +G.sub.x x+G.sub.y y+G.sub.z z, and the direction of the field is usually thought to be along the z-axis. This description, however, is not exactly correct. As long as a linear magnetic field gradient is applied, the overall magnetic field is nutated away from the z-axis and its amplitude exhibits higher-order spatial dependencies (x.sup.2, y.sup.2, z.sup.2, x.sup.2 z, . . .). These phenomena are a direct consequence of the Maxwell equations which require that the overall magnetic field satisfy the following two conditions: .gradient..multidot.B=0 and .gradient..times.B=0 (assuming the true and displacement current densities are zero). The higher-order magnetic fields, referred to as "Maxwell terms" (or Maxwell fields), represent a fundamental physics effect, and are not related to eddy currents or imperfection in hardware design and manufacture. Although Maxwell terms have been known for at least a decade, their effect on imaging has been largely ignored because of their negligible effect under conventional imaging conditions.

In axial echo planar imaging (EPI) carried out in a horizontal field magnet, it has been observed that the image in an off-center slice (z.noteq.0) can shift along the phase-encoded direction. The amount of the shift (.DELTA.y) is proportional to the square of the slice location z (.DELTA.y.varies.z.sup.2). This parabolic shift can cause image misregistration problems in neuro-functional imaging, where activation maps obtained from EPI images are overlaid onto a high resolution non-EPI image to establish the correlation between brain function and anatomy. Even in non-functional imaging, the shift can confuse clinical diagnosis and therapeutic planning, especially when the axial 2D slices are reformatted into arbitrary planes. Clearly, the source of this shift needs to be identified, and methods to remove the shift need to be developed.

The present invention is a method for reducing artifacts in off-center EPI images caused by Maxwell terms produced by the readout gradient. More specifically, a phase shift (.phi.) and a frequency shift (.DELTA.f) are calculated for each NMR echo signal in an EPI pulse sequence based on the imaging gradient waveforms, the main field, and the slice location. These calculated values are either used to control the phase and frequency of the rf reference signal applied to the receiver as these NMR echo signals are received, or applied to the acquired signals (k-space data) during image reconstruction to remove the artifacts.

In axial echo planar imaging (EPI) carried out in a horizontal field magnet (e.g. a superconducting magnet), an image in an off-center slice can shift along the phase-encoded direction. The amount of the shift (.DELTA.y) is proportional to the square of the slice location z (.DELTA.y.varies.z.sup.2). This phenomenon, which we call parabolic shift, is caused by the z.sup.2 -Maxwell term arising mainly from the EPI readout gradient. The parabolic shift, although most frequently encountered in axial images in a horizontal field magnet, can also be observed in off-center sagittal and coronal slices. For certain magnets whose magnetic field corresponds to the anterior/posterior direction of the patient, the parabolic shift is most visible in the coronal planes. In addition to the parabolic shift, Maxwell terms can also produce a Nyquist ghost whose intensity is proportional to the difference of two images shifted in opposite directions.

We have developed two techniques which remove the parabolic shift as well as the Nyquist ghost caused by the Maxwell term. In the first technique, frequency and phase errors caused by the Maxwell terms are calculated analytically on an individual-slice and an individual-echo basis. The calculated errors are then compensated by adjusting the receiver frequency and phase during data acquisition. In the second technique, the frequency and phase errors are first converted to two linear phase shifts in the k-space data, one in the readout direction and the other in the phase-encoding direction. The linear phase shifts are removed in the image reconstruction process after data acquisition.

FIG. 1 is a block diagram of a known MRI system which is modified to employ the present invention;

FIG. 2 is an electrical block diagram of a known transceiver which forms part of the MRI system of FIG. 1;

FIG. 3 is a graphic representation of readout gradient used in an EPI pulse sequence;

FIG. 4 is a graphic representation of an EPI pulse sequence to which the present invention applies; and

FIG. 5 is a graphic representation of the changes made to the pulse sequence of FIG. 4 to practice the first preferred method of the present invention.

The Maxwell terms are basically the higher-order spatial gradients (second order, third order, etc.) produced by the linear magnetic field gradients (x, y, and z-gradients). These terms can be directly derived from the Maxwell equations. According to the Maxwell equations, a magnetic field B must satisfy the following two conditions: ##EQU1## where .gradient. is the derivative operator (.gradient..ident.i.differential./.differential.x+j.differential./.differe ntial.y+k.differential./.differential.z), E is the electric field, J is the current density, and .mu..sub.0 and .epsilon..sub.0 are the magnetic permeability and the electric permeability of the free space, respectively. If there is no displacement current and the electric field is static, equation 1b reduces to:

From equations 1a and 1c, we obtain: ##EQU2## The above 4 equations 2 and 3a-c contain a total of 9 partial derivatives, among which only 5 are independent. Our next task is to select these five independent variables. Recognizing that ##EQU3## (G.sub.x, G.sub.y, and G.sub.z are the linear gradients), we can readily choose G.sub.x, G.sub.y, and G.sub.z as the first three independent variables. For a radially symmetric G.sub.z -field in cylindrical coordinates .differential.B.sub.x /.differential.x and .differential.B.sub.y /.differential.y should be identical. However, to cover a more general case, we choose a dimensionless symmetry parameter .alpha. as the fourth independent variable: ##EQU4## The last independent variable can be conveniently chosen as (based on equation 3a): ##EQU5## At this point, all the partial derivatives described in equations 2 and 3 can be expressed using the 5 independent variables G.sub.x, G.sub.y, G.sub.z, .alpha., and g: ##EQU6## With all the terms, the overall magnetic field becomes:

where, to first order, ##EQU7## The above equations have two important implications. First, the B.sub.0 -field is no longer aligned along the z-axis because of the transverse fields B.sub.x and B.sub.y. Second, the amplitude of the B.sub.0 field is not simply given by B=B.sub.0 +G.sub.x x+G.sub.y y+G.sub.z z, but instead by ##EQU8## (B.sub.0 +G.sub.x x+G.sub.y y+G.sub.z z merely represents the z-component of the overall field). If we perform three sequential Taylor series expansions on equation 9 with respect to x, y, and z, respectively, it can be seen that the magnetic field not only has its regular zeroth and first order spatial dependencies, but also shows higher-order spatial components. The result of the Taylor expansion to the second order is given by equation 10: ##EQU9## For most gradient systems used in MRI, we have g=0, and .alpha..apprxeq.1/2 (due to the cylindrical symmetry). Under these conditions, equation 10 is simplified to: ##EQU10## Equations 10 and 11 indicate that, whenever a linear magnetic field gradient is applied, higher-order gradient fields will be produced to satisfy the Maxwell equations. These higher-order gradient fields are referred to as the "Maxwell terms", or "Maxwell fields."With the inclusion of the Maxwell terms, the 2D NMR signal equation becomes: ##EQU11## where B.sub.M is the Maxwell magnetic field and .phi..sub.M is the associated phase error. As implied by equations 12a-c, the Maxwell phase error depends on the details of each pulse sequence. In some pulse sequences, the phase error can be negligible or zero and thus causes no image degradation. In most other sequences, a non-negligible phase error is produced, giving various image quality problems such as distortion, ghosting, image shift, shading, blurring, and intensity reduction.

In axial EPI carried out in a horizontal field magnet, the readout gradient G.sub.ro can be applied along either x or y-axis. For a general case, G.sub.ro can be expressed as: ##EQU12## where i and j are unit vectors in the G.sub.x and G.sub.y directions, respectively. Note that equation 13a also accounts for oblique scans, although we rarely use oblique scans in the axial plane. According to equation (11), the z.sup.2 -Maxwell term produced by the G.sub.x and/or G.sub.y gradients is given by: ##EQU13## where B.sub.0 is the main magnetic field. Combining equations 13b and 14, we obtain: ##EQU14## The Maxwell field B.sub.M vanishes at the central slice (z=0), provided the slice is infinitesimally thin. For off-center slices, however, B.sub.M increases parabolically with respect to the slice location z. Such a quadratic magnetic field has two consequences.

First, the Larmor resonance frequency of the spins is shifted away from the pre-set receiver demodulation frequency (Eq. 16), leading to a spatial shift along the readout direction (Eq. 17): ##EQU15## Since the readout gradient is bipolar, the spatial shift .delta..sub.ro alternates between the odd and even echoes, producing a Nyquist ghost which is the difference between two shifted images in the opposite directions (i.e., .+-..delta..sub.ro) along the readout direction. This effect is due to the use of the bipolar readout gradient. If we use G.sub.ro =2.2G/cm, z=24 cm and B.sub.0 =1.5T, the maximum ghost-to-image intensity ratio is found to be 0.5%.

Second, the z.sup.2 -Maxwell term introduces a cumulative phase shift at each echo acquisition throughout the echo train. The phase at the center of echo m is given as: ##EQU16## where t.sub.m is the time at the center of the mth echo, and t.sub.s is the time when the readout gradient (including pre-phasing gradient) is first activated. For a periodic readout gradient G.sub.ro, equation 18 can be expressed as: ##EQU17## t.sub.esp is the inter-echo spacing, and .phi..sub.p (z) is the Maxwell phase introduced by the pre-phasing readout gradient. For a trapezoidal readout gradient train illustrated in FIG. 3, .lambda.(z) can be explicitly expressed as: ##EQU18## where g.sub.ro, t.sub.1 and t.sub.2 are defined in FIG. 3. In single-shot EPI, consecutive m values correspond to consecutive k-space lines along the phase-encoding direction. Thus, the phase in equation 19 represents a linear phase in the k-space, which gives a shift .delta..sub.pe in the image domain according to the shift theorem of Fourier transform: ##EQU19## where L.sub.pe is the field-of-view (FOV) in the phase-encoding direction. If we choose g.sub.ro =2.2G/cm, z=24 cm, B.sub.0 =1.5T, t.sub.1 =184 .mu.s, and t.sub.2 =512 .mu.s (corresponding to a receiver bandwidth of .+-.125 kHz for 128 readout data points), we find that .delta..sub.pe =0.251L.sub.pe, a shift of more than a quarter of the FOV! For multi-shot EPI with N shots, equation 20 becomes: ##EQU20## Even in this case, the shift along the phase-encoding direction can be substantial.

To summarize, the z.sup.2 -Maxwell term produced by the readout gradient in axial EPI scans can cause a Nyquist ghost and an image shift along the phase-encoding direction. Although the level of the ghost could be negligible, the image shift is substantial. Since the amount of the shift is quadratically proportional to the slice location and the readout gradient, and inversely proportional to the static magnetic field, increasing the slice offset or the readout gradient, or decreasing the main magnetic field can aggravate the parabolic shift problem.

To remove the parabolic shift, as well as the Nyquist ghost, we have developed two techniques, a compensation technique during data acquisition and a post-acquisition correction technique during image reconstruction.

In the first technique the frequency and phase errors are first calculated from equations 16-19 for each echo m and for each slice located at z. These errors are compensated by dynamically adjusting the receiver frequency and phase on an individual-echo, individual-slice basis while the data acquisition proceeds. The frequency and phase compensation values are: .DELTA.f.sub.comp (z)=-.DELTA.f(z), and .phi..sub.comp (m,z)=-.phi.(m,z), respectively, where .DELTA.f(z) and .phi.(m,z) are respectively given by equations 16 and 19.

In the second technique, the echo signals are first acquired in the presence of the frequency and phase errors. After data acquisition, a linear phase shift .psi..sub.ro (z) given by equation 22 is applied along the readout direction to correct for the frequency error (equation 16), and another linear phase shift .psi..sub.pe (z)=-.lambda.(z) is applied along the phase-encoding direction to remove the phase error (equation 19). ##EQU21## (L.sub.ro is the FOV in the readout direction). For multi-shot EPI, equation 22 remains the same, but the linear phase shift along the phase-encoding direction must be done with respect to a group of k-space lines, instead of individual ones. It is interesting to note that the polarity of the linear phase shift along the readout direction alternates with the readout gradient polarity, as indicated by equation 22.

The above methods to remove the parabolic shift and the Nyquist ghost assume that the z.sup.2 -Maxwell term mainly arises from the EPI readout gradient. This is generally a good approximation since the phase-encoding gradient blip is significantly smaller than its readout counterpart and is only active for a brief period of time. In cases where the Maxwell term due to the phase-encoding gradient becomes significant, such as multi-shot EPI with a large number of shots, equations 15 and 19a should be respectively modified as follows: ##EQU22## where G.sub.pe is the phase-encoding gradient waveform. Although we did not consider the phase-encoding gradient in the above discussion, inclusion of this gradient is straightforward and only involves minor changes.

Although the Maxwell term induced parabolic shifts and ghosts are most evident in axial slices for a horizontal field magnet, similar problems can also be observed in sagittal and coronal images. For example, in a coronal image (xz-plane) acquired using G.sub.z as the readout gradient, the corresponding Maxwell term becomes: ##EQU23## The first term in equation 25 causes in-plane image distortion (R. M. Weiskoff, et al., Magnetic Resonance in Medicine, Vol. 29, 796-803, 1993), and the second term generates a quadratic phase along the slice-selection direction (y-direction), resulting in parabolic shifts and ghosts similar to those discussed for axial images. All the theoretical analyses and correction methods presented for axial images are equally applicable to coronal and sagittal images, except that the coefficient of the Maxwell term is four times smaller. It is interesting to note that when the readout gradient is chosen along x-axis in a coronal image, no quadratic Maxwell field along the y-direction is produced. Therefore, the image does not experience parabolic shift and ghosting. The same argument is also true for sagittal images with the readout direction along the y-axis.

In the analysis of the Maxwell terms, we have assumed a horizontal field magnet configuration. When a vertical field magnet is used, the physical z-axis is changed to the anterior-posterior direction of the patient. Therefore, the coronal images, instead of axial ones, will exhibit the most prominent parabolic shift and ghosting. The same principles to correct for the artifacts presented for axial slices in horizontal field magnet are equally applicable to this case with only minor notational changes.

Although we have employed an EPI pulse sequence with a trapezoidal readout gradient as an example in the above discussions, it should be understood that the present invention can be practiced in other EPI pulse sequences, and other sequences employing bipolar readout gradients for multiple echo acquisitions. For example, the trapezoidal readout gradient in FIGS. 3 and 4 can be replaced by a sinusoidal gradient. Even for trapezoidal gradients, data acquisition is not limited to the flat-top and can be acquired while the gradient ramps up or down. In these cases, equations (16) and (19a) can still be used to calculate the frequency and phase errors with the understanding that .DELTA.f(z) in equation (16) now becomes time-dependent: .DELTA.f(z,t). Another variation of the EPI pulse sequence is the so called "skip-echo" EPI where only odd or even echoes are acquired when the readout gradient lobes have the same polarity. All the above equations derived for correcting the Maxwell-field induced artifacts are still valid for skip-echo EPI since we have defined t.sub.esp as the inter-echo spacing (i.e., the time interval between two consecutively acquired echoes), instead of the time interval between the positive and negative gradient lobes. Echo planar imaging may also be combined with the fast spin echo technique to acquire multiple gradient echoes within multiple spin echoes. The correction methods, presented in equations (15)-(24), are equally applicable to this imaging technique. Description Of The Preferred Embodiment

Referring first to FIG. 1, there is shown the major components of a preferred MRI system which incorporates the present invention. The operation of the system is controlled from an operator console 100 which includes a keyboard and control panel 102 and a display 104. The console 100 communicates through a link 116 with a separate computer system 107 that enables an operator to control the production and display of images on the screen 104. The computer system 107 includes a number of modules which communicate with each other through a backplane. These include an image processor module 106, a CPU module 108 and a memory module 113, known in the art as a frame buffer for storing image data. The computer system 107 is linked to a disk storage 111 and a tape drive 112 for storage of image data and programs, and it communicates with a separate system control 122 through a high speed serial link 115.

The system control 122 includes a set of modules connected together by a backplane. These include a CPU module 119 and a pulse generator module 121 which connects to the operator console 100 through a serial link 125. It is through this link 125 that the system control 122 receives commands from the operator which indicate the scan sequence that is to be performed. The pulse generator module 121 operates the system components to carry out the desired scan sequence. It produces data which indicates the timing, amplitude and shape of the RF pulses which are to be produced, and the timing of and length of the data acquisition window. The pulse generator module 121 connects to a set of gradient amplifiers 127, to indicate the timing and shape of the gradient pulses to be produced during the scan. The pulse generator module 121 also receives patient data from a physiological acquisition controller 129 that receives signals from a number of different sensors connected to the patient, such as ECG signals from electrodes or respiratory signals from a bellows. And finally, the pulse generator module 121 connects to a scan room interface circuit 133 which receives signals from various sensors associated with the condition of the patient and the magnet system. It is also through the scan room interface circuit 133 that a patient positioning system 134 receives commands to move the patient to the desired position for the scan.

The gradient waveforms produced by the pulse generator module 121 are applied to a gradient amplifier system 127 comprised of G.sub.x, G.sub.y and G.sub.z amplifiers. Each gradient amplifier excites a corresponding gradient coil in an assembly generally designated 139 to produce the magnetic field gradients used for spatially encoding acquired signals. The gradient coil assembly 139 forms part of a magnet assembly 141 which includes a polarizing magnet 140 and a whole-body RF coil 152. A transceiver module 150 in the system control 122 produces pulses which are amplified by an RF amplifier 151 and coupled to the RF coil 152 by a transmit/receive switch 154. The resulting signals radiated by the excited nuclei in the patient may be sensed by the same RF coil 152 and coupled through the transmit/receive switch 154 to a preamplifier 153. The amplified NMR signals are demodulated, filtered, and digitized in the receiver section of the transceiver 150. The transmit/receive switch 154 is controlled by a signal from the pulse generator module 121 to electrically connect the RF amplifier 151 to the coil 152 during the transmit mode and to connect the preamplifier 153 during the receive mode. The transmit/receive switch 154 also enables a separate RF coil (for example, a head coil or surface coil) to be used in either the transmit or receive mode.

The NMR signals picked up by the RF coil 152 are digitized by the transceiver module 150 and transferred to a memory module 160 in the system control 122. When the scan is completed and an entire array of data has been acquired in the memory module 160, an array processor 161 operates to Fourier transform the data into an image data set. This image data set is conveyed through the serial link 115 to the computer system 107 where it is stored in the disk memory 111. In response to commands received from the operator console 100, this image data set may be archived on the tape drive 112, or it may be further processed by the image processor 106 and conveyed to the operator console 100 and presented on the display 104.

Referring particularly to FIGS. 1 and 2, the transceiver 150 produces the RF excitation field B1 through power amplifier 151 at a coil 152A and receives the resulting signal induced in a coil 152B. As indicated above, the coils 152A and B may be separate as shown in FIG. 2, or they may be a single wholebody coil as shown in FIG. 1. The base, or carrier, frequency of the RF excitation field is produced under control of a frequency synthesizer 200 which receives a set of digital signals from the CPU module 119 and pulse generator module 121. These digital signals indicate the frequency and phase of the RF carrier signal produced at an output 201. The commanded RF carrier is applied to a modulator and up converter 202 where its amplitude is modulated in response to a signal R(t) also received from the pulse generator module 121. The signal R(t) defines the envelope of the RF excitation pulse to be produced and is produced in the module 121 by sequentially reading out a series of stored digital values. These stored digital values may, in turn, be changed from the operator console 100 to enable any desired RF pulse envelope to be produced.

The magnitude of the RF excitation pulse produced at output 205 is attenuated by an exciter attenuator circuit 206 which receives a digital command, from the backplane 118. The attenuated RF excitation pulses are applied to the power amplifier 151 that drives the RF coil 152A. For a more detailed description of this portion of the transceiver 122, reference is made to U.S. Pat. No. 4,952,877.

Referring still to FIG. 1 and 2 the signal produced by the subject is picked up by the receiver coil 152B and applied through the preamplifier 153 to the input of another receiver amplifier 207. The receiver amplifier 207 further amplifies the signal by an amount determined by a digital attenuation signal received from the backplane 118.

The received signal is at or around the Larmor frequency, and this high frequency signal is down converted in a two step process by a down converter 208 which first mixes the NMR signal with the carrier signal on line 201 and then mixes the resulting difference signal with the 2.5 MHz reference signal on line 204. The down converted NMR signal is applied to the input of an analog-to-digital (A/D) converter 209 which samples and digitizes the analog signal and applies it to a digital detector and signal processor 210 which produces 16-bit in-phase (I) values and 16-bit quadrature (Q) values corresponding to the received signal. The resulting stream of digitized I and Q values of the received signal are output through backplane 118 to the memory module 160 where they are employed to reconstruct an image.

The EPI pulse sequence employed in the preferred embodiment of the invention is illustrated in FIG. 4. A 90.degree. RF excitation pulse 250 is applied in the presence of a slice-selection gradient pulse 251 to produce transverse magnetization in a slice. The excited spins are rephased by a negative lobe 252 on the slice-selection gradient and then a time interval elapses before a 180.degree. RF refocusing pulse 260 is applied in the presence of a slice-selection gradient pulse 262, surrounded by crusher gradients. A total of N.sub.y (e.g., N.sub.y =128) separate NMR echo signals, indicated at 253 are acquired during the EPI pulse sequence.

The NMR echo signals 253 are gradient recalled echoes produced by the application of an oscillating readout gradient 255. The readout sequence is started with a pre-phasing readout gradient lobe 256 and the echo signals 253 are produced as the readout gradient oscillates between positive and negative values. A total of N.sub.x (e.g., N.sub.x =128) samples are taken of each NMR echo signal 253 during each readout gradient pulse 255. The successive N.sub.y NMR echo signals 253 are separately phase encoded by a series of phase encoding gradient pulses 258. A pre-phasing phase-encoding lobe 259 occurs before the echo signals acquired to position the central view (k.sub.y =0) at the desired echo time (TE). Subsequent phase encoding pulses 258 occur as the readout gradient pulses 255 switch polarity, and they step the phase encoding monotonically upward through k.sub.y -space.

At the completion of the EPI pulse sequence, therefore, N.sub.x separate frequency encoded samples of N.sub.y separately phase encoded NMR echo signals 253 have been acquired. After time reversing every other echo, this N.sub.x .times.N.sub.y element array of complex numbers is Fourier transformed along both of its dimensions (k.sub.y and k.sub.x) to produce an image data set that indicates the NMR signal magnitude along each of its two dimensions (x and y).

Referring still to FIG. 4, during a conventional EPI pulse sequence the receiver reference signal used to demodulate the NMR signals 253 is set to a fixed frequency as indicated at 264 and its phase accumulates at a fixed rate as indicated at 266. As shown in FIG. 2, the 2.5 MHz reference signal as well as the 250 kHz sampling signal and the 5, 10 and 60 MHz reference signals are produced by a reference frequency generator 203 from a common 20 MHz master clock signal.

One way the present invention is implemented is by controlling the frequency and phase of the 2.5 MHz reference signal on line 201 to correct for the frequency shift .DELTA.f(z) given in equation (16) and the cumulative phase shift given in equation (19). Referring to FIGS. 4 and 5, after the RF pulses 250 and 260 are produced during the EPI pulse sequence, the frequency and phase of the 2.5 MHz reference signal on line 201 is controlled by command signals from the pulse generator module 121 to make the corrective changes as the readout gradient 255 is played out. The amount of this correction depends on the location of the slice along the z-axis, and hence the amount of the corrections .DELTA.f(z) and .phi.(m,z) will be different for each slice (1 through n) in a multi-slice acquisition. This is indicated in FIG. 5 by the different levels of the receiver frequency at 268 and the different rates of phase accumulation at 270.

The other way the present invention can be implemented is by first calculating the two linear phase shifts .psi..sub.ro (z) and .psi..sub.pe (z) for individual slices, as discussed earlier. The linear phase shifts are then respectively applied to the readout and the phase-encoding directions of the complex k-space data. After the phase correction, a 2D Fourier transform on the k-space data produces an image where the parabolic shifts and ghosts caused by the Maxwell terms are removed.