1 Introduction
Quantitative susceptibility mapping (QSM) [13] is a novel imaging technique that visualizes the magnetic susceptibility distribution from the measured field data associated with magnetization induced in the body by an MR scanner. The magnetic susceptibility is an intrinsic property of the material which relates and the magnetic field through [45]. As physiological and/or pathological processes alter tissues’ magnetic susceptibilities, QSM has been widely applied in biomedical image analysis [45]. Applications include demyelination, inflammation, and iron overload in multiple sclerosis [8], neurodegeneration and iron overload in Alzheimer’s disease [1], Huntington’s disease [50], changes in metabolic oxygen consumption [25], hemorrhage including microhemorrhage and blood degradation [28], bone mineralization [14], drug delivery using magnetic nanocarriers [34].
QSM uses the phase data of a complex gradient echo (GRE) signal as the phase linearly increases with respect to the field perturbation induced by the magnetic susceptibility distribution in an MR scanner [52]. More concretely, assume that an object is placed in an MR scanner with the main static magnetic field where is a positive constant. Then, for any , the observed complex GRE signal at an echo time is modeled as
(1.1) 
where is the proton gyromagnetic ratio, is the total field induced by the susceptibility distribution in an MR scanner, and is the coil sensitivity dependent phase offset. The magnitude image in (1.1) is proportional to the proton density [52], and the phase in is written as
(1.2) 
Based on the observations , QSM aims at visualizing the susceptibility distribution in the region of interests (ROI) which occupies the water and brain tissues. Note that the ROI can be readily determined by (and thus by ) as whenever [30, 44, 52]. The standard QSM consists of the following four steps: offset correction, phase unwrapping, background field removal and dipole inversion (see Fig. 1 for the overview of the process). The first three steps extract the local field that is contained in the total field : the offset correction removes/corrects from to obtain (the offset corrected phase) lying in
; the phase unwrapping removes the artificial jumps in the offset corrected phase when estimating the total field
; the background field removal eliminates the field induced by the susceptibility outside such as skulls and nasal cavity. Interested readers may refer to [24, 44, 52] and references therein for more details.Given the local field , the dipole inversion recovers the susceptibility distribution in by solving the following convolution relation [31, 32, 33]:
(1.3) 
where denotes the principal value [48] of the singular integral with the kernel :
In the frequency domain, Eq. 1.3 reads
(1.4) 
where is the Fourier transform of and by the definition of [13, 24]. From (1.4), it is easy to see that recovering the susceptibility distribution is illposed as on the critical manifold . This illposedness leads to the streaking artifacts unless the data satisfies a proper compatibility condition [9].
1.1 Existing QSM Reconstruction Methods
In the literature, various QSM reconstruction methods have been explored to deal with the illposed nature of the inverse problem Eq. 1.4. Early attempts mainly focus on the direct methods based on the modification of Eq. 1.4 near [27]. One benchmark method, called the truncated Kspace division (TKD) [47], finds the approximate solution to Eq. 1.4 via:
(1.5) 
with a threshold level . Another method recovers via solving the following Tikhonov regularization [29]:
(1.6) 
where and denotes the forward operator that is obtained by discretizing the kernel . Recently, some other direct methods are proposed, e.g. the iterative susceptibility weighted imaging and susceptibility mapping [49], the analytic continuation [39] and so on. Even though these direct methods are simple to implement, they can introduce additional artifacts due to the modification of near in the frequency domain [9, 27, 40].
In recent years, the regularization based methods have been proposed and show the superior performance over the direct method [27, 51]. Mathematically, it is formulated as solving the minimization problem:
(1.7) 
where denotes the data fidelity term and is the regularization term which mostly promotes the sparse approximation of
under some linear transformation such as total variation and wavelet frames. According to the choices of
, the regularization based methods can be classified into the
integral approaches and the differential approaches [27]. The most widely used integral approaches are based on the convolution relation Eq. 1.3. For example, when the data is corrupted by a white Gaussian noise. Even though the integral approach is capable of suppressing streaking artifacts, it is empirically reported in [27] that the reconstructed image can contain the shadow artifacts in the region of piecewise constant susceptibility. The differential approaches are based on the following partial differential equation (PDE)(1.8) 
which is derived from the Maxwell’s equation [23, 45]. In this case, one typical fidelity term is by considering as a measurement. Compared with the integral approach, the differential approach is able to restore susceptibility image with less shadow artifacts. However, the noise in the data can be amplified by , which leads to the streaking artifacts [52]. In [27], the differential approach is implemented by incorporating the spherical mean value (SMV) filter with a radius [30] into the integral approach:
(1.9) 
Since the implementation of causes the erosion of according to the choice of , the loss of anatomical information near is inevitable at the cost of the shadow artifact removal [27].
1.2 Motivations and Contributions of Our Approach
Even though the equations Eqs. 1.8 and 1.3 are known to be equivalent [9, 27, 40], it is observed that the local field defined as Eq. 1.3 is a particular solution of the PDE Eq. 1.8. Whenever the data acquisition is based on the PDE Eq. 1.8, the measured local field data will be written as the superposition of in Eq. 1.3 and the ambiguity of , which will be referred as the harmonic incompatibility. Therefore, there is a need to identify/remove the harmonic incompatibility from the measured local field data for better reconstruction results as it is smooth, analytic and satisfies the mean value property in an open set [19], which are different from the noise properties.
It is noted that the background field removal aims at obtaining the local field via solving a Poisson equation with certain boundary condition as the background field is harmonic in [30, 43, 52, 57]. In this case, the measured local field is represented by the Green’s function associated with the boundary condition [9]. Thus, it is inevitable that contains the incompatibility associated with the imposed boundary condition. In this paper, we investigate the incompatibility of the local filed data in QSM and establish that this incompatibility consists of two harmonic functions inside and outside respectively, and its (distributional) Laplacian defines a surface measure on (see Theorem 2.1 for details and Figs. 5, 4, 3 and 2 for illustrations). Therefore, we can establish a new forward model in QSM by taking this harmonic incompatibility into account.
Based on this discovery, we impose a constraint on harmonic incompatibility term in susceptibility reconstruction model. Since our theoretical results suggest that the incompatibility is harmonic except on , one straightforward approach is to penalize its (discrete) Laplacian on points . However, it is in general difficult to explicitly model this harmonic incompatibility and/or to directly impose its property into the susceptibility reconstruction model due to the complicated geometries of human brains and the limited spatial resolution in real MRI data. Instead, we impose the sparse regularization of the incompatibility as the support of its Laplacian is small compared to the size of image. Combing it with traditional regularization on the susceptibility image, we propose a novel regularization based QSM model by imposing additional constraints on the incompatibility term. Within the new model, we can suppress the incompatibility other than the noise, achieving the whole brain imaging with less artifacts together with the regularization of susceptibility image. Experiments on both brain phantom and vivo MR data consistently show the advantages of the proposed HIRE model which achieves the stateoftheart performance. Besides, our experiments suggest that tight frame regularization of the susceptibility image can avoid the constant offset [27] and lead to efficient computation.
1.3 Organization of Paper
In Section 2, we introduce our HIRE model for whole brain susceptibility imaging. More precisely, we first briefly review the biophysics forward model of QSM in Section 2.1, and characterize the harmonic incompatibility in the local field data in Section 2.2. Based on the characterization, we introduce the proposed HIRE model in Section 2.3, followed by an alternating minimization algorithm in Section 2.4. In Section 3, we present experimental results for both brain phantom and in vivo MR data, and the concluding remarks are given in Section 4.
2 Harmonic Incompatibility Removal (HIRE) Model for Whole Brain Imaging
2.1 Preliminaries on Biophysics of QSM
In an MRI scanner with the main static magnetic field where is a positive constant, objects gain a magnetization . This magnetization generates a macroscopic field satisfying the following magnetostatic Maxwell’s equation [23, 45]
(2.1) 
where is the vacuum permittivity. Since the MRI signal is generated by the microscopic field experienced by the spins of water protons [27], we use the following Lorenz sphere correction model [23]:
(2.2) 
to relate and .
Note that since is generated by field, we have . Moreover, since we consider the linear magnetic materials with , can be approximated as
(2.3) 
Finally, we introduce the total field as
(2.4) 
where denotes the third component of .
Combining Eqs. 2.4, 2.3, 2.2 and 2.1 and taking the third component into account only, we obtain the following relation between and in the frequency domain:
(2.5) 
which gives
(2.6) 
Then for a given susceptibility distribution (in ), the general solution which is bounded everywhere in is expressed as
(2.7) 
where is some constant, and .
In MRI, the phase of a complex GRE MR signal is linear with respect to the total field in Eq. 2.7 [52], and the constant is determined by the coil sensitivity of an MR scanner as the coil sensitivity dependent phase offset is in general assumed to be a constant [24, 44]. However, since we can remove it during the phase estimation from the multi echo GRE signal [12], we assume that and
(2.8) 
in the rest of this paper. Note that defined as above is induced by the susceptibility distribution in the entire space, which is different from in Eq. 1.3.
2.2 Characterization of Harmonic Incompatibility in Local Field Data
In QSM, the total field is obtained from the phase data of a complex GRE MR signal [44, 52]. In fact, if the information of is available over the entire space, then we can directly solve inverse problem from the knowledge of without the background field removal step. However, since the GRE signal is not available outside , the information of is available only inside . Moreover, even if is compactly supported, the support of may not necessarily coincide with that of , which inevitably leads to the information loss outside [44, 52].
Since the total field depends on the susceptibility distribution throughout the entire space [44], it consists of the background field induced from the susceptibility outside , which is of no interest, and the local field by the susceptibility in which we aim to visualize. Since the substantial susceptibility sources are usually located outside which makes the background field dominant in compared to the local field , we need to remove the background field from the (incomplete) total field prior to the dipole inversion [44, 52].
In the literature, given that the background field is harmonic in [52, 57], the background field removal methods take the form of the following Poisson’s equation in [57]:
(2.10) 
Under this setting, we present Theorem 2.1 which characterizes the relation between Eq. 2.10 and the PDE Eq. 2.6, and the measured local field obtained by solving Eq. 2.10 contains an incompatibility which consists of two harmonic functions both inside and outside due to the imposed boundary condition.
Theorem 2.1

There exists such that
(2.11) for , and satisfies
(2.12) for , where and denote the restriction of in and respectively, and
denotes the outward unit normal vector of
. 
Moreover, we have
(2.13) whenever in . Hence, in , and on in this case.
Proof
Since , the governing equation in Eq. 2.10 becomes
Let denote the Green’s function in :
where for each , the corrector function satisfies
Note that since for and , we have for . Consequently, we have
(2.14) 
Then the solution to Eq. 2.10 is represented as
where we used the fact that , and is defined as
for . Then we can see that satisfies
(2.15) 
where the first equation of Eq. 2.15 comes from Eq. 2.14. Here, is induced by the information of only in :
(2.16) 
with
being the characteristic function of
.Based on the fact that for , we define
Hence, we obtain Eq. 2.11, and we can further see that and satisfy
(2.17)  
(2.18)  
(2.19) 
respectively, where Eq. 2.18 comes from Eq. 2.16; , i.e. in , and in .
To prove Eq. 2.12, let , and we consider
Using Eqs. 2.19 and 2.17 and the Green’s identity (e.g. [19]), we have
Similarly using Eqs. 2.19 and 2.18, we also have
where the second equality comes from the fact that we need to compute the inward normal derivatives on . Hence, combining these two equalities, we obtain Eq. 2.12.
To prove 2, we assume on the contrary that there exists such that
for some open and connected set such that and . Choose such that is contained in , where denotes an open ball centered at with radius . Then since in and in , it follows that in by the analyticity of in . Since this means that in , together with the fact that is harmonic in , we have in by the analyticity of in . Since on , we have in . Hence, in , and thus, in , which is a contradiction.
Remark
From the proof of Theorem 2.1, the incompatibility in Eq. 2.11 is from the boundary condition of Eq. 2.10 which is not related to the regularity of . More precisely, inside is generated by the information of the unknown true local field on so that the boundary condition of Eq. 2.10 is matched. In addition, it is obvious that outside is due to the information loss outside .
Remark
Notice that is a “wave type” differential operator (by considering as the time variable). Indeed, the proof of Eq. 2.13 tells us that if in , such has a wave type structure in regardless of its regularity, whereas the susceptibility of human brain does not have such a wave type structure [9]. Hence in QSM, it follows that defined as Eq. 2.12 is a nonvanishing surface measure on , i.e. in , but on .
We present Figs. 5, 4, 3 and 2 to illustrate Theorem 2.1 by using the SheppLogan phantom (Figs. 3 and 2) and the brain phantom (Figs. 5 and 4). Using the limited total field in Figs. (d)d, (d)d, (d)d and (d)d which are derived from Eq. 2.8 by placing strong susceptibilities outside , we solve Eq. 2.10 using multigrid (MG) based the finite difference method [57] to obtain the measured local field in Figs. (f)f, (f)f, (f)f and (f)f which are used for the susceptibility reconstruction. We also display the true local field obtained from Eq. 2.16 in Figs. (e)e, (e)e, (e)e and (e)e for the comparison with the measured . Finally, and are displayed in Figs. (g)g, (g)g, (g)g and (g)g and Figs. (h)h, (h)h, (h)h and (h)h respectively for better illustrations. Compared to the SheppLogan phantom, the brain phantom shows the artifacts as shown in Fig. (h)h. There are two possible reasons of the artifacts. Firstly, since the boundary of human brain is more complicated than the SheppLogan phantom, the erroneous boundary values may have affected the background field removal in the case of brain phantom, as pointed out in [57]. Secondly, unlike the SheppLogan phantom with isotropic spatial resolution , the brain phantom has an anisotropic spatial resolution of . As pointed out in [20], the MG method has a limitation that errors in certain directions ( direction in our case) are not smoothed by standard relaxation and as a consequence, it is inappropriate to coarsen in these directions, which may lead to artifacts in Fig. (h)h along the direction. Since the real spatial resolution of phase data is not necessarily isotropic, an efficient and effective numerical solver of Eq. 2.10 need to be investigated in the future, which is beyond the scope of this paper at this point.
2.3 Proposed HIRE Susceptibility Reconstruction Model
We begin with introducing some notation. Let denote the set of indices of grids, and let denote the set of indices corresponding to the ROI. Denote to be the indices of the boundary of ROI . Finally, the space of real valued functions defined on is denoted as .
Let be the (noisy) measured local field data obtained from Eq. 2.10, which satisfies in . From the viewpoint of Theorem 2.1, we can model it as
where denotes the discretization of the forward operator in Eq. 2.11. Here, is the unknown true susceptibility image supported in , is the incompatibility arising from solving Eq. 2.10, and is the additive noise.
We observe that in the discrete setting, Eq. 2.12 in Theorem 2.1 can be understood as
(2.20) 
with the discrete Laplacian , as is supported on . However, it is in general difficult to directly impose Eq. 2.20 into the susceptibility reconstruction model (e.g. penalizing ) for the following reasons: 1) the estimation of always contains error due to the complicated geometry of the human brain; 2) the real MRI data may not exactly satisfy Eq. 2.20 due to its spatial resolution [23]; 3) the discretization can introduce the error on the boundary of . However, it is a fact that the support of is small compared to , i.e. is sparse. Consequently, we penalize the for the incompatibility term . Although the does not necessarily satisfy the harmonic constraints on , it is a good relaxation approach when considering the error source of the forward model in QSM. In addition, motivated by the successful results on the wavelet frame based image restoration (e.g. [3, 4]), we assume the sparse approximation of under a given wavelet transformation , and propose our HIRE model as follows:
(2.21) 
where with the SNR weight which is estimated from the MRI [32, 38]. Here, is the isotropic norm of the wavelet frame coefficients [4] defined as
(2.22) 
(See Appendix A for the brief introduction on the wavelet frames.)
There are many variational regularizations for the susceptibility image including total variation (TV) [42], total generalized variation (TGV) [2, 55], and weighted TV for morphological consistency [26]. However, since , the subproblem is a rank deficient system matrix, when using the alternating direction method of multipliers (ADMM) methods [18] to solve the regularization model. As a consequence, we may need additional prior information such as the zero susceptibility value in the cerebrospinal fluid region [36] for the stable reconstruction. In contrast, by using the tight frame regularization, the system matrix of subproblem has a full column rank, which can lead to the computational efficiency over the existing variational methods.
Remark
For better understanding of our HIRE model, we temporarily assume that and consider
If , our model reduces to the integral approach model:
In addition, if we fix , our model reduces to the fidelity version of the following differential approach model:
as discretizes the PDE Eq. 2.6 in the sense of [9, Proposition A.1.].
From Remark , we can see that our model considers the incompatibility and noise separately, thereby providing a more precise forward model for QSM. This is because is obtained from the Poisson’s equation Eq. 2.10 and it inevitably contains the harmonic incompatibility related to the imposed boundary condition, as described in Theorem 2.1. Even though more rigorous theoretical analysis is needed, we can somehow explain the effect of harmonic incompatibility in this manner; since the standard arguments on the harmonic functions (e.g. [19]) tell us that is smooth and satisfies the mean value property except on , it has slow variations on this region. As a consequence, it mostly affects the low frequency components in compared to the noise which mainly affects the high frequency components. Together with the fact that the critical manifold forms a conic manifold in the frequency domain, the harmonic incompatibility in mainly leads to the loss of in low frequency components.
As empirically observed in [27], the incompatibility in low frequency components of leads to the shadow artifacts in the reconstructed image, while that in high frequency components leads to the streaking artifacts. Therefore, the simultaneous consideration on the incompatibilities in both components is crucial for better susceptibility imaging. The integral approach does not take the harmonic incompatibility in into account, which may not be capable of suppressing the incompatibility in low frequency components of , and leads to the shadow artifacts in the reconstructed images. The differential approach can be viewed as a preconditioned integral approach since the harmonic incompatibility in has been removed in advance. However, the noise in can be amplified by at the cost of harmonic incompatibility removal, and this leads to the streaking artifacts propagating from the noise in final image [52]. In contrast, the HIRE model takes the form of integral approach which explicitly considers the incompatibility other than the noise by incorporating its sparsity under . By doing so, we expect that the HIRE model can suppress both the noise (cause of streaking artifacts) and the harmonic incompatibility (cause of shadow artifacts), so that we can achieve the whole brain imaging with less artifacts.
We would like to mention that the formulation of HIRE model is not limited to Eq. 2.21. In fact, we can use the nonlinear fidelity term to further compensate the errors in phase unwrapping, which will be more coincident with the GRE signal model [27, 35]. However, we will not discuss the details on such nonlinear variants as this is beyond the scope of this paper. We will focus on Eq. 2.21 throughout this paper.
2.4 Numerical Algorithm
In the literature, there are numerous algorithms which can solve the proposed HIRE model Eq. 2.21. In this paper, we adopt the split Bregman algorithm given in [38] in the framework of ADMM [18] as we can convert Eq. 2.21 into several subproblems which can be solved efficiently. More precisely, let , , , and . Then Eq. 2.21 is reformulated as follows:
Under this reformulation, we summarize the split Bregman algorithm for Eq. 2.21 in Algorithm 1.
(2.23)  
(2.24) 
(2.25)  
(2.26)  
(2.27)  
(2.28) 
(2.29)  
(2.30)  
(2.31)  
(2.32) 
It is easy to see that each subproblem has a closed form solution. The solutions to Eqs. 2.24 and 2.23 can be written as
(2.33)  
(2.34) 
Since we use the periodic boundary conditions, both Eqs. 2.34 and 2.33 can be easily solved by using the fast Fourier transform. In addition, the solutions to Eqs. 2.26 and 2.25 are expressed in terms of the soft thresholding:
(2.35)  
(2.36) 
Here, is the isotropic soft thresholding in [4]: given defined as