Quantitative Susceptibility Mapping (QSM) strives to estimate the underlying tissue magnetic susceptibility from gradient echo MRI phase data. In order to reconstruct the susceptibility map an ill-conditioned dipole inversion has to be solved, after performing phase unwrapping and background field estimations. Recently, a number of regularization algorithms (e.g. MEDI, HEIDI)^1-5 have been proposed to address this issue. While MEDI¹ and HEIDI² produce high quality and to a certain degree quantitatively accurate susceptibility maps^1,2,6, they suffer from reconstruction artifacts and long reconstruction times (~20 min)^1,2,5. On the other hand, the L1-Split-Bregman (SB) method⁷ offers fast reconstruction times (< 5 min), but suffers from over-smoothing that could render visual inspection of some ROIs infeasible. Herein, we propose a novel two-step inversion algorithm with faster reconstruction times than SB, while producing qualitatively appealing susceptibility maps with a lower root-mean-square error (RMSE) and a higher coefficient of determination (R²) than MEDI and HEIDI when compared to COSMOS.

In the interest of facilitating a comparison of the proposed method with the previously mentioned state-of-the-art algorithms, reconstruction and analysis was carried out using the dataset of Wang and Liu's summary paper⁵, which the authors generously posted online at http://weill.cornell.edu/mri/pages/qsmreview.html.
The tissue magnetic susceptibility $$$\chi$$$ can be calculated from the measured field perturbations $$$\phi$$$ by inverting: $$$F^{-1} D F \chi = \phi$$$, where $$$F$$$ means Fourier Transform, and $$$D$$$ is the unit dipole. In order to address the ill-conditioned k-space regions of the unit dipole during this inversion, we propose an incremental inversion scheme^2-4. First the k-space is divided into a well-conditioned and an ill-conditioned region. Next, a Gaussian filter ($$$\sigma = 0.55$$$) is applied to the field map, as the ill-conditioned inversion is highly susceptible to noise. Then, the proposed method solves the well-conditioned region using a LSMR solver⁸ to estimate $$$DF \chi = F \phi$$$, exiting the iterations early to avoid streaking artifacts.
Once the well-conditioned region has been reconstructed the following total variation model will be solved to fill in the ill-conditioned k-space region:
$$\chi^* = argmin_{\chi} ||\chi||_{TV} + \frac{\mu}{2} ||MF\chi - MF\chi_{LSMR}||_2^2,$$
where $$$M$$$ is a binary mask according to $$$|D| < \delta$$$, and $$$||.||_{TV}$$$ is the anisotropic TV norm. The minimization was solved using an augmented lagrangian method utilizing the alternating direction method (ADM) adopted from Chan et al.⁹.

For the reconstruction of the well-conditioned k-space region the stopping criteria was empirically determined to be 4 iterations. For the augmented lagrangian method all parameters, except the regularization parameter $$$\mu$$$, were set to the default values presented in the original publication ($$$\rho_r = 2, \gamma = 0.7$$$)⁹, as they performed robustly in practice. The regularization parameter $$$\mu$$$ was empirically determined to be $$$\mu = 2500$$$, and the threshold parameter $$$\delta = 0.15$$$. The stopping criteria was set to $$$\epsilon = 10^{-3}$$$.

Fig. 1 shows the reconstructed susceptibility maps. The RMSE, restricted to voxels inside the ROI, of the proposed method with respect to COSMOS is 2.42x10^-2 ppm. Compared to MEDI (3.07x10^-2 ppm) and HEIDI (2.98x10^-2 ppm), the proposed method yields a slight improvement. Further, the R² value with respect to COSMOS of the proposed method is 0.61, showing again a slight improvement compared to MEDI (0.59) and HEIDI (0.55).

The reconstruction of the well-conditioned k-space region took 1.9 sec, and solving the minimization took 29.7 sec for a total reconstruction time of 31.6 sec. Compared to MEDI (1008 sec)⁵ and HEIDI (715 sec)⁵ the proposed method is more than 20x faster, while producing qualitatively comparable susceptibility maps. The matrix size was 256x256x146.

As demonstrated, the proposed inversion algorithm outperforms state-of-the-art algorithms in terms of reconstruction speed and produces qualitatively comparable susceptibility maps. Nevertheless, further validation is needed in order to determine the extent of the quantitative accuracy. The immense improvement in reconstruction speed is in part due to the fact that no a priori information of the underlying susceptibility distribution is used. Another advantage of omitting a priori information is the reduction of reconstruction artifacts in the ill-conditioned k-space region as seen in Fig. 2. On the other hand, not incorporating a priori information will result in underestimated susceptibility values, due to the missing information along the two conical surfaces of the unit dipole. However, while the proposed method omitted the use of any a priori information, it can be trivially modified to incorporate it.