Cite this Article
Chen Lin, Zheng Zhi-Wei, Bao Li-Jun, Fang Jin-Sheng, Yang Tian-He, Cai Shu-Hui, Cai Cong-Bo. Weighted total variation using split Bregman fast quantitative susceptibility mapping reconstruction method. Chinese Physics B, 2018, 27(8): 088701  
Permissions
Weighted total variation using split Bregman fast quantitative susceptibility mapping reconstruction method
Chen Lin1, Zheng Zhi-Wei1, Bao Li-Jun1, Fang Jin-Sheng1, Yang Tian-He2, Cai Shu-Hui1, †, Cai Cong-Bo3, ‡
Department of Electronic Science, Fujian Provincial Key Laboratory of Plasma and Magnetic Resonance, Xiamen University, Xiamen 361005, China
Magnetic Resonance Center, Zhongshan Hospital, Medical College of Xiamen University, Xiamen 361004, China
Department of Communication Engineering, Xiamen University, Xiamen 361005, China

 

† Corresponding author. E-mail: shcai@xmu.edu.cn cbcai@xmu.edu.cn

Project supported by the National Natural Science Foundation of China (Grant Nos. 11474236, 81671674, and 11775184) and the Science and Technology Project of Fujian Province, China (Grant No. 2016Y0078).

Abstract

An ill-posed inverse problem in quantitative susceptibility mapping (QSM) is usually solved using a regularization and optimization solver, which is time consuming considering the three-dimensional volume data. However, in clinical diagnosis, it is necessary to reconstruct a susceptibility map efficiently with an appropriate method. Here, a modified QSM reconstruction method called weighted total variation using split Bregman (WTVSB) is proposed. It reconstructs the susceptibility map with fast computational speed and effective artifact suppression by incorporating noise-suppressed data weighting with split Bregman iteration. The noise-suppressed data weighting is determined using the Laplacian of the calculated local field, which can prevent the noise and errors in field maps from spreading into the susceptibility inversion. The split Bregman iteration accelerates the solution of the L1-regularized reconstruction model by utilizing a preconditioned conjugate gradient solver. In an experiment, the proposed reconstruction method is compared with truncated k-space division (TKD), morphology enabled dipole inversion (MEDI), total variation using the split Bregman (TVSB) method for numerical simulation, phantom and in vivo human brain data evaluated by root mean square error and mean structure similarity. Experimental results demonstrate that our proposed method can achieve better balance between accuracy and efficiency of QSM reconstruction than conventional methods, and thus facilitating clinical applications of QSM.

PACS: 87.19.lf;33.15.Kr;02.30.Zz
Keyword:quantitative susceptibility mapping;ill-posed inverse problem;noise-suppressed data weighting;split Bregman iteration
1. Introduction

Nuclear magnetic resonance (NMR) has now reached a remarkable degree of maturity and become a widely available technique in medical and biological applications.[1–3] Magnetic susceptibility serves as an intrinsic physical property of a substance that can provide tissue contrast mechanism in proton magnetic resonance imaging (MRI).[4,5] Quantitative susceptibility mapping (QSM) based on magnetic resonance phase images has become a promising technique for quantifying this property.[6] The QSM has the ability to accurately evaluate the strong susceptibility sources including deoxyhemoglobin and hemoglobin degradation products,[7] investigate iron deposition to differentiate neurodegenerative patients and healthy control,[8–10] and quantify changes of demyelination and iron levels during the formation of multiple sclerosis lesions.[11] As a result, it provides a valuable insight into the diagnosis of specific diseases, such as cerebral hemorrhage, neurological diseases, and multiple sclerosis.

Reconstruction of the susceptibility map is an ill-posed inversion problem due to the singularity around the magic angle of the dipole convolution kernel for the phase to susceptibility transformation.[12] Various strategies have been developed to solve this ill-conditioned problem. Calculation of susceptibility through multiple orientation sampling (COSMOS), which is referred to as the gold standard when compared with the QSM algorithms, may be impractical in clinical applications because of inconvenient data acquisition in different orientations.[13] In regard to single-orientation reconstruction methods, a fast and simple method is truncated k-space division (TKD),[14] which applies a threshold to the dipole kernel. Despite easy implementation, threshold selection is usually a trade-off between the level of streaking artifact and the under-estimation of susceptibility. Another choice of single-orientation methods is to introduce a penalty term and solve the regularization problem by using an optimization algorithm.[15–21] For example, morphology enabled dipole inversion (MEDI),[18,19] where the priori anatomical information from the magnitude map is utilized, has proved to be an effective method of reconstructing susceptibility.[22] However, the long reconstruction time is a major issue for MEDI and other similar regularization methods, which may restrict their clinical applications. The total variation using split Bregman (TVSB) proposed by Bilgic et al. accelerates the L1-regularized susceptibility reconstruction.[23,24] Nevertheless, this method abandons the data weighting matrix in the data fidelity term to employ split Bregman iteration, while data weighting has proved to be significant in artifact suppression,[25] especially around the region of strong susceptibility variation like blood vessels.

In this work, we propose a weighted total variation using split Bregman (WTVSB) to accelerate the L1-regularized susceptibility reconstruction. By combining split Bregman iteration with noise-suppressed data weighting, WTVSB realizes the balance between quality and efficiency in QSM reconstruction. The performance of WTVSB is compared with those of TKD, MEDI and TVSB using numerical simulation, gadolinium phantom, and in vivo human brain data evaluated in terms of root mean square error (RMSE) and mean structure similarity (MSSIM). The influence of parameters in noise-suppressed data weighting on the quality of susceptibility map is also measured. This work may promote the potential applications of QSM in clinical diagnosis and treatment.

2. Theory
2.1. Magnetic field and susceptibility

The local magnetic field can be represented as the convolution of the dipole kernel with the susceptibility distribution, as well as multiplication between the susceptibility and dipole kernel in the Fourier domain, which can be generally written as

where denotes the susceptibility kernel in the Fourier domain, which is defined as the magnetic field generated by a unit dipole; kx, ky, and kz are coordinates in three directions; B(k) and χ(k) respectively represent the measured local field and susceptibility distribution in the Fourier domain; symbol · denotes point-wise multiplication. Zeroes on or near a cone in k-space at a magic angle 54.7° make the direct inversion from magnetic field to susceptibility invalid and lead to severe streaking artifacts.[12]

Regularization constitutes a common approach to process an ill-posed inverse problem. Regularized QSM methods often make use of the observation that the edges in the susceptibility maps are nearly the same as those in magnitude images obtained in the same subject, and their inconsistency can be considered to be sparse. To satisfy this sparsity, we apply a weighted L1 minimization that penalizes voxels on the susceptibility map, which are inconsistent with the magnitude image. Accordingly, the susceptibility reconstruction is formulated as a constrained minimization problem:

Here G is the three-dimensional gradient operator on the susceptibility distribution χ, and Wg denotes the binary gradient weighting diagonal matrix. FD is the matrix representing the convolution with a unit dipole field: FD = F−1DF, b is the calculated local field map, Wn represents the data weighting to account for the field errors, and ε is the tolerance error.

The solution of the constrained convex optimization problem in Eq. (2) coincides with the solution of the unconstrained Lagrangian problem in the following equation with an appropriately chosen regularization parameter λ

where χ* is the estimated susceptibility for the minimization problem.

Many algorithms have been proposed to solve this optimization problem, such as a nonlinear conjugate gradient method,[16,21,23] and lagged diffusivity fixed point method.[17,26] The time-consuming process restricts these solvers’ clinical applications. Therefore, a valid solution with an accelerated process is expected.

2.2. Weighted split Bregman based fast reconstruction method

In this study, we propose a new reconstruction method which combines split Bregman iteration and noise-suppressed weighting in the data fidelity term. The split Bregman iteration is achieved by introducing an auxiliary variable y. The noise-suppressed data weighting matrix Wn is used to suppress the spread of noise and errors from calculated field to susceptibility map. The susceptibility reconstruction is formulated as

The constrained objective is transformed into an unconstrained optimization problem and then solved by two-phase split Bregman iteration:
where μ is an introduced regularization parameter. The variable η adds the mismatch to compensate for y = WgGχ in the unconstrained problem shown in Eq. (5). Then equation (5) can be solved by iterative minimization with respect to χ and y:
Note that the gradient operator in the direction i (i = x, y, z) can be represented asGi = F−1EiF, where Ei is the differential operator in Fourier domain. The diagonal element of Ei is , where j = 0, 1,…, Ni − 1, and Ni denotes the matrix size along the direction i. The noise-suppressed weighting Wn is derived from setting thresholds to the Laplacian of the calculated local field to deal with different levels of field unwrapping errors. It serves as a mask that is zero in voxels with large field changes, one in voxels with small field changes and a transition value between two thresholds,
where ∇2bmax and ∇2bmin are the thresholds for the Laplacian of the magnetic field, which measures its smooth, and Wn is designed because noise and unwrapping errors are prone to exist in the regions of sharp field changes. The same formulation was also utilized to improve the QSM quality of the least-squares QR (LSQR) algorithm.[27] The structural weighting Wg is the binary mask with one and zero in voxels with small gradients and large gradients respectively,
where m is the three-dimensional magnitude image, ∇im denotes the i-dimensional gradient of the magnitude image, and r is the threshold according to the noise level of image m.

With simultaneously adding and subtracting a temporary term (I − Wn)F−1DFχ (I is a unit matrix with the same size as Wn), equation (7) becomes

which can be simplified into
In Eq. (12), if the temporary term (I − Wn)F−1DFχ can be predicted, the minimization problem can be solved effectively. Actually, what needs to be predicted here is a part of a voxel in a field, and we propose to predict the unknown χ in (I − Wn)F−1DFχ with the previous step result χt and update it after each iteration. The reason for making such a prediction lies in the facts as follows. On one hand, most elements in the diagonal matrix (I − Wn) are equal to or close to 0, except in the place where the data weighting value is small. This means that good field points always remain unchanged and only those noisy field points are predicted and updated in iterations. On the other hand, with smoothness constraints and priori structure imposed on the penalty term to obtain new χt in each iteration, the field forward-calculated by χt is expected to be less noisy and more reliable than the previous one, and thus making the substitution of those noisy field points effective. After each iteration, a corrected noisy part is obtained and combined with the unchanged part to form a new whole field to be used in the next iteration.

Fig. 1. Flow chart of the WTVSB algorithm.

The optimal condition for Eq. (12) can be found by taking the gradient and setting it to be zero:

When setting the thresholds to be the gradients of the magnitude image, 30% voxels within the brain mask are considered to be edges,[22] and the matrix Wg is equal to the same size unit matrix I except for approximate 5% elements. According to the theory, we can make an approximation and equation (13) can be efficiently solved by a preconditioned conjugate gradient (PCG) solver with the pre-conditioner (D2 + μE2). As the algorithm is solved iteratively, the difference between χt and χt+1 decreases progressively, so that the problem will be solved within 1% tolerance after a couple of steps.

Equation (8) can be quickly solved by soft thresholding operation:

To demonstrate our method, the flow chart of algorithm implementation is given below.

3. Methods

The performance of the WTVSB method was evaluated on a numerical model with pre-known susceptibility, gadolinium phantom, and in vivo data in comparison with TKD,[14] MEDI,[18,19] and TVSB.[23] The RMSE and MSSIM with respect to the references (known susceptibility of numerical model and COSMOS for gadolinium phantom and in vivo human brain data) were used for quantitative performance evaluation. The datasets of gadolinium phantom were downloaded from website (http://weill.cornell.edu/mri/pages/qsmreview.html). All data processing was implemented in MATLAB 2014a using a workstation with an Intel Xeon processor E5-2630 and 8 GB of RAM.

3.1. Numerical susceptibility model

A three-dimensional numerical brain model with a matrix size of 256 × 256 × 150 was created based on an actual susceptibility map followed by manual segmentation of different anatomical regions. The background area and cerebrospinal fluid regions were set to be zero susceptibility. In cerebral regions, the white matter region was set to be −0.05 ppm, cortical gray matter to 0.04 ppm, veins to 0.35 ppm, and the caudate nucleus, putamen, globus pallidus, red nucleus, substantia nigra, and dentate nucleus to 0.08, 0.1, 0.19, 0.14, 0.16, and 0.09 ppm, respectively. Normally distributed noise with a standard deviation of 3% susceptibility values was added to mimic random physiological susceptibility variations. The field perturbation of the model was generated by fast forward field computation. In addition, the GRE magnitude image was simulated according to the corresponding T1, and spin density values. Magnitude and phase were combined to create a complex signal. The noise was then added into real and imaginary parts of the signal separately and the signal-to-noise ratio (SNR) was set to be 30 dB.

3.2. Gadolinium phantom

A phantom was utilized for experimental validation, the same as the one carried out by Wang et al.[5] Five spherical balloons were filled with gadolinium solutions of various concentrations, which were immersed in a cylindrical container filled with 2% agarose gel. The susceptibility values in five balloons were 0.8, 0.4, 0.2, 0.1, and 0.05 ppm respectively. This phantom was imaged on a 3 T scanner (HDx, GE healthcare, Waukesha, WI) with a multi-echo gradient echo sequence. The imaging parameters were as follows: repetition time (TR) = 70 ms, echo time (TE) = 5 ms ∼ 40 ms, echo spacing (ΔTE) = 5 ms, bandwidth = 480 Hz/pixel, and flip angle (FA) = 15°. The resolution was 1 mm×1 mm×1 mm with a matrix size of 130 × 130 × 86.

3.3. In vivo human brain data

In vivo human data were acquired at the F M Kirby Research Center in the Kennedy Institute and Johns Hopkins University Hospital. Institutional review board (IRB) approval and written informed consent were obtained following a complete description of the study. A healthy volunteer was acquired on a 7 T Philips Healthcare MRI equipped with a 32-channel Novamedical head receive coil. A three-dimensional multi-echo gradient recalled echo (GRE) sequence was utilized. The scan parameters were set to be as follows: TR = 45 ms, TE1 = 2 ms, ΔTE = 2 ms, field of view (FOV) = 220 mm × 220 mm × 110 mm, matrix size = 224 × 224 × 110, the number of echoes = 9, flip angle = 9°, and SENSE factor = 2.5 × 2 for the phase encoding directions. Fat suppression was accomplished by using a water-selective ProSet 121 pulse. Data in multiple head orientations were acquired to calculate susceptibility maps with COSMOS,[13] which were used as references for the adjustment of the regularization parameters and standard values for comparing with the results from different single-orientation QSM reconstruction methods. The volunteer’s head was at four positions: normal supine position, titled to the subject’s right shoulder, titled to the subject’s left shoulder, and titled to the subject’s back. The rotation angle for each orientation varied randomly in a range from 5° to 22° with respect to the main B0 direction.

3.4. Data processing and parameter selection

For the gadolinium phantom, the field was first estimated from the multi-echo datasets with a nonlinear fitting algorithm,[28] and then a magnitude image guided spatial phase unwrapping[29] was implemented to address the frequency aliasing on the field map. A method based on projection onto the dipole field was used to eliminate the background field.[30] For in vivo human brain data, phase wraps were removed by a Laplacian-based phase unwrapping method.[31] The frequency shift of each voxel was calculated by least squares fitting of the linear slope with the phase as a function of TEs using multiple echoes. The initial phase could be estimated by the linear fitting of the phase along with time and we utilized the initial phase to exclude some voxels which had unreliable phase measurement. To remove the background field mainly caused by the susceptibility variations around the air-tissue interfaces, a modified SHARP method (radius = 5 mm, regularization parameter = 0.06) was utilized.[32] Specifically, the radius of the spherical mean filter decreased as it approached to the brain boundary.[21] Different background removal methods were used to test whether or not WTVSB is suitable for inversion from the local field processed by different methods to the tissue susceptibility. We utilized the brain extraction tool (BET) in FMRIB Software Library (FSL, University of Oxford) (FSL BET) to generate the binary mask with the magnitude image. The mean normalized GRE magnitude image at the sixth echo was used to generate the weighting matrix, which performed appropriate tissue contrast at 7 T scanner. Serving as the gold standard susceptibility reference, COSMOS based susceptibility maps were calculated by using the LSQR algorithm with a relative convergence tolerance of 10−5.

4. Results
4.1. Parameter determination

The optimal regularization parameters for MEDI and WTVSB were selected according to RMSEs and MSSIMs of the reconstructed susceptibility maps with different values relative to the numerical susceptibility model, and the corresponding COSMOS maps of the gadolinium phantom and human brain, respectively. Regularization parameter λ adjusts the weight of data fidelity. The decrease of λ leads to strong constraint on the data fidelity but streaking artifacts may not be suppressed effectively, and vice visa. So the choice of the values of regularization parameter λ in MEDI and WTVSB achieves a trade-off between data fidelity and streaking artifacts. Figure 2 shows that RMSE and MSSIM change with λ from 0.0005 to 0.08 for the in vivo data. We determine parameters considering RMSE, MSSIM, and visual perception. The optimal parameter values are 0.0017 for MEDI and 0.0011 for WTVSB and TVSB, respectively. The optimal parameters for numerical susceptibility model and gadolinium phantom were set in the same way and illustrated in the following sections. It has been reported that the regularization parameter μ cannot change the solution to the optimization problem, but it will influence the convergence speed,[24] so we test the convergence behavior of WTVSB with different μ values using the in vivo data. As shown in Fig. 3, the optimal μ value is 0.07. This value also proves to result in fast convergence for the numerical susceptibility model and gadolinium phantom.

Fig. 2. (color online) Selection of regularization parameters for MEDI and WTVSB according to RMSEs and MSSIMs using in vivo human brain data. Optimal parameters are shown with black dots. Unequal interval of the horizontal coordinate is used for better comparison.
Fig. 3. (color online) Convergence behavior of WTVSB for different values of μ using in vivo human brain data.
4.2. Numerical susceptibility model

The simulated local field map, the Laplacian of the local field, and the corresponding noise-suppressed data weighting are described with two typical slices in Fig. 4. It can be seen that some unwrapping errors exist around veins and some edges in the regions of the simulated field map (red arrows). These errors are subsequently amplified and spread into the susceptibility inversion, leading to severe streaking artifacts in the results of TKD and TVSB methods. The noise-suppressed data weighting utilized in the WTVSB method is illustrated in Fig. 4, which availably indicates the positions of errors in the field map.

Fig. 4. (color online) Simulated local field map, Laplacian of the local field, and corresponding noise-suppressed weighting in the numerical susceptibility model.

In Fig. 5(a), the susceptibility maps of the brain model reconstructed by the TKD, MEDI, TVSB, and WTVSB methods are compared on two different slices accompanied with their error maps relative to the true susceptibility maps. For the results from TKD, there exist serious streaking artifacts in several regions (see the error maps). The results from the TVSB method also significantly deviate from the references. The QSM maps generated by MEDI and WTVSB are close to the true ones overall, with an improvement on the streaking artifact suppression. The results imply that data weighting is important in reconstructing the susceptibility map, especially when there are obvious errors in the calculated field. The susceptibility profiles of the line drawn on the two slices validate the accuracies of different reconstruction algorithms in Fig. 5(b). The susceptibility values of TKD and TVSB fluctuate widely around the true values, while the deviations decrease for the values of MEDI and WTVSB (the regions pointed by arrows). The results of quantitative comparison and the reconstruction time among different methods are illustrated in Table 1. Note that even with similar results, the processing time of the proposed WTVSB is only 90 s, much less than that of MEDI (699 s). The threshold in TKD is set to be 0.10 empirically and the regularization parameters are 0.0052 for MEDI, and 0.0032 for TVSB and WTVSB according to the way described in the previous Subsection 4.1.

Fig. 5. (color online) (a) Comparison of susceptibility maps of two slices of numerical susceptibility models reconstructed by TKD, MEDI, TVSB, and WTVSB and their error maps. Error maps of TKD and TVSB present obvious remnants, which are suppressed in MEDI and WTVSB results. (b) Susceptibility values along the line drawn on two slices, where arrows indicate obvious differences due to different reconstructed methods.
Table 1.

Quantitative comparison of different methods in accuracy and reconstruction time.

.
4.3. Gadolinium phantom

The reconstructed results from COSMOS, TKD, MEDI, TVSB, and WTVSB for the gadolinium phantom are shown in Fig. 6. Two zoomed-in local regions in the gadolinium solution are shown in black rectangles for more detail. Although with a relatively large threshold (0.15), there are substantial streaking artifacts in the result of TKD. The TVSB method partly reduces the artifacts, but strong artifacts still appear around the balloons with strong susceptibility. The WTVSB effectively suppresses the artifacts with reconstruction time only increasing by 5 s relative to TVSB. Similar suppression is exerted by the MEDI, but the reconstruction time is eight times longer than that of WTVSB. The regularization parameters for the TVSB, WTVSB, and MEDI are 0.0019, 0.0019, and 0.0033 respectively.

Fig. 6. Comparison of susceptibility maps of two slices of the gadolinium phantom reconstructed by COSMOS, TKD, MEDI, TVSB, and WTVSB.
4.4. In vivo human brain data

Figure 7(a) presents reconstructed maps of a human brain in axial, coronal, and sagittal directions together with the corresponding error maps. The zoomed-in views in black rectangles and the susceptibility values of the line drawn on slices are also exhibited. Compared with the COSMOS map, the TKD result again suffers streaking artifacts, which are obvious in the error maps. The noise and artifacts are suppressed in the MEDI result but the time consumption for reconstruction is heavy. The TVSB method validly reduces the time consumption but loses partial performance in reducing artifacts, which is obvious around the veins and boundary as shown in the zoomed-in region and error map in sagittal view. By contrast, the WTVSB reconstructs the susceptibility map with few artifacts and preserves the fine structures, such as tissue variations and blood vessels. Moreover, the reconstruction does not suffer heavy time consumption. Compared with the results from the TKD and TVSB, the susceptibility values from MEDI and WTVSB are close to those from COSMOS, especially in the regions indicated by the arrows in Fig. 7(b). The threshold for TKD is set to be 0.15 empirically.

Fig. 7. (color online) (a) Comparison of susceptibility maps reconstructed by TKD, MEDI, TVSB, and WTVSB of in vivo human brain data of three views (axial view, 34th slice; coronal view, 106th slice; sagittal view, 113th slice). (b) Susceptibility values of the line drawn on slices, where arrows indicate obvious differences due to different reconstructed methods.

Figures 8(a) and 8(b) show the selected ROI contours of deep gray matter (caudate nucleus (CN), putamen (PT), globus pallidum (GP), substantia nigra (SN), and red nucleus (RN)) in the susceptibility map. The mean susceptibility values and their standard deviations of the five ROIs reconstructed by different algorithms are presented in Fig. 8(c). The susceptibility values calculated using MEDI and WTVSB are closest to those estimated from COSMOS among four algorithms.

Fig. 8. (color online) ((a) and (b)) ROI contours of deep gray matter in susceptibility map; (c) comparison of the mean susceptibility values and standard deviations of ROIs in deep gray matter reconstructed by TKD, MEDI, TVSB, and WTVSB.

Quantitative comparisons of accuracy and reconstruction time among all methods are listed in Table 1. Among the four reconstruction methods, the TKD is the fastest, but it often has the highest RMSE and the lowest MSSIM. As an L1-regularization method, MEDI proves to have a good performance in accuracy but with long reconstruction time. The TVSB effectively accelerates the reconstruction and keeps better accuracy than the TKD. However, as we can see in previous subsections, the results of the TVSB fail to present fine boundaries of structures around strong susceptibility variation regions, which are also reflected from higher RMSE and lower MSSIM values. As shown in Table 1, our WTVSB method can achieve a favorable balance between the computational efficiency and reconstruction accuracy.

4.5. Influence of thresholds ∇2φmax and ∇2φmin

As shown in Fig. 9(a), the curves of RMSE and MSSIM are flat, which suggests that the WTVSB method is insensitive to the change of ∇2φmin. It may lie in the fact that those lower ∇2φ are close to 0 and have an ignorable influence on the calculation of Wn. In Fig. 9(b), when the value of ∇2φmax increases from 0.99 to 0.999, the curves of RMSE and MSSIM are relatively stable. However, when the value is very close to 1, the curves change dramatically. It can be analyzed from Eq. (9) that overlarge ∇2φmax will cause Wn to approach to a unit matrix and lose its efficacy in artifact suppression, which is consistent with the scenario in Fig. 9(b).

Fig. 9. Influences of thresholds ∇2φmin and ∇2φmax on WTVSB results of in vivo human brain data. (a) Changes of RMSE and MMSIM with ∇2φmin in a range of [0.5, 0.9]. (b) Changes of RMSE and MMSIM with ∇2φmax in a range of [0.99, 1].
5. Discussion

In this study, an improved fast L1-regularized QSM reconstruction method is proposed to suppress the artifacts and minimize noise amplification in the susceptibility map. In this method, noise-suppressed data weighting is used to restrain the spread of errors from magnetic field to susceptibility and is appropriately combined with the split Bregman iterations, which results in the acceleration of susceptibility reconstruction. The qualitative and quantitative analyses indicate that the QSM map calculated from WTVSB is comparable to MEDI and superior to TVSB and TKD, and the reconstruction speed is much faster than those of MEDI and the conventional L1-regularization method.

Due to the structure property of the dipole kernel, the noise and errors in the field map will introduce severe streaking artifacts in the reconstructed susceptibility map. The data weighting is utilized to identify these bad points and suppresses the spread of noise and errors in the process of susceptibility inversion. Although a data weighting term exists in both MEDI and WTVSB methods, they play different roles in the implementation. On one hand, the structures of weighting matrix are different. The data weighting in MEDI serves as a normalized diagonal matrix whose diagonal elements are determined by the corresponding magnitudes. This weighting mask can be adjusted after each iteration to further reduce the weightings of bad points.[28] Whereas in the process of WTVSB reconstruction, based on the observation that the noise and errors in the field map usually exist in regions with sharp field changes, the weighting is determined by the Laplacian of the calculated local field. Most of diagonal elements in the resulting diagonal matrix are computed as 1 except for those positions with noisy field points. On the other hand, the weighting effects between them are distinct. Each diagonal element in the matrix of MEDI represents the reliability of corresponding field point, thus different field points would have different influences during the inversion of QSM. By contrast, the field points in the WTVSB method with low weighting are updated in each iteration to obtain a more accurate field. In other words, the effect of data weighting is converted into the update of field data. Anyway, these two weighting schemes both perform effective artifact suppression in the process of QSM reconstruction. A more concrete comparison is beyond the scope of this work.

Split Bergman iteration, which obviously accelerates L1-regularized reconstruction, has also been applied in MRI and CT image reconstructions.[33,34] In the original TVSB method, the introduction of preconditioner D2 + μE2 makes reconstruction quickly solved, because the valid approximation can lead to fast convergence of PCG. However, this approximation is unavailable if the original magnitude-calculated data weighting is incorporated, thus the preconditioner is no longer effective. Our WTVSB employs a new data weighting and proper conversion to make PCG feasible as done previously. Figure 3 shows the convergence behavior of WTVSB algorithm (μ = 0.07) using in vivo human brain data. The relative error is 6.89% at the 4th iteration and it can be reduced to 0.98% at the 18th iteration, which demonstrates a fast convergence rate of WTVSB. All simulation and experimental QSM results demonstrate that WTVSB can substantially suppress the artifacts with computational time comparable to TVSB. This improvement may facilitate the clinical application of regularized susceptibility mapping. The increased computational efficiency of WTVSB compared with MEDI stems from two sources: to reach the same convergence criterion, WTVSB needs approximately 25 times fewer iterations than MEDI, and at each iteration, WTVSB uses FFT operation, whereas MEDI needs FFT and spatial differencing operations.

For TV regularized methods in QSM reconstruction, a precise gradient weighting helps to mitigate the over-smoothness of edges and improve the accuracy of susceptibility map.[22] In this study, like MEDI, the gradient weighting is determined by the magnitude based on the structural consistency between the magnitude and susceptibility. Phase information always serves as another choice for priori gradient in previous studies. A piece-wise gradient weighting from an initial susceptibility obtained by a fast QSM method[27] is also developed and proves to be effective. However, in the cases where noise and errors in field map are abrupt, the obtained initial susceptibility is usually polluted by streaking artifacts, which affect the final susceptibility value in the reconstruction. So it seems infeasible to obtain valid priori gradient using the fast-calculated initial susceptibility. Advanced image processing algorithms may be helpful to obtain a better priori gradient and further improve the quality of WTVSB reconstruction. Nevertheless, WTVSB is introduced to suppress the artifacts in QSM maps as well as to keep the time efficiency, so it may be improper to bring in more complicated processing.

A fault for WTVSB is that it needs a pair of appropriate thresholds to adjust the contribution of data weighting. According to the in vivo human brain experiments detailed in Section 4, the WTVSB never shows great dependence on the choice of ∇2φmin, but it is sensitive to the choice of ∇2φmax. The value of ∇2φmax at 0.999 yields favorable reconstruction results in this work. However, a field with a different level of errors may need a different optimal ∇2φmax. Thus an automatic threshold selection method would be beneficial to further improving the feasibility of WTVSB. Moreover, it should be mentioned that for WTVSB and other currently commonly used QSM methods, a simple dipole model is used to relate the measured field to susceptibility distribution, which is comprehensible but may be insufficient. The accuracy is limited by the assumption that susceptibility is isotropic in cerebral tissues. A more sophisticated model, such as susceptibility tensor model,[35] may help to compensate for this inadequacy and further improve the accuracy of susceptibility map, but it will also increase the difficulty in clinical applications.

6. Conclusions

In this work, a noise-suppressed data weighting determined by the Laplacian of the calculated local field is proposed to help to restrain the spread of noise and errors in susceptibility inversion, which is suited to be incorporated with split Bregman iterations for the acceleration of L1-regularized susceptibility reconstruction. Through appropriate combination of data weighting and split Bregman iterations, the weighted total variation using split Bregman (WTVSB) QSM method effectively facilitates the suppression of streaking artifacts in susceptibility map and retain the efficiency of reconstruction. The balance between computational efficiency and reconstruction accuracy is reasonably realized, which may promote the clinical diagnostic applications of QSM.

Reference
[1] Jiang B Luo F Ding Y Sun P Zhang X Ling J Li C Mao X Yang D Tang C Liu M 2013 Anal. Chem. 85 2523
[2] Zhong K Ernst T Buchthal S Speck O Anderson L Chang L 2011 NeuroImage 55 1068
[3] Ni H Zhou P Zeng P Huang X Liu H Ning X 2015 Chin. Phys. 24 070502
[4] Li J Chang S Liu T Wang Q Cui D Chen X Jin M Wang B Pei M Wisnieff C Spincecmaille P Zhang M Wang Y 2012 Magn. Reson. Med. 68 1563
[5] Wang Y Liu T 2015 Magn. Reson. Med. 73 82
[6] Langkammer C Schweser F Shmueli K Kames C Li X Li G Milovic C Kim J Wei H Bredies K Buch S Guo Y Liu Z Meineke J Rauscher A Marques J Bilgic B 2017 Magn. Reson. Med. 10.1002/mrm
[7] Liu T Surapaneni K Lou M Cheng L Spincemaille P Wang Y 2012 Radiology 262 269
[8] Langkammer C Schweser F Krebs N Deistung A Goessler W Scheurer E Sommer K Reishofer G Yen K Fazekas F 2012 NeuroImage 62 1593
[9] Murakami Y Kakeda S Watanabe K Ueda I Ogasawara A Moriya J Ide S Futatsuya K Sato T Okada K 2015 Am. J. Neuroradiol. 36 1102
[10] van Bergen J Hua J Unschuld P Lim I Jones C Margolis R Ross C van Zijl P Li X 2016 Am. J. Neuroradiol. 37 789
[11] Chen W Gauthier S A Gupta A Comunale J Liu T Wang S Pei M Pitt D Wang Y 2014 Radiology 271 183
[12] Marques J Bowtell R 2005 Concepts Magn. Reson. Part. B 25 65
[13] Liu T Spincemaille P de Rochefort L Kressler B Wang Y 2009 Magn. Reson. Med. 61 196
[14] Shmueli K de Zwart J A van Gelderen P Li T Q Dodd S J Duyn J H 2009 Magn. Reson. Med. 62 1510
[15] Bao L Li X Cai C Chen Z van Zijl P C 2016 IEEE Trans. Med. Imaging 35 2040
[16] Bilgic B Pfefferbaum A Rohlfing T Sullivan E V Adalsteinsson E 2012 NeuroImage 59 2625
[17] de Rochefort L Liu T Kressler B Liu J Spincemaille P Lebon V Wu J Wang Y 2010 Magn. Reson. Med. 63 194
[18] Liu J Liu T de Rochefort L Ledoux J Khalidov I Chen W Tsiouris A J Wisnieff C Spincemaille P Prince M R 2012 NeuroImage 59 2560
[19] Liu T Liu J de Rochefort L Spincemaille P Khalidov I Ledoux J R Wang Y 2011 Magn. Reson. Med. 66 777
[20] Schweser F Sommer K Deistung A Reichenbach J R 2012 NeuroImage 62 2083
[21] Wu B Li W Guidon A Liu C 2012 Magn. Reson. Med. 67 137
[22] Liu T Xu W Spincemaille P Avestimehr A S Wang Y 2012 IEEE Trans. Med. Imaging 31 816
[23] Bilgic B Fan A P Polimeni J R Cauley S F Bianciardi M Adalsteinsson E Wald L L Setsompop K 2014 Magn. Reson. Med. 72 1444
[24] Goldstein T Osher S 2009 SIAM J. Imaging Sci. 2 323
[25] Wang S Liu T Chen W Spincemaille P Wisnieff C Tsiouris A J Zhu W Pan C Zhao L Wang Y 2013 IEEE Trans. Biomed. Eng. 60 3441
[26] Vogel C R Oman M E 1996 SIAM J. Sci. Comput. 17 227
[27] Li W Wang N Yu F Han H Cao W Romero R Tantiwongkosi B Duong T Q Liu C 2015 NeuroImage 108 111
[28] Liu T Wisnieff C Lou M Chen W Spincemaille P Wang Y 2013 Magn. Reson. Med. 69 467
[29] Cusack R Papadakis N 2002 NeuroImage 16 754
[30] Liu T Khalidov I de Rochefort L Spincemaille P Liu J Tsiouris A J Wang Y 2011 NMR Biomed. 24 1129
[31] Li W Wu B Liu C 2011 NeuroImage 55 1645
[32] Schweser F Deistung A Lehr B W Reichenbach J R 2011 NeuroImage 54 2789
[33] Fang S Wu W Ying K Guo H 2013 Acta Phys. Sin. 62 048702 in Chinese
[34] Mao B Chen X Xiao D Fan S Teng Y Kang Y 2014 Acta Phys. Sin. 63 138701 in Chinese
[35] Liu C 2010 Magn. Reson. Med. 63 1471