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 November 1.
Published in final edited form as:
PMCID: PMC2962531
NIHMSID: NIHMS239940

Quantifying Effective Magnetic Moments of Narrow Cylindrical Objects in MRI

Abstract

A new procedure for accurately measuring effective magnetic moments of long cylinders is presented. Partial volume, dephasing, and phase aliasing effects are naturally included and overcome in our approach. Images from a typical gradient echo sequence at one single echo time are usually sufficient to quantify the effective magnetic moment of a cylindrical-like object. Only pixels in the neighborhood of the object are needed. Our approach can accurately quantify the magnetic moments and distinguish subpixel changes of cross sections between cylindrical objects. Uncertainties of our procedure are studied through the error propagation method. Images acquired with different parameters are used to test the robustness of our method. Alternate approaches and their limitations to extract magnet moments of objects with different orientations are also discussed. Our method has the potential to be applied to any long object whose cross section is close to a disk.

Keywords: CISSCO, magnetic moment, MRI, quantification, subpixel, susceptibility

1 Introduction

Changes of blood oxygenation cause a change in the magnetic susceptibility of venous blood. This is the basis of functional magnetic resonance imaging (fMRI) (e.g., see [1]). For this reason, there has been a strong interest to quantify the magnetic susceptibility of venous blood. As a forward problem, researchers in MRI usually model a blood vessel as an infinitely long cylinder and calculate its MR signal with a known susceptibility of the vessel obtained from ex vivo studies [2]. In order to avoid the partial volume effect in MRI but to obtain accurate results, an ex vivo study usually analyzes the phase information both inside and outside a cylindrical test tube whose diameter occupies at least 10 image pixels [2]. If the wall thickness of the test tube can be neglected, this approach can accurately recover the susceptibility difference between the materials inside and outside the cylinder. On the other hand, the in vivo measurement of the blood susceptibility has also been performed on vessels with large enough diameters (e.g., [3]). Nonetheless, the quantification of susceptibility of a cylindrical object with a small diameter has been a challenging problem and has not yet been tackled. Until recently, we [4] and Sedlacik et al. [5] have attempted to solve this inverse problem with two independent approaches. The former requires a priori knowledge of the cylinder radius while the latter fits several variables to MR images from multiple long echo times, which can become impractical in in vivo studies. In this paper, we will demonstrate a general approach in the quantification of the effective magnetic moments of narrow cylindrical objects from only one echo time of a gradient echo image. This new approach relies only on the known imaging parameters but not on the object size, which is usually not known in advance. From the magnetostatic theory [6], in the first order approximation, the magnetic moment per length of a long but narrow object is equal to that of an infinitely long cylinder. This concept makes our approach described in this paper more applicable to long objects whose cross sections may not be exactly a disk. Practically, our approach may be applied to clinical MR images and industrial problems [7].

In this paper, we will also show our improved method to identify the center of the cylindrical object. This improved method is robust as we will prove mathematically. Uncertainties of our methods including both the thermal noise from an object and the systematic uncertainty due to discrete pixels will be studied from the error propagation method [8]. By minimizing the uncertainty due to the thermal noise, we will derive criteria of how to effectively evaluate the magnetic moment of an object from images. Certain technical procedures that have been discussed in detail in [4] will be briefly mentioned here but will not be repeated. In the following sections, we will provide the general formulas but focus on null-signal objects perpendicular to the MRI main field in our phantom studies. Results from different imaging parameters will be shown. Objects with different orientations and limitations of our approach will be discussed toward the end of this paper.

2 Material and Methods

The basic concept of our method is to add the complex MR signal of each voxel around a long cylindrical object of interest. As we model each object as an infinitely long cylinder, we only need to examine signals on a plane whose normal is parallel to the axis of the cylindrical object. The cross section of the cylinder on the plane is a disk of which the circumference is a circle. Thus, these terms will also be used in this paper. Theoretically, the MR signal around the infinitely long cylinder can be easily modeled. The total signal within a co-axial cylinder can be formulated and can be compared to the actual signal obtained from images, such that the effective magnetic moment can be solved. However, as different types of noise exist in images and the actual center of the object requires a subpixel determination, the challenge is to identify a set of robust procedures to fully solve this problem.

2.1 MR signal of a cylindrical object

Applying an external magnetic field on an infinitely long cylindrical object with an absolute susceptibility χi embedded in an environment of susceptibility χo will induce a magnetic field ΔB both inside and outside the object [6]. If the distortion of the object can be neglected in images, after the correction of the Lorentzian sphere term and frequency adjustment by the MR systems, the phase value of the complex MR signal from a gradient echo sequence outside the cylinder is [1]

φout=γΔBTE=γΔχ2a2ρ2B0TEcos2ψsin2θ
(1)

where γ is the proton gyromagnetic ratio (2π · 42.58 MHz/T), TE is the echo time, Δχ [equivalent] χiχo, a is the radius of the cylindrical object, ρ is the perpendicular distance measured from the axis of the cylinder, B0 is the main field of the MR system, ψ is the polar angle associated with ρ, and θ is the angle between the main field direction and the axis of the cylinder. Figure 1 shows the coordinate system used in Eq. 1.

Figure 1
A schematic drawing shows the coordinate systems used in Eq. 1. Note that the axis of the cylinder is chosen to be parallel to the z-axis only for the derivation of Eq. 5. After deriving Eq. 5, only angle θ is relevant for further discussions. ...

In the context of this paper, it is Δχ that we wish to measure from experiments. Thus, the term susceptibility in this paper refers to Δχ rather than the absolute susceptibility of the object. In addition, we have adopted SI units. From Eq. 1, we define the maximal (or minimal, depending on the sign of Δχ) phase value g as

g0.5γΔχB0TE
(2)

and the effective magnetic moment (or “magnetic moment, hereafter) as

pga2
(3)

Under the same consideration, the induced phase value inside the cylinder is [1]

φin=γΔχ6(3cos2θ1)B0TE=g3(13cos2θ)
(4)

If we consider a pseudo cylinder with radius R that is co-axial with the cylindrical object (see Fig. 2), for g ≠ 0, the overall MR signal of the pseudo cylinder from a gradient echo sequence is [4]

Figure 2
A schematic drawing shows the cross section of a cylindrical object with radius a, enclosed by a co-axial pseudo cylinder whose radius is R. The total MR signal within the pseudo cylinder can be formulated.

S=ρ0aRdρρ02πdψeiφout+πa2ρ0,ceiφin=πρ0pp/R2gdxx2J0(xsin2θ)+πa2ρ0,ceiφin
(5)

where [ell] is an arbitrary length of the cylindrical object and can be the slice thickness of the image, J0 is the zeroth order Bessel function, ρ0,c is the effective spin density of the cylindrical object, and ρ0 is the effective spin density of the tissue around the object. Both ρ0,c and ρ0 are constants but depend on imaging and tissue parameters such as relaxation times. Although Eq. 5 is derived for an infinitely long cylinder, practically, it is the MR signal per unit length of the pseudo cylinder (S/[ell]) that is needed in our analysis. Thus, an image slice whose normal component is parallel to the cylinder axis should be used in all quantifications. On that slice, the cross section of the cylindrical object appears as a disk. The same is true for the co-axial pseudo cylinder. The orientation of the cylinder, θ, can be estimated from the coordinates of the two ends of the cylindrical object in images. Therefore, the unknowns in Eq. 5 are reduced to g, p, ρ0, and ρ0,c.

2.2 Determining the center of the cylindrical object

Equation 5 is valid only when the center of the disk (i.e., the cross section of the cylindrical object) is identical to the center of the pseudo disk. Thus the first step is to identify the center of the disk. When a cylindrical object itself has no MR signal, i.e., ρ0,c = 0, we can identify the center of the object by minimizing the imaginary part of the signal (Eq. 5) within the pseudo disk [4].

When a cylindrical object has a non-zero ρ0,c, even if the pseudo circle is replaced by two pseudo concentric circles such that the imaginary part from the annular ring is theoretically zero, the above approach still fails. This is because when the object has an imaginary part of the MR signal, its point spread function leads to an additional imaginary signal in each pixel outside the object. On the other hand, although the “subvoxel shift approach by [5] is valid, it requires many Fourier transformations of an image and thus longer computing time. We offer two alternate approaches below.

The first alternate approach is to maximize the real part of the signal from an annular ring region where the disk is completely inside the smaller pseudo circle. Assume that the annular ring region is formed by two pseudo circles with radii R1 and R2 and a < R1 < R2. The center of the object is located at coordinates (x0, y0) with x02+y02+aR1. The MR signal from the annular ring is

S=ρ0R1R2dρρ02πdψexp{iga2ρ2cos2ψsin2θ}
(6)

where the relations between ρ, ρ, ψ, and ψ in the integrand are given by ρcosψ = x0 + ρcosψ and ρsinψ = y0 + ρsinψ. One can analytically prove that at x0 = y0 = 0, Sx0=Sy0=0, and 2Sx0y0=0. Appendix A proves that the real part of the signal is maximum at x0 = y0 = 0 under certain conditions.

The second approach is to minimize the real part of the signal within one pseudo circle with radius R. With the parameters defined above and with x02+y02<a, the MR signal of the pseudo disk is

S=ρ002πdψr(ψ)Rdρρexp{iga2ρ2cos2ψsin2θ}+πa2ρ0,ceiφin
(7)

where the relations ρcosψ = x0+ ρcosψ and ρsinψ = y0+ ρsinψ still hold. In addition, r(ψ) is the minimal value of ρ and satisfies r(ψ) cosψ = x0 + a cosψ and r(ψ) sinψ = y0 + a sinψ. This leads to

(r(ψ)cosψx0)2+(r(ψ)sinψy0)2=a2
(8)

and

r(ψ)=(x0cosψ+y0sinψ)+(x0cosψ+y0sinψ)2+(a2x02y02)=x02+y02cos(ψψ0)+a2(x02+y02)sin2(ψψ0)
(9)

where cosψ0x0/x02+y02 and sinψ0y0/x02+y02.

Again, one can analytically prove that at x0 = y0 = 0, Sx0=Sy0=0, and 2Sx0y0=0. Appendix B proves that the real part of the signal is minimum at x0 = y0 = 0.

All these approaches will fail when θ = 0. In this special case, no phase value is outside the object (see Eq. 1). The MR signal becomes:

S=ρ0(Aπa2)+πa2ρ0,ce2ig/3
(10)

where A is an arbitrarily shaped area that encloses the disk. Only the size of the object can be determined in this case.

2.3 Magnetic moment of the cylindrical object

As shown in Fig. 3, if three arbitrary concentric circles with radii R1, R2, and R3 are chosen, then three complex signals S1, S2, and S3 can be calculated from MR images as described in [4]. In order to improve the accuracy of our approach, each image pixel is further divided into 100 subpixels in the calculations of the complex signals Si [4]. From Eq. 5, we obtain

Figure 3
(a) Magnitude and (b) its associated phase image at an echo time of 5 ms show an air cylinder in the gel phantom. (c) Magnitude and (d) phase image at an echo time of 20 ms show the same air cylinder in the gel phantom. Although the actual radius of the ...

(S1S2)p/R22p/R32dxx2J0(xsin2θ)=(S2S3)p/R12p/R22dxx2J0(xsin2θ)
(11)

where the effective magnetic moment, p, becomes the only unknown in the equation. Note that both ±p satisfy Eq. 11 but the correct sign of p may be determined from MR phase images with the help of Eq. 1. From Eq. 1, it is clear that psin2θ/Ri2 is the maximal (or minimal) phase value at the circumference of the i-th circle (where i = 1, 2, 3). Therefore, if any Ri is chosen larger than the phase aliasing area, then psin2θ/Ri2 will be always less than π and the solution of |p sin2 θ| can be numerically searched between 0 and πRmin2 with a Van Wijngaarden-Dekker-Brent method [9], where Rmin is the smallest radius among three circles Ri. With the choices of psin2θ/Ri2 in the following subsection, the solution of p2 sin4 θ in Eq. 11 is unique. This is because ∫ dx/x2 J0(x) can be well approximated by −1/xx/4 in our consideration, which leads to a monotonic function of p2 sin4 θ in Eq. 11. For example, when psin2θ/R12 is 0.5 and psin2θ/R22 is 1, the deviation between the exact integral and its approximation is 0.5%. When psin2θ/R22 is 1 and psin2θ/R32 varies from 2.5 to 3.0, the deviation increases monotonically from 7% to 12%. Thus, the initial guess of |p sin2 θ| may be chosen as πRmin2/2. After p is solved, the effective spin density ρ0 can be solved from the signal of an annular ring, e.g., S1S2.

The above discussion is generally valid except when the angle θ is close to zero. At θ = 0, no magnetic moment appears in the signal, which is given by Eq. 10.

2.4 Uncertainty of the magnetic moment

It is important to study the uncertainty of the method when it is applied to actual images. In addition, in order to determine the optimal choices of three radii for the quantification of the magnetic moment, the study of the uncertainty can provide some insights. Thermal noise due to the presence of an object and systematic noise due to discrete pixels in images always exist. They lead to the uncertainty of p. By defining p [equivalent] p sin2 θ and rewriting Eq. 11, we can derive the uncertainty of p through error propagation [8]:

δpp=δpp=1D(δ(S2S3))2h122+(δ(S1S2))2h232
(12)

where hij is defined as

hijp/Ri2p/Rj2dxx2J0(x)
(13)

with i, j = 1, 2, 3 and

D(S1S2)J0(φ3)φ3+(S2S3)J0(φ1)φ1+(S3S1)J0(φ2)φ2=πρ0p(h12J0(φ3)φ3+h23J0(φ1)φ1+h31J0(φ2)φ2)
(14)

Here φip/Ri2 is chosen to be between 0 and π for all i. In the derivation of Eq. 12, we have already assumed that R3 < R2 < R1 such that the uncertainty from the annular ring region between R2 and R3 is uncorrelated with the uncertainty from the area between R1 and R2. The uncertainty from each annular ring consists of the thermal noise and the systematic noise. Both noise sources are also uncorrelated. The former can be approximated by σΔxΔyπRi2Rj2 where σ is the standard deviation of the thermal noise in the image and ΔxΔy is the image in-plane resolution [4]. The systematic noise is represented by εij calculated from δ(SiSj) [equivalent] εij|SiSj| when the thermal noise is neglected. The uncertainty due to the systematic noise δ(SiSj) can be calculated if the center of the cylindrical object, the radius of the object, and the susceptibility difference are all known. Some examples were shown in [4, 10].

Equation 12 can be rewritten as

δpp=ΔxΔyπpSNR2(|1φ21φ3|h122+|1φ11φ2|h232)+(ε122+ε232)h122h232|h12J0(φ3)φ3+h23J0(φ1)φ1+h31J0(φ2)φ2|
(15)

where SNR [equivalent] ρ0/σ is the signal to noise ratio of the magnitude image. In order to determine the optimal combination of [var phi]1, [var phi]2, and [var phi]3, for simplicity, we assume ε12 = ε23 = 0 in Eq. 15 and numerically search for the minimum of δp/p with 0 < [var phi]iπ for i = 1, 2, 3.

2.5 Resolving the susceptibility when the object has no spin density

The susceptibility can be calculated from the known effective magnetic moment, if the radius of the object can be determined. Two approaches are considered to determine the radius of the object when the object has no spin density. The first approach is to utilize a spin echo sequence. The second approach is to perform a gradient echo sequence at a short echo time.

If a spin echo sequence is used to image the cross section of the object, the MR signal from an arbitrarily uniform area A that includes the object is:

SSE=ρ0,SE(Aπa2)
(16)

where ρ0,SE is the effective spin density (a constant) of the spin echo images. We have neglected the dephasing effect during the sampling time. With two arbitrary areas A1 and A2 such that A2 is completely enclosed by A1, ρ0,SE and the cross section of the object πa2 can be determined. The uncertainty of the cross section can be derived from the error propagation method.

δ(πa2)=ΔxΔySNRA2+(A2πa2)2A1A2
(17)

Note that SNR in Eq. 17 refers to the signal to noise ratio of the spin echo image. The noise between A2 and A1A2 is uncorrelated. In order to minimize the uncertainty, it is obvious from Eq. 17 that the area A1A2 should be chosen as large as possible before heterogeneity is encountered in the images, while the area A2 should be chosen as small as possible but larger than πa2. In fact, A2 has to be larger than the area of the distorted object in the image so the cross section of the object can be estimated more accurately.

In the gradient echo approach, an additional pseudo circle with radius R may be used to solve Eq. 5 as g becomes the only unknown after both p and ρ0 are found. However, Fig. 4a demonstrates that the integral in Eq. 5 oscillates and quickly approaches an asymptotic value when |g| or echo time increases (with a fixed p). This is due to the J0 function and 1/x2 in the integrand. In order to obtain a unique solution of g from Eq. 5, the echo time TE needs to be chosen short enough such that |g| is at least less than 2.4, which is the first root of the J0(x) function. However, in order to distinguish the solution from the asymptotic value of the integral with the presence of noise in images, the value of |g| may have to be much less than 2.4. The radius R also needs to be large enough such that |p/R2| is much less than 2.4. We will show some examples in the Results section.

Figure 4
Demonstration of the theoretical signals due to long echo times. For simplicity, no noise has been added in any of these simulations. (a) A plot of the integral in Eq. 5 as a function of g with ρ0[ell] = 10 units, a = 1 mm, TE = 5 ms, ρ ...

2.6 Phantom studies

Images from three different gel phantoms are analyzed in this paper. The first set of images come from the phantom images in [4]. A narrow but long air cylinder with a diameter of roughly 1.6 mm appears in the images. Images acquired from a coronal 3D gradient echo sequence (Fig. 3) and a spin echo sequence are available for our re-analyses. The imaging parameters of the gradient echo sequence were: TE = 5 ms and 20 ms, TR = 50 ms, flip angle = 15°, read bandwidth = 390 Hz per pixel (Hz/pixel), resolution = 1 mm ×1 mm ×1 mm, and fields of view = 256 mm ×128 mm ×64 mm. The imaging parameters of the spin echo sequence were: TE = 8.4 ms, TR = 400 ms, read bandwidth = 130 Hz/pixel, resolution ≈ 0.55 mm ×0.55 mm, fields of view = 140 mm ×140 mm, slice thickness = 1.5 mm, 18 slices, and slice spacing = 1.95 mm.

The second set of images is from a phantom that consisted of a hollow straw in the gel. As one end of the straw was sealed, the hollow straw created another air cylinder in the gel phantom. The diameter of the straw was 5.24 ± 0.01 mm. This phantom was prepared with 107 g of the gel powder in 1200 ml distilled water. The straw was placed vertically in the phantom when the gel was in its liquid form. Coronal images of a multiple echo 3D gradient echo sequence and a spin echo sequence were acquired. The imaging parameters of the gradient echo sequence were: TE = 5 ms, 10 ms, and 15 ms, TR = 25 ms, flip angle = 15°, read bandwidth = 610 Hz/pixel for TE = 5 ms and 220 Hz/pixel for the other two echo times, resolution = 1 mm ×1 mm ×1 mm, and fields of view = 256 mm ×256 mm ×96 mm. The imaging parameters of the spin echo sequence were TE = 15 ms, TR = 550 ms, read bandwidths = 90 Hz/pixel and 590 Hz/pixel, resolution = 1 mm ×1 mm, fields of view = 256 mm ×256 mm, slice thickness = 2 mm, 20 slices, and slice spacing = 2 mm. Both air cylinders were placed perpendicular to the main field of the MR system. All images were acquired from a 1.5 T Siemens Sonata system. The phase images were all filtered by a high pass filter [11] that has a minimal effect on our study here [4]. A high pass filter will also remove the eddy current effects and the constant background phase in a given set of phase images.

The distortion of objects in images (Fig. 5) due to the susceptibility effect was studied by varying the read bandwidth. These studies were preformed on another straw phantom made on a different day. The read bandwidths of the gradient echo sequence at TE = 10 ms were varied from 70 Hz/pixel, 220 Hz/pixel, 350 Hz/pixel, to 700 Hz/pixel. Other imaging parameters remained the same as those listed in the above paragraph.

Figure 5
Gradient echo gel images of an empty straw perpendicular to the main field with different read bandwidths are shown. The echo time of these images is 10 ms. (a) Magnitude and (b) its associated phase image were acquired with a bandwidth of 70 Hz/pixel. ...

2.7 Image simulations

In order to calculate the correct systematic noise (or uncertainty) due to discrete pixels, both image parameters and measured magnetic moments are used. In addition, the susceptibility difference between the air and gel is assumed to be 9.4 ppm [12]. The simulations and the addition of Gaussian noise were described in [4]. It is worth mentioning again that proper rotations of data after Fourier transformation are required. This is due to different definitions of the Fourier transformation in different computer software, but we follow the definition provided in [1].

In all simulations, we assume that θ = π/2, Δχ= 9.4 ppm, and B0 = 1.5 T. The image resolution is set to 1 mm ×1 mm (or 1 pixel ×1 pixel). Other parameters are provided with corresponding simulated result in the next Section.

3 Results

3.1 Center of the object

The one circle approach appears to be the easiest method in the identification of the object center. The radius of the pseudo circle cannot be too small or too large. In the former case, the thermal noise and discrete pixels will contribute significant uncertainty to the process. We would recommend a radius of at least three pixels. In the latter case, the overall MR signal will be dominated from the pixels without magnetic moment information and become insensitive to identifying the center. Equation 36 in Appendix C is consistent with these considerations. Based on the phase images, a circle whose circumference intersects with phase values around ±2 radians along the vertical and horizontal axes is a reasonable choice while the phase of 2.6 radians would be the optimal choice. In the approach of using two concentric circles, we suggest using radii that differ by at least one pixel such that enough pixels are used in the analysis. Based on our experience in simulation and phantom studies, the center of the object identified by any of our three procedures often differs from the actual center by 0.1 to 0.3 pixel. In addition, all the results from the one circle approach agree with the theoretical criterion (Eq. 36 in Appendix C). We have also tried the approach by Sedlacik et al. [5] on some of our simulated images and have identified the same centers determined by our one-circle method. However, as the “subvoxel shift approach by [5] involves a 2D Fourier transformation of every possible center, while our method only requires a sum of complex numbers from the neighboring subpixels, our method takes less than half the computing time of their approach.

3.2 Optimal choice of the radii combination

Simulations of Eq. 15 reveal that the uncertainty of the magnetic moment will be at a minimum when the three phase values are smallest, roughly 1 radian, and 2.5–3 radians. This result means that one circle should be chosen as large as possible while the other two circles should be chosen when their radii are roughly the lengths from the center of the object to the desired phase values. All circles should be outside the phase aliasing region. Due to the left handed setting by MRI manufacturers, note that the signs in all phase images of Fig. 3 are opposite to the signs in Eq. 1 (right handed system). Practically, as the uncertainty in the phase image is inversely proportional to SNR in the magnitude image [1], the radius of the largest circle can be chosen where the phase value is not smaller than 0.1 radian. Furthermore, as the MR signal is discretized and the center of the object may be different from the center of a pixel, the phase value measured from MR images may not be exactly equal to any of the above desired phase values and may not be symmetric around the object. In addition, phase values close to ±π radians may not be available from the phase profiles due to the partial volume effect. Simulated phase profiles shown in Fig. 6 reflect these practical problems.

Figure 6
An infinitely long cylinder with radius 0.8 pixel is simulated without the thermal noise. The center of the object is purposely shifted to the 128.9th pixel in simulations. Our method identified a center that is only 0.1 pixel away from the simulated ...

In any case, the phase profiles are only used as references for choosing the radii of the three circles for the measurements of the magnetic moment. It is not necessary to identify the exact phase values for our measurements. In fact, any choice of three circles can lead to a solution of the magnetic moment except that the uncertainty may be larger than desired (see below). Nevertheless, in order to better estimate the uncertainty (Eq. 15) due to each radius used in our analyses, we have averaged four radii and their corresponding phase values taken from a vertical and a horizontal phase profile through the center of the object.

3.3 Measurements of the magnetic moment

Table 1a lists the magnetic moments obtained from our previous data [4]. As the radii of the three circles are chosen as described above, the uncertainties calculated from Eq. 15 are less than 10%. As a comparison, when the radii of the largest circles are reduced, the uncertainties of the measured magnetic moments are increased but are still within 16% (Table 1b).

Table 1
(a) Magnetic moments of six slices from the existing images at echo times of 5 ms and 20 ms. The first column lists the slice number, with the smaller number representing the top slice. The second and the seventh columns list the values of the three radii ...

In theory, the effective magnetic moment is proportional to the echo time. The results of the magnetic moment obtained from different echo times clearly show that relationship within uncertainties in Table 1. With the known radius of the object 0.8 mm and an assumed Δχ= 9.4 ppm, the theoretical magnetic moments at TE = 5 ms and 20 ms are 6.04 radian·mm2 and 24.14 radian·mm2, respectively. That the magnetic moment at the bottom slice is smaller than that at the top is likely due to the collapsing of the air cylinder in the gel phantom.

Two slices of the second phantom (straw phantom) at three echo times are analyzed and their results are shown in Table 2. As the diameter of the straw is 5.24 ± 0.01 mm, the theoretical magnetic moments at TE = 5 ms, 10 ms and 15 ms are 64.7 radian·mm2, 129.5 radian·mm2, and 194.2 radian·mm2, respectively. Even though the cross section of the straw is slightly deformed during the solidification of the gel solution, the measured magnetic moments are still in good agreement with the theoretical values and between themselves at different echo times. However, when the calculated uncertainty is less than 1%, one will need to consider the inaccuracy due to pure numerical algorithms which can be as large as 0.3% [4].

Table 2
Effective magnetic moments of the air straw measured at two different slices. The notation of each column is identical to that described in Table 1. The SNRs of images at TE = 10 ms and 15 ms are 23:1. The SNR of images at TE = 5 ms is 15:1. The measured ...

Table 3 shows magnetic moments of the straw at four different read bandwidths from the same slice. All the measured magnetic moments agree with each other within uncertainties. Therefore, our measurements of the magnetic moment from a gradient echo sequence do not seem to be affected by the distortion effect even at a relatively low bandwidth. The measured moments in both Tables 2 and and33 are slightly different from the theoretical values as the presence of the gel powder used in the phantoms can shift the susceptibility away from that of pure water. As we choose the radii outside the phase aliasing regions, we also effectively avoid the distortion effect.

Table 3
Study of the distortion effect in gradient echo images (Fig. 5). The theoretical value is 129.5 radian·mm2. The results indicate that distortion has minimal effect to our method.

3.4 Estimations of object volumes from spin echo images

When the object has no spin density, the cross section of the object and its uncertainty are calculated based on Eq. 16 and Eq. 17. Three slices of the spin echo images from our previous study are analyzed as they correspond to slice 10, 19, and 37 of the gradient echo images. The measured cross sections are shown in Table 4. The result from slice 19 agrees well with the theoretical value of the object cross section, 2.01 mm2. In general, the results from Table 4 also indicate a collapse at the bottom of the phantom and an enlargement of the object at the top of the phantom.

Table 4
Three cross sections Ameasured of the narrow air cylinder are analyzed from the spin echo images. The areas of A1 and A2 are shown in the number of pixels while the values of Ameasured and Apredicted are shown in units of mm2. The uncertainties of the ...

For the straw phantom, one slice of the spin echo images at two different read bandwidths is analyzed. The SNR is 33:1 for images acquired with bandwidth 90 Hz/pixel and is 11:1 for images with bandwidth 590 Hz/pixel. The measured cross sections are 21.4 ± 0.3 mm2 and 22.2 ± 0.8 mm2 at bandwidth 90 Hz/pixel and 590 Hz/pixel, respectively. In the former case, A1 = 1148 mm2 and A2 = 120 mm2, while in the latter case A1 = 821 mm2 and A2 = 69 mm2. The theoretical value of the cross section is 21.6 mm2 if the cross section is a perfect disk. As the straw has slightly deformed inside the gel phantom, its cross section may not be a perfect disk but has a smaller cross section than that of the perfect disk, given the same circumference. These results show good agreement between the measurements and the theoretical value within uncertainties predicted by Eq. 17. The results also indicate that the geometric distortion does not significantly affect the quantification.

3.5 Resolving susceptibility from gradient echo simulations with no spin density of the object

In addition to our efforts in the previous subsection, we also explore the possibility of quantifying susceptibility from gradient echo images. We investigate this possibility with simulations. If all variables and parameters are known, we can calculate the MR signal in each case. Such an MR signal is a number and is displayed as a horizontal straight line in Fig. 4a or 4b. On the other hand, if the magnetic moment p is known but the susceptibility and object radius are not known, then we can simulate the MR signal as a function of g, as shown as a curve in Fig. 4a or 4b. The intersections of the curve and the line indicate the possible solutions of the susceptibility. The result shown in Fig. 4a indicates that the minimal value of |Δχ| (with a fixed p) can be extracted at a given echo time. In this example, the minimal value of |Δχ| is roughly 1.5 ppm, much less than 9.4 ppm used for the image simulations. If the echo time is shortened, as shown in Fig. 4b, then g or susceptibility is likely to be uniquely determined. If the thermal noise is added as one arbitrary unit (compared to ρ0[ell], which is 10 units) and an image resolution is assumed as 1 mm2, then the uncertainty of the signal S within a radius R (3 mm) is roughly 5.3 units ( πR2/ΔxΔy). This means that the signal S is within one standard deviation of the asymptotic value shown in Fig. 4a when g is larger than roughly 3.5. On the other hand, with a much shorter echo time (which may be impractical), the susceptibility and radius of the object can be resolved from Eq. 5 and Fig. 4b. However, this result implies that the spin echo approach presented in the above subsection seems better.

The examples shown in Fig. 4a and 4b indicate that one cannot accurately determine the susceptibility and object size through a curve fitting of data from long echo times, where g becomes too large. If sufficient noise had been included in the simulations, both Fig. 4c and 4d demonstrate that the volume fraction and susceptibility each can only be roughly determined at orders of magnitude through curve fitting. However, when an object has enough MR signal, the curve fitting method appears applicable to determine the susceptibility and volume of the object [5]. This is because Eq. 4 only contains the susceptibility information (i.e., g), rather than the magnetic moment (i.e., p).

4 Discussion

Three common problems in MRI have been overcome in our method of obtaining the magnetic moment: partial volume effect, dephasing effect (signal loss shown in Fig. 6a and 6c), and the phase aliasing effect (shown in Fig. 6d). In our examples, the partial volume effect is overshadowed by the dephasing effect. These are practical reasons why our method applied on small objects is much better than the conventional least squares fitting method or the method of measuring the relaxation time T2. Some results of how bad the least squares fitting method has performed are shown in [4].

As shown in Fig. 6, even when the object has no signal, the actual center of the object may not be in the pixel with the lowest magnitude signal. This fact implies that the center of the object has to be determined based on the information outside the object, as suggested by our procedures here or by Sedlacik et al. [5]. However, the approach by [5] will take much longer time than ours.

After the center of the object is identified, two phase profiles through the object center along the vertical and horizontal directions are used for the selections of the radii of three circles. All three circles should be larger than the phase aliasing area around the object but the smallest circle should be as small as possible (i.e., as close to the aliased phase as possible). However, as shown in Fig. 6, the maximum phase value may not be close to ±π radians but a choice of a slightly larger radius of the smallest circle does not seem to increase the uncertainty much. Due to the discrete pixels, a function can be used to interpolate the phase values such that more accurate phase values can be used for the selections of three radii at non-integer pixels. For this purpose, we have tried a linear function and a function based on Eq. 1 and we have found no difference of which function to use. This is clearly supported by the fact that Eq. 11 does not rely on the knowledge of phase values, which only serve as references of how to choose the radii. The radii of non-integer pixels are an important feature in our method. For example, if we consider two radii at 2.3 and 2.6 pixels, both radii refer to the same pixel along the vertical and horizontal directions. However, as the signal inside a circle can change even with a slight change of its radius, Eq. 11 offers a much more accurate way to solve the magnetic moment. Practically, a uniform region may appear only around the neighborhood of the object. Even with non-optimal choices of radii, our method appears to be effective as most of the uncertainties listed in Table 1b are still within 10%.

The uncertainty δp/p derived from Eq. 15 is inversely proportional to SNR and p (or TE). This is clearly shown in Table 1. Although Eq. 15 implies that the longer the echo time, the smaller the uncertainty, it is important to remember that the SNR of the material around the object will also reduce when the echo time increases. Thus, if the echo time is longer than the relaxation time T2 of the surrounding material, the uncertainty will likely increase.

The systematic uncertainty and the uncertainty due to the thermal noise in the measurements of magnetic moments are comparable in most of our phantom studies. Although all sources of uncertainties are considered uncorrelated in Eq. 15, actually, the systematic uncertainties from two annular rings S1S2 and S2S3 may be partially correlated as they share the same boundary (2πR2). However, we have neglected that correlation in this research.

Given the uncertainties and magnetic moments shown in Table 1, the internal collapsing of the air cylinder in the first phantom seems the only explanation. The measurements of the cross sections from the spin echo images also support this finding (Table 4). This collapsing is not obvious from the visual inspection of images. Assuming Δχ= 9.4 ppm, the radius of the air cylinder appears to change from 0.8 mm (slice 19) to 0.67 mm (slice 46). With an image resolution of 1 mm ×1 mm, such a subpixel change of the object radius can be clearly distinguished from our magnetic moment measurements.

The presence of the object susceptibility can lead to distortion artifacts and dephasing effects in images. They occur in all images but in the gradient echo images the dephasing effect at the echo dominates the artifacts unless the read bandwidth is very low. A usual way to minimize the distortion artifacts is to increase the read bandwidth in a sequence. However this change will lead to a lower SNR in images and thus a larger uncertainty (Eq. 15 and Table 4). Our phantom studies indicate that the geometric distortion does not seem to affect the measurements of magnetic moment as shown in Table 3. In the spin echo images, as the dephasing effect due to the susceptibility difference can be enhanced during the sampling period, the measured volume of the object should be treated as the maximal volume. This is supported by the results shown in Table 4. The geometric distortion also does not seem to affect the volume measurements of the objects from the spin echo images.

Although the absolute susceptibility of gel is assumed as the same as the absolute susceptibility of water, −9 ppm, in our analyses, the actual susceptibility of the gel may be slightly different than that of water.

4.1 Resolving susceptibility and volume when the object has no spin density

Further examination of Eq. 17 indicates some interesting aspects and the limitation of the spin echo approach. If the cross section of the object occupies more than one pixel in the image, then the uncertainty estimated from Eq. 17 does not depend on the image resolution but it depends on the fields of view. This is because SNR is proportional to LxLyΔxΔy where Lx and Ly are the fields of view along the two orthogonal directions. However, if the cross section of the object is within one pixel, then the smallest A2 that can be chosen is ΔxΔy. If the second term in Eq. 17 can be neglected, then the uncertainty is proportional to ΔxΔy/LxLy. This means that a spin echo image with a high resolution can determine the size of the object more accurately, even at the expense of a low SNR or a long imaging time. On the other hand, if the object is less than one pixel after the image has been acquired, the uncertainty of Eq. 17 in this case becomes ΔxΔy/SNR. This means that if the volume fraction of the object in a pixel (πa2/(ΔxΔy)) is less than 1/SNR, the estimated uncertainty of the cross section is more than 100%.

Practically, a spin echo sequence is routinely performed in each clinical diagnosis. As the magnetic moment can be determined from one gradient echo sequence described in Section 2.3, knowing the radius of the object will automatically lead to the value of susceptibility. The only exception is when θ = 0. In this special case, only the size of the object can be determined. When a gradient echo image and a spin echo image are analyzed, no registration of the images is required in our approach.

In the gradient echo approach, the susceptibility and volume cannot be resolved individually when the susceptibility is larger than a certain value with a fixed magnetic moment (see Fig. 4a). This is consistent with the theory of electromagnetism, as in the far field the magnetic moment of an object is the leading term [6]. In general, our experience shows that the gradient echo approach for determining object susceptibility is not favorable, as the desired echo time to minimize distortions and signal losses may need to be very short and thus require a very high read bandwidth. For this reason, the spin echo approach is a better method.

5 Conclusion

Our method can accurately quantify the effective magnetic moment of a narrow cylinder from MR images. It has overcome three general problems in MRI and is applicable to small objects in a locally uniform background. However, the application of our method is not limited to small objects. The effective magnetic moment of an object can be obtained from standard gradient echo images while the volume of the object may be derived from typical spin echo images. We have demonstrated the feasibility of our approaches through phantom studies. In addition, a phantom study has revealed that a subpixel change of the object volume can be distinguished from the change of magnetic moment using our method in MRI.

Acknowledgments

This work was supported in part by the Department of Radiology at Wayne State University and NIH HL062983-04A2. The authors would like to thank Mr. Zahid Latif, R.T., and Mr. Yang Xuan for technical assistance in phantom imaging. The authors would also like to acknowledge the use of Mathematica. The author Y.-C. N. Cheng would like to thank Dr. William Condit for his stimulating discussions.

Appendix

A Maximizing the real signal from an annular ring to identify the object center

We first define β [equivalent] ρ2 = (ρcosψx0)2+(ρsinψy0)2 and α [equivalent] (ρcosψx0)2− (ρsinψy0)2 such that

cos2ψ=(ρcosψx0)2(ρsinψy0)2ρ2αβ
(18)

We also define

feipcos2ψ/ρ2=eipα/β2
(19)

such that the MR signal S is

S=ρ0R1R2dρρ02πdψf
(20)

The first derivatives of f with respect to x0 and y0 are

fx0=ipf2β3(x0ρcosψ)[3(ρsinψy0)2(ρcosψx0)2]
(21)
fy0=ipf2β3(y0ρsinψ)[3(ρcosψx0)2(ρsinψy0)2]
(22)

The second derivatives of f with respect to x0 and y0 at x0 = y0 = 0 are

2fx02|x0=y0=0=2p2f1ρ6(1+cos6ψ)6ipf1ρ4cos4ψ
(23)
2fy02|x0=y0=0=2p2f1ρ6(1cos6ψ)+6ipf1ρ4cos4ψ
(24)

As

02πdψf=2πJ0(p/ρ2)
(25)

02πdψf·cos4ψ=2πJ2(p/ρ2)
(26)

02πdψf·cos6ψ=2πiJ3(p/ρ2)
(27)

it becomes clear that

Im2Sx02|x0=y0=0=Im2Sy02|x0=y0=0
(28)

and

Re2Sx02|x0=y0=0=Re2Sy02|x0=y0=0=2πρ0(pR22J1(pR22)pR12J1(pR12))
(29)

As R2 > R1 and the maximum of xJ1(x) occurs at x ≈ 2.4, the first root of J0(x), proper choices of R1 and R2 can lead to a negative value of Eq. 29. Therefore, the object center can be identified by maximizing the real part of the signal from an annular ring.

B Minimizing the real signal of a circle to identify the object center

The calculations used in this Section require the results derived in the Appendix A and Eq. 9. At x0 = y0 = 0, r(ψ) = a and the relevant MR signal is

S=ρ002πdψaRdρρf
(30)

which is identical to Eq. 20 if the upper and lower limits of the ρ integral are properly replaced.

The second derivatives of Eq. 7 are the results derived from Appendix A with proper changes of the integral limits, plus the following term:

2Sx02|extra=x0(Sx0|extra)ρ002πdψr(ψ)r(ψ)x0fx0|ρ=r(ψ)
(31)

where

Sx0|extra=ρ002πdψr(ψ)r(ψ)x0f|ρ=r(ψ)
(32)

The second derivative with respect to y0 has the identical format with x0 replaced by y0.

After tedious derivations, we obtain

Im2Sx02|x0=y0=0=Im2Sy02|x0=y0=0
(33)

and

Re2Sx02|x0=y0=0=Re2Sy02|x0=y0=0=2πρ0(pR2)J1(pR2)>0
(34)

when 0 < |p/R2| < 3.8, the second root of J1(x). That the result of Eq. 34 only depends on R rather than both R and a indicates that this procedure of finding the object center is still valid even if the center of the circle R is outside the object. From a different point of view, when an object (such as a nanoparticle) is very small, one can identify the center of the object by replacing the small object by a much larger object but with the identical magnetic moment p. This is because when the value of |g| is larger than roughly 5, one cannot unambiguously determine the volume of the object based on the gradient echo signal (see Fig. 4a).

C Accuracy of the object center

Even though the center of the object can be determined by the method described in the Appendix B, the center cannot be accurately located due to the presence of noise in images. Consider the Taylor expansion of the MR signal

S(x0,y0)=Sx0=y0=0+12x022Sx02|x0=y0=0+12y022Sy02|x0=y0=0+(higherorderterms)
(35)

where the first order terms vanish, as described in the main text. The center of the object cannot be accurately determined if the sum of the second order terms is less than the noise term and even if the higher order terms are neglected. Using the result of Eq. 34 and considering only the thermal noise, σΔxΔyπR2, the following criterion can be established

x02+y02ΔxΔyπR2ΔxΔyπSNR(pR2)J1(pR2)=pΔxΔyπSNR(pR2)3/2J1(pR2)
(36)

The maximum of x1.5J1(x) occurs at x ≈ 2.6 and the maximum is roughly 2. At x = 1.5 and 2, x1.5 J1(x) = 1.02 and 1.63, respectively. With a given magnetic moment p in an image, Eq. 36 clearly defines how accurate the center can be determined.

References

1. Haacke EM, Brown RW, Thompson MR, Venkatesan R. Magnetic Resonance Imaging: Physical Principles and Sequence Design. John Wiley & Sons; 1999.
2. Weisskoff RM, Kiihne S. MRI susceptometry: image-based measurement of absolute susceptibility of MR contrast agent and human blood. Magn Reson Med. 1992;24:375–383. [PubMed]
3. Lai S, Hopkins AL, Haacke EM, Li D, Wasserman BA, Buckley P, Friedman L, Meltzer H, Hedera P, Friedland R. Identification of vascular structures as a major source of signal contrast in high resolution 2D and 3D functional activation imaging of the motor cortex at 1.5 T: Preliminary results. Magn Reson Med. 1993;30:387–392. [PubMed]
4. Cheng YCN, Hsieh CY, Neelavalli J, Liu Q, Dawood MS, Haacke EM. A complex sum method of quantifying susceptibilities in cylindrical objects: The first step toward quantitative diagnosis of small objects in MRI. Magn Reson Imaging. 2007;25(8):1171–1180. [PubMed]
5. Sedlacik J, Rauscher A, Reichenbach JR. Obtaining blood oxygenation levels from MR signal behavior in the presence of single venous vessels. Magn Reson Med. 2007;58:1035–1044. [PubMed]
6. Jackson JD. Classical Electrodynamics. New York: John Wiley & Sons; 1999.
7. Robson P, Hall L. Identifying particles in industrial systems using MRI susceptibility artifacts. American Institute of Chemical Engineers. 2005;51(6):1633–1640.
8. Bevington PR, Robinson DK. Data Reduction and Error Analysis for the Physical Sciences. New York: McGraw-Hill; 1992.
9. Press WH, Flannery BP, Teukolsky SA, Vetterling WT. Numerical Recipes in C. New York: Cambridge University Press; 1992.
10. Hsieh C-Y, Cheng Y-CN, Neelavalli J, Haacke EM. An improved approach of quantifying magnetic moments of small in-vivo objects. Proceedings of the International Society for Magnetic Resonance in Medicine; 2007. p. 2596.
11. Haacke EM, Xu Y, Cheng YCN, Reichenbach J. Susceptibility weighted imaging (SWI) Magn Reson Med. 2004;52(3):612–618. [PubMed]
12. Lide DR, editor. Handbook of Chemistry and Physics. 87. New York: CRC Press; 2006–2007.