Current-density imaging using ultra-low-field MRI with zero-field encoding

Magnetic Resonance Imaging 32 (2014) 766–770

Contents lists available at ScienceDirect

Magnetic Resonance Imaging journal homepage: www.mrijournal.com

Current-density imaging using ultra-low-field MRI with zero-field encoding Panu T. Vesanen a,⁎, Jaakko O. Nieminen a, Koos C.J. Zevenhoven a, Yi-Cheng Hsu a, b, Risto J. Ilmoniemi a a b

Department of Biomedical Engineering and Computational Science, Aalto University School of Science, P.O. Box 12200, FI-00076 AALTO, Finland Department of Mathematics, National Taiwan University, Taipei, TW-10617, Taiwan

a r t i c l e

i n f o

Article history: Received 14 February 2013 Revised 4 August 2013 Accepted 19 January 2014 Keywords: Ultra-low-field MRI Current-density imaging Zero field Rotation-free

a b s t r a c t Electric current density can be measured noninvasively with magnetic resonance imaging (MRI). Determining all three components of the current density, however, requires physical rotation of the sample or current injection from several directions when done with conventional methods. However, the emerging technology of ultra-low-field (ULF) MRI, in which the signal encoding and acquisition is conducted at a microtesla-range magnetic field, offers new possibilities. The low applied magnetic fields can even be switched off completely within the pulse sequence, increasing the flexibility of the available sequences. In this article, we present a ULFMRI sequence designed for obtaining all three components of a current-density pattern without the need of sample rotations. The sequence consists of three steps: prepolarization of the sample, signal encoding in the current-density-associated magnetic field without applying any MRI fields, and spatial encoding in a microtesla-range field using any standard ULF-MRI sequence. The performance of the method is evaluated by numerical simulations. The method may find applications, e.g., in noninvasive conductivity imaging of tissue. © 2014 Elsevier Inc. All rights reserved.

1. Introduction Electric current density can be measured noninvasively with magnetic resonance imaging (MRI). Methods for imaging of both static and alternating magnetic fields have been presented [1–8]. The methods are based on measuring the phase changes that the magnetic field BJ, generated by the current density J, produces on the spin system. Typically, the B0 field is orders of magnitude stronger than BJ so that only the component of BJ along B0 can be detected. The measurement of the remaining components typically requires physical rotation of the sample. Current-density-imaging (CDI) methods without the need to rotate the object have been presented as well. Radio-frequency currentdensity imaging allows the measurement of all three components of a current density varying at the Larmor frequency [4,5]. A related method is electric-properties tomography [9], which allows determining the conductivity and permittivity of a sample without rotation. Finally, methods for calculating all the three components of the current density based on only a single magnetic field component have also been presented [10,11]. However, the methods require current injection from several different directions, and the reconstruction algorithms are vulnerable to measurement noise. Ultra-low-field (ULF) MRI is based on signal encoding in a microtesla-range magnetic field and signal detection using superconducting quantum interference devices (SQUIDs) [12–14]. To ⁎ Corresponding author. Tel.: +358 405074216. E-mail address: [email protected] (P.T. Vesanen). http://dx.doi.org/10.1016/j.mri.2014.01.012 0730-725X/© 2014 Elsevier Inc. All rights reserved.

increase the signal-to-noise ratio (SNR) of the measurement, the sample is typically prepolarized [15] in a stronger magnetic field of 10–100 mT. To suppress external magnetic field noise, ULF-MRI devices are typically operated inside a magnetically shielded room. In the context of current-density imaging, a clear advantage of ULF-MRI is the ability to rapidly switch off the B0 field. As will be shown, in the presence of zero external field, information about all three components of BJ can be encoded in the magnetization direction, and no sample rotations are needed. An important application of current-density imaging is the determination of the sample’s conductivity [10,11,16–19]. Accurate conductivity information, e.g., about the human brain is important in several functional neuroimaging methods. In magnetoencephalography (MEG) and especially in electroencephalography (EEG), estimating the location and amplitude of brain activity benefits from accurate knowledge of the conductivity. Because the conductivity of tissue depends on frequency, the existing methods that give conductivity at a megahertz-range Larmor frequency are not useful for assisting in MEG and EEG source estimation. Another possible application of current-density imaging might be direct neuronal imaging (DNI) [20–28], which aims to use MRI to directly measure the currents produced by brain activity. In this paper, we present a method to image static magnetic fields and current densities with ULF MRI. The approach is designed to reveal all three components of the current density and the associated magnetic field without the need to rotate the object. In the following, we refer to the approach as zero-field-encoded current-density imaging (Z-CDI). The performance and stability of the method is evaluated using

