Fast quantitative susceptibility reconstruction via total field inversion with L0 norm approximation
Shuhui Cai1, Li Zhang1, Congbo Cai1, and Zhong Chen1

1Department of Electronic Science, Xiamen University, Xiamen, China


Quantitative susceptibility mapping (QSM) is a meaningful MRI technique owing to its unique relation to actual physical tissue magnetic properties. The reconstruction of QSM is usually decomposed into three sub-problems which are solved independently. Here, we propose a fast reconstruction method named as fast TFI based on total field inversion. It accelerates the total field inversion by using specially selected preconditioner and the advanced solution of weighted L0 regularization. Results from gadolinium phantom and in vivo data verified that the new method has good performance.


Compared to conventional magnitude imaging, QSM has the advantage that it can provide a better contrast in gray matter and white matter based on their differences in the magnetic susceptibility. Although many QSM reconstruction methods have been proposed, some problems still exist, such as the existence of noise, residual artifacts and accumulated errors. In this work, a fast TFI method based on the total field was proposed. Through the variable splitting method and hard thresholding algorithm, the proposed method is more effective than conventional methods.


The whole susceptibility distribution $$$\chi$$$ of a brain is conventionally comprised of background susceptibility $$$\chi _{H}$$$ and local susceptibility $$$\chi_{L}$$$. In this work, we estimate $$$\chi_{H}$$$ and $$$\chi_{L}$$$ jointly using the fast TFI:

$$w_{opt}=\underset{w}{\arg \min }\left \| M(F^{H}DFL_{\chi}^{-1}w-b)\right \|^{2}+\lambda\left \| W_{g}GL_{\chi}^{-1}w \right \| _{0} \; with\: G=\left [ G_{x};G_{y};G_{z}\right ]$$

where the susceptibility distribution can be calculated by $$$\chi=L_{\chi}^{-1}w_{opt}$$$, $$$M$$$ is a binary mask of the reliable phase inside ROI, $$$F$$$ is the discrete Fourier transform operator, $$$D$$$ is the dipole kernel in k-space, $$$b$$$ is the total magnetic field shift, $$$\lambda$$$ is a regularization parameter, $$$W_{g}$$$ is the weighting matrix, and $$$G$$$ is the gradient operator in r-space. When prior information concerning the susceptibility distribution $$$\chi$$$ is available, this derivation process helps to introduce a reasonable regularization operator $$$L_{\chi}^{-1}$$$, which greatly speeds up the convergence of the algorithm. In this work, $$$L_{\chi}^{-1}$$$ is constructed as follow:

$$L_{\chi }^{-1}=\left\{\begin{matrix}1,M\cap M_{R_{2}^{*}} \\ Q,otherwise\end{matrix}\right.\: M_{R_{2}^{*}}=binary(R_{2}^{*}< th)$$

where $$$Q$$$ and $$$th$$$ are constant. The matrix $$$L_{\chi}^{-1}$$$ is structured such that the difference of the matrix $$$\Gamma_{\chi }=L_{\chi}^{-1}(L_{\chi}^{-1})^{T}$$$ between voxels inside and outside ROI is roughly equivalent to the difference between the local and the background susceptibility distributions. After introducing an auxiliary variable $$$A$$$, a residual term $$$\eta$$$ and a regularization parameter $$$\mu$$$, we can obtain

$$\left\{\begin{matrix}w_{t+1}=\underset{w}{\arg \min }\left \|M(F^{H}DFL_{\chi }^{-1}w-b) \right \|_{2}^{2}+\mu \left \| A_{t}-W_{g}GL_{\chi }^{-1}w-\eta _{t} \right \|_{2}^{2}\\ A_{t+1}=\underset{A}{\arg \min }\: \mu \left \|A- W_{g}GL_{\chi }^{-1}w_{t+1}-\eta _{t} \right \|_{2}^{2}+\lambda \left \|A \right \|_{0}\\ \eta _{t+1}=W_{g}GL_{\chi }^{-1}w_{t+1}+\eta _{t}-A_{t+1}\\ \chi _{opt}=L_{\chi}^{-1}w_{opt}\end{matrix}\right.$$

Gadolinium phantom data, human healthy brain data, and pathological brain data were utilized to demonstrate the performance of the proposed method. The phase was unwrapped using a region growing algorithm. For three-step QSM methods, we removed the background field by using SHARP.


Figure 1 shows the phantom QSMs reconstructed using six different methods. The linear regressions between the estimated (y) and reference (x) balloon susceptibilities were y = 0.847x+0.026 (R2=0.998) for COSMOS, y = 0.866x+0.017 (R2=0.999) for MEDI, y = 0.771x+0.011 (R2=0.996) for iLSQR, y = 0.829x+0.046 (R2=0.993) for WL1, y = 0.909x+0.016 (R2=0.999) for TFI, and y = 0.919x+0.019 (R2=0.996) for fast TFI.

The three plane views of the QSM results from the Cornell healthy brain are displayed in Figure 2. The fast TFI method is quantitatively compared with three three-step QSM methods (MEDI1, iLSQR2, and WL13) and a one-step QSM method (SSQSM4 and TFI5). The COSMOS6 result is presented as a ground truth. The average reconstruction time was 1394 s for MEDI, 152 s for iLSQR, 88 s for WL1, 819 s for SSQSM, 986 s for TFI, and 471 s for fast TFI. SSQSM, TFI, and fast TFI are all one-step QSM methods, but fast TFI is fastest among them. Compared to the results from the other methods, the fast TFI result reveals excellent fidelity to the COSMOS map for all three plane views. The susceptibility value from TFI in globus pallidus deviates from the ground truth. In particular, the coronal and sagittal views of the WL1 and SSQSM results depict streaking artifacts that deteriorate the image quality.

Figure 3 shows the QSM results of a male patient with cerebral hematoma. It is obvious that the reconstructed results of all methods other than fast TFI obviously fail to restore the susceptibility distribution due to violent artifacts around the hemorrhagic area. The proposed fast TFI, benefited from the advanced solution and the prior known matrix $$$L_{\chi}^{-1}$$$, is effective in identifying the injury and successfully reconstructs the distribution.


In this paper, we proposed a fast QSM reconstruction method under the TFI frame, named fast TFI. The proposed method performs well on phantom, healthy and pathological brain data with a short reconstruction time. In comparison with three-step QSM methods, our experimental results verify the special advantages of the fast TFI in protecting sample details, eliminating artifacts and restoring lesion areas. Moreover, a set of fixed parameters is used in QSM reconstruction, which would greatly promote the clinical applications of QSM.


This work was supported in part by the National Natural Science Foundation of China under Grant 11775184, and Science and Technology Project of Fujian Province of China under Grant 2016Y0078.


1. Liu J, Liu T, de Rochefort L, et al. Morphology enabled dipole inversion for quantitative susceptibility mapping using structural consistency between the magnitude image and the susceptibility map. Neuroimage, 2012;59(3):2560-2568.

2. Li W, Wang N, Yu F, et al. A method for estimating and removing streaking artifacts in quantitative susceptibility mapping. Neuroimage, 2015;108:111-122.

3. Bilgic B, Fan AP, Polimeni JR, et al. Fast quantitative susceptibility mapping with L1-regularization and automatic parameter selection. Magn Reson Med, 2013;72(5):1444-1459.

4. Chatnuntawech I, McDaniel P, Cauley SF, et al. Single-step quantitative susceptibility mapping with variational penalties. NMR Biomed, 2017;30(4):e3570.

5. Liu Z, Kee Y, Zhou D, et al. Preconditioned total field inversion (TFI) method for quantitative susceptibility mapping. Magn Reson Med, 2016;78(1):303-315.

6. Liu T, Spincemaille P, de Rochefort L, et al. Calculation of susceptibility through multiple orientation sampling (COSMOS): a method for conditioning the Inverse problem from measured magnetic field map to susceptibility source image in MRI. Magn Reson Med, 2009;61(1):196-204.


Figure 1 QSM results of the gadolinium phantom reconstructed with COSMOS (A), MEDI (B), iLSQR (C), WL1 (D), TFI (E) and fast TFI (F) methods.

Figure 2 Coronal, sagittal, and axial views of the QSM results of the Cornell healthy brain reconstructed with COSMOS (A), MEDI (B), iLSQR (C), WL1 (D), SSQSM (E), TFI (F), and fast TFI (G) methods.

Figure 3 Results of the patient with cerebral hematoma. (A) Amplitude image, (B) R2* map; (C) total field map; (D-H) QSM results reconstructed with MEDI (D), iLSQR (E), WL1 (F), SSQSM (G), TFI (H), and fast TFI (I) methods.

Proc. Intl. Soc. Mag. Reson. Med. 27 (2019)