Correction of magnetic resonance images

[0001] This invention relates to magnetic resonance imaging (MRI). In particular, the invention is directed to the use of spherical harmonic deconvolution methods to compensate for the effects on MRI images of (i) non-linear gradients (hereinafter referred to as the “non-linear problem”), (ii) misalignments of the symmetry points of one or more of the MR system's field generating devices, i.e., the main (B₀) magnet, the X-gradient magnet, the Y-gradient magnet, and the Z-gradient magnet (hereinafter referred to as the “non-coincident center problem”), (iii) deviations from ideal geometric alignment of the field generating devices (hereinafter referred to as the “axes misalignment problem”), and (iv) winding and/or design errors of the field generating devices (hereinafter referred to as the “winding error problem”).

[0002] The quality of images in magnetic resonance imaging is dependent upon the homogeneity of the respective instrumental components used to generate the image. Much attention has been given to the design and development of magnets having the necessary homogeneity over the imaging volume. Variations in magnetic field homogeneity manifest themselves as image distortions and intensity variation from the true values. Radio frequency (rf) coil inhomogeneity results in non-uniform excitation of the nuclear spins resulting in intensity variations, while deviation from linearity of the gradient field from application of the pulsed field gradients induces image translation and distortions.

[0003] More specifically, the quality of a magnetic resonance image is dependent upon the accuracy by which physical position is spatially encoded. As MRI data is now used routinely for stereotaxy, longitudinal studies of atrophy and functional studies, ensuring images have no distortion and inhomogeneity is critical. The principal machine dependant sources of this inhomogeneity previously recognized in the art are eddy currents, gradient non-linearity, B₀ and B₁ inhomogeneity. See Tanner S F, Finnigan D J, Khoo V S, Mayles P, Dearnaley D P, Leach M O. Radiotherapy planning of the pelvis using distortion corrected MR images: the removal of system distortions. Phys Med Biol. 2000 August; 45(8):2117-32.

[0004] As discussed in detail below, the present invention in accordance with certain of its aspects is directed to an analytical approach for calculating and removing the effects of non-linear gradients. Non-linear gradients have been addressed because of the recent interest in short-bore high-speed gradients.

[0005] In particular, there has been recent significant interest in the development of gradients having rise times less than 200 μseconds. Although peripheral nerve stimulation is a limiting feature of short rise times, such gradients have found use in high-speed echo planar imaging (EPI) of the heart and diffusion tensor imaging of the brain. To achieve short rise times and to avoid peripheral nerve stimulation, gradient designers have restricted the length of the gradients and limited the number of turns. These limitations, although suitable for the implementation of pulse sequences having the desired speed, have the undesired consequence of increased non-linearity in the gradient fields. The result is distorted images.

[0006] More particularly, non-linear pulsed field gradients induce image distortions due to incorrect spatial encoding of the signal. If we assume field gradients are linear, it follows that k-space is sampled linearly and thus a fast Fourier transform (FFT) is suitable for reconstruction. However, any deviation from linearity in the gradients results in non-linear data sampling and subsequent errors in image spatial encoding. A non-linear Fourier transformation would allow this data to be correctly transformed to an image. Unfortunately, non-linear Fourier transformation greatly increases computation time by N/log₂N, relative to a FFT, making real-time image generation computationally prohibitive.

[0007] As discussed fully below, the present invention in accordance with certain of its aspects addresses this problem and provides a general analytical solution to correct image distortions induced by gradient non-linearity which: (1) is applicable to any gradient configuration, (2) is robust, and (3) importantly, maintains the FFT as the method for image reconstruction.

[0008] Mapping and correcting of MRI distortions has been previously considered in the art and attempts have been made to correct the image distortion induced by the non-linearity of the gradients produced by the gradient coils. Two early papers by Schad et al discuss “pincushion” effects seen with 2D phantoms. See Schad L, Lott S, Schmitt F, Sturm V, Lorenz W J. Correction of spatial distortion in MR imaging: a prerequisite for accurate stereotaxy. J Comput Assist Tomogr. 1987 May-June; 11(3):499-505; and Schad L R, Ehricke H H, Wowra B, Layer G, Engenhart R, Kauczor H U, Zabel H J, Brix G, Lorenz W J. Correction of spatial distortion in magnetic resonance angiography for radiosurgical treatment planning of cerebral arteriovenous malformations. Magn Reson Imaging. 1992; 10(4):609-21.

[0009] Subsequent schemes include:

• • [0010] (1) using a stereotaxic frame as a reference marker (see Maurer C R Jr, Aboutanos G B, Dawant B M, Gadamsetty S, Margolin R A, Maciunas R J, Fitzpatrick J M. Effect of geometrical distortion correction in MR on image registration accuracy. J Comput Assist Tomogr. 1996 July-August; 20(4):666-79; and Moerland M A, Beersma R, Bhagwandien R, Wijrdeman H K, Bakker C J. Analysis and correction of geometric distortions in 1.5 T magnetic resonance images for use in radiotherapy treatment planning. Phys Med Biol. 1995 October; 40(10):1651-4);
  • [0011] (2) comparing phantom images from CT and MRI (see Fransson A, Andreo P, Potter R. Aspects of MR image distortions in radiotherapy treatment planning. Strahlenther Onkol. 2001 February; 177(2):59-73; and Yu C, Petrovich Z, Apuzzo M L, Luxton G. An image fusion study of the geometric accuracy of magnetic resonance imaging with the Leksell stereotactic localization system. J Appl Clin Med Phys. 2001 Winter; 2(1):42-50); and
  • [0012] (3) imaging specific MRI phantoms that typically consist of an array of tubes filled with a suitable contrast agent (see Bakker C J, Moerland M A, Bhagwandien R, Beersma R. Analysis of machine-dependent and object-induced geometric distortion in 2DFT MR imaging. Magn Reson Imaging. 1992; 10(4):597-608; Mizowaki T, Nagata Y, Okajima K, Murata R, Yamamoto M, Kokubo M, Hiraoka M, Abe M. Development of an MR simulator: experimental verification of geometric distortion and clinical application. Radiology. 1996 June; 199(3):855-60; Woo J H, Kim Y S, Kim S I. The correction of MR images distortion with phantom studies. Stud Health Technol Inform. 1999; 62:388-9.; and Tanner et al., supra).

[0013] A later Mizowaki et al paper (Mizowaki T, Nagata Y, Okajima K, Kokubo M, Negoro Y, Araki N, Hiraoka M. Reproducibility of geometric distortion in magnetic resonance imaging based on phantom studies. Radiother Oncol. 2000 November; 57(2):237-42) alludes to an important problem with phantom studies; after assessing the reproducibility of their own phantom based correction, the authors conclude that such correction may be only applicable in limited situations.

[0014] Field mapping is another approach to determining distortion via induced phase shift (see Chang H, Fitzpatrick J M. A technique for accurate magnetic resonance imaging in the presence of field inhomogeneities. IEEE Trans. Med. Imaging, 11: 319-329, 1992; Chen N K, Wyrwicz A M. Optimized distortion correction technique for echo planar imaging. Magn Reson Med. 2001 March; 45(3):525-8; and Cusack R, Brett M, Osswald K. An evaluation of the use of magnetic field maps to undistort echo-planar images. Neuroimage. 2003 January; 18(1):127-42). In particular, the Chang et al. approach relies on the acquisition of two images using different gradient strengths. From the induced phase shift across the difference image a distortion map of the gradient field can be calculated and the images corrected. Other correction schemes have been proposed for particular problems including a post-processing method for correcting distortions induced when using bi-planar gradients (see Liu H, An efficient geometric image distortion correction method for a biplanar planar gradient coil. Magnetic Resonance Materials in Physics, Biology and Medicine, 10: 75-79, 2000). Regardless, these approaches are of limited use in attempting to discern, model or correct distortion due to part of the MRI system in isolation, specifically, the gradients.

[0015] In a recent paper (see Langlois S, Desvignes M, Constans J M, Revenu M. MRI geometric distortion: a simple approach to correcting the effects of non-linear gradient fields. J Magn Reson Imaging. 1999 June; 9(6):821-31), a method has been suggested for the correction of both the intensity variations and geometric distortions induced by non-linear gradients based on treating gradient coils as either an opposed Helmoltz pair for the Z-gradient and a Golay arrangement for the transverse gradients. However, this approach treats the non-linear component of the gradient field as a constant for the intensity correction and demonstrates correction of the non-linear induced distortion by correcting a large phantom (a cube of approximately 20 cm length). It is also assumed that for the chosen gradient geometry only second order terms are important. Whilst this approach is useful for gradients with a significant length to diameter ratio (>2), it is of limited use with current high-speed gradients with significant higher order gradient field impurities.

[0016] In a recent patent (see Krieg R, Schreck O. U.S. Pat. No. 6,501,273, issued Dec. 31, 2002), it has been suggested that spherical harmonics may be useful in developing a general robust method, but none is given. Another patent application (see Zhu Y, U.S. Patent Application Publication No. US 2002/0093334, dated Jul. 18, 2002) has suggested a method to reconstruct k-space based upon an assumed non-linearity. The method proposed involves approximately 10¹⁴ triple integrations and thus is of limited practical value. Indeed, this patent publication highlights one of the primary deficiencies in the art—namely, the failure to realize that any method, which is of practical value, must use the FFT as its basis.

[0017] From the foregoing, it can be seen that there exists an established need in the art for a general method for correcting image distortions due to gradient non-linearity. In particular, a need exists for a method that is based on a rigorous approach to the problem allowing clinicians confidence in their image interpretation.

[0018] This invention in accordance with certain of its aspects provides such a general method to correct image distortions induced by gradient non-linearity. Significantly, the method is applicable to any gradient configuration.