P.T. Vesanen et al. / Magnetic Resonance Imaging 32 (2014) 766–770

767

numerical simulations. In an accompanying paper [29], we present an alternative method, called adiabatic current-density imaging (A-CDI), for achieving the same goal. 2. Methods ULF-MRI encoding fields are typically in the microtesla range, which allows switching of all the applied fields within the pulse sequence. For example, if a current density J causing a magnetic field BJ is applied when all the MRI fields have been turned off, the precession axis and frequency of the pre-polarized magnetization are set only by BJ. Thus, information about BJ will be encoded in the direction of the magnetization; this information can be extracted from the phase and magnitude of the MRI signal. The following description of the method concentrates on finding the magnetic field BJ; the current density J can then be obtained simply by J¼

1 ∇  BJ ; μ0

ð1Þ

where μ0 is the permeability of vacuum. Specifically, consider a single voxel with a homogeneous magnetic field BJ. The proton magnetization dynamics in the voxel can be described by the Bloch equation, dm ¼ γm  BJ ; dt

ð2Þ

where h m is the iT magnetization, γ is the gyromagnetic ratio, and BJ ¼ Bx By Bz . This can be written in matrix form as dm ¼ Am; dt

ð3Þ

where 2

3 2 0 mx m ¼ 4 my 5; and A ¼ γ 4 −Bz By mz

Bz 0 −Bx

3 −By Bx 5: 0

ð4Þ

Consider a pre-polarized phase-contrast sequence in which Bp ¼ Bp ex is first applied to create an x-directional magnetization m0 = [m0 0 0]T, after which only BJ is applied during a time τ. For now, suppose an ideal case, where no magnetic fields other than BJ and the imaging fields are present. In this ideal case, the slopes of the magnetic field ramps do not affect the spin dynamics. Thus, the magnetization after these events can be written as a solution of Eq. (3): Aτ

m ¼ e m0 ¼ Φm0 ;

ð5Þ

where Φ = e Aτ is a rotation matrix that describes the precession around BJ for a duration τ. The problem of finding BJ by an MRI experiment is now transformed into finding the elements Φij of the rotation matrix Φ. With a main field B0 = B0ez applied after the sequence described above, the elements Φ11 and Φ21 corresponding to mx and my are obtained from the real and imaginary parts of the detected signal. Thus, finding all the elements of Φ requires additional sequences. To accomplish this, consider the sequences depicted in Fig. 1 and described by m1 ¼ Φm0 ; m2 ¼ ΦRz ðπ=2Þm0 ; m3 ¼ ΦRy ð−π=2Þm0 ;

ð6Þ and

ð7Þ ð8Þ

where Rz(π/2) and Ry(− π/2) are rotation matrices describing π/2 excitations around z and y axes that can be accomplished by pulsing,

Fig. 1. The sequence diagrams for Z-CDI corresponding to Eqs. (6–8). Three sequences are required to calculate the three components of BJ. Each of these sequences can be divided in four phases. The first phase consists of the spin polarization to establish m0. In the second phase, the spin system is prepared by different π/2 pulses to obtain sensitivity for all components of BJ. The third phase is the zero-field encoding period of a duration τ, corresponding to Φ. In the fourth phase, any spatial encoding sequence can be applied to obtain images. The phase and magnitude of each resulting voxel contain information from which BJ can be uniquely reconstructed.

e.g., the B0 and B1 = B1ey fields, respectively [30]. As discussed above, the sequence described by Eq. (6) provides elements Φ11 and Φ21. Similarly, the Rz(π/2) pulse in Eq. (7) provides a y-directional magnetization before applying Φ so that elements Φ12 and Φ22 can be obtained. Finally, with the sequence described by Eq. (8), elements Φ13 and Φ23 can be measured. Having measured the first two rows of Φ, the elements of the third row can be calculated as the cross product of the vectors formed by the first and second row, as a rotation matrix is always orthogonal. However, in a realistic ULF-MRI measurement, noise from the SQUIDs and their conductive surroundings contaminates the elements of Φ so that the first and second row of Φ are not exactly orthogonal. In this case, it is useful to first apply the Löwdin transformation [31], which orthogonalizes and normalizes the two vectors with a result that is the closest to the original vectors in a least-squares sense. After determining Φ, it remains to deduce BJ from this information. First, the rotation angle ϕ of Φ can be calculated from the trace of Φ, 3 X i¼1

Φii ¼ 1 þ 2cosϕ:

ð9Þ

768

P.T. Vesanen et al. / Magnetic Resonance Imaging 32 (2014) 766–770

Second, matrix A can be found from Φ by a matrix logarithm, which for a rotation matrix is given by Aτ ¼

 ϕ  T Φ−Φ ; 2sinϕ

ð10Þ

as can be shown geometrically or by substitution. Finally, BJ can be found by simply picking the right elements of A according to Eq. (4). In practice, also the spin density m0 of the sample, effects of relaxation, and the sensitivity profile of the receiver coil are unknown. Thus, a fourth sequence with no applied BJ is required. Measuring such a calibration image is useful also in the case that non-idealities cause systematic errors in the image phase. Another practical constraint is that ϕ should be constrained to − π ≤ ϕ ≤ π to avoid phase wrapping. This can be accomplished by choosing τ accordingly. This constraint can be relaxed if phase unwrapping is carefully considered in the reconstruction. To evaluate the performance of Z-CDI, the sequences depicted in Fig. 1 were numerically simulated with a full time-domain Blochequation solver using Matlab (Mathworks, Inc., Natick, MA, USA). For simplicity, we decided to simulate the spatial encoding by a 3D gradient-echo sequence. We assumed ideal homogeneous fields, i.e., Bp ¼ Bp ex , B0 = B0ez, and B1 = B1ey. The gradient fields were purely z-directional and each of them varied linearly along one of the Cartesian coordinate axes. Frequency encoding was applied in the z direction. The field amplitudes were Bp = 100 mT and B0 = 50 μT. In the phase-encoded x and y directions, 20 k-space lines were acquired. The time to echo was 50 ms and the polarizing time was 1 s. The digital sample in our simulations shown in Fig. 2 was a volume of 20 × 20 × 20 cubic voxels positioned at the origin. The voxel side length was 4 mm. We assumed the spin density of water and relaxation times of T1 = T2 = 200 ms uniformly across the phantom. For simplicity, the current density J was assumed to be that of three long 24-mm-thick tubes with circular cross sections passing through the origin along the x, y, and z axes, each carrying a current of 12 mA, which corresponds to a constant current density of 27 A/m 2. This is an approximation of a conductor geometry where the tubes are filled with conducting liquid and their surfaces are made of insulating material. The magnetic field BJ generated by J was

Fig. 2. Simulation geometry. The orange tubes represent the region of the simulated current flow. Simulation results in Figs. 3 and 4 are plotted along the planes indicated in green.

calculated using the Biot–Savart law and varied from 0 nT at the center of the phantom to around 420 nT near the surfaces of the tubes. The simulated signals were sampled with a receiver having homogeneous sensitivity to the x component of the precessing magnetization. To avoid committing the inverse crime, we added small random deviations in the positions of the dipoles representing the voxels: the dipoles were deviated in all three Cartesian directions with the maximum deviation being 5% of the voxel side length. To simulate a realistic ULF-MRI measurement, we added white Gaussian noise to the simulated signals. The noise standard deviation σ was chosen such that the signal-to-noise ratio (SNR) of the simulated data matched that of an experimental ULF-MRI data set [14], which was acquired using a resolution similar to that in our simulated data. Specifically, using the definition

SNR ¼

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u∫jsðt Þj2 dt t

1 ¼ σ ∫σ 2 dt

sffiffiffiffiffiffiffi sT s ; Ns

ð11Þ

