Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2962550

Formats

Article sections

Authors

Related links

Phys Med Biol. Author manuscript; available in PMC 2010 October 22.

Published in final edited form as:

Published online 2009 January 30. doi: 10.1088/0031-9155/54/5/005

PMCID: PMC2962550

NIHMSID: NIHMS239939

Corresponding author: Y.-C. N. Cheng, MRI Center/Concourse Research, Harper University Hospital, 3990 John R Street, Detroit MI, 48201, Tel: 313-966-8220, Fax: 313-745-9182, Email: ude.enyaw@61cxy

The publisher's final edited version of this article is available at Phys Med Biol

See other articles in PMC that cite the published article.

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.

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.

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

$$\overrightarrow{B}(\overrightarrow{r})=\frac{{\mu}_{0}}{4\pi}{\int}_{V}{d}^{3}r\left\{\frac{3\overrightarrow{M}(\overrightarrow{r})\xb7(\overrightarrow{r}-\overrightarrow{r})}{{\mid \overrightarrow{r}-\overrightarrow{r}\mid}^{5}}(\overrightarrow{r}-\overrightarrow{r})-\frac{\overrightarrow{M}(\overrightarrow{r})}{{\mid \overrightarrow{r}-\overrightarrow{r}\mid}^{3}}\right\}$$

(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 *B _{z}*(

$${B}_{z}(\overrightarrow{r})=\frac{{\mu}_{0}}{4\pi}{\int}_{V}{d}^{3}r\left\{\frac{3{M}_{z}(\overrightarrow{r}){(z-z)}^{2}}{{\mid \overrightarrow{r}-\overrightarrow{r}\mid}^{5}}-\frac{{M}_{z}(\overrightarrow{r})}{{\mid \overrightarrow{r}-\overrightarrow{r}\mid}^{3}}\right\}\equiv {\mu}_{0}{\int}_{V}{d}^{3}r{M}_{z}(\overrightarrow{r}){G}_{z,3D}(\overrightarrow{r}-\overrightarrow{r})$$

(2)

$${B}_{z}(\overrightarrow{k})\equiv \mathcal{F}\{{B}_{z}(\overrightarrow{r})\}={\mu}_{0}{M}_{z}(\overrightarrow{k}){G}_{c,3D}(\overrightarrow{k})$$

(3)

where *M _{z}* (

$${G}_{z,3D}(\overrightarrow{r})\equiv \frac{1}{4\pi}\xb7\frac{3{z}^{2}-{r}^{2}}{{r}^{5}}\equiv \frac{1}{4\pi}\xb7\frac{3{cos}^{2}{\theta}_{3D}-1}{{r}^{3}}$$

(4)

where *θ*_{3}* _{D}* is the azimuthal angle in the spherical coordinate system and

$${G}_{c,3D}(\overrightarrow{k})\equiv \mathcal{F}\{{G}_{z,3D}(\overrightarrow{r})\}=\frac{1}{3}-\frac{{k}_{z}^{2}}{{k}^{2}}$$

(5)

where
${k}^{2}\equiv {k}_{x}^{2}+{k}_{y}^{2}+{k}_{z}^{2}$ and the notation *G _{c,}*

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 *$\stackrel{\u20d7}{M}$* (*$\stackrel{\u20d7}{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.

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 (*$\stackrel{\u20d7}{H}$*) on a 2D plane due to the scalar potential of a dipole *$\stackrel{\u20d7}{m}$* may be written as (also see Eq. (13) in [11])

$$\overrightarrow{H}(\overrightarrow{x})=-\frac{1}{2\pi}\overrightarrow{\nabla}\left(\frac{\overrightarrow{m}\xb7\overrightarrow{x}}{{\rho}^{2}}\right)=\frac{1}{2\pi}\left(\frac{2\overrightarrow{x}(\overrightarrow{m}\xb7\overrightarrow{x})}{{\rho}^{4}}-\frac{\overrightarrow{m}}{{\rho}^{2}}\right)$$

(6)

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

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

$$\overrightarrow{B}(\overrightarrow{x})=\frac{{\mu}_{0}}{2\pi}\int {d}^{2}x\left\{\frac{2\overrightarrow{M}(\overrightarrow{x})\xb7(\overrightarrow{x}-\overrightarrow{x})}{{\mid \overrightarrow{x}-\overrightarrow{x}\mid}^{4}}(\overrightarrow{x}-\overrightarrow{x})-\frac{\overrightarrow{M}(\overrightarrow{x})}{{\mid \overrightarrow{x}-\overrightarrow{x}\mid}^{2}}\right\}+{\mu}_{0}\overrightarrow{M}(\overrightarrow{x})$$

(7)

It is obvious that *$\stackrel{\u20d7}{M}$* (*$\stackrel{\u20d7}{x}$*) in the last term can be replaced by ∫*d*^{2}*x$\stackrel{\u20d7}{M}$* (*$\stackrel{\u20d7}{x}$*)*δ*^{(2)}(*$\stackrel{\u20d7}{x}$ − $\stackrel{\u20d7}{x}$*) 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.

$$\overrightarrow{B}(\overrightarrow{x})\simeq \frac{{\mu}_{0}}{2\pi}\int {d}^{2}x\left\{\frac{2\overrightarrow{M}(\overrightarrow{x})\xb7(\overrightarrow{x}-\overrightarrow{x})}{{\mid \overrightarrow{x}-\overrightarrow{x}\mid}^{4}}(\overrightarrow{x}-\overrightarrow{x})-\frac{\overrightarrow{M}(\overrightarrow{x})}{{\mid \overrightarrow{x}-\overrightarrow{x}\mid}^{2}}\right\}+\left(1-\frac{1}{2}{sin}^{2}\theta \right){\mu}_{0}\overrightarrow{M}(\overrightarrow{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

$$\overrightarrow{B}(\overrightarrow{x})\simeq \frac{{\mu}_{0}}{2\pi}\int {d}^{2}x\left\{\frac{2\overrightarrow{M}(\overrightarrow{x})\xb7(\overrightarrow{x}-\overrightarrow{x})}{{\mid \overrightarrow{x}-\overrightarrow{x}\mid}^{4}}(\overrightarrow{x}-\overrightarrow{x})-\frac{\overrightarrow{M}(\overrightarrow{x})}{{\mid \overrightarrow{x}-\overrightarrow{x}\mid}^{2}}\right\}+\left(\frac{1}{3}-\frac{1}{2}{sin}^{2}\theta \right){\mu}_{0}\overrightarrow{M}(\overrightarrow{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,

$$\begin{array}{l}{B}_{\omega}(\overrightarrow{x})=\frac{{\mu}_{0}{sin}^{2}\theta}{2\pi}\int {d}^{2}x\left\{\frac{2{(z-z)}^{2}{M}_{\omega}(\overrightarrow{x})}{{\mid \overrightarrow{x}-\overrightarrow{x}\mid}^{4}}-\frac{{M}_{\omega}(\overrightarrow{x})}{{\mid \overrightarrow{x}-\overrightarrow{x}\mid}^{2}}\right\}+\left(\frac{1}{3}-\frac{1}{2}{sin}^{2}\theta \right){\mu}_{0}{M}_{\omega}(\overrightarrow{x})\\ \equiv {\mu}_{0}\int {d}^{2}x{M}_{\omega}(\overrightarrow{x}){G}_{\omega ,2D}(\overrightarrow{x}-\overrightarrow{x})\end{array}$$

(10)

where

$${G}_{\omega ,2D}(\overrightarrow{x})\equiv \frac{{sin}^{2}\theta}{2\pi}\xb7\frac{{z}^{2}-{y}^{2}}{{({z}^{2}+{y}^{2})}^{2}}+\frac{3{cos}^{2}\theta -1}{6}{\delta}^{(2)}(\overrightarrow{x})\equiv \frac{{sin}^{2}\theta}{2\pi}\xb7\frac{2{cos}^{2}\phi -1}{{\rho}^{2}}+\frac{3{cos}^{2}\theta -1}{6}{\delta}^{(2)}(\overrightarrow{x})$$

(11)

where *ρ*^{2} *z*^{2} + *y*^{2} and is the polar angle between *z* and *y* axes. Here, the *y*–*z* 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 sin^{2} *θ* in Eq. (10) as discussed in Appendix B. Again, in the derivation of Eq. (11), the term, (2 cos^{2} − 1)*/ρ*^{2}, is already assumed to be zero at *ρ* = 0. The Fourier transformation of *G _{ω}*

$${G}_{c,2D}(\overrightarrow{k})\equiv \mathcal{F}\{{G}_{\omega ,2D}(\overrightarrow{x})\}=\frac{{sin}^{2}\theta}{2}\xb7\frac{{k}_{y}^{2}-{k}_{z}^{2}}{{k}_{y}^{2}+{k}_{z}^{2}}+\frac{1}{6}(3{cos}^{2}\theta -1)=\frac{1}{3}-{sin}^{2}\theta \frac{{k}_{z}^{2}}{{k}_{y}^{2}+{k}_{z}^{2}}$$

(12)

Similar to the 3D case, Eq. (12) is valid only when
${k}_{y}^{2}+{k}_{z}^{2}$ is not equal to zero or infinity. When *$\stackrel{\u20d7}{k}$* = 0, *G _{c,}*

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.

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), *G _{z,}*

For a non-ferromagnetic object, its magnetization can be expressed in terms of its magnetic susceptibility and an external field *B*_{0}.

$${\mu}_{0}\overrightarrow{M}(\overrightarrow{r})\simeq \chi (\overrightarrow{r}){\overrightarrow{B}}_{0}$$

(13)

If an object has a constant susceptibility *χ*, then the Fourier transformation of its magnetization along the main field direction, *M _{z}* (

$${B}_{z}(\overrightarrow{r})={\mathcal{F}}^{-1}\{\mathcal{F}\{\chi (\overrightarrow{r})\}\xb7G(\overrightarrow{k})\}{B}_{0}$$

(14)

where *G*(*$\stackrel{\u20d7}{k}$*) represents either the 2D or 3D discrete Greens function or the continuous Greens functions. We evaluate the 3D and 2D Greens functions and *M _{z}*(

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.

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=\sum _{i=1}^{n}{\left(\frac{{B}_{i}-\chi {B}_{0}{g}_{i}}{\delta {B}_{i}}\right)}^{2}=\sum _{i=1}^{n}{\left(\frac{{\phi}_{i}-\gamma {T}_{E}\chi {B}_{0}{g}_{i}}{\delta {\phi}_{i}}\right)}^{2}=\sum _{i=1}^{n}{\text{SNR}}_{i}^{2}{({\phi}_{i}-\gamma {T}_{E}\chi {B}_{0}{g}_{i})}^{2}$$

(15)

where *n* is the total number of voxels used in the fit, *B _{i}* is the measured magnetic field from MR images at each voxel

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

$$\chi {B}_{0}=\frac{{\displaystyle \sum _{i=1}^{n}}{B}_{i}{g}_{i}{\text{SNR}}_{i}^{2}}{{\displaystyle \sum _{i=1}^{n}}{g}_{i}^{2}{\text{SNR}}_{i}^{2}}$$

(16)

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

$$\delta \chi =\frac{1}{\gamma {T}_{E}{B}_{0}\sqrt{{\displaystyle \sum _{i=1}^{n}}{g}_{i}^{2}{\text{SNR}}_{i}^{2}}}$$

(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, Δ*T _{E}*, 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 _{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=\sum _{i=1}^{n}{\text{SNR}}_{i}^{2}{({\phi}_{i}-{\phi}_{0}-\gamma {T}_{E}\chi {B}_{0}{g}_{i})}^{2}$$

(18)

Similar to the above derivation, both values of *χ* and _{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].

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 _{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 _{0}. The actual phase value due to the object at the *i*-th voxel is * _{i}* −

$$\mid {\phi}_{i}-{\phi}_{0}-\gamma {T}_{E}\chi {B}_{0}{g}_{i}\mid \phantom{\rule{0.16667em}{0ex}}\le p/{\text{SNR}}_{i}$$

(19)

where SNR* _{i}* is the signal to noise ratio of the

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 NiSO_{4} 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.

Images of the water phantom. (a) Transverse magnitude image at *T*_{E} = 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: *T _{R}* 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

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

$$\begin{array}{l}{\mu}_{0}{M}_{z}(\overrightarrow{k})=\int {d}^{3}r\chi {B}_{0}\mathrm{\Theta}(a-r){e}^{-i2\pi \overrightarrow{k}\xb7\overrightarrow{r}}\\ =4\pi \chi {B}_{0}{\int}_{0}^{a}dr\phantom{\rule{0.16667em}{0ex}}{r}^{2}{j}_{0}(2\pi kr)\\ =\frac{2\chi {B}_{0}{a}^{2}}{k}{j}_{1}(2\pi ka)\end{array}$$

(20)

where Θ is the step function, and *j*_{0} and *j*_{1} are spherical Bessel functions. Based on Eqs. (5) and (29), the induced magnetic field in the spatial domain is

$$\begin{array}{l}{B}_{z}(\overrightarrow{r})={\mathcal{F}}^{-1}\{{\mu}_{0}{M}_{z}(\overrightarrow{k}){G}_{c,3D}(\overrightarrow{k})\}\\ =\frac{8\pi}{3}\chi {B}_{0}{a}^{2}(3{cos}^{2}{\theta}_{3D}-1){\int}_{0}^{\infty}dk\phantom{\rule{0.16667em}{0ex}}{kj}_{1}(2\pi ka){j}_{2}(2\pi kr)\\ =\{\begin{array}{ll}\frac{1}{3}\chi {B}_{0}(3{cos}^{2}{\theta}_{3D}-1){\left(\frac{a}{r}\right)}^{3}\hfill & \text{when}\phantom{\rule{0.16667em}{0ex}}r>a\hfill \\ 0\hfill & \text{when}\phantom{\rule{0.16667em}{0ex}}r<a\hfill \end{array}\end{array}$$

(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 *R _{i}* and

$${\mu}_{0}{M}_{z}(\overrightarrow{r})=({\chi}_{i}-{\chi}_{e}){B}_{0}\mathrm{\Theta}({R}_{i}-r)+{\chi}_{e}{B}_{0}\mathrm{\Theta}({R}_{e}-r)$$

(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

$${B}_{z}(\overrightarrow{r})=\{\begin{array}{ll}0\hfill & \text{when}\phantom{\rule{0.16667em}{0ex}}r<{R}_{i}\hfill \\ \frac{1}{3}({\chi}_{i}-{\chi}_{e}){B}_{0}(3{cos}^{2}{\theta}_{3D}-1)\phantom{\rule{0.16667em}{0ex}}{\left(\frac{{R}_{i}}{r}\right)}^{3}\hfill & \text{when}\phantom{\rule{0.16667em}{0ex}}{R}_{i}<r<{R}_{e}\hfill \\ \frac{1}{3}\left[({\chi}_{i}-{\chi}_{e})\phantom{\rule{0.16667em}{0ex}}{\left(\frac{{R}_{i}}{r}\right)}^{3}+{\chi}_{e}{\left(\frac{{R}_{e}}{r}\right)}^{3}\right]\phantom{\rule{0.16667em}{0ex}}{B}_{0}(3{cos}^{2}{\theta}_{3D}-1)\hfill & \text{when}\phantom{\rule{0.16667em}{0ex}}{R}_{e}<r\hfill \end{array}$$

(23)

This result agrees with solutions from solving the Laplace equation in the magnetostatics case when higher order terms of *χ _{i}* and

Another example is to consider a sphere with radius *R _{i}* and constant susceptibility

$${\mu}_{0}{M}_{z}(\overrightarrow{r})=({\chi}_{i}-{\chi}_{e}){B}_{0}\mathrm{\Theta}({R}_{i}-r)+{\chi}_{e}{B}_{0}\mathrm{\Theta}({R}_{e}-\sqrt{{x}^{2}+{y}^{2}})$$

(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

$${B}_{z}(\overrightarrow{r})=\{\begin{array}{ll}\frac{1}{3}{\chi}_{e}{B}_{0}\hfill & \text{when}\phantom{\rule{0.16667em}{0ex}}r<{R}_{i}\hfill \\ \frac{1}{3}({\chi}_{i}-{\chi}_{e}){B}_{0}(3{cos}^{2}{\theta}_{3D}-1)\phantom{\rule{0.16667em}{0ex}}{\left(\frac{{R}_{i}}{r}\right)}^{3}+\frac{1}{3}{\chi}_{e}{B}_{0}\hfill & \text{when}\phantom{\rule{0.16667em}{0ex}}{R}_{i}<r<{R}_{e}\hfill \\ \frac{1}{3}({\chi}_{i}-{\chi}_{e})\phantom{\rule{0.16667em}{0ex}}{B}_{0}(3{cos}^{2}{\theta}_{3D}-1)\phantom{\rule{0.16667em}{0ex}}{\left(\frac{{R}_{i}}{r}\right)}^{3}\hfill & \text{when}\phantom{\rule{0.16667em}{0ex}}{R}_{e}<r\hfill \end{array}$$

(25)

Even if radius *R _{e}* is increased to a size that is much larger than

In this subsection, we assume an infinitely long cylinder with radius *a* and susceptibility *χ* at an angle *θ* with the main field *B*_{0} 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.

$${\mu}_{0}{M}_{\omega}(\overrightarrow{x})=\chi {B}_{0}\mathrm{\Theta}(a-\sqrt{{z}^{2}+{y}^{2}})$$

(26)

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

$${\mu}_{0}{M}_{\omega}(\overrightarrow{k})=\mathcal{F}\{{\mu}_{0}{M}_{w}(\overrightarrow{x})\}=2\pi \chi {B}_{0}{\int}_{0}^{a}d\rho \phantom{\rule{0.16667em}{0ex}}\rho {J}_{0}(2\pi \sqrt{{k}_{z}^{2}+{k}_{y}^{2}}\rho )=\chi {B}_{0}a\frac{{J}_{1}(2\pi \sqrt{{k}_{z}^{2}+{k}_{y}^{2}}a)}{\sqrt{{k}_{z}^{2}+{k}_{y}^{2}}}$$

(27)

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

$$\begin{array}{l}{B}_{\omega}(\overrightarrow{x})={\mathcal{F}}^{-1}\{{\mu}_{0}{M}_{\omega}(\overrightarrow{k}){G}_{c,2D}(\overrightarrow{k})\}\\ =\chi {B}_{0}a\frac{{sin}^{2}\theta}{2}\xb7\frac{{z}^{2}-{y}^{2}}{{z}^{2}+{y}^{2}}{\int}_{0}^{\infty}{\mathit{dkJ}}_{1}(ka){J}_{2}(k\sqrt{{z}^{2}+{y}^{2}})+\frac{\chi {B}_{0}}{6}(3{cos}^{2}\theta -1)a{\int}_{0}^{\infty}{\mathit{dkJ}}_{1}(ka){J}_{0}(k\sqrt{{z}^{2}+{y}^{2}})\\ =\{\begin{array}{ll}\chi {B}_{0}{sin}^{2}\theta \frac{{a}^{2}}{2({z}^{2}+{y}^{2})}\xb7\frac{{z}^{2}-{y}^{2}}{{z}^{2}+{y}^{2}}\hfill & \text{outside}\phantom{\rule{0.16667em}{0ex}}\text{the}\phantom{\rule{0.16667em}{0ex}}\text{cylinder}\hfill \\ \frac{\chi {B}_{0}}{6}(3{cos}^{2}\theta -1)\hfill & \text{inside}\phantom{\rule{0.16667em}{0ex}}\text{the}\phantom{\rule{0.16667em}{0ex}}\text{cylinder}\hfill \end{array}\end{array}$$

(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].

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 *G _{c}*) and the discrete Greens functions (labeled as

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.

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 **...**

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].

Susceptibility error in percentage versus the ratio of the diameter of the object to the field of view. The notations *G*_{c} and *G*_{d} represent the continuous and discrete Greens functions, respectively. All pixels outside the objects are used in the quantifications. **...**

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 128^{3} 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).

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.

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 128^{3} 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 256^{3} 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 256^{3} 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.

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, _{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 _{0} of roughly 0.75 radian (Table 2). Nonetheless, whether one starts _{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.

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

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

$${e}^{i\overrightarrow{k}\xb7\overrightarrow{r}}\equiv 4\pi \sum _{\ell =0}^{\infty}{i}^{\ell}{j}_{\ell}(kr)\sum _{m=-\ell}^{\ell}{Y}_{\ell m}^{\ast}(\theta ,\phi ){Y}_{\ell m}({\theta}_{k},{\phi}_{k})$$

(29)

where *θ _{k}* and

The Fourier transformation of the Greens function, *G _{z,}*

$$\begin{array}{l}{G}_{c,3D}(\overrightarrow{k})\equiv \mathcal{F}\{{G}_{z,3D}(\overrightarrow{r})\}\\ \equiv \int {d}^{3}{rG}_{z,3D}(\overrightarrow{r}){e}^{-i2\pi \overrightarrow{k}\xb7\overrightarrow{r}}\\ =(1-3{cos}^{2}{\theta}_{k}){\int}_{0}^{\infty}dr\frac{1}{r}{j}_{2}(2\pi kr)\\ =\frac{1}{3}-{cos}^{2}{\theta}_{k}\\ \equiv \frac{1}{3}-\frac{{k}_{z}^{2}}{{k}^{2}}\end{array}$$

(30)

where cos *θ _{k}*

Throughout this paper, we choose exp{−*i*2*π$\stackrel{\u20d7}{k}$*·*$\stackrel{\u20d7}{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.

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 sin^{2} *θ* 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*

$$\begin{array}{l}{\int}_{\rho <R}\widehat{m}\xb7\overrightarrow{B}(\overrightarrow{x})dS=-\frac{{\mu}_{0}}{2\pi}{\int}_{\rho <R}dS\widehat{m}\xb7\overrightarrow{\nabla}\left(\frac{\overrightarrow{m}\xb7\overrightarrow{x}}{{\rho}^{2}}\right)+{\int}_{\rho <R}dS{\mu}_{0}\widehat{m}\xb7\overrightarrow{m}{\delta}^{(2)}(\overrightarrow{x})\\ =-\frac{{\mu}_{0}m}{2\pi}\oint \frac{{(\widehat{m}\xb7\widehat{n})}^{2}}{R}(Rd\phi )+{\mu}_{0}m\\ =-\frac{{\mu}_{0}m{sin}^{2}\theta}{2\pi}{\int}_{0}^{2\pi}d\phi \phantom{\rule{0.16667em}{0ex}}{cos}^{2}(\phi -{\phi}_{m})+{\mu}_{0}m\\ =-\frac{{\mu}_{0}m{sin}^{2}\theta}{2}+{\mu}_{0}m\\ =\left(1-\frac{{sin}^{2}\theta}{2}\right){\mu}_{0}m\end{array}$$

(31)

Note that in this derivation we have defined the unit vector *$\stackrel{\u20d7}{x}$*/*R* = cos *$\widehat{x}$* + sin *ŷ* where the *x* and *y* axes define the plane of the disk and sin *θ* cos * _{m}$\widehat{x}$* + sin

The continuous Greens function *G _{c,}*

$${J}_{m}(x)=\frac{1}{2\pi {i}^{m}}{\int}_{0}^{2\pi}d\phi \phantom{\rule{0.16667em}{0ex}}{e}^{\mathit{ix}cos\phi -im\phi}$$

(32)

Therefore,

$$\begin{array}{l}{G}_{c,2D}(\overrightarrow{k})\equiv \frac{{sin}^{2}\theta}{2\pi}\int {d}^{2}x\phantom{\rule{0.16667em}{0ex}}{e}^{-i2\pi \overrightarrow{k}\xb7\overrightarrow{x}}\frac{2{cos}^{2}\phi -1}{{\rho}^{2}}+\frac{1}{6}(3{cos}^{2}\theta -1)\\ =-{sin}^{2}\theta cos2{\theta}_{k,2D}{\int}_{0}^{\infty}d\rho \frac{1}{\rho}{J}_{2}(2\pi k\rho )+\frac{1}{6}(3{cos}^{2}\theta -1)\\ =\frac{{sin}^{2}\theta}{2}\xb7\frac{{k}_{y}^{2}-{k}_{z}^{2}}{{k}_{y}^{2}+{k}_{z}^{2}}+\frac{1}{6}(3{cos}^{2}\theta -1)\end{array}$$

(33)

where
$k\equiv \sqrt{{k}_{z}^{2}+{k}_{y}^{2}}$ in this derivation and
$cos2{\theta}_{k,2D}\equiv ({k}_{z}^{2}-{k}_{y}^{2})/({k}_{z}^{2}+{k}_{y}^{2})$. This result is valid only when *k* is not equal to zero or infinity.

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 ^{3}He. 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.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |