Quantitative susceptibility mapping (QSM) allows the generation of three-dimensional maps showing the variation of the relative magnetic susceptibility within the human body [1-3]. These maps, which are derived from phase images [4] produced using simple gradient echo imaging, carry useful diagnostic information, particularly relating to variation the amount of paramagnetic iron in different tissues [5]. They also allow the differentiation of diamagnetic calcifications from paramagnetic blood products, and the identification of strongly myelinated regions in the brain. In order to produce susceptibility maps it is necessary to carry out a number of processing steps, first to convert the wrapped phase measurements into a map of the field variation inside the region of interest, second to separate out the field perturbation generated by the tissue in the region of interest from that produced by external sources of background field variation and finally to calculate the susceptibility map from the field perturbation. Each of these steps will be briefly described here, following a short discussion of the relationships between tissue magnetic properties and measured signals.

The volume magnetic susceptibility, χ, which is the property that is measured in susceptibility mapping, characterizes the magnetization per unit volume, M, induced in a material when it is exposed to a magnetic field, H, via the expression, M = χH. in many materials, χ, is a scalar quantity so that the induced magnetization is parallel to the applied magnetic field, and in general in biological materials χ is very small (the susceptibility of most tissues in the body being of the order of -9 x 10^-6). The magnetic flux density, B, is in turn related to M and H and thus to the susceptibility by B = μ₀(H+M) = μ₀(1+χ) H.

When an object, such as the human body is exposed to the strong uniform magnetic flux density, B₀ $$$ \widehat{\bf{z}} $$$, of an MRI magnet, it becomes weakly magnetized and the induced magnetization generates a small spatially varying magnetic field perturbation, ΔB(r), which disturbs the homogeneity of the applied field. The resulting field inhomogeneity, broadens lines in magnetic resonance spectroscopy, and if large enough can cause distortion and signal drop-out in magnetic resonance images.

The form of the field perturbation depends on the spatial variation of magnetic susceptibility in the object, described by χ(r), with large field variations generally occurring close to boundaries between regions of different magnetic susceptibility, where the normal to the boundary is parallel to B₀. This, for example, gives rise to the large field perturbations produced in the brain above the air-tissue boundaries in the air-filled sinuses in the skull. We can find the form of ΔB(r) by summing up the dipolar fields, D(r), generated by all of the induced magnetization, M_z(r′) = (B₀/µ₀) χ(r′), induced in the object, yielding

$$\frac{\Delta B(\textbf r)}{B_{0}}=\frac{1}{4\pi}\int_{}^{}\chi\otimes D \space dV'=\frac{1}{4\pi}\int_{}^{}\chi(\textbf r')\frac{3cos^{2}\theta-1}{|\textbf r'-\textbf r'|^{3}} \space dV \space Eq. [1] $$

where $$$ \theta $$$ is the angle between the field direction and the vector (r - r′) linking the magnetization element at r’ to the point where the field is measured, r. This expression is cumbersome to use, since it requires a three-dimensional integral to be evaluated for each position where we want to find the field, but fortunately it can be simplified by three dimensional Fourier transformation. With the application of the convolution theorem this gives

$$\frac{\Delta B(\textbf k)}{B_{0}}=\chi (\textbf k) \times D(\textbf k) \space =\chi (\textbf k)\left( \frac{1}{3}-\frac{{k_{z}}^2}{k^{2}}\right) \space Eq. [2]$$

where k_z is the component of the k-vector along the field direction [6,7]. This expression allows rapid evaluation of the field perturbation produced by any structure, through application of a single forward and a single inverse 3D Fourier transformation. It can be usefully applied to simulating the field perturbations generated by many structures when placed inside an MR scanner and importantly it incorporates the sphere of Lorentz correction [8] needed to estimate the actual field offsets experienced by the nuclear magnetization.

In the presence of a field offset, ΔB, the Larmor frequency is changed by an amount, Δω=γΔB, where γ is the magnetogyric ratio. The phase of the complex signal measured from a voxel experiencing such a frequency offset at echo time, TE, is consequently ∅ =-ΔωTE =-γ ΔB TE, which means the field perturbation can be measured directly from the (unwrapped) phase, by dividing by -(γ TE).

Equation [3] can also easily be rearranged to give the following

$$ \chi (\textbf k)=\frac{\Delta B(\textbf k)}{B_{0}\left( \frac{1}{3}-\frac{{k_{z}}^2}{k^{2}}\right)} \space Eq. [3] $$

which forms the basis of susceptibility mapping, by allowing the susceptibility to be calculated from measurements of the field variation derived from scaled phase maps [5,6].

Although Equation [3] potentially provides a straightforward route to calculating a susceptibility map, by Fourier transformation of the measured phase, followed by division by the function, D(k), and inverse Fourier transformation, there are some significant complications which have to be addressed.

I) Since phase values in the maps generated at the scanner are wrapped into the interval -π to π, phase maps generally need to be "unwrapped" in order to produce correct phase values which can be used to estimate, ΔB(r).

II) The measured phase maps do not only show the effects of field perturbations due to the susceptibility distribution in the imaging region of interest: they are also sensitive to sources of field variation lying outside the ROI (e.g. due to currents in shim coils or magnetic susceptibility in other body regions). It is necessary to filter out the effects of these externally generated fields before the susceptibility map can be calculated.

III) The function, D(k), tends to zero over two conical surfaces in k-space (where $$$ \frac{{k_z^2}}{k^{2}}=\frac{1}{3} $$$) which subtend the magic angle (54.7^o) with respect to the field direction. It is necessary therefore to take account of the fact that χ(k) is not defined at these surfaces and also that k-space noise contributions in the measured ΔB(k) can be greatly amplified by division by D(k).

IV) The values of ΔB(r) can only be measured in regions of the object that generate an MR signal - e.g. in the brain, but not in the skull or air around in the head - so we do not have access to the full field pattern needed to calculate ΔB(k).

Fortunately, as detailed below, methods of addressing all of these issues have been developed in recent years.

Phase Unwrapping Phase wrapping is a common problem in MRI and other imaging methods, and a range of different approaches for resolving phase ambiguities via phase unwrapping have consequently been developed [2]. These include iterative methods that operate directly on the three dimensional images, as well as approaches that work on the Fourier transform of the phase or complex data data (k-space methods). One of these k-space methods which exploits the properties of the Laplacian of the phase has been most widely applied to phase unwrapping for QSM [11], because it is fast and fits nicely into the general framework of susceptibility map calculation [12].

Background Field Removal Magnetic fields (ΔB_z-ext) generated by sources outside the region of interest have different spatial properties compared with those produced by the susceptibility distribution inside the ROI (ΔB_z-int). In particular fields from external sources obey Laplace's equation $$$ \triangledown^{2}$$$B_z-ext, while those due to internal sources do not [13,14]. This opens up a range of different approaches for filtering out the background fields from external sources, including identifying the distribution of dipole fields from sources outside the ROI which best fits the field from inside the ROI and then subtracting this from the measured fields [19,22]. The most widely used method, known as SHARP (Sophisticated Harmonic Artifact reduction for Phase data), relies on the harmonic property of functions that obey Laplace's equation, which means that convolution with any radially symmetric function leaves them unchanged [13]. Consequently, convolving the measured field with a small sphere shape and subtracting the resulting field from the measurement, eliminates ΔB_z-ext. Subsequent deconvolution with the sphere shape function provides a map of ΔB_z-int alone. This approach to background field removal is fast because the convolution and deconvolution can be implemented as multiplications in k-space. Problems with SHARP and other approaches to background field removal arise very near to the external boundary of the ROI, but these are being addressed in advanced approaches which vary the size of the spherical convolution kernel or extend the field variation outside the ROI [15,16].

Susceptibility Map Calculation The simplest method of dealing with the divide-by-zero problem inherent in calculating the susceptibility from the measured field using Eq. [4] is to threshold the value of D(k) so that in the region where |D(k)| < ε , is set equal to sgn[D(k)] × ε [17-19]. The choice of threshold in this TKD (thresholded k-space division) approach is made as a compromise between avoiding noise amplification (by choosing too low a value of ε) and not significantly underestimating susceptibility values (by choosing too high a value of ε). The TKD approach to susceptibility map calculation is fast to carry out and yields surprisingly good results for such a simple technique [17].

Alternatively, the divide-by-zero issue can be avoided by combining field data acquired with the object of interest oriented at multiple angles to the field [19-20]. This is possible because object rotation rotates the conical surfaces with respect to χ(k), and hence the combination of data from multiple orientations can yields estimates of χ(k) at all locations without divide-by-zero problems. The resulting COSMOS (Calculation of Susceptibility by using Multiple Orientation Sampling) approach yields optimal susceptibility maps, but the requirement for rotation of the organ of interest is obviously not practical in many circumstances. New approaches for rapid acquisition of high resolution phase maps will help in this regard [21].

To account for (IV) the susceptibility map van be formed for single or multiple orientation data by finding the susceptibility values in the ROI that produce the best fit to the measured field variation [22,23]. In applying this approach, it is often helpful to apply some regularization to stabilize the fit. A range of different approaches are available [1], with some using additional information from the simultaneously measured magnitude image data or from the phase data itself. The morphology enabled dipole inversion (MEDI) approach [24,25], for example, penalizes the presence of boundaries in the generated susceptibility map which do not overly boundaries in the magnitude data, thus promoting a similar morphological structure in the two image types. Other approaches use regularization to force smoothness within regions where the phase is relatively constant, or to minimize the total norm of the susceptibility map [1].