where s(t) is the signal as a function of time, s a vector containing the sampled data points, and Ns the number of acquired samples, the SNR of the experimental data set in Ref. [14] was 8.3. Using the same definition, the noise standard deviation σ for the simulation was adjusted to obtain equal SNR for the simulated data set. The computational time of reconstructing the current density from the simulated MRI data on a quad-core 2.83-GHz desktop computer was only a fraction of a second. 3. Results and discussion Fig. 3 shows the simulated magnitude and phase images of the sample obtained using the sequences described by Eqs. (6–8). In the images produced by the sequences in Eqs. (6) and (7), the magnetization before the zero-field encoding period is prepared to be x- and y-directional, respectively. In these cases, the z component of BJ gives rise to the phase variations in the images. On the other hand, in images produced by the sequence in Eq. (8), the magnetization before the zero-field encoding period is prepared to be z-directional. In this case, the intensity values obtained from the regions with weak BJ are small, as indicated by the corresponding magnitude image. Here, x and y components of BJ give rise to the observed phase and intensity variations. Fig. 4 shows the current-density reconstruction from the image data shown in Fig. 3. It can be seen that the correct current density is reconstructed with no visible bias in the current direction or amplitude. Furthermore, noise is evenly distributed, indicating the stability of the method. In fact, numerical simulations using sequences described by Eqs. (6–8) suggest that, for Gaussian noise added to the data, the noise is Gaussian also in the reconstructed components of BJ. The sensitivity to a weak current density pattern can be increased by using a long zero-field evolution period τ. However, as τ approaches the relaxation times, the SNR starts to decrease because of magnetization dephasing. The value of τ can be optimized to maximize the SNR when the relaxation times of the sample are known. The measurement time of the proposed method is equal to that of conventional CDI [1–3], when equal SNR is assumed and the calibration phase is ignored. However, because of the ULF setting required to apply Z-CDI, SNR of Z-CDI measurements is inevitably lower than that of conventional CDI measurements. For an SNR comparison between ULF MRI and conventional MRI, see Ref. [32]. In the presented simulation, it was assumed that no magnetic fields other than BJ and the imaging fields were present. However, in practice a fraction of Earth’s magnetic field Brem exists even inside magnetically shielded rooms. Furthermore, pulsing of the polarizing

P.T. Vesanen et al. / Magnetic Resonance Imaging 32 (2014) 766–770

769

Fig. 3. Magnitude and phase images of the sample in the z = 20 mm plane obtained using the sequences described by Eqs. (6–8). The data are normalized such that in a noise-free case, the maximum magnitude of the signal over the whole sample volume corresponds to unity.

field induces eddy currents in the conductive structures of the scanner and shielded room, creating unwanted magnetic field transients during the measurement. The presence of unwanted magnetic fields affect the measurement in several ways. First, the ramps of all the pulses have to be sufficiently steep to reach nonadiabatic manipulation of the spin system. Reaching the nonadiabatic condition for the ramps of the polarizing field is technically challenging but has, nevertheless, been demonstrated in the ULFMRI literature [13]. Second, the presented formalism allows only the measurement of the total external field, i.e., the superposition of BJ and the unwanted Brem. However, it is possible to distinguish BJ from

Brem by first imaging only Brem without applying BJ at all. A subsequent measurement of the total field Brem + BJ allows the recovery of BJ by subtracting the results of the two measurements. During this study, we were unable to experimentally verify the proposed method because of hardware limitations in rapidly switching the low-noise currents in the B0 and frequency-encoding gradient coils. However, ULF-MRI measurements with sufficiently fast field switching, using second-order-gradiometric sensors, have previously been reported [13]. In principle, Z-CDI can be used not only with ULF MRI but also with any MRI system capable of B0 field cycling. For example, existing pre-polarization-based systems using

Fig. 4. Current density reconstruction from the data in Fig. 3 at the planes indicated in Fig. 2. The dashed lines indicate the regions of the current flow; inside these regions, the normalized current density should be 1, which corresponds to 27 A/m2, and outside them 0.

770

P.T. Vesanen et al. / Magnetic Resonance Imaging 32 (2014) 766–770