[0019] In addition to the problem of non-linearity, in the course of developing the analytical techniques described below, it was discovered that the field generating devices of a typical MRI system can have (i) substantially displaced symmetry points, i.e., the symmetry points do not share a common center (the non-coincident center problem), and (ii) substantial deviations from ideal geometric alignment, i.e., the individual Cartesian coordinate systems associated with the individual field generating devices, in addition to having displaced centers, can also be rotated relative to each other and/or relative to the overall Cartesian coordinate system of the MRI system (the axes misalignment problem). Such translational and rotational misalignments in and of themselves produce distortions in MR images. Winding and/or design errors for the field generating devices are another source of distortion (the winding error problem). In accordance with others of its aspects, the present invention also addresses the problems associated with these sources of distortion.

[0020] This invention provides methods for correcting distortions of MR images due to gradient non-linearities in a nuclear magnetic resonance (NMR) tomography apparatus and/or distortions due to translational and/or rotational misalignments of the field generating devices of such apparatus and/or distortions due to errors in the winding and/or design of the field generating devices. The apparatus typically comprises a magnet (with shim coils as required) for generating a substantially constant main (B₀) magnetic field in a volume of interest, as well as gradient coils for generating gradient magnetic fields, as is known in the art. The method comprises the initial step of deconvolving the main magnetic field. Preferably, the magnetic field is deconvolved using a bases set of spherical harmonics, but other suitable bases sets can be used. The deconvolution is performed using an origin defined by the MRI apparatus, e.g., the center of the image region.

[0021] The gradient coils are then preferably excited individually, i.e., one at a time, using a constant current. Typically, a low current is applied in order not to burn out the coils. For each constantly excited gradient coil, the magnetic field is deconvolved again preferably using the same origin. As the gradient field is intended to be linear, for a perfect gradient magnet having its symmetry point and axes coincident with the origin and its associated axes, a single harmonic coefficient should result from the deconvolution, e.g., the a(1,1) coefficient for the X-gradient magnet, the b(1,1) coefficient for the Y-gradient magnet, and the a(1,0) coefficient for the Z-gradient magnet (hereinafter referred to as the “primary coefficients”). The harmonics of the main magnet (optionally shimmed as necessary) obtained from the first deconvolution are then subtracted, in order to obtain the effect of just the gradient and its impurities and/or alignment errors. The result is a set of spherical harmonic coefficients for the gradient coil, i.e., the “a” and “b” coefficients for the coil discussed below and referred to in the claims. The presence of coefficients with substantial values other than the primary coefficients is indicative of non-linearities in the gradient magnets and/or non-coincident symmetry points, and/or rotated axes, and/or winding or design errors.

[0022] Although it is preferred to obtain spherical harmonic coefficients for each gradient coil separately, in some cases, the symmetry of the system will permit the use of spherical harmonic coefficients obtained for one of the gradient coils for another coil. For example, for B₀ directed along the z-axis of an x,y,z-coordinate system and for x-gradient, y-gradient, and z-gradient coils having symmetry points at the same or effectively the same location as (i) the main (B₀) coil (after any shimming), (ii) one another, and (iii) the center of the image region, and for well-aligned axes for the field generating devices and low levels of winding/design errors, measurements may only be needed for the z-gradient coil and one of the x and y coils, with the spherical harmonic coefficients for the measured x or y coil being used for the non-measured coil. Similarly, it is most preferred to obtain spherical harmonic coefficients for each individual MRI machine as installed and shimmed. However, generic coefficients obtained for a representative sample of a particular series or model of MRI machine may be used if desired. Along these same lines, theoretical coefficients, as opposed to measured coefficients, can be used, but are not preferred.

[0023] As a general rule, the more measurements made, e.g., measurements for each of the x, y, and z gradient coils, and the more the measurements are for a specific machine, e.g., for an individual machine after installation and shimming, the better the correction of gradient-based distortion. However, notwithstanding this general rule, the invention will still be of benefit even if the sets of spherical harmonic coefficients used are not the optimum coefficients. In the end, all that is required to apply the corrective techniques of the invention is a set of spherical harmonic coefficients, however obtained, for each of the gradient magnets used to produce the MRI image which is to be corrected. In practice, these sets will be stored in a look-up table, e.g., on the hard drive of a computer system, for use in performing the corrective methods of the invention.

[0024] In its typical embodiment, the invention comprises acquiring time domain data using gradients which exhibit one or more of the non-linear, non-coincident center, axes misalignment, or winding error problems and then transforming the uncorrected time domain data using a fast Fourier transformation to provide a distorted image in image space. As the origin or cause of the distortion is known (namely, the non-primary harmonics calculated above), they are then used to correct the image space image.

[0025] More particularly, in accordance with a first aspect (see generally, the left column of FIG. 8), the invention provides a method for forming an image of an object using a magnetic resonance imaging system which (i) comprises X, Y and Z gradient magnets and (ii) defines an origin of an (x, y, z) coordinate system (the “x,y,z-coordinate system”), said X, Y, and Z gradient magnets producing magnetic fields having gradients (G_X, G_Y, G_Z), respectively, said method comprising:

[0026] (A) acquiring a magnetic resonance signal data set (the “MR data set”) for the object in three dimensions in data space using the X, Y, and Z gradient magnets;

[0027] (B) producing an uncorrected image function ρ(x,y,z) from the MR data set by a process which comprises performing a three dimensional fast Fourier transform on said data set; and

[0028] (C) producing a corrected image function ρ(x′,y′,z′) from the uncorrected image function ρ(x,y,z),where the coordinates (x′,y′,z′) represent a corrected position in the x,y,z-coordinate system for the value of the image function at (x,y,z), said (x′,y′,z′) coordinates being determined using an iterative mapping which employs coefficients of spherical harmonic expansions for the magnetic fields produced by the X, Y, and Z gradient magnets (the “spherical harmonic coefficients”).

[0029] Preferably, the iterative mapping is of the form of equation (21) below and employs a “quitting” parameter of the type set forth in equation (22). A suitable value for the quitting parameter is, for example, 10⁻³, which is sufficiently small to provide a correct image appearance. Generally, the smaller the quitting parameter, the longer the correction takes to perform. For example, about 3 iterations are needed to achieve the 10⁻³ criterion with a computation time of about 2 seconds. Increasing the computation time to, for example, 2.5 seconds would easily achieve a lower value, e.g., 10⁻⁴. In general, between two and four iterations are preferred as a reasonable compromise between computation time and satisfactory image correction. Other quitting parameters and iteration numbers can, of course, be used in the practice of this aspect of the invention or other aspects which involve iterations and/or quitting parameters.

[0030] In accordance with a second aspect (see generally, the middle column of FIG. 8), the invention provides a method for forming an image of an object using a magnetic resonance imaging system which (i) comprises X, Y and Z gradient magnets and (ii) defines an origin of an (x, y, z) coordinate system (the “x,y,z-coordinate system”), said X, Y, and Z gradient magnets producing magnetic fields having gradients (G_X, G_Y, G_Z), respectively, said method comprising:

[0031] (A) acquiring a magnetic resonance signal data set (the “MR data set”) for the object in three dimensions in data space using the X, Y, and Z gradient magnets;

[0032] (B) producing an uncorrected image function ρ(x,y,z) from the MR data set by a process which comprises performing a three dimensional fast Fourier transform on said data set; and

[0033] (C) producing a corrected image function ρ(x′,y′,z′) from the uncorrected image function ρ(x,y,z),where the coordinates (x′,y′,z′) represent a corrected position in the x,y,z-coordinate system for the value of the image function at (x,y,z),said (x′,y′,z′) coordinates being determined by solving a system of non-linear equations so that an error function f(x′,y′,z′) satisfies the relationship:
∥f(x′,y′,z′)∥≦ε,
where ε is a predetermined constant;

[0034] wherein the system of non-linear equations employs coefficients of spherical harmonic expansions for the magnetic fields produced by the X, Y, and Z gradient magnets (the “spherical harmonic coefficients”).

[0035] Preferably, the system of non-linear equations is solved iteratively (e.g., using a Levenberg-Marquardt non-linear curve fitting technique or a non-linear minimization technique) using (x₀′=x, y₀′=y, z₀′=z) as the initial (x′,y′,z′) values and an error function and system of non-linear equations of the form: $f\left(x^{\prime },y^{\prime },z^{\prime }\right)=\{\begin{array}{c}{f}_X\left(x^{\prime },y^{\prime },z^{\prime }\right)=B_X^{LN}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)-x{G}_X^L\cong 0\\ f_Y\left(x^{\prime },y^{\prime },z^{\prime }\right)=B_Y^{LN}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)-y{G}_Y^L\cong 0\\ f_Z\left(x^{\prime },y^{\prime },z^{\prime }\right)=B_Z^{LN}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)-z{G}_Z^L\cong 0\end{array}$
where:

[0036] (a) (x′, y′, z′) values and (r′, θ′, φ′) values are related by: $\hspace{1em}\{\begin{array}{c}{r}^{\prime }=\sqrt{x^{\mathrm{\prime 2}}+y^{\mathrm{\prime 2}}+z^{\mathrm{\prime 2}}}\\ {\theta }^{\prime }={\tan }^{-1}\left(\frac{\sqrt{x^{\mathrm{\prime 2}}+y^{\mathrm{\prime 2}}}}{z^{\prime }}\right)\\ {\varphi }^{\prime }={\tan }^{-1}\left(\frac{y^{\prime }}{x^{\prime }}\right)\end{array}$

[0037] (b) the a_V(n,m) and b_V(n,m) used in the expression for B_V(n,m)(r′,θ′,φ′) (see below) are the spherical harmonic coefficients, where the subscript V is one of X, Y, and Z;

[0038] (c) (G_X^L=a_X(1,1), G_Y^L=b_Y(1,1), G_Z^L=a_Z(1,0)) represent the linear parts of the gradients (G_X, G_Y, G_Z) of the magnetic fields produced by the X, Y, and Z gradient magnets, respectively; and $\{\begin{array}{c}{B}_X^{LN}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)=\sum _{n=0}\sum _m{B}_{X\left(n,m\right)}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)\\ B_Y^{LN}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)=\sum _{n=0}\sum _m{B}_{Y\left(n,m\right)}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)\\ B_Z^{LN}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)=\sum _{n=0}\sum _m{B}_{Z\left(n,m\right)}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)\end{array}\hspace{1em}$
represent the linear and non-linear parts of the magnetic fields produced by the X, Y, and Z gradient magnets, where:

[0039] (i) m and n are the indices of the summations and satisfy the relationship m≦n; and

[0040] (ii) B_V(n,m)(r′,θ′,φ′) is given by:
B_V(n,m)(r′,θ′,φ′)=r′ⁿ[a_V(n,m)cos(mφ′)+b_V(n,m)sin(mφ′)]P_(n,m)(cos θ′)
where the P_(n,m)(cos θ′) are associated Legendre functions and the subscript V is one of X, Y, and Z.

[0041] It should be noted that if each of the gradient magnets and the main (B₀) magnet (after any shimming) are effectively free of the non-coincident center problem, the axes misalignment problem, and the winding error problem (hereinafter referred to as the “base case”), then the summation in subparagraph (d) can begin with n=1 rather than n=0 if desired, since the lower order terms will be equal to, or effectively equal to, zero, i.e., their magnitudes will be zero or numerically small relative to other terms.

[0042] In these preferred embodiments, E preferably satisfies the relationship:
ε=R_es×10⁻³
where R_es is the resolution of the image given by:
R_es=A_image/N_pix,
where A_image is the area of the image (e.g., a two-dimensional image on a display screen) and N_pix is total number of pixels making up the image.

[0043] As with the first aspect of the invention, larger or smaller values of ε can be used in the practice of the invention, with larger values leading to reduced computation time and smaller values to increased time.

[0044] In accordance with the foregoing first and second aspects of the invention, the data space preferably comprises one of:

• • [0045] (a) k-space;
  • [0046] (b) time space;
  • [0047] (c) phase space;
  • [0048] (d) frequency space;
  • [0049] (e) time space along at least one coordinate direction and phase space along at least one other coordinate direction; or
  • [0050] (f) time space along at least one coordinate direction and frequency space along at least one other coordinate direction.

[0051] The first and second aspects of the invention have been described above in terms of correcting 3D imaging where the imaging coordinate system corresponds to the x,y,z-coordinate system. These aspects of the invention can also be used In the case of arbitrary 3D imaging, i.e., 3D imaging where the imaging coordinates are oriented in an arbitrary way with respect to the x,y,z-coordinate system. Such arbitrary 3D imaging can be of benefit in, for example, avoiding artifacts in brain images, specifically, in avoiding artifacts associated with the presence of air in nasal structures. The arbitrary 3D imaging is achieved in practice by selecting values for the X, Y, and Z gradient magnets which cause the imaging volume to have associated therewith a Cartesian coordinate system which is rotated relative to the x,y,z-coordinate system.

[0052] The relevant geometry is shown in FIG. 1A, where (ζ₁,ζ₂,ζ₃) represent the arbitrary coordinates, with ζ₁ and ζ₂ being perpendicular to one another and intersecting at O, and ζ₃ being an axis that passes through O and is perpendicular to the plane defined by ζ₁ and ζ₂. In this figure, θ(0≦θ≦π) is the angle between Oζ₃ and OZ, ψ (0≦ψ<2π) is the angle between OA and OX, and φ (0≦φ<2π) is the angle between OA and Oζ₁, where OA is the intersection of the plane Oζ₁ζ₂ and the plane OXY.

[0053] The coordinate transformation between the x,y,z-coordinate system and the (ζ₁,ζ₂,ζ₃) coordinate system can be written as: $r=A\zeta \mathrm{and}\zeta =A^{-1}r$ $\mathrm{where}$ $r={\left[x,y,z\right]}^T,\zeta ={\left[{\zeta }_1,{\zeta }_2,{\zeta }_3\right]}^T,A=\left[\begin{array}{ccc}{c}_{11}& c_{12}& c_{13}\\ c_{21}& c_{22}& c_{23}\\ c_{31}& c_{32}& c_{33}\end{array}\right],\mathrm{and}$ $c_{11}=\cos \mathrm{\psi cos\varphi }-\cos \mathrm{\theta sin}\mathrm{\psi sin}\varphi ,c_{21}=-\cos \mathrm{\psi sin}\varphi -\cos \mathrm{\theta sin}\mathrm{\psi cos}\varphi ,c_{31}=\sin \mathrm{\theta sin}\psi ,c_{12}=\sin \mathrm{\psi cos}\varphi -\cos \mathrm{\theta cos}\mathrm{\psi sin}\varphi ,c_{22}=-\sin \mathrm{\psi sin}\varphi +\cos \mathrm{\theta cos}\mathrm{\psi cos}\varphi ,c_{32}=-\sin \mathrm{\theta cos}\psi ,c_{32}=-\sin \mathrm{\theta cos}\psi ,c_{23}=\sin \mathrm{\theta cos}\varphi ;\mathrm{and}$ $c_{33}=\cos \theta$
and
C₁₁=cos ψ cos φ−cos θ sin ψ sin φ,
C₂₁=−cos ψ cos φ−cos θ sin ψ sin φ,
C₃₁=sin φ sin ψ,
C₁₂=sin ψ cos φ−cos θ cos ψ sin φ,
C₂₂=−sin ψ sin φ+cos θ cos ψ cos φ,
C₃₂=−sin θ cos ψ,
C₃₂=−sin θ cos ψ,
C₃₂=sin θ cos φ;
C₃₂=cos θ,
are the directional cosines.

[0054] The equations r=Aζ and ζ=A⁻¹r provide a complete relationship between the coordinate system r and the coordinate system ζ.

[0055] For such an arbitrarily oriented 3D image, the correction is performed by:

[0056] (1) taking an arbitrarily orientated image in the ζ coordinate system;

[0057] (2) transforming from ζ coordinate system to the r coordinate system using the r=Aζ; transformation;

[0058] (3) correcting the image in the r coordinate system using the correction methods of either the first or second aspects of the invention; and

[0059] (4) transforming back from the r coordinate system to the ζ coordinate system using the ζ=A⁻¹r transformation.

[0060] Thereafter, the corrected image can be displayed in the ζ coordinate system.

[0061] In accordance with a third aspect (see generally, the right column of FIG. 8), the invention provides a method for forming an image of an object using a magnetic resonance imaging system which (i) comprises X, Y and Z gradient magnets and (ii) defines an origin of an (x, y, z) coordinate system (the “x,y,z-coordinate system”), said X, Y, and Z gradient magnets producing magnetic fields having gradients (G_X, G_Y, G_Z), respectively, said method comprising:

[0062] (A) acquiring a magnetic resonance signal data set (the “MR data set”) for the object in three dimensions in data space using the X, Y, and Z gradient magnets, where the MR data set is of the form s(m_XΔt_X, m_YΔt_Y, m_ZΔt_Z) where:

[0063] (i) m_X, m_Y and m_Z are indices for the data points in the x, y, and z directions, respectively; and

[0064] (ii) Δt_X, Δt_Y and Δt_Z are equivalent time steps in the x, y, and z directions, respectively;

[0065] (B) producing an uncorrected image function ρ(l_XΔω_X^L,l_YΔψ_Y^L,l_ZΔω_Z^L) from the MR data set by a process which comprises performing a three dimensional fast Fourier transform on said data set, said uncorrected image function being in frequency space at a set of frequency points given by (l_XΔω_X^L, l_YΔω_Y^L, l_ZΔω_Z^L) where:

[0066] (i) l_X, l_Y and l_Z are indices for the frequency points in the ω_X^L, ω_Y^L and ω_Z^L directions, respectively;

[0067] (ii) Δω_X^L, Δω_Y^L and Δω_Z^L are frequency steps in the ω_X^L, ω_Y^L and ω_Z^L directions, respectively;

[0068] (iii) ω_X^L=γxG_X^L, ω_Y^L=γyG_Y^L, ω_Z^L=γzG_Z^L, where:

• • [0069] (a) γ is a constant; and
  • [0070] (b) G_X^L, G_Y^L, and G_Z^L represent the linear parts of the gradients (G_X, G_Y, G_Z) of the magnetic fields produced by the X, Y, and Z gradient magnets, respectively;

[0071] (iv) the ranges of l_X, l_Y, and l_Z are selected so that l_X·Δω_X^L, l_Y·Δω_Y^L, and l_Z·Δω_Z^L respectively span the frequency ranges that would arise from the X, Y, and Z gradient magnets if those magnets produced purely linear gradient fields in the x, y, and z directions of the x,y,z-coordinate system (e.g., the ranges of l_X, l_Y, and l_Z are selected based on the number of image points in each of the x, y, and z directions); and

[0072] (v) the fast Fourier transform is of the form: $\begin{array}{c}\rho \left(l_X{\mathrm{\Delta \omega }}_X^L,l_Y{\mathrm{\Delta \omega }}_Y^L,l_Z{\mathrm{\Delta \omega }}_Z^L\right)=\sum _{m_X}\sum _{m_Y}\sum _{m_Z}s\left(m_X\Delta t_X,m_Y\Delta t_Y,m_Z\Delta t_Z\right)\\ e^{i\left(l_X{\mathrm{\Delta \omega }}_X^L·m_X\Delta t_X+l_Y{\mathrm{\Delta \omega }}_Y^L·m_Y\Delta t_Y+l_Z{\mathrm{\Delta \omega }}_Z^L·m_Z\Delta t_Z\right)}\\ \Delta t_X\Delta t_Y\Delta t_Z\end{array}$

[0073] (C) producing a corrected image function ρ(x′,y′,z′) at (x′,y′,z′) in the x,y,z-coordinate system from the uncorrected image function ρ(l_XΔω_X^L,l_YΔω_Y^L,l_ZΔω_Z^L) using an interpolation function L in frequency space of the form: $\rho \left({\omega }_X^N,{\omega }_Y^N,{\omega }_Z^N\right)=L\left({\omega }_X^N,{\omega }_Y^N,{\omega }_Z^N,\rho \left(l_X{\mathrm{\Delta \omega }}_X^L,l_Y{\mathrm{\Delta \omega }}_Y^L,l_Z{\mathrm{\Delta \omega }}_Z^L\right)\right)$
where:

[0074] (i) ω_X^N, ω_Y^N, and ω_Z^N are each a function of (x′,y′,z′); and

[0075] (ii) ρ(ω_X^N, ω_Y^N, ω_Z^N) equals the desired ρ(x′,y′,z′);

[0076] wherein G_X^L, G_Z^L, G_Z^L, ω_X^N, ω_Y^N, and ω_Z^N are determined using coefficients of spherical harmonic expansions for the magnetic fields produced by the X, Y, and Z gradient magnets (the “spherical harmonic coefficients”).

[0077] The quantities Δt_X, Δt_Y and Δt_Z are referred to above as “equivalent time steps” for the following reasons. To encode a spatial position a gradient is used. If it is the read gradient, the time issue does not arise as time is determined by the sampling rate of the analog to digital converter (ADC). But this only encodes one direction, e.g., the y-direction. One must then ask how the other directions are encoded, e.g., the x-direction.

[0078] In the original spin warp experiment proposed by Ernst in 1976, the second gradient (the phase gradient) (and the third gradient) were applied at constant gradient strength and 256 time increments of each gradient were acquired. If the gradient strength is G and the time increment t, this means that 256 applications of the x gradient with gradient strength G but with time increasing from 0·t·G to 255·t·G were acquired, each being followed by a constant read gradient to give the xy encoding. The phase shift for a point with x-coordinate x was thus x·G·n·t, where n is the step in time. Multiplying this by γ gives the phase shift in units of radians.

[0079] But the same result is achieved if the gradient is stepped for a fixed on time, call this T. If g is the gradient step size, then the phase shift after the n^th incrementation of a constant gradient step is n·g·T·x. That is, n·g·x·T=x·G·n·t. Thus, the same phase encoding results whether time is incremented for a constant gradient or a gradient is incremented for constant time. It is in this sense that the phrase “equivalent time” is used in the above description of the third aspect of the invention and in the claims.

[0080] Certain preferred embodiments of this aspect of the invention are characterized by the following:

[0081] (a) a_V(n,m) and b_V(n,m) are the spherical harmonic coefficients, where the subscript V is one of X, Y, and Z;

[0082] (b) G_X^L=a_X(1,1), G_Y^L=b_Y(1,1), and G_Z^L=a_Z(1,0); and

[0083] (c) ω_X^N, ω_Y^N, and ω_Z^N are given by: $\hspace{1em}\{\begin{array}{c}{\omega }_X^N=\gamma \left(x^{\prime }+{\eta }_X\left(x^{\prime },y^{\prime },z^{\prime }\right)\right)G_X^L\\ {\omega }_Y^N=\gamma \left(y^{\prime }+{\eta }_Y\left(x^{\prime },y^{\prime },z^{\prime }\right)\right)G_Y^L\\ {\omega }_Z^N=\gamma \left(z^{\prime }+{\eta }_Z\left(x^{\prime },y^{\prime },z^{\prime }\right)\right)G_Z^L\end{array}$
where η_X(x′,y′,z′), η_Y(x′,y′,z′) and η_Z(x′,y′,z′) are given by: $\{\begin{array}{c}{\eta }_X\left(x^{\prime },y^{\prime },z^{\prime }\right)=\frac{B_X^N\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)}{G_X^L}\\ {\eta }_Y\left(x^{\prime },y^{\prime },z^{\prime }\right)=\frac{B_Y^N\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)}{G_Y^L}\\ {\eta }_Z\left(x^{\prime },y^{\prime },z^{\prime }\right)=\frac{B_Z^N\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)}{G_Z^L}\end{array}\mathrm{where}\{\begin{array}{c}{B}_X^N\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)=B_{X\left(0,0\right)}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)+B_{X\left(1,0\right)}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)+\\ \sum _{n=2}\sum _m{B}_{X\left(n,m\right)}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)\\ B_Y^N\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)=B_{Y\left(0,0\right)}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)+B_{Y\left(1,0\right)}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)+\\ \sum _{n=2}\sum _m{B}_{Y\left(n,m\right)}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)\\ B_Z^N\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)=B_{Z\left(0,0\right)}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)+B_{Z\left(1,1\right)}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)+\\ \sum _{n=2}\sum _m{B}_{Z\left(n,m\right)}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)\end{array}$
represent the non-linear parts of the magnetic fields produced by the X, Y, and Z gradient magnets, where:

[0084] (i) m and n are the indices of the summations and satisfy the relationship m≦n;

[0085] (ii) B_V(n,m)(r′,θ′,φ′) is given by:
B_V(n,m)(r′,θ′,φ′)=r′ⁿ[a_V(n,m)cos(mφ′)+b_V(n,m) sin(mφ′)]P_(n,m)(cos θ′)
where the P_(n,m)(cos θ′) are associated Legendre functions and the subscript V is one of X, Y, and Z; and

[0086] (iii) (x′, y′, z′) values and (r′, θ′, φ′) values are related by: $\hspace{1em}\{\begin{array}{c}{r}^{\prime }=\sqrt{x^{\mathrm{\prime 2}}+y^{\mathrm{\prime 2}}+z^{\mathrm{\prime 2}}}\\ {\theta }^{\prime }={\tan }^{-1}\left(\frac{\sqrt{x^{\mathrm{\prime 2}}+y^{\mathrm{\prime 2}}}}{z^{\prime }}\right)\\ {\varphi }^{\prime }={\tan }^{-1}\left(\frac{y^{\prime }}{x^{\prime }}\right)\end{array}.$

[0087] For the base case, the first two terms of each of the expressions set forth above for the non-linear parts of the magnetic fields produced by the X, Y, and Z gradient magnets are zero or effectively zero.

[0088] In other preferred embodiments, the interpolation function L is a sinc function such that: $L\left({\omega }_X^N,{\omega }_Y^N,{\omega }_Z^N,\rho \left(l_X{\mathrm{\Delta \omega }}_X^L,l_Y{\mathrm{\Delta \omega }}_Y^L,l_Z{\mathrm{\Delta \omega }}_Z^L\right)\right)=\sum _{l_Z}\sum _{l_Y}\sum _{l_X}\rho \left(l_X{\mathrm{\Delta \omega }}_X^L,l_Y{\mathrm{\Delta \omega }}_Y^L,l_Z{\mathrm{\Delta \omega }}_Z^L\right)\sin c\left(\frac{{\omega }_X^N}{{\mathrm{\Delta \omega }}_X^L}-l_X,\frac{{\omega }_Y^N}{{\mathrm{\Delta \omega }}_Y^L}-l_Y,\frac{{\omega }_Z^N}{{\mathrm{\Delta \omega }}_Z^L}-l_Z\right)$

[0089] In still further preferred embodiments, the interpolation function L is of the form: $L\left({\omega }_X^N,{\omega }_Y^N,{\omega }_Z^N,\rho \left(l_X{\mathrm{\Delta \omega }}_X^L,l_Y{\mathrm{\Delta \omega }}_Y^L,l_Z{\mathrm{\Delta \omega }}_Z^L\right)\right)=\sum _{l_Z}\sum _{l_Y}\sum _{l_X}{C}_{\left(l_X,l_Y,l_Z\right)}{\psi }_{X\left(l_X\right)}\left({\omega }_X^N\right){\psi }_{Y\left(l_Y\right)}\left({\omega }_Y^N\right){\psi }_{Z\left(l_Z\right)}\left({\omega }_Z^N\right)$
where:

[0090] (a) ψ_X(lX), ψ_Y(lY) and ψ_Z(lZ) are orthogonal bases functions; and

[0091] (b) the C_(lX,lY,lZ) are coefficients which satisfy the matrix equation:
AC=ρ
where the vector C is of the form:
C=[C_(1,1,1) C_(2,1,1) . . . C_(lX,lY,lZ)]^T
the vector ρ is of the form:
ρ=[ρ(ω_X(1)^L,ω_Y(1)^L,ω_Z(1)^L) ρ(ω_X(2)^L,ω_Y(1)^L,ω_Z(1)^L) . . . ρ(ω_X(lX)^L,ω_Y(lY)^L,ω_Z(lZ)^L)]^T
and the matrix A is of the form: $\left[\begin{array}{cccc}{\beta }_{\left(1,1\right)}& {\beta }_{\left(1,2\right)}& \dots & {\beta }_{\left(1,q\right)}\\ {\beta }_{\left(2,1\right)}& {\beta }_{\left(2,2\right)}& \dots & {\beta }_{\left(2,q\right)}\\ ⋮& ⋮& ⋰& ⋮\\ {\beta }_{\left(q,1\right)}& {\beta }_{\left(q,2\right)}& \dots & {\beta }_{\left(q,q\right)}\end{array}\right]$ $\mathrm{where}$ ${\beta }_{\left(1,1\right)}={\psi }_{X\left(1\right)}\left({\omega }_{X_{\left(1\right)}}^L\right){\psi }_{Y\left(1\right)}\left({\omega }_{Y_{\left(1\right)}}^L\right){\psi }_{Z\left(1\right)}\left({\omega }_{Z\left(1\right)}^L\right),{\beta }_{\left(1,2\right)}={\psi }_{X\left(2\right)}\left({\omega }_{X_{\left(1\right)}}^L\right){\psi }_{Y\left(1\right)}\left({\omega }_{Y_{\left(1\right)}}^L\right){\psi }_{Z\left(1\right)}\left({\omega }_{Z\left(1\right)}^L\right),⋮$ ${\beta }_{\left(1,q\right)}={\psi }_{X\left(l_X\right)}\left({\omega }_{X_{\left(1\right)}}^L\right){\psi }_{Y\left(l_Y\right)}\left({\omega }_{Y_{\left(1\right)}}^L\right){\psi }_{Z\left(l_Z\right)}\left({\omega }_{Z\left(1\right)}^L\right);$ ${\beta }_{\left(2,1\right)}={\psi }_{X\left(1\right)}\left({\omega }_{X\left(2\right)}^L\right){\psi }_{Y\left(1\right)}\left({\omega }_{Y_{\left(1\right)}}^L\right){\psi }_{Z\left(1\right)}\left({\omega }_{Z_{\left(1\right)}}^L\right),{\beta }_{\left(2,2\right)}={\psi }_{X\left(2\right)}\left({\omega }_{X\left(2\right)}^L\right){\psi }_{Y\left(1\right)}\left({\omega }_{Y_{\left(1\right)}}^L\right){\psi }_{Z\left(1\right)}\left({\omega }_{Z_{\left(1\right)}}^L\right),⋮$ ${\beta }_{\left(2,q\right)}={\psi }_{X\left(l_X\right)}\left({\omega }_{X\left(2\right)}^L\right){\psi }_{Y\left(l_Y\right)}\left({\omega }_{Y_{\left(1\right)}}^L\right){\psi }_{Z\left(l_Z\right)}\left({\omega }_{Z_{\left(1\right)}}^L\right);$ $⋮$ ${\beta }_{\left(q,1\right)}={\psi }_{X\left(1\right)}\left({\omega }_{X\left(l_X\right)}^L\right){\psi }_{Y\left(1\right)}\left({\omega }_{Y_{\left(l_Y\right)}}^L\right){\psi }_{Z\left(1\right)}\left({\omega }_{Z\left(l_Z\right)}^L\right),{\beta }_{\left(q,2\right)}={\psi }_{X\left(2\right)}\left({\omega }_{X\left(l_X\right)}^L\right){\psi }_{Y\left(1\right)}\left({\omega }_{Y_{\left(l_Y\right)}}^L\right){\psi }_{Z\left(1\right)}\left({\omega }_{Z\left(l_Z\right)}^L\right),{\beta }_{\left(q,q\right)}={\psi }_{X\left(l_X\right)}\left({\omega }_{X\left(l_X\right)}^L\right){\psi }_{Y\left(l_Y\right)}\left({\omega }_{Y_{\left(l_Y\right)}}^L\right){\psi }_{Z\left(l_Z\right)}\left({\omega }_{Z\left(l_Z\right)}^L\right),$
where q≦[l_X×l_Y×l_Z]_max.

[0092] In connection with these embodiments:

[0093] (a) the interpolation preferably comprises Lagrange's interpolation, Newton interpolation, finite difference interpolation, Hermite interpolation, spline interpolation, or combinations thereof; and/or

[0094] (b) the interpolation is applied to:

• • [0095] (i) all of the data points obtained by performing the three dimensional fast Fourier transform on the MR data set; or
  • [0096] (ii) a subset of the data points obtained by performing the three dimensional fast Fourier transform on the MR data set, said subset being in the vicinity of a region of interest of the object being imaged.

[0097] As with the first and second aspects of the invention, the third aspect of the invention can be employed to achieve corrected images for arbitrary 3D imaging. In terms of frequencies, the transformations are as follows:
ω_(X,Y,Z)=Aω_{(ζ1,ζ2,ζ3)}, and
ω_{(ζ1,ζ2,ζ3)}=A⁻¹ω_(X,Y,Z),
where

[0098] ω_(X,Y,Z)=[ω_X,ω_Y,ω_Z]^T,

[0099] ω_{(ζ1,ζ2,ζ3)}=[ω_{ζ1,ωζ2,ωζ3}]^T, and

[0100] A is the same as above.

[0101] The correction for this case is performed by:

[0102] (1) taking an arbitrarily orientated image in terms of the frequencies ω_{(ζ1,ζ2,ζ3)};

[0103] (2) transforming from the ω_{(ζ1,ζ2,ζ3)} frequencies to the ω_(X,Y,Z) frequencies using the ω_(X,Y,Z)=Aω_{(ζ1,ζ2,ζ3)} transformation;

[0104] (3) correcting the image in the ω_(X,Y,Z) frequency space using the correction methods of the third aspect of the invention; and

[0105] (4) transforming from the ω_(X,Y,Z) frequencies back to the ω_{(ζ1,ζ2,ζ3)} frequencies using the ω_{(ζ1,ζ2,ζ3)}=A⁻¹ω_(X,Y,Z) transformation.

[0106] Thereafter, the corrected image can be displayed in the ζ coordinate system.

[0107] In accordance with each of the first, second, and third aspects of the invention, ρ(x′,y′,z′) is preferably displayed as a corrected image. Even more preferably, as illustrated in FIG. 9, an image function ρ(x″,y″,z″) having uniform pixel dimensions is produced from ρ(x′,y′,z′) by performing an interpolation of the form:
ρ(x″,y″,z″)=H(x″,y″,z″,ρ(x′,y′,z′))
where H is an interpolation function and (x″,y″,z″) are uniformly distributed point coordinates in the x,y,z-coordinate system. Examples of suitable interpolation procedures include: Lagrange interpolation, Newton interpolation, finite differences interpolation, Hermite interpolation, spline interpolation, or combinations thereof. Once computed, the image function ρ(x″,y″,z″) is preferably displayed as a corrected image. In the case of case of arbitrary 3D imaging, the interpolation can be performed on the corrected image either before or after the back transformation step.

[0108] The foregoing three aspects of the invention relate to three dimensional imaging. The following three aspects relate to two dimensional slice imaging. In overview, the fourth, fifth, and sixth aspects employ correction techniques which generally correspond to those of the first, second, and third aspects, respectively, but in two dimensions, rather than three. FIG. 1B shows the coordinate nomenclature employed in the discussion of these two-dimensional aspects of the invention.

[0109] As in the first three aspects of the invention, the magnetic resonance imaging system for the fourth, fifth, and sixth aspects of the invention defines an origin of an (x, y, z) coordinate system (the “x,y,z-coordinate system”) and comprises X, Y and Z gradient magnets which produce magnetic fields having gradients (G_X, G_Y, G_Z), respectively.

[0110] Based on the geometry of FIG. 1B, the slice select gradient for an oblique slice can be defined as:
G₃=(p_X3G_X, p_Y3G_Y, p_Z3G_Z)
where (p_X3, p_Y3, p_Z3) are gradient strength coefficients of the X, Y, and Z gradient magnets (G_X, G_Y, G_Z), respectively, |G₃|=√{square root over ((p_X3G_X)²+(p_Y3G_Y)²+(p_Z3G_Z)²)}{square root over ((p_X3G_X)²+(p_Y3G_Y)²+(p_Z3G_Z)²)}{square root over ((p_X3G_X)²+(p_Y3G_Y)²+(p_Z3G_Z)²)} is the slice selection gradient amplitude, and the axis direction of that gradient is: ${\hat{\zeta }}_3=\left(\sin {\theta }_3\cos {\varphi }_3,\sin {\theta }_3\sin {\varphi }_3,\cos {\theta }_3\right)$ $\mathrm{where}$ ${\theta }_3={\tan }^{-1}\frac{\sqrt{{\left(p_{\mathrm{X3}}{G}_X\right)}^2+{\left(p_{\mathrm{Y3}}{G}_Y\right)}^2}}{p_{\mathrm{Z3}}{G}_Z},\mathrm{and}$ ${\varphi }_3={\tan }^{-1}\left(\frac{p_{\mathrm{Y3}}{G}_Y}{p_{\mathrm{X3}}{G}_X}\right).$

[0111] Transverse, coronal, and sagittal slices then correspond to:

[0112] Transverse:
(p_X3, p_Y3, p_Z3)=(0, 0, 1)
 G₃=(0, 0, G_Z)
|G₃|=|G_Z|
{circumflex over (ζ)}₃={circumflex over (z)};

[0113] Coronal:
(p_X3, p_Y3, p_Z3)=(0, 1, 0)
G₃=(0, G_Y, 0)
|G₃|=|G_Y|
{circumflex over (ζ)}₃=ŷ;

[0114] Sagittal:
(p_X3, p_Y3, p_Z3)=(1, 0, 0)
G₃=(G_X, 0, 0)
|G₃|=|G_X|
{circumflex over (ζ)}₃={circumflex over (x)}.

[0115] Again based on the geometry of FIG. 1B, the phase encoding gradient for an oblique slice can be defined as:
G₂=(p_X2G_X, p_Y2G_Y, p_Z2G_Z) such that {circumflex over (ζ)}₃·{circumflex over (ζ)}₂=0,
where (p_X2, p_Y2, p_Z2) are gradient strength coefficients of the X, Y and Z gradient magnets (G_X, G_Y, G_Z), respectively, |G₂|=√{square root over ((p_X2G_X)²+(p_Y2G_Y)²+(p_Z2G_Z)²)}{square root over ((p_X2G_X)²+(p_Y2G_Y)²+(p_Z2G_Z)²)}{square root over ((p_X2G_X)²+(p_Y2G_Y)²+(p_Z2G_Z)²)} is the phase encoding gradient amplitude, and the axis direction of that gradient is: ${\hat{\zeta }}_2=\left(\sin {\theta }_2\cos {\varphi }_2,\sin {\theta }_2\sin {\varphi }_2,\cos {\theta }_2\right)$ $\mathrm{where}$ ${\theta }_2={\tan }^{-1}\frac{\sqrt{{\left(p_{\mathrm{X2}}{G}_X\right)}^2+{\left(p_{\mathrm{Y2}}{G}_Y\right)}^2}}{p_{\mathrm{Z2}}{G}_Z},\mathrm{and}{\varphi }_2={\tan }^{-1}\left(\frac{p_{\mathrm{Y2}}{G}_Y}{p_{\mathrm{X2}}{G}_X}\right).$

[0116] Finally, based on FIG. 1B, the read gradient for an oblique slice can be defined as:
G₁=(p_X1G_X, p_Y1G_Y, p_Z1G_Z), such that {circumflex over (ζ)}₃·{circumflex over (ζ)}₁=0,
where (p_X1, P_Y1, p_Z1) are gradient strength coefficients of the X, Y and Z gradient magnets (G_X, G_Y, G_Z), respectively, |G₁|=√{square root over ((p_X1G_X)²+(p_Y1G_Y)²+(p_Z1G_Z)²)}{square root over ((p_X1G_X)²+(p_Y1G_Y)²+(p_Z1G_Z)²)}{square root over ((p_X1G_X)²+(p_Y1G_Y)²+(p_Z1G_Z)²)} is the read gradient amplitude, and the axis direction of that gradient is ${\hat{\zeta }}_1=\left(\sin {\theta }_1\cos {\varphi }_1,\sin {\theta }_1\sin {\varphi }_1,\cos {\theta }_1\right)$ $\mathrm{where}$ ${\theta }_1={\tan }^{-1}\frac{\sqrt{{\left(p_{\mathrm{X1}}{G}_X\right)}^2+{\left(p_{\mathrm{Y1}}{G}_Y\right)}^2}}{p_{\mathrm{Z1}}{G}_Z},\mathrm{and}{\varphi }_1={\tan }^{-1}\left(\frac{p_{\mathrm{Y1}}{G}_Y}{p_{\mathrm{X1}}{G}_X}\right).$

[0117] The coordinate system (ζ₁, ζ₂, ζ₃) defined in this way can be seen to be a linear combination of the X, Y and Z gradient magnets (G_X, G_Y, G_Z) in the coordinate system (x, y, z). The new coordinate system (ζ₁, ζ₂, ζ₃) has the same origin as the (x, y, z) coordinate system and the transformations between the coordinate systems can be written as:
r=r₀+A·(ζ−ζ₀), and
ζ=ζ₀+A⁻¹(r−r₀),
where A is a slice rotation matrix given by: $A=\left[\begin{array}{ccc}\sin {\theta }_1\cos {\varphi }_1& \sin {\theta }_1\sin {\varphi }_1& \cos {\theta }_1\\ \sin {\theta }_2\cos {\varphi }_2& \sin {\theta }_2\sin {\varphi }_2& \cos {\theta }_2\\ \sin {\theta }_3\cos {\varphi }_3& \sin {\theta }_3\sin {\varphi }_3& \cos {\theta }_3\end{array}\right]$
A⁻¹ is the inverse of the matrix, and
r=(x,y,z)^T,
r₀=(x₀,y₀,z₀)^T,
ζ=(ζ₁,ζ₂,ζ₃)^T, and
ζ₀=(0,0,ζ₃)^T.

[0118] In this transformation, r₀=(x₀,y₀,z₀)^T is the slice position in the (x, y, z) coordinate system with:
(x₀,y₀,z₀)=(ζ₃ sin θ₃ cos φ₃, ζ₃ sin θ₃ sin φ₃, ζ₃ cos θ₃),
and (0,0,ζ₃) is the slice position along the slice axis {circumflex over (ζ)}₃ in the (ζ₁, ζ₂, ζ₃) coordinate system, with:
ζ₃=√{square root over (x₀²+y₀²+z₀²)}.

[0119] As discussed above, the linear parts of the gradient fields of the X, Y, and Z gradient magnets in the x,y,z-coordinate system are given by (G_X^L=a_X(1,1), G_Y^L=b_Y(1,1), Z_X^L=a_Z(1,0)). The linear parts of the G₁ and G₂ gradients (also referred to herein and in the claims as, respectively, the “G₁ gradient” and the “G₂ gradient,” or collectively as the “gradients (G₁, G₂)”) can be written as:
G₁^L=G₁^L·ζ₁=p_X1a_X(1,1)+p_Y1a_Y(1,1)+p_Z1a_Z(1,0)
G₂^L=G₂^L·ζ₂=p_X2a_X(1,1)+p_Y2a_Y(1,1)+p_Z2a_Z(1,0).

[0120] In terms of the foregoing nomenclature, the invention in accordance with a fourth aspect provides a method for forming an image of an object using a magnetic resonance imaging system which (i) comprises X, Y and Z gradient magnets and (ii) defines an origin of an (x, y, z) coordinate system (the “x,y,z-coordinate system”), said X, Y, and Z gradient magnets producing magnetic fields having gradients (G_X, G_Y, G_Z), respectively, said method comprising:

[0121] (A) acquiring a magnetic resonance signal data set (the “MR data set”) for a slice (e.g., a transverse, coronal, sagittal, or oblique slice) of the object in two dimensions in data space using a linear combination of the X, Y, and Z gradient magnets;

[0122] (B) producing an uncorrected image function ρ(ζ₁,ζ₂) from the MR data set by a process which comprises performing a two dimensional fast Fourier transform on said data set, where ζ₁ and ζ₂ are read and phase encoding coordinates, respectively, of a (ζ₁,ζ₂,ζ₃) coordinate system as defined above;

[0123] (C) producing a corrected image function ρ(ζ₁′,ζ₂′) from the uncorrected image function ρ(ζ₁,ζ₂), where the coordinates (ζ₁′,ζ₂′) represent a corrected position in the (ζ₁,ζ₂,ζ₃) coordinate system for the value of the image function at (ζ₁,ζ₂), said (ζ₁′,ζ₂′) coordinates being determined using an iterative mapping which employs coefficients of spherical harmonic expansions for the magnetic fields produced by the X, Y, and Z gradient magnets (the “spherical harmonic coefficients”).

[0124] In certain preferred embodiments of this aspect of the invention, the iterative mapping is of the form: $\{\begin{array}{c}{\zeta }_1^{\prime }={\zeta }_1-{\eta }_1\left({\zeta }_1^{\prime },{\zeta }_2^{\prime }\right)\\ {\zeta }_2^{\prime }={\zeta }_2-{\eta }_2\left({\zeta }_1^{\prime },{\zeta }_2^{\prime }\right)\end{array}\mathrm{where}\{\begin{array}{c}{\zeta }_{1\left(1\right)}^{\prime }={\zeta }_1-{\eta }_1\left({\zeta }_1,{\zeta }_2\right)\\ {\zeta }_{2\left(2\right)}^{\prime }={\zeta }_2-{\eta }_2\left({\zeta }_1,{\zeta }_2\right)\end{array},\{\begin{array}{c}{\zeta }_{1\left(2\right)}^{\prime }={\zeta }_1-{\eta }_1\left({\zeta }_{1\left(1\right)}^{\prime },{\zeta }_{2\left(1\right)}^{\prime }\right)\\ {\zeta }_{2\left(2\right)}^{\prime }={\zeta }_2-{\eta }_2\left({\zeta }_{1\left(1\right)}^{\prime },{\zeta }_{2\left(1\right)}^{\prime }\right)\end{array},\dots ,\{\begin{array}{c}{\zeta }_{1\left(n\right)}^{\prime }={\zeta }_1-{\eta }_1\left({\zeta }_{1\left(n-1\right)}^{\prime },{\zeta }_{2\left(n-1\right)}^{\prime }\right)\\ {\zeta }_{2\left(n\right)}^{\prime }={\zeta }_2-{\eta }_2\left({\zeta }_{1\left(n-1\right)}^{\prime },{\zeta }_{2\left(n-1\right)}^{\prime }\right)\end{array}$
where n is the iteration number and η_l(ζ₁, ζ₂)=η₁(ζ₁, ζ₂, ζ₃) and η₂(ζ₁, ζ₂)=η₂(ζ₁, ζ₂, ζ₃) are given by: $\hspace{1em}\{\begin{array}{c}{\eta }_1\left({\zeta }_1,{\zeta }_2,{\zeta }_3\right)=\frac{B_1^N\left(r,\theta ,\varphi \right)}{G_1^L·{\hat{\zeta }}_1}\\ {\eta }_2\left({\zeta }_1,{\zeta }_2,{\zeta }_3\right)=\frac{B_2^N\left(r,\theta ,\varphi \right)}{G_2^L·{\hat{\zeta }}_2}\end{array}$
where

[0125] (a) r, θ and φ are defined relative to the same origin as the ζ₁, ζ₂ and ζ₃ coordinate system and are given by: $\left(i\right)r=r_0+A·\left(\zeta -{\zeta }_0\right);\mathrm{and}\left(ii\right)\left(r=\sqrt{x^2+y^2+z^2},\theta ={\tan }^{-1}\left(\frac{\sqrt{x^2+y^2}}{z}\right),\varphi ={\tan }^{-1}\left(\frac{y}{x}\right)\right);$

[0126] (b) (G₁^L, G₂^L) represent the linear part of the gradients (G₁, G₂) of the magnetic fields produced by the X, Y, and Z gradient magnets; and

[0127] (c) $\begin{array}{cc}\hspace{1em}\{\begin{array}{c}\begin{array}{c}{B}_1^N\left(r,\theta ,\varphi \right)=p_{\mathrm{X1}}(B_{X\left(0,0\right)}\left(r,\theta ,\varphi \right)+B_{X\left(1,0\right)}\left(r,\theta ,\varphi \right)+\\ \sum _{n=2}\sum _m{B}_{X\left(n,m\right)}\left(r,\theta ,\varphi \right))+p_{\mathrm{Y1}}(B_{Y\left(0,0\right)}\left(r,\theta ,\varphi \right)+\\ B_{Y\left(1,0\right)}\left(r,\theta ,\varphi \right)+\sum _{n=2}\sum _m{B}_{Y\left(n,m\right)}\left(r,\theta ,\varphi \right))+\\ p_{\mathrm{Z1}}\left(B_{Z\left(0,0\right)}\left(r,\theta ,\varphi \right)+B_{Z\left(1,1\right)}\left(r,\theta ,\varphi \right)+\sum _{n=2}\sum _m{B}_{Z\left(n,m\right)}\left(r,\theta ,\varphi \right)\right)\end{array}\\ \begin{array}{c}{B}_2^N\left(r,\theta ,\varphi \right)=p_{\mathrm{X2}}(B_{X\left(0,0\right)}\left(r,\theta ,\varphi \right)+B_{X\left(1,0\right)}\left(r,\theta ,\varphi \right)+\\ \sum _{n=2}\sum _m{B}_{X\left(n,m\right)}\left(r,\theta ,\varphi \right))+p_{\mathrm{Y2}}(B_{Y\left(0,0\right)}\left(r,\theta ,\varphi \right)+B_{Y\left(1,0\right)}\left(r,\theta ,\varphi \right)+\\ \sum _{n=2}\sum _m{B}_{Y\left(n,m\right)}\left(r,\theta ,\varphi \right))+p_{\mathrm{Z2}}(B_{Z\left(0,0\right)}\left(r,\theta ,\varphi \right)+B_{Z\left(1,1\right)}\left(r,\theta ,\varphi \right)+\\ \sum _{n=2}\sum _m{B}_{Z\left(n,m\right)}\left(r,\theta ,\varphi \right))\end{array}\end{array}& \end{array}$
represent the non-linear parts of the magnetic fields produced by the X, Y, and Z gradient magnets; where

[0128] (i) m and n are the indices of the summations and satisfy the relationship m≦n;

[0129] (ii) the B_{V(n m)}(r,θ,φ) are given by:
B_V(n,m)(r,θ,φ)=rⁿ[a_V(n,m)cos(mφ)+b_V(n,m)sin(mφ)]P_(n,m)(cos θ),
where the P_(n,m)(cos θ) are the associated Legendre functions and the subscript V is one of X, Y, and Z; and

[0130] (iii) a_V(n,m) and b_V(n,m) are the spherical harmonic coefficients, where again, the subscript V is one of X, Y, and Z.

[0131] For the base case, the magnitudes of the X(0,0), X(1,0), Y(0,0), Y(1,0), and Z(0,0) field components are zero or effectively zero.

[0132] Preferably, the iterative mapping is continued until the following condition is satisfied: $\frac{\sqrt{{\left({\zeta }_{1\left(n\right)}^{\prime }-{\zeta }_{1\left(n-1\right)}^{\prime }\right)}^2+{\left({\zeta }_{2\left(n\right)}^{\prime }-{\zeta }_{2\left(n-1\right)}^{\prime }\right)}^2}}{\sqrt{{\left({\zeta }_{1\left(n\right)}^{\prime }-{\zeta }_1\right)}^2}+{\left({\zeta }_{2\left(n\right)}^{\prime }-{\zeta }_2\right)}^2}\le \varepsilon$
where, for example, ε=10⁻³. In such a case, the number of iterations is typically between 2 and 4.

[0133] In accordance with a fifth aspect, the invention provides a method for forming an image of an object using a magnetic resonance imaging system which (i) comprises X, Y and Z gradient magnets and (ii) defines an origin of an (x, y, z) coordinate system (the “x,y,z-coordinate system”), said X, Y, and Z gradient magnets producing magnetic fields having gradients (G_X, G_Y, G_Z), respectively, said method comprising:

[0134] (A) acquiring a magnetic resonance signal data set (the “MR data set”) for a slice (e.g., a transverse, coronal, sagittal, or oblique slice) of the object in two dimensions in data space using a linear combination of the X, Y, and Z gradient magnets;

[0135] (B) producing an uncorrected image function ρ(ζ₁,ζ₂) from the MR data set by a process which comprises performing a two dimensional fast Fourier transform on said data set, where ζ₁ and ζ₂ are read and phase encoding coordinates, respectively, of a (ζ₁,ζ₂,ζ₃) coordinate system as defined above;

[0136] (C) producing a corrected image function ρ(ζ₁′,ζ₂′) from the uncorrected image function ρ(ζ₁,ζ₂), where the coordinates (ζ₁′,ζ₂′) represent a corrected position in the (ζ₁,ζ₂,ζ₃) coordinate system for the value of the image function at (ζ₁,ζ₂), said (ζ₁′,ζ₂′) coordinates being determined by solving a system of non-linear equations so that an error function f (ζ₁′,ζ₂′) satisfies the relationship:
∥f((ζ₁′,ζ′)∥≦ε
where ε is a predetermined constant;

[0137] wherein the system of non-linear equations employs coefficients of spherical harmonic expansions for the magnetic fields produced by the X, Y, and Z gradient magnets (the “spherical harmonic coefficients”).

[0138] Preferably, the system of non-linear equations is solved iteratively (e.g., using a Levenberg-Marquardt non-linear curve fitting technique or a non-linear minimization technique) using ζ′₁₍₀₎=ζ₁ and ζ′₂₍₀₎=ζ₂ as the initial (ζ′₁,ζ′₂) values and an error function and system of non-linear equations of the form: $f\left({\zeta }_1^{\prime },{\zeta }_2^{\prime }\right)=\{\begin{array}{c}{f}_1\left({\zeta }_1^{\prime },{\zeta }_2^{\prime },{\zeta }_3^{\prime }\right)=B_1^{LN}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)-G_1^L·{\zeta }_1\cong 0\\ f_2\left({\zeta }_1^{\prime },{\zeta }_2^{\prime },{\zeta }_3^{\prime }\right)=B_2^{LN}\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)-G_2^L·{\zeta }_2\cong 0\end{array}$
where:

[0139] (a) (r′, θ′, φ′) and (ζ₁′, ζ₂′, ζ₃′) values are related by: $\left(r^{\prime }=\sqrt{x^{\mathrm{\prime 2}}+y^{\mathrm{\prime 2}}+z^{\mathrm{\prime 2}}},{\theta }^{\prime }={\tan }^{-1}\left(\frac{\sqrt{x^{\mathrm{\prime 2}}+y^{\mathrm{\prime 2}}}}{z^{\prime }}\right),{\varphi }^{\prime }={\tan }^{-1}\left(\frac{y^{\prime }}{x^{\prime }}\right)\right)$ $\mathrm{and}$ $r^{\prime }=r_0^{\prime }+A·\left({\zeta }^{\prime }-{\zeta }_0^{\prime }\right)$

[0140] (b) (G₁^L, G₂^L) represent the linear part of the gradients (G₁, G₂) of the magnetic fields produced by the X, Y, and Z gradient magnets; and

[0141] (c) $\{\begin{array}{c}\begin{array}{c}{B}_1^{\mathrm{LN}}\left(r,\theta ,\varphi \right)=p_{\mathrm{X1}}\left(\sum _{n=0}\sum _m{B}_{X\left(n,m\right)}\left(r,\theta ,\varphi \right)\right)+\\ p_{\mathrm{Y1}}\left(\sum _{n=0}\sum _m{B}_{Y\left(n,m\right)}\left(r,\theta ,\varphi \right)\right)+p_{\mathrm{Z1}}\left(\sum _{n=0}\sum _m{B}_{Z\left(n,m\right)}\left(r,\theta ,\varphi \right)\right)\end{array}\\ \begin{array}{c}{B}_2^{\mathrm{LN}}\left(r,\theta ,\varphi \right)=p_{\mathrm{X2}}\left(\sum _{n=0}\sum _m{B}_{X\left(n,m\right)}\left(r,\theta ,\varphi \right)\right)+\\ p_{\mathrm{Y2}}\left(\sum _{n=0}\sum _m{B}_{Y\left(n,m\right)}\left(r,\theta ,\varphi \right)\right)+p_{\mathrm{Z2}}\left(\sum _{n=0}\sum _m{B}_{Z\left(n,m\right)}\left(r,\theta ,\varphi \right)\right)\end{array}\end{array}\hspace{1em}$
represent the linear and non-linear parts of the magnetic fields produced by the X, Y, and Z gradient magnets, where

[0142] (i) m and n are the indices of the summations and satisfy the relationship m≦n;

[0143] (ii) the B_V(n,m)(r′,θ′,φ′) are given by:
B_V(n,m)(r′,θ′,φ′)=r′ⁿ[a_V(n,m)cos(mφ′)+b_V(n,m)sin(mφ′)]P_(n,m)(cos θ′)
where the P_(n,m)(cos θ′) are the associated Legendre functions and the subscript V is one of X, Y, and Z; and

[0144] (iii) a_V(n,m) and b_V(n,m) are the spherical harmonic coefficients, where again, the subscript V is one of X, Y, and Z.

[0145] For the base case, the summations in subparagraph (c) can begin with n=1 rather than n=0 if desired.

[0146] As in the second aspect of the invention, ε preferably satisfies the relationship:
ε=R_es×10⁻³
where R_es is the resolution of the image given by:
R_es=A_image/N_pix,
where A_image is the area of the image and N_pix is total number of pixels making up the image.

[0147] In accordance with a sixth aspect, the invention provides a method for forming an image of an object using a magnetic resonance imaging system which (i) comprises X, Y and Z gradient magnets and (ii) defines an origin of an (x, y, z) coordinate system (the “x,y,z-coordinate system”), said X, Y, and Z gradient magnets producing magnetic fields having gradients (G_X, G_Y, G_Z), respectively, said method comprising:

[0148] (A) acquiring a magnetic resonance signal data set (the “MR data set”) for a slice (e.g., a transverse, coronal, sagittal, or oblique slice) of the object in two dimensions in data space using a linear combination of the X, Y, and Z gradient magnets, where the MR data set is of the form s(m₁Δt₁, m₂Δt₂) where:

[0149] (i) m₁ and m₂ are respectively indices for the k-space data points in the {circumflex over (ζ)}₁ and {circumflex over (ζ)}₂ directions defined above; and

[0150] (ii) Δt₁ and Δt₂ are respectively equivalent time steps in k-space for the {circumflex over (ζ)}₁ and {circumflex over (ζ)}₂ directions defined above, where equivalent time has the meaning discussed above in connection with the third aspect of the invention;

[0151] (B) producing an uncorrected image function ρ(l₁Δω₁^L,l₂Δω₂^L) from the MR data set by a process which comprises performing a two dimensional fast Fourier transform on said data set, said uncorrected image function being in frequency space at a set of frequency points given by (l₁Δω_l^L, L₂Δω₂^L) where:

[0152] (i) l₁ and l₂ are indices for the frequency points in the ω₁^L and ω₂^L directions, respectively;

[0153] (ii) Δω₁^L and Δω₂^L are frequency steps in the ω₁^L and ω₂^L directions, respectively;

[0154] (iii) ω₁^L=γG₁^L·ζ₁, ω₂^L=γG₂^L·ζ₂, where:

• • [0155] (a) γ is a constant; and
  • [0156] (b) (G₁^L, G₂^L) represent the linear parts of the gradients (G₁, G₂) of the magnetic fields produced by the X, Y, and Z gradient magnets;

[0157] (iv) the range of l₁ is selected so that l₁·Δω₁^L spans the frequency range that would arise from the G₁ gradient if that gradient were purely linear in the {circumflex over (ζ)}₁ direction and the range of l₁ is selected so that l₁·Δω₁^L spans the frequency range that would arise from the G₂ gradient if that gradient were purely linear field in the {circumflex over (ζ)}₂ direction; and

[0158] (v) the fast Fourier transform is of the form: $\rho \left(l_1{\mathrm{\Delta \omega }}_1^L,l_2{\mathrm{\Delta \omega }}_2^L\right)=\sum _{m_2}\sum _{m_1}s\left(m_1\Delta t_1,m_2\Delta t_2\right)e^{i\left(l_1{\mathrm{\Delta \omega }}_1^L·m_1\Delta t_1+l_2\Delta {\omega }_2^L·m_2\Delta t_2\right)\Delta t_1\Delta t_2}$

[0159] (C) producing a corrected image function ρ(ζ₁′,ζ₂′) in frequency space at (ζ₁′, ζ₂′, ζ₃′) in the (ζ₁,ζ₂,ζ₃) coordinate system from the uncorrected image function ρ(l₁Δω₁^L,l₂Δω₂^L) using an interpolation function L in frequency space of the form:
ρ(ω₁^N(ζ₁′,ζ₂′,ζ₃′),ω₂^N(ζ₁′,ζ₂′,ζ₃′))=L(ω₁^N,ω₂^N,ρ(l₁Δω₁^L,l₂Δω₂^L))
where:

[0160] (i) ρ(ω₁^N(ζ₁′,ζ₂′,ζ₃′),ω₂^N(ζ₁′,ζ₂′,ζ₃′)) is the desired ρ(ζ₁′,ζ₂′); and

[0161] (ii) G₁^L,G₂^L,ω₁^N(ω₁^N(ζ₁′,ζ₂′,ζ₃′), and ω₂^N(ζ₁′,ζ₂′,ζ₃′) are determined using coefficients of spherical harmonic expansions for the magnetic fields produced by the X, Y, and Z gradient magnets (the “spherical harmonic coefficients”).

[0162] Certain preferred embodiments of this aspect of the invention are characterized by the following:

[0163] (a) a_V(n,m) and b_V(n,m) are the spherical harmonic coefficients, where the subscript V is one of X, Y, and Z;

[0164] (b) G₁^L=p_X1a_X(1,1)+p_Y1a_Y(1,1)+p_Z1a_Z(1,0), and
G₂^L=p_X2a_X(1,1)+p_Y2a_Y(1,1)+P_Z2a_Z(1,0); and

[0165] (c) ω₁^N(ζ₁′, ζ₂′, ζ₃′) and ω₂^N(ζ₁′, ζ₂′, ζ₃′) are given by: $\hspace{1em}\{\begin{array}{c}{\omega }_1^N\left({\zeta }_1^{\prime },{\zeta }_2^{\prime },{\zeta }_3^{\prime }\right)=\gamma \left({\zeta }_1^{\prime }+{\eta }_1\left({\zeta }_1^{\prime },{\zeta }_2^{\prime },{\zeta }_3^{\prime }\right)\right)\left(G_1^L·{\hat{\zeta }}_1\right)\\ {\omega }_2^N\left({\zeta }_1^{\prime },{\zeta }_2^{\prime },{\zeta }_3^{\prime }\right)=\gamma \left({\zeta }_2^{\prime }+{\eta }_2\left({\zeta }_1^{\prime },{\zeta }_2^{\prime },{\zeta }_3^{\prime }\right)\right)\left(G_2^L·{\hat{\zeta }}_2\right)\end{array}$
where η₁(ζ₁′,ζ₂′,ζ₃′) and η₂(ζ₁′,ζ₂′,ζ₃′) are given by: $\{\begin{array}{c}{\eta }_1\left({\zeta }_1^{\prime },{\zeta }_2^{\prime },{\zeta }_3^{\prime }\right)=\frac{B_1^N\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)}{G_1^L·{\hat{\zeta }}_1}\\ {\eta }_2\left({\zeta }_1^{\prime },{\zeta }_2^{\prime },{\zeta }_3^{\prime }\right)=\frac{B_2^N\left(r^{\prime },{\theta }^{\prime },{\varphi }^{\prime }\right)}{G_2^L·{\hat{\zeta }}_2}\end{array}\mathrm{where}\{\begin{array}{c}\begin{array}{c}{B}_1^N\left(r,\theta ,\varphi \right)=p_{\mathrm{X1}}\left(B_{X\left(0,0\right)}\left(r,\theta ,\varphi \right)+B_{X\left(1,0\right)}\left(r,\theta ,\varphi \right)+\sum _{n=2}\sum _m{B}_{X\left(n,m\right)}\left(r,\theta ,\varphi \right)\right)+\\ p_{\mathrm{Y1}}\left(B_{Y\left(0,0\right)}\left(r,\theta ,\varphi \right)+B_{Y\left(1,0\right)}\left(r,\theta ,\varphi \right)+\sum _{n=2}\sum _m{B}_{Y\left(n,m\right)}\left(r,\theta ,\varphi \right)\right)+\\ p_{\mathrm{Z1}}\left(B_{Z\left(0,0\right)}\left(r,\theta ,\varphi \right)+B_{Z\left(1,1\right)}\left(r,\theta ,\varphi \right)+\sum _{n=2}\sum _m{B}_{Z\left(n,m\right)}\left(r,\theta ,\varphi \right)\right)\end{array}\\ \begin{array}{c}{B}_2^N\left(r,\theta ,\varphi \right)=p_{\mathrm{X2}}\left(B_{X\left(0,0\right)}\left(r,\theta ,\varphi \right)+B_{X\left(1,0\right)}\left(r,\theta ,\varphi \right)+\sum _{n=2}\sum _m{B}_{X\left(n,m\right)}\left(r,\theta ,\varphi \right)\right)+\\ p_{\mathrm{Y2}}\left(B_{Y\left(0,0\right)}\left(r,\theta ,\varphi \right)+B_{Y\left(1,0\right)}\left(r,\theta ,\varphi \right)+\sum _{n=2}\sum _m{B}_{Y\left(n,m\right)}\left(r,\theta ,\varphi \right)\right)+\\ p_{\mathrm{Z2}}\left(B_{Z\left(0,0\right)}\left(r,\theta ,\varphi \right)+B_{Z\left(1,1\right)}\left(r,\theta ,\varphi \right)+\sum _{n=2}\sum _m{B}_{Z\left(n,m\right)}\left(r,\theta ,\varphi \right)\right)\end{array}\end{array}$
represent the non-linear parts of the magnetic fields produced by the X, Y, and Z gradient magnets, where: