PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Phys Med Biol. Author manuscript; available in PMC 2010 October 22.
Published in final edited form as:
PMCID: PMC2962550
NIHMSID: NIHMS239939

Limitations of Calculating Field Distributions and Magnetic Susceptibilities in MRI using a Fourier Based Method

Abstract

A discrete Fourier based method for calculating field distributions and local magnetic susceptibility in MRI is carefully studied. Simulations suggest that the method based on discrete Greens functions in both 2D and 3D spaces has less error than the method based on continuous Greens functions. The 2D field calculations require the correction of the “Lorentz disk, which is similar to the Lorentz sphere term in the 3D case. A standard least squares fit is proposed for the extraction of susceptibility for a single object from MR images. Simulations and a phantom study confirm both the discrete method and the feasibility of the least squares fit approach. Finding accurate susceptibility values of local structures in the brain from MR images may be possible with this approach in the future.

Keywords: Fourier transformation, Greens function, Lorentz disk, Lorentz sphere, MRI, susceptibility

1 Introduction

The ability to measure the magnetic susceptibility from magnetic resonance imaging (MRI) data could prove very useful in a number of in vivo applications. The basal ganglia appear to increase their iron content in diseases such as Parkinsons, Huntingtons and multiple sclerosis [1, 2]. Knowing the magnetic susceptibility makes it possible to remove unwanted phase effects in gradient echo images in methods such as susceptibility weighted imaging [3]. Further, using the geometry and the susceptibility makes it possible to remove the phase (field) effects caused by the presence of air/tissue interfaces and even to predict and correct geometric distortion in echo planar images [4]. Recently a Fourier based method has been introduced to find local magnetic fields [4, 5, 6, 7] and it is similar to the methods presented in [8, 9]. If sufficient accuracy can be achieved with this method, it can be easily implemented in the form of a software tool and is computationally fast.

Initial investigations [4, 6, 7] suggest that the method used in forward calculations of the magnetic field seems promising with tolerable errors if the field-of-view (FOV) is more than 50% of the object size (i.e., length, width, or height). In this paper, we systematically study the source of the uncertainty from the method itself. We present an improved method based on the use of discrete Greens functions that can be applied on both 2D and 3D objects. We demonstrate the field distributions from different Greens functions for a cylinder and a sphere. Through a least squares fit approach, we quantify the susceptibilities from simulations of a spherical shell and a water phantom. In our approach, we use an iterative algorithm for the quantification of susceptibility in the presence of noise and aliasing in MRI phase maps.

2 Theory and Methods

2.1 3D induced field distribution and its Greens function

The induced magnetic field due to a magnetization distribution, M (r), is given by (see Eq. (5.64) of [10] or Eq. (3) of [7])

B(r)=μ04πVd3r{3M(r)·(rr)rr5(rr)M(r)rr3}
(1)

Equation (1) indicates that the induced field distribution can be expressed as a convolution between the magnetization distribution and the Greens function. As Deville et al. [5] pointed out, Eq. (1) may be easily calculated in the Fourier domain (k-space domain) and then Fourier transformed back to the spatial domain. In addition, as discussed in [7] or implied by Eq. (5.64) of [10], the Lorentz sphere correction has been included in Eq. (1).

If the main field is along the z-axis, then only the z-component of the magnetic field is important in most MRI research. The z-component of the induced magnetic field Bz(r) and its Fourier pair in k-space are

Bz(r)=μ04πVd3r{3Mz(r)(zz)2rr5Mz(r)rr3}μ0Vd3rMz(r)Gz,3D(rr)
(2)

Bz(k)F{Bz(r)}=μ0Mz(k)Gc,3D(k)
(3)

where Mz (k) is the Fourier transformation of Mz (r) and the Greens function Gz,3D(r) is defined as

Gz,3D(r)14π·3z2r2r514π·3cos2θ3D1r3
(4)

where θ3D is the azimuthal angle in the spherical coordinate system and r2 [equivalent] x2 + y2 + z2. The Fourier transformation of Gz,3D(r) is

Gc,3D(k)F{Gz,3D(r)}=13kz2k2
(5)

where k2kx2+ky2+kz2 and the notation Gc,3D denotes the continuous Fourier transformation of the 3D Greens function. Equation (5) is valid when k is not equal to zero or infinity. When k = 0, the derivation of Gc,3D(k) shown in Eq. (30) (in the Appendix) indicates that Gc,3D(0) = 0 as long as the integral is calculated over the entire space larger than a small sphere with radius ε. In fact, Eq. (30) at k = 0 is equivalent to the calculation of the Lorentz sphere term [10]. Because the Lorentz sphere correction has been included in Eq. (1), it is consistent to assign zero to Gc,3D(0) for future numerical calculations [4, 7].

Although Eq. (5) may be derived from other methods, as pointed out by Koch et al. [4], derivations involving derivatives require special attention at the boundary of M (r). Therefore, we believe the derivation through Fourier transformation shown in [7] is the most appropriate approach. Our derivation is also shown in Appendix A.

2.2 2D induced field distributions and its Greens function

For an infinitely long object, with a uniform cross section along its length, its induced magnetic field can be calculated from a 2D Greens function. Similar to the electric field distribution due to an electric dipole presented in [10], the magnetic field (H) on a 2D plane due to the scalar potential of a dipole m may be written as (also see Eq. (13) in [11])