induction-coil detection [33–36] could be suitable for Z-CDI with relatively minor changes in the hardware. If Z-CDI is applied using a device without a magnetically shielded room, additional coils for compensating the magnetic field of the Earth are needed. The low magnetic fields used in ULF MRI allow flexibility in the coil setup. In addition to having B0 in all three directions, one could imagine a system with all the nine different Cartesian linear gradients of the magnetic field. With such a system, current-density imaging could be performed exactly as is done conventionally [1–3]: to image the x component of BJ, the B0 and gradient fields in the x direction would be turned on. Similarly, the corresponding y- and z-directional B0 and encoding fields would subsequently be used to measure the respective components of BJ. The A-CDI method presented in the accompanying paper [29] differs from Z-CDI in several ways. First, the physical principles of the methods are different; while Z-CDI is based on zero-field encoding, A-CDI uses adiabatic ramps to encode the direction of the magnetization. Second, the reconstruction of Z-CDI is calculated voxel-wise, while that of A-CDI is calculated globally based on Maxwell’s Equations. 3. Conclusion We have presented a method for imaging magnetic fields and current densities using ultra-low-field MRI. Specifically, the magnetization is allowed to dephase in the unknown magnetic field without a B0 field present. By suitably preparing the magnetization, the method gives complete 3D magnetic field and current-density information without a need to rotate the object. The method may be useful, e.g., for conductivity imaging. Acknowledgment The research leading to these results has received funding from the European Community’s Seventh Framework Programme (FP7/ 2007–2013) under grant agreement No. 200859, from the Emil Aaltonen Foundation, the Academy of Finland, the International Doctoral Programme in Biomedical Engineering and Medical Physics (iBioMEP), and the Finnish Cultural Foundation. References [1] Joy MLG, Scott GC, Henkelman RM. In-vivo detection of applied electric currents by magnetic resonance imaging. Magn Reson Imaging 1989;7:89–94. [2] Pešikan P, Joy MLG, Scott GC, Henkelman RM. Two-dimensional current density imaging. IEEE Trans Instrum Meas 1990;39:1048–53. [3] Scott GC, Joy MLG, Armstrong RL, Henkelman RM. Measurement of nonuniform current density by magnetic resonance. IEEE Trans Med Imaging 1991;10: 362–74. [4] Scott GC, Joy MLG, Armstrong RL, Henkelman RM. Electromagnetic considerations for RF current density imaging. IEEE Trans Med Imaging 1995;14:515–24. [5] Scott GC, Joy MLG, Armstrong RL, Henkelman RM. Rotating frame RF current density imaging. Magn Reson Med 1995;33:355–69. [6] Ider YZ, Muftuler LT. Measurement of AC magnetic field distribution using magnetic resonance imaging. IEEE Trans Med Imaging 1997;16:617–22. [7] Mikac U, Demšar F, Beravs K, Serša I. Magnetic resonance imaging of alternating electric currents. Magn Reson Imaging 2001;19:845–56. [8] Wang D, DeMonte TP, Ma W, Joy MLG, Nachman AI. Multislice radio-frequency current density imaging. IEEE Trans Med Imaging 2009;28:1083–92. [9] Voigt T, Katscher U, Doessel O. Quantitative conductivity and permittivity imaging of the human brain using electric properties tomography. Magn Reson Med 2011;66:456–66.

[10] Oh SH, Lee BI, Woo EJ, Lee SY, Cho MH, Kwon O, et al. Conductivity and current density image reconstruction using harmonic Bz algorithm in magnetic resonance electrical impedance tomography. Phys Med Biol 2003;48:3101–16. [11] Seo JK, Yoon J-R, Woo EJ, Kwon O. Reconstruction of conductivity and current density images using only one component of magnetic field measurements. IEEE Trans Biomed Eng 2003;50:1121–4. [12] McDermott R, Lee SK, ten Haken B, Trabesinger AH, Pines A, Clarke J. Microtesla MRI with a superconducting quantum interference device. Proc Natl Acad Sci USA 2004;101:7857–61. [13] Zotev VS, Matlashov AN, Volegov PL, Savukov IM, Espy MA, Mosher JC, et al. Microtesla MRI of the human brain combined with MEG. J Magn Reson 2008;194:115–20. [14] Vesanen PT, Nieminen JO, Zevenhoven KCJ, Dabek J, Parkkonen LT, Zhdanov AV, et al. Hybrid ultra-low-field MRI and MEG based on a commercial whole-head neuromagnetometer. Magn Reson Med 2013;69:1795–804. [15] Packard M, Varian R. Free nuclear induction in the Earth’s magnetic field. Phys Rev 1954;93:941. [16] Lee J-Y. A reconstruction formula and uniqueness of conductivity in MREIT using two internal current distributions. Inverse Probl 2004;20:847–58. [17] Nachman A, Tamasan A, Timonov A. Conductivity imaging with a single measurement of boundary and interior data. Inverse Probl 2007;23:2551–63. [18] Woo EJ, Seo JK. Magnetic resonance electrical impedance tomography (MREIT) for high-resolution conductivity imaging. Physiol Meas 2008;29:R1-26. [19] Hasanov KF, Ma AW, Nachman AI, Joy MLG. Current density impedance imaging. IEEE Trans Med Imaging 2004;27:1301–9. [20] Bodurka J, Jesmanowicz A, Hyde JS, Xu H, Estkowski L, Li S-J. Current-induced magnetic resonance phase imaging. J Magn Reson 1999;137:265–71. [21] Bodurka J, Bandettini PA. Toward direct mapping of neuronal activity: MRI detection of ultraweak, transient magnetic field changes. Magn Reson Med 2002;47:1052–8. [22] Kraus Jr RH, Espy MA, Volegov PL, Matlachov AN, Mosher JC, Urbaitis AV, et al. Toward SQUID-based direct measurement of neural currents by nuclear magnetic resonance. IEEE Trans Appl Supercond 2007;17:854–7. [23] Kraus Jr RH, Volegov P, Matlachov A, Espy M. Toward direct neural current imaging by resonant mechanisms at ultra-low field. NeuroImage 2008;39: 310–7. [24] Burghoff M, Albrecht HH, Hartwig S, Hilschenz I, Körber R, Sander Thömmes T, et al. SQUID system for MEG and low field magnetic resonance. Metrol Meas Syst 2009;16:371–5. [25] Burghoff M, Hartwig S, Cassará A, Trahms L. Approaches to detect neuronal currents using a DC-mechanism by means of low field magnetic resonance. In: World Congress on Medical Physics and Biomedical Engineering, September 7 – 12, 2009, Munich, Germany, pp. 31–33, 2009 [26] Burghoff M, Albrecht HH, Hartwig S, Hilschenz I, Körber R, Höfner N, et al. On the feasibility of neurocurrent imaging by low-field nuclear magnetic resonance. Appl Phys Lett 2010;96:233701. [27] Höfner N, Albrecht H-H, Cassará AM, Curio G, Hartwig S, Haueisen J, et al. Are brain currents detectable by means of low-field NMR? A phantom study. Magn Reson Imaging 2011;29:1365–73. [28] Körber R, Nieminen JO, Höfner N, Jazbinšek V, Scheer H-J, Kim K, et al. An advanced phantom study assessing the feasibility of neuronal current imaging by ultra-low-field NMR. J Magn Reson 2013;237:182–90. [29] Nieminen JO, Zevenhoven KCJ, Vesanen PT, Hsu Y-C, Ilmoniemi RJ. Currentdensity imaging using ultra-low-field MRI with adiabatic pulses. Magn Reson Imaging 2014;32:54–9. [30] Friedman LJ, Wennberg AKM, Ytterboe SN, Bozler HM. Direct detection of lowfrequency NMR using a dc SQUID. Rev Sci Instrum 1986;57:410–3. [31] Löwdin P-O. On the non-orthogonality problem connected with the use of atomic wave functions in the theory of molecules and crystals. J Chem Phys 1950;18:365–75. [32] Myers W, Slichter D, Hatridge M, Busch S, Mössle M, McDermott R, et al. Calculated signal-to-noise ratio of MRI detected with SQUIDs and Faraday detectors in fields from 10 μT to 1.5 T. J Magn Reson 2007;186:182–92. [33] Venook RD, Matter NI, Ramachandran M, Ungersma SE, Gold GE, Giori NJ, et al. Prepolarized magnetic resonance imaging around metal orthopedic implants. Magn Reson Med 2006;56:177–86. [34] Matter NI, Scott GC, Grafendorfer T, Macovski A, Conolly SM. Rapid polarizing field cycling in magnetic resonance imaging. IEEE Trans Med Imaging 2006;25: 84–93. [35] Gilbert KM, Handler WB, Scholl TJ, Odegaard JW, Chronik BA. Design of fieldcycled magnetic resonance systems for small animal imaging. Phys Med Biol 2006;51:2825–41. [36] Halse ME, Coy A, Dykstra R, Eccles C, Hunter M, Ward R, et al. A practical and flexible implementation of 3D MRI in the Earth’s magnetic field. J Magn Reson 2006;182:75–83.