H(x)=12π(m·xρ2)=12π(2x(m·x)ρ4mρ2)
(6)

whose spatial distribution agrees with that shown in [12] (p. 751).

Recall that the B field is equal to μ0(H + M). After replacing the magnetic dipole moment by its magnetization distribution, the induced magnetic field on a 2D plane becomes

B(x)=μ02πd2x{2M(x)·(xx)xx4(xx)M(x)xx2}+μ0M(x)
(7)

It is obvious that M (x) in the last term can be replaced by ∫d2xM (x)δ(2)(xx) for all later discussions.

However, neither Eq. (6) nor Eq. (7) has been properly defined at the origin, where the magnetic dipole moment is located. This term has been calculated in Appendix B. Similar to Eq. (5.64) of [10], the corrected magnetic field embedded in the 3D space is now the combination of Eq. (6) and Eq. (31) from Appendix B.

B(x)μ02πd2x{2M(x)·(xx)xx4(xx)M(x)xx2}+(112sin2θ)μ0M(x)
(8)

where θ is the angle between the main field direction and the axis of the infinitely long object. In the derivation of Eq. (8), we have automatically assumed that the integrand is zero when x = x. Furthermore, as we are interested in the magnetic field measured in the 3D space, the Lorentz sphere correction needs to be included in Eq. (8). The magnetic field now becomes

B(x)μ02πd2x{2M(x)·(xx)xx4(xx)M(x)xx2}+(1312sin2θ)μ0M(x)
(9)

If the main field is along the ω-direction, only the ω-component of the induced magnetic field is important in our discussions here. Hence,

Bω(x)=μ0sin2θ2πd2x{2(zz)2Mω(x)xx4Mω(x)xx2}+(1312sin2θ)μ0Mω(x)μ0d2xMω(x)Gω,2D(xx)
(10)

where

Gω,2D(x)sin2θ2π·z2y2(z2+y2)2+3cos2θ16δ(2)(x)sin2θ2π·2cos2φ1ρ2+3cos2θ16δ(2)(x)
(11)

where ρ2 [equivalent] z2 + y2 and [var phi] is the polar angle between z and y axes. Here, the yz plane is chosen to be perpendicular to the axis of the infinitely long object and the plane contains the normal cross section of the object. We have properly inserted the factor sin2 θ in Eq. (10) as discussed in Appendix B. Again, in the derivation of Eq. (11), the term, (2 cos2 [var phi] − 1)2, is already assumed to be zero at ρ = 0. The Fourier transformation of Gω, 2D(x), as derived in Appendix C, is

Gc,2D(k)F{Gω,2D(x)}=sin2θ2·ky2kz2ky2+kz2+16(3cos2θ1)=13sin2θkz2ky2+kz2
(12)

Similar to the 3D case, Eq. (12) is valid only when ky2+kz2 is not equal to zero or infinity. When k = 0, Gc,2D(k = 0) = (3 cos2 θ − 1)/6, as the integral of the first term in Gω,2D (x) is zero. This is because that term is (assumed) zero at ρ = 0 and the integral of (2 cos2 [var phi] − 1)/ρ2 is also zero as long as the integral is calculated outside a small circle with radius ε. Furthermore, when the infinitely long object is parallel to the x direction in the 3D space with the main field parallel to the z-direction, Eq. (12) is equal to Eq. (5) with kx = 0.

In the rest of the paper, we will consider that the main field is parallel to the z-direction for the 3D case and ω-direction for the 2D case. We will label the induced field and magnetization with subscript z for any general discussions below. However, because our main interest here is the measured magnetic field from MR images, the consideration of the field direction is not important at all.

2.3 Discrete Greens functions in k-space

The continuous Greens functions were derived based on an infinite field of view. However, all MR images have finite fields-of-view and are in a discrete form. In order to obtain consistent results, both a discrete magnetization and a discrete Greens function should be used in k-space. The discrete Greens functions in k-space in the 3D and 2D cases can be numerically calculated from the discrete Fourier transformation of the spatial Greens functions in Eqs. (4) and (11), respectively, with a given matrix size. The spatial Greens functions are discretized and evaluated in units of the image resolutions. Nonetheless, the values at the origin of the Greens functions require some discussion. In the 3D case, because the Lorentz sphere correction is included in Eq. (4), Gz,3D = 0 at r = 0. In the 2D case, the 3D Lorentz sphere is corrected with the presence of the 2D Lorentz disk. Therefore, at x= 0, Gω,2D = (3 cos2 θ − 1)/6.

2.4 Field calculations using the Greens functions as a forward problem

For a non-ferromagnetic object, its magnetization can be expressed in terms of its magnetic susceptibility and an external field B0.

μ0M(r)χ(r)B0
(13)

If an object has a constant susceptibility χ, then the Fourier transformation of its magnetization along the main field direction, Mz (k), is simply the Fourier transformation of the geometry of the object multiplied with χB00. As shown in Eq. (3), μ0Mz (k) should be multiplied with the Greens function in k-space and the total result should be inverse Fourier transformed back to the spatial domain. The induced magnetic field along the main field direction will then be given by

Bz(r)=F1{F{χ(r)}·G(k)}B0
(14)

where G(k) represents either the 2D or 3D discrete Greens function or the continuous Greens functions. We evaluate the 3D and 2D Greens functions and Mz(r) at integer coordinates, i.e., Δx = Δy = Δz = 1 unit. Although one can evaluate these functions at a different spacing between grid points, any choice of spacing has to be consistent between the evaluation of Greens functions and geometry of the object (i.e., Mz (r)).

When the induced field distribution of an infinitely long object is calculated from the 2D Greens functions, the geometry of the object should be taken as the cross section perpendicular to the object axis. The induced field should be calculated based on the coordinate systems on the plane whose normal vector is parallel to the object axis. For example, if an infinitely long cylinder intersects at an angle θ with the main field direction, regardless of the value of θ, the geometry used in calculating the magnetic field distribution should be a disk rather than an ellipse.

For the forward calculations, we have simulated the magnetic fields of a sphere and a cylinder with the continuous and discrete Greens function. As the field decreases when a voxel is away from the object, the field difference between the calculated value (Eq. (14)) and the theoretical value is quoted as error and is presented in percentage. However, when the field inside a sphere is calculated, its value is shown as in theory the magnetic field is zero inside the sphere after the Lorentz sphere correction. In order to minimize the numerical error due to the unsmoothed surface of the object [4], the radius of any simulated object in this paper is at least 16 pixels with a fixed field of view of 256 pixels along each dimension. The ratio of the object diameter to the field of view is used as an indication of how large an object is in the simulations.

2.5 Quantifying magnetic susceptibility from a least squares fit

Because the Greens functions in Eqs. (5) and (12) contain null values, it is not appropriate to determine susceptibility through Eq. (3) by inverting the Greens functions. Instead, for an object with a constant susceptibility χ, the following goodness-of-fit least squares function can be established based on Eq. (14) [13]. A similar concept was briefly mentioned in [6].

f=i=1n(BiχB0giδBi)2=i=1n(φiγTEχB0giδφi)2=i=1nSNRi2(φiγTEχB0gi)2
(15)

where n is the total number of voxels used in the fit, Bi is the measured magnetic field from MR images at each voxel i, gi is the induced field per unit susceptibility and per unit main field calculated based on Eq. (14), δBi is the uncertainty of the measured field, [var phi]i [equivalent]γTE Bi is the phase value at each voxel, γ is the gyromagnetic ratio, TE is the echo time, and SNRi is the signal-to-noise ratio at each voxel from the magnitude image. Note that gi depends on the geometry of the object, field of view, and the Greens function.

By minimizing function f in Eq. (15) with respect to χ, i.e., [partial differential]f/[partial differential]χ = 0, we obtain

χB0=i=1nBigiSNRi2i=1ngi2SNRi2
(16)

The uncertainty of χ can be found through the error propagation method [13].

δχ=1γTEB0i=1ngi2SNRi2
(17)

Equations (16) and (17) show that the uncertainty of susceptibility measurement can be reduced when the SNR is increased.

It is important to note that Eqs. (16) and (17) will fail if phase aliasing is not properly unwrapped. The phase aliasing may be removed by complex dividing two phase images acquired at two different echo times with a short difference, ΔTE, between the echo times. However, due to rapid field change within a voxel, usually phase aliasing would exist at voxels with low SNRs. For this reason we exclude those voxels from our analyses. Furthermore, voxels outside imaged objects containing no sufficient SNR in the magnitude images are also removed. However, it can be seen from Eqs. (16) and (17) that voxels with low SNR used in the fitting analysis will not significantly affect the susceptibility quantification or its uncertainty, provided that enough voxels with sufficient SNR has been used in the analysis.

When acquiring images, a uniform phase [var phi]0 can exist due to central frequency adjustment by the scanner or rf excitation in the pulse sequence. Therefore, even for imaging one single object, practically, Eq. (15) should be modified to read

f=i=1nSNRi2(φiφ0γTEχB0gi)2
(18)

Similar to the above derivation, both values of χ and [var phi]0 can be determined by minimizing the function f through the standard least squares fit method [13]. The solutions and their associated standard deviations (i.e., uncertainties) derived through the error propagation method are shown in [13].

2.6 Quantification of susceptibility from simulations

All simulations use MATLAB (The Mathworks Inc., Natick, MA) on a Windows XP platform and with a 1.54 GHz AMD Turion Mobile processor and with 2 GB RAM memory. Unless otherwise mentioned, in all simulations, the input main field is 1 Tesla and the input susceptibility of an object is 1 ppm in SI Units. The simulated objects are also assumed to be surrounded by the vacuum. Equation (14) is used to simulate magnetic fields based on discrete and continuous Greens functions as well as object geometries, which are spheres and cylinders at 14 different radii from 16 pixels to 120 pixels with an incremental size of 8 pixels. The field of view is always fixed at 256 pixels along each dimension.

With no white noise included in the simulations, Eq. (16) can be used to quantify the susceptibility with a constant SNR in all pixels. The difference between the calculated and input susceptibility basically indicates the accuracy of the methods (or Greens functions) themselves. Only pixels outside the objects are used in the quantification of susceptibilities, as the actual field inside a sphere is zero. One set of quantifications uses all pixels outside the object. The other set uses pixels from a spherical shell or an annular ring, whose locations of pixels satisfy R + 1 < r ≤ R + 5 where R is the radius of the object.

In a simulation more close to practical imaging situations, we simulate magnitude and phase images with white noise added [12]. An SNR of 5:1 is assigned and the object is a spherical shell defined by two concentric spheres in this simulation. The inner sphere has a radius of 16 pixels and the outer sphere has a radius of 96 pixels. The field of view is again 256 pixels. The input susceptibility is 10 ppm in the shell and zero outside the shell and inside the inner sphere. The main field strength is 1 T and the echo time is 4 ms such that some phase aliasing occurs around the surface of the inner sphere. A constant background phase value [var phi]0 = 1 radian is also given as an input parameter in the simulation.

Only the discrete Greens function is used in the quantification of susceptibility from this shell simulation. Through Eq. (18), the least squares fit approach allows us to determine both the susceptibility and the constant background phase [var phi]0. The actual phase value due to the object at the i-th voxel is [var phi]i[var phi]0 where [var phi]i is the simulated phase value. In order to remove the aliased phase points, the following criterion is used:

φiφ0γTEχB0gip/SNRi
(19)

where SNRi is the signal to noise ratio of the i-th voxel in the magnitude images and its inverse is the standard deviation of the phase in the same voxel. We choose p as an integer and it is usually one or two such that 68% or 95% of the non-aliased phase values are included in the susceptibility quantification. Due to this criterion which helps to remove phase aliased voxels, an iterative procedure is used for the quantification as both χ and [var phi]0 are unknowns. A flow chart of the iterative procedure is shown in Fig. 1. When the susceptibility value changes less than 0.1% (i.e., ε = 0.001 in Fig. 1), the iterative procedure will stop. In our program, we often choose the initial values of [var phi]0 and χ as 0 radian and 7 ppm, respectively, compared to the input values of 1 radian and 10 ppm, respectively.

Figure 1
A flow chart showing the iterative procedure of solving the susceptibility through the least squares fit method.

2.7 Phantom study

One object with a somewhat complicated geometry was considered for our phantom study. As shown in Fig. 2, a CD case-top made of polypropylene was filled with water and was scanned on a 1.5 T Siemens Sonata MRI system. The CD case-top had a wall thickness of 1.3 mm, height of 180 mm, and a diameter of 120 mm. It had a hollow cylindrical post in the middle of one of its flat ends. The post had a varying diameter along its length, from 12 mm at its base to 9 mm at its apex. Prior to the scans, a standard spherical phantom (with a diameter of 170 mm) containing the NiSO4 solution for the system quality control was used for shimming. The full width at half maximum frequency was reduced to 6 Hz at the end of the shimming process. All the first and second order shim current values were noted and were subsequently used for imaging the water phantom.

Figure 2
Images of the water phantom. (a) Transverse magnitude image at TE = 6.58 ms. (b) Its associated phase image. (c) The phase image after applying the complex division method. Voxels inside the dashed box shown in (a) were used for susceptibility quantification. ...

The water phantom was imaged transversely twice with echo times 6.58 ms and 9.58 ms of a 3D gradient echo sequence and with the following parameters: TR 15 ms, flip angle 6°, voxel size 0.78 mm × 0.78 mm × 0.78 mm, matrix size 256 × 256 × 192 and read bandwidth 610 Hz per pixel. The images were later placed in a matrix of 256 × 256 × 256 with the additional slices filled with zeros. The data acquired at the echo time 6.58 ms were complex divided into data acquired at the echo time 9.58 ms. A set of images at an effective echo time 3 ms was created and its phase values were used for the susceptibility quantification. The SNR of the complex divided images was roughly 5:1. The complex images acquired from the echo time 6.58 ms were used to obtain the geometry of the water phantom with a complex thresholding algorithm [14]. The iterative procedure shown in Fig. 1 was again used for the susceptibility quantification. The initial values of [var phi]0 and χ used in the iterative program were 0 radian and 6 ppm, respectively.

3 Results

3.1 Spherical case: analytical solutions

We first assume a sphere with radius a and susceptibility χ embedded in vacuum. Using Eq. (29), its magnetization in k-space is

μ0Mz(k)=d3rχB0Θ(ar)ei2πk·r=4πχB00adrr2j0(2πkr)=2χB0a2kj1(2πka)
(20)

where Θ is the step function, and j0 and j1 are spherical Bessel functions. Based on Eqs. (5) and (29), the induced magnetic field in the spatial domain is

Bz(r)=F1{μ0Mz(k)Gc,3D(k)}=8π3χB0a2(3cos2θ3D1)0dkkj1(2πka)j2(2πkr)={13χB0(3cos2θ3D1)(ar)3whenr>a0whenr<a
(21)

The last step requires the identity from [15]. The solution agrees with the result shown in [12] (see p. 755).

If two concentric spheres with radii Ri and Re (Re > Ri) and constant susceptibilities χi and χe are embedded in vacuum, then the overall magnetization in the spatial domain can be written as

μ0Mz(r)=(χiχe)B0Θ(Rir)+χeB0Θ(Rer)
(22)

Because the Fourier transformation is a linear operation, based on the above derivations and results from Eq. (21), the induced magnetic field due to the two concentric spheres can be easily written down

Bz(r)={0whenr<Ri13(χiχe)B0(3cos2θ3D1)(Rir)3whenRi<r<Re13[(χiχe)(Rir)3+χe(Rer)3]B0(3cos2θ3D1)whenRe<r
(23)

This result agrees with solutions from solving the Laplace equation in the magnetostatics case when higher order terms of χi and χe are neglected.

Another example is to consider a sphere with radius Ri and constant susceptibility χi and an infinitely long cylinder with radius Re (Re > Ri) and constant susceptibility χe parallel to the main field along the z-direction. Outside the cylinder it is vacuum. The origin of the sphere is located on the axis of the cylinder. The magnetization of the entire system is

μ0Mz(r)=(χiχe)B0Θ(Rir)+χeB0Θ(Rex2+y2)
(24)

Following the derivations of Eq. (21) and with the results shown in the next subsection, the induced magnetic field can be readily written down

Bz(r)={13χeB0whenr<Ri13(χiχe)B0(3cos2θ3D1)(Rir)3+13χeB0whenRi<r<Re13(χiχe)B0(3cos2θ3D1)(Rir)3whenRe<r
(25)

Even if radius Re is increased to a size that is much larger than Ri in both Eqs. (23) and (25), it is clear that a uniform induced field, χeB0/3, only exists in Eq. (25) within a region less than radius Re. Thus, when calculating the magnetic field distribution, one needs to include the geometry of the largest object in the magnetic field. Outside the largest object it should be vacuum or air with essentially only the main field. One cannot simply adjust the background field to be χeB0/3 in solving the field distributions, as incorrectly mentioned in [7].

3.2 Infinitely long cylindrical case: analytical solution

In this subsection, we assume an infinitely long cylinder with radius a and susceptibility χ at an angle θ with the main field B0 which is taken to be parallel to the ω-direction. As discussed in Sec. 2.4, the geometry of the object should be treated as a disk and the calculations should be performed in its own cross sectional coordinate system.

μ0Mω(x)=χB0Θ(az2+y2)
(26)

With the help of Eq. (32), its 2D Fourier transformation is

μ0Mω(k)=F{μ0Mw(x)}=2πχB00adρρJ0(2πkz2+ky2ρ)=χB0aJ1(2πkz2+ky2a)kz2+ky2
(27)

In combination with the 2D continuous Greens function shown in Eq. (12), the induced magnetic field in the spatial domain is

Bω(x)=F1{μ0Mω(k)Gc,2D(k)}=χB0asin2θ2·z2y2z2+y20dkJ1(ka)J2(kz2+y2)+χB06(3cos2θ1)a0dkJ1(ka)J0(kz2+y2)={χB0sin2θa22(z2+y2)·z2y2z2+y2outsidethecylinderχB06(3cos2θ1)insidethecylinder
(28)

The last step requires an identity from [15]. The coordinates y and z are defined by the plane perpendicular to the axis of the cylinder. This solution agrees with the result from [12] (p. 755). The solution of the θ = π/2 case was also derived in [4].

3.3 Simulations

Figure 3 shows the magnetic field profiles for a sphere and a cylinder based on the theoretical predictions for the continuous Greens functions (labeled as Gc) and the discrete Greens functions (labeled as Gd). When judging carefully, one can see the ringing effect of the field profiles across the entire field of view calculated from the continuous Greens functions. Actually this is the blurring effect due to the singularity (or discontinuity) of the continuous Greens functions at the origin (r = 0). This blurring effect is not observed with the use of the discrete Greens functions as every step involved in such a method is under the discrete Fourier transformation. No truncation or finite sampling is introduced in k-space when the discrete Greens functions are applied. Compared to the continuous Greens function, it is also clear that the fields calculated from the discrete Greens function agree well with the theoretical values especially at distances further away from the object.

Figure 3
Magnetic field profiles of (a) a sphere and (b) a cylinder along the main field direction. The cylinder is perpendicular to the main field. (c) The difference of fields in percentage between the theoretical value from a sphere and the calculated values ...

Another way to assess the accuracy due to different Greens function is to compare the actual predictions with the theoretically known field values (Table 1). In this comparison, the FOV is fixed but the diameter of either the spherical or cylindrical object is changed from 32 pixels by steps of 32 pixels. The pixels at the corners of the FOV whose locations are more than 128 pixels from the center of the FOV are excluded in Table 1. Only the numbers of pixels outside the object that lie within a given uncertainty are listed. The results from Table 1 clearly show that the field calculations from the discrete Greens function are in better agreement with the theoretical values than the calculations from the continuous Greens function. Further, the true magnetic field should be zero inside the sphere (Eq. (21)). The root-mean-squared field inside a sphere as a function of its diameter is shown in Fig. 4. It is again seen that the results based on the discrete Greens function are better than the results from the continuous Greens function.

Figure 4
The root-mean-squared magnetic field divided by the main field inside a sphere as a function of the object size. The horizontal axis is displayed as the ratio of the object diameter to the field of view. The field of view is fixed at 256 pixels. Voxels ...
Table 1
This table lists the number of pixels outside a given object within a specified uncertainty of the magnetic field distribution. Voxels that are more than 128 pixels away from the center of the field of view are not counted in this table. The first column ...

In the quantification of susceptibility, both Fig. 5 and Fig. 6 show that the method using the discrete Greens function is better than the method using the continuous Greens function. When all the pixels outside the object and the discrete Greens function are used for susceptibility quantification, no significant error is produced by the method itself if the diameter of the sphere or cylinder studied is less than 40% of the field of view (see Fig. 5). If only certain pixels in spherical shells or annular rings are used, as shown in Fig. 6, the method using the discrete Greens functions produces no significant error for objects with diameters up to almost 60% of the field of view. The results shown in Table 1, Fig. 3, Fig. 5, and Fig. 6 imply that voxels close to the object are more important for an accurate quantification with the Fourier based method than voxels away from the object. This is consistent with the fact that fields in voxels close to the edge of the field of view are prone to errors [7].

Figure 5
Susceptibility error in percentage versus the ratio of the diameter of the object to the field of view. The notations Gc and Gd represent the continuous and discrete Greens functions, respectively. All pixels outside the objects are used in the quantifications. ...
Figure 6
Susceptibility error in percentage versus the ratio of the diameter of the object to the field of view. This figure is similar to Fig. 5 but only certain pixels outside each object are used in the quantifications, as described in the text. (a) for spheres ...

In the simulation of the spherical shell with noise included, several results are observed. First, if we choose all voxels inside the shell for the susceptibility quantification, then we end up with a susceptibility that agrees with the result shown in Fig. 5. Thus, for a more accurate quantification, we only use the central 1283 voxels of the shell for the susceptibility quantification. Second, when the initial choice of the susceptibility is between 7 and 10 ppm, the final answers from the iterative procedure do not seem to be affected. Third, whether one or two standard deviations (i.e., p = 1 or 2) are used in Eq. (19) does not seem to affect the final answers either. The results are shown in Table 2a. As the chi-square per number of points in each calculation is close to one, each result is from a good fit. Although in general, the results in Table 2a are in very good agreement with the input values, the differences between the calculated values and the input values for p = 1 and 2 are due to the Fourier based method itself (or the Greens function).

Table 2
This table lists the results from (a) the simulated shell and (b) the phantom studies. The first column lists the number of standard deviations used in the computer program. The second column lists the number of iterations. The third column lists the ...

As the number of voxels in the quantification is large, the statistical uncertainty is small in the shell calculation. However, as the exact location or geometry of an object may not be perfectly defined in MRI magnitude images due to the partial volume effect or the signal loss, the uncertainty due to the object geometry can be a concern. For that, we purposely calculated the discrete Greens function based on the radius of the inner sphere as 17 pixels. The calculated susceptibility through the least squares fit deviated from the input value by roughly 18%. This agrees with the volume increase of the inner sphere as the magnetic field distribution is proportional to the magnetic moment of the object, which is proportional to the product of the volume and susceptibility.

3.4 Phantom results

Most of the results of the phantom study are similar to the results of the shell simulation. In order to remove any significant error due to the method itself, we again restrict ourselves to analyze the central 1283 voxels of the magnitude and complex divided images (see the dashed box shown in Fig. 2a). The results are shown in Table 2b. The calculated susceptibilities are within 3% of the theoretical susceptibility of water relative to air, which is −9.4 ppm in SI units [16]. The results indicate that the statistical uncertainties are smaller than the actual uncertainties which are probably due to the non-prefect definition of the object geometry.

In either the simulations or the phantom studies, at least two 2563 matrices are saved and renewed for data processing in the computer memory. As shown in Fig. 1, fast Fourier transformation of the geometry or of the Greens function is performed outside the iteration loop. Fast Fourier transformation of one 2563 matrix takes roughly two minutes. An iteration of susceptibility quantification also takes roughly two minutes. For example, the total time of the entire procedure shown in Fig. 1 with ten iterations is roughly 22 minutes.

4 Discussion and Conclusions

As known in image reconstruction (see Ch. 12 of [12]), sampling in k-space domain will lead to a periodic reconstructed object in the spatial domain. As a result, the magnetic field distribution is also aliased [7]. This is the main reason why all results based on the continuous Greens function are worse than those from the discrete Greens function. Ideally, the magnetic field calculated based on a discrete Greens function should not produce any error. However, as the discrete Greens function is generated by finite sampling the original Greens function in the spatial domain, it is also inevitable that the discrete Greens function used in this method are also aliased (albeit to a lesser degree) in k-space. Generally, the discrete Greens function will start to fail when the object size increases and becomes comparable to the field of view, as shown in Table 1 and Figs. 4, ,5,5, and and6.6. The problem due to the finite sampling can be alleviated by enlarging the field of view relative to the object size [4, 7]. This is consistent with the results from all of our simulations. However, the immediate problem is that the computing time will increase dramatically when a larger matrix size is used in calculations or practically when collecting data with a pre-determined resolution. Despite these limitations, as long as the discrete Greens function is used, the method can be used accurately even when the object size (i.e., diameter) is as large as 60% of the field of view.

When the field of view is enlarged with a fixed image resolution, the resolution in k-space is increased. This obviously implies that enough points around the k = 0 point are needed for the reconstruction of accurate magnetic field in the spatial domain. This is also due to the discontinue value at k = 0. On the other hand, increasing the spatial resolution does not improve the method. Increasing the spatial resolution at a fixed FOV is the same as increasing the window function in k-space. A window function that consists of less than 10 points can lead to a significant blurring effect but this is often not the case in MRI. However, an improved image resolution with sufficient SNR will help in the definition of the object geometry and in the measurement of local phase values near the edge of the object.

The numbers listed in Table 1 may be used to explain the results shown in Fig. 6 and to estimate the size of a spherical shell used for the quantification of susceptibility, when a given uncertainty is desired. For example, in the third column of Table 1a, when the radius of the spherical object is 96 pixels, a total of 781,096 voxels have uncertainties within 5% of the exact field values. As shown in Fig. 6a, if voxels in a shell are between radius 101 pixels and 97 pixels and are used for the susceptibility quantification, then the calculated susceptibility is roughly 6% deviated from its true value. The number of voxels in the shell is roughly 490,000 voxels, less than the above number of voxels. The difference (or excess) in the number of voxels indicates that more voxels can be included in the susceptibility quantification such that the error of the susceptibility will remain the same. The 781,096 voxels imply that the radius of the larger sphere can be increased to 103 pixels rather than 101 pixels. The calculated susceptibility from the case of the larger shell (between 97 and 103 pixels along the radius) is also roughly 6% deviated from the true susceptibility value. The method with the discrete Greens function basically allows us to include more points for analysis without penalty.

The background phase value, [var phi]0, may be estimated from the phase value at the k = 0 point and be used as an initial value in the iterative program. In theory, this is true if an object has a uniform zero phase value or certain symmetry in its complex images. We have confirmed this conjecture from our simulated shell. In addition, the phase is roughly 0.54 radian at the k = 0 point from the phantom images while our program found a [var phi]0 of roughly 0.75 radian (Table 2). Nonetheless, whether one starts [var phi]0 with a zero estimate or the central k-space estimate seems to make little difference in our least squares fit algorithm. In order to further obtain a better set of initial values, the value p can be assigned to three first and then reduced to either one or two in the program. The program seems to converge faster with a higher value of p. However, the higher the value of p, the more the noisy data are included. Although the 2D discrete Greens function seems to produce slightly larger error than that from the 3D Greens functions (see Figs. 5 and and6),6), the 2D method only requires one slice for analysis. Therefore one can afford computing time and enlarge the field of view. When the Fourier based method is utilized, the overall geometry of the object is better to be included within the field of view. Although it may be reasonable to neglect the field distribution due to a boundary that is relatively far away from the region of interest [4], one should not adjust the background field by a constant without care as stated in results.

Voxels containing low SNR in the magnitude images will lead to high uncertainty in the quantification of object susceptibility. In order to remove phase aliasing, phase unwrapping or complex division method may be applied in combination with the Fourier based method. Alternatively, our iterative procedure can automatically exclude the phase aliasing voxels. This is better than to select non-aliasing voxels manually as we want to keep voxels close to the object which have high phase to noise ratios. Induced eddy currents or shimming performed on the scanner can change the phase values in images. The latter problem is minimized by shimming a spherical phantom for quality control prior to scans.

Despite the fact that the statistical uncertainty of each susceptibility value shown in Table 2a is very small, the statistical uncertainty is not large enough to explain the difference between the quantified susceptibility and the true susceptibility. This means that the uncertainty of the method itself (see, for example, Figs. 5 and and6)6) rather than the statistical uncertainty from white noise is a major source of the uncertainty. Fundamentally, the Fourier based method will completely fail when the object size is equal to the image resolution [17, 18]. Practically, the geometry of small objects and their induced field distributions cannot be accurately prescribed without sufficient resolution. The problem is that the magnetic field will decrease quickly so that there will not be enough voxels to accurately estimate the susceptibility. In addition, due to the partial volume effect or the signal loss, the definition of object geometry can also significantly affect the answer. As a result, the uncertainty of susceptibility quantification will increase. These factors require further careful studies of this method.

In summary, we have shown that discrete Greens functions used in the Fourier based method will lead to much better results in estimating local magnetic fields from a given object of interest and that a least squares fit approach based on this method can be applied to quantify susceptibility accurately.

Acknowledgments

This work was supported by the Department of Radiology at Wayne State University and by NIH HL062983-04A2.

Appendix

A Derivation of the continuous 3D Greens function in k-space

Equation (5) can be easily derived with the following identity provided in [10]

eik·r4π=0ij(kr)m=Ym(θ,φ)Ym(θk,φk)
(29)

where θk and [var phi]k are the angles used in the spherical coordinate system.

The Fourier transformation of the Greens function, Gz,3D(r), is

Gc,3D(k)F{Gz,3D(r)}d3rGz,3D(r)ei2πk·r=(13cos2θk)0dr1rj2(2πkr)=13cos2θk13kz2k2
(30)

where cos θk [equivalent] kz/k. This derivation is valid only when k is not equal to zero or infinity.

Throughout this paper, we choose exp{−i2πk·r} as the Fourier transformation. Although it is well understood that such a choice is not unique, it is acceptable as long as its inverse Fourier transformation is consistently defined.

B Calculation of the 2D Lorentz disk

Equation (6) was established on a 2D plane. If the direction of a dipole is always assumed to be parallel to the main field and if the dipole is embedded in 3D space, due to the inner product in Eq. (6), only the transverse component of the dipole on the 2D plane will generate a nonzero magnetic field. If only the field produced by Eq. (6) parallel to the main field is of interest, only the projected component along the main field direction is needed. For these two reasons, an additional factor sin2 θ is added in front of Eq. (10) and is also seen in the first term of Eq. (31) below.

Similar to discussions under Sec. 5.6 in [10] and Eq. (6), the magnetic field at the dipole with a magnitude of m can be calculated from averaging the field within an infinitesimal disk with a radius R

ρ<Rm^·B(x)dS=μ02πρ<RdSm^·(m·xρ2)+ρ<RdSμ0m^·mδ(2)(x)=μ0m2π(m^·n^)2R(Rdφ)+μ0m=μ0msin2θ2π02πdφcos2(φφm)+μ0m=μ0msin2θ2+μ0m=(1sin2θ2)μ0m
(31)

Note that in this derivation we have defined the unit vector [n with circumflex] [equivalent] x/R = cos [var phi]x̂ + sin [var phi]ŷ where the x and y axes define the plane of the disk and m [equivalent] sin θ cos [var phi]mx̂ + sin θ sin [var phi]mŷ + cos θz. The component perpendicular to the plane of the disk (i.e., the component parallel to the axis of the infinitely long object) has been neglected for all discussions relevant to 2D calculations. The result of Eq. (31) can be considered as the Lorentz disk and it is the induced magnetic field inside an infinitely long cylinder (see Eq. (25.35) of [12]). The Lorentz sphere in the 3D case can also be calculated with a similar derivation of Eq. (31).

C Derivation of the continuous 2D Greens function in k-space

The continuous Greens function Gc,2D(k) in k-space can be derived with the help of the following definition of the Bessel functions.

Jm(x)=12πim02πdφeixcosφimφ
(32)

Therefore,

Gc,2D(k)sin2θ2πd2xei2πk·x2cos2φ1ρ2+16(3cos2θ1)=sin2θcos2θk,2D0dρ1ρJ2(2πkρ)+16(3cos2θ1)=sin2θ2·ky2kz2ky2+kz2+16(3cos2θ1)
(33)

where kkz2+ky2 in this derivation and cos2θk,2D(kz2ky2)/(kz2+ky2). This result is valid only when k is not equal to zero or infinity.

References

1. Haacke EM, Cheng YCN, House MJ, Liu Q, Neelavalli J, Ogg RJ, Khan A, Ayaz M, Kirsch W, Obenaus A. Imaging iron stores in the brain using magnetic resonance imaging. Magn Reson Imaging. 2005;23(1):1–25. [PubMed]
2. Haacke EM, Ayaz M, Khan A, Manova ES, Krishnamurthy B, Gollapalli L, Ciulla C, Kim I, Petersen F, Kirsch W. Establishing a baseline phase behavior in magnetic resonance imaging to determine normal vs. abnormal iron content in the brain. J Magn Reson Imaging. 2007;26(2):256–264. [PubMed]
3. Haacke EM, Xu Y, Cheng YCN, Reichenbach J. Susceptibility weighted imaging (SWI) Magn Reson Med. 2004;52(3):612–618. [PubMed]
4. Koch KM, Papademetris X, Rothman DL, de Graaf RA. Rapid calculations of susceptibility-induced magnetostatic field perturbations for in vivo magnetic resonance. Phys Med Biol. 2006;51(24):6381–6402. [PubMed]
5. Deville G, Bernier M, Delrieux JM. NMR multiple echoes observed in solid 3He. Phys Rev B. 1979;19:5666–5688.
6. Salomir R, de Senneville BD, Moonen CT. A fast calculation method for magnetic field inhomogeneity due to an arbitrary distribution of bulk susceptibility. Concepts Magn Reson Part B (Magn Reson Engineering) 2003;198:26–34.
7. Marques JP, Bowtell R. Application of a Fourier-based method for rapid calculation of field inhomogeneity due to spatial variation of magnetic susceptibility. Concepts Magn Reson Part B (Magn Reson Engineering) 2005;25B(1):65–78.
8. Jenkinson M, Wilson JL, Jezzard P. Perturbation method for magnetic field calculations of nonconductive objects. Magn Reson Med. 2004;52:471–477. [PubMed]
9. Yoder DA, Zhao YS, Paschal CB, Fitzpatrick JM. MRI simulator with object-specific field map calculations. Magn Reson Imaging. 2004;22:315–328. [PubMed]
10. Jackson JD. Classical Electrodynamics. 3. Hoboken: John Wiley and Sons; 1999.
11. Sen PN, Axelrod S. Inhomogeneity in local magnetic field due to susceptibility contrast. J Appl Phys. 1999;86(8):4548–4554.
12. Haacke EM, Brown RW, Thompson MR, Venkatesan R. Magnetic resonance imaging: physical principles and sequence design. New York: John Wiley and Sons; 1999.
13. Bevington PR, Robinson DK. Data Reduction and Error Analysis for the Physical Sciences. New York: McGraw-Hill; 1992.
14. Pandian DS, Ciulla C, Haacke EM, Jiang J, Ayaz M. Complex threshold method for identifying pixels that contain predominantly noise in magnetic resonance images. J Magn Reson Imaging. 2008;28:727–735. [PubMed]
15. Gradshteyn IS, Ryzhik IM. Table of Integrals, Series, and Products. 5. San Diego: Academic Press; 1994.
16. Lide DR, editor. Handbook of Chemistry and Physics. 87. New York: CRC Press; pp. 2006–2007.
17. Sepulveda NG, Thomas IM, Wikswo JP. Magnetic susceptibility tomography for three dimensional imaging of diamagnetic and paramagnetic objects. IEEE Trans Magn. 1994;30(6):5062–5069.
18. Tan S, Ma YP, Thomas IM, Wikswo JP. Reconstruction of two-dimensional magnetization and susceptibility distributions from the magnetic field of soft magnetic materials. IEEE Trans Magn. 1996;32(1):230–234.