1961

Primal-Dual Implementation for Quantitative Susceptibility Mapping (QSM)
Youngwook Kee1, Kofi Deh1, Alexey Dimov1,2, Pascal Spincemaille1, and Yi Wang1,2

1Weill Cornell Medical College, New York, NY, United States, 2Cornell University, Ithaca, NY, United States

Synopsis

We investigate the computational aspects of the prior term in the field-to-susceptibility inversion problem for QSM. Providing a spatially continuous formulation of the problem, we analyze 1) its Euler-Lagrange equation that appears degeneracy and 2) the Gauss-Newton conjugate gradient (GNCG) algorithm that employs numerical conditioning. We propose a primal-dual (PD) formulation that avoids such degeneracy and use the Chambolle-Pock algorithm to solve this alternative formulation; thus numerical conditioning is not required. The two methods were tested and validated on numerical/gadolinium phantoms and ex-vivo/in-vivo MRI data. The PD formulation with the Chambolle-Pock algorithm was faster and more accurate than GNCG.

Purpose

We propose a primal-dual (PD) implementation for QSM than the classical Gauss-Newton conjugate gradient (GNCG) algorithm. It turns out that PD is not only faster but is also more accurate than GNCG. This is because GNCG employs numerical conditioning (hence modifies the prior term) which also affects the spectrum of eigenvalues which is closely related to the speed of convergence.

Theory

1) Energy Functional

In a spatially continuous setting, the problem of morphology-enabled dipole inversion (MEDI) (1) can be expressed as follows:

$$ \text{Find } \hat \chi \text{ such that } \hat \chi \in \operatorname{Argmin}\left\{\chi \in BV(\Omega; \mathbb{R}): \int_\Omega |M\nabla \chi|_1 + \frac{\lambda}{2}\int_\Omega |w(d*\chi-b)|^2\right\}$$

Here, the susceptibility distribution $$$\chi:(\Omega\subset\mathbb{R}^3)\to \mathbb{R}$$$ is estimated from the measured local field $$$b:\Omega\to\mathbb{R}$$$. The space of bounded variation, $$$BV(\Omega;\mathbb{R})$$$ , allows the existence of jumps modeling sharp edges. The map $$$w:L^1\to L^1$$$ compensates for the non-uniform phase noise (SNR weighting), and the edge indicator map $$$M:[L^1]^3 \to [L^1]^3$$$ is derived from the corresponding magnitude image.

2) GNCG

GNCG assumes that $$$\chi$$$ is differential, which allows us to derive the Euler-Lagrange (EL) equation associated with the energy functional as follows:

$$ \sum_{\ell=1}^3 \frac{\partial}{\partial r_\ell}\left( \frac{M_{r_\ell(r)}\partial_{r_\ell}\chi(r)}{|M_{r_\ell(r)}\partial_{r_\ell}\chi(r)|} \right) - \lambda (w(d*))^H(w(d*\chi-b))=0, \quad \forall r=(r_1, r_2, r_3) \in int\,\Omega, $$

where $$$(w(d*))^H$$$ is the adjoint of $$$w(d*)$$$. The EL equation, however, exhibits degeneracy when $$$|M_{r_\ell(r)}\partial_{r_\ell}\chi(r)|=0$$$ so that numerical conditioning is required for computational purposes: $$$|M_{r_\ell(r)}\partial_{r_\ell}\chi(r)|$$$ is replaced with $$$\sqrt{|M_{r_\ell(r)}\partial_{r_\ell}\chi(r)|^2 +\epsilon}$$$ for $$$\ell=1,2,3$$$ and $$$\epsilon>0$$$. As a consequence, this conditioning deviates from the definition of total variation, which leads to solving a slightly different objective.

3) PD

We can derive a generalized definition of anisotropic weighted total variation as follows:

$$\int_\Omega |M\nabla\chi|_1 = \sup \left\{ -\int_\Omega \chi(r)\operatorname{div}(M^H(r)\xi(r))dr : \xi \in C_c^1(\Omega;\mathbb{R}^3), ||\xi(r)||_\infty \leq 1, \forall r \in \Omega \right\}$$,

where $$$C_c^1$$$ is the set of continuously differentiable functions with compact support. $$$\xi$$$ is the so-called dual variable with the constraint of pointwise maximum norm. Note that the discontinuities are allowed in $$$\xi$$$ in the above dual formulation. Additional dualization for the data fidelity term (2) leads us to the following primal-dual formulation:

$$ \hat \chi \in \operatorname{Argmin} \left\{ \operatorname{sup} \left\{ (\xi, \nu) \in C_c^1 \times L^2 : -\int_\Omega \chi \operatorname{div}(M^H\xi) + \langle w(d*\chi-b), \nu \rangle_{L^2} - \iota_P(\xi) - \frac{1}{2\lambda} ||\nu||_2^2 \right\} \right\}$$,

where the indicator function $$$\iota(\xi)$$$ of a set $$$P\subset C_c^1$$$ is defined as $$$\iota_P(\xi)=0$$$ if $$$\xi \in P$$$, $$$+\infty$$$ otherwise. Here, the convex set $$$P=\{\xi\in C_c^1(\Omega;\mathbb{R}^3) : ||\xi(r)||_\infty \leq 1, \forall r \in \Omega\}$$$. This saddle-point problem can be solved by the Chambolle-Pock algorithm (3). Note that GNCG now can be viewed as the primal approach.

Methods

A simulated brain, in vivo brain MRI data in 4 healthy subjects, and ex vivo brain blocks were used for experimental validation of the two methods (GNCG and PD) using as reference the ground truth when available and COSMOS (4) for the healthy subjects and brain blocks. One set of the in vivo brain MRI data was downloaded at http://weill.cornell.edu/mri/pages/qsm.html. The rest of the in vivo data were extracted from (5). For phase preprocessing, see (5) for details. For the excised brain blocks, see the caption of Fig. 2 for details.

Results

Convergence behavior of the two methods on a simulated brain is seen in Fig. 1.

Figs. 2 and 3 show the solutions obtained by the GNCG and PD algorithms and results of quantitative analysis for ex vivo brain tissues and in vivo brain data, respectively.

Discussion

The accuracy of a minimizing sequence in GNCG can be improved by decreasing $$$\epsilon$$$ from $$$10^{-6}$$$ to $$$10^{-8}$$$, but at the cost of significantly slowing down the convergence rate of the algorithm. This is because a linear system solved by the nested CG loop becomes ill-conditioned. The rate of convergence of the CG iteration is determined by the location of the eigenvalues (spectrum) of a linear system (6), and the value of the conditioning parameter $$$\epsilon$$$ affects the location of the spectrum. From a computational perspective, therefore, the algorithm would become impractical if $$$\epsilon$$$ is set to a very small number $$$\sim 10^{-10}$$$ to achieve better accuracy. This makes GNCG less preferable than PD when it comes to non-smooth convex optimization.

Conclusion

We have shown the influence of specific implementation of the regularization term in QSM as it relates to convergence rate and accuracy. The slower convergence rate and accuracy associated with the numerical conditioning required for the standard GNCG algorithm is substantially improved using a primal-dual formulation that no longer requires the modification of the regularization term in QSM.

Acknowledgements

This research was supported by R01NS072370, R01NS090464, and R01NS095562.

References

(1) Liu J, Liu T, de Rochefort L, Ledoux J, Khalidov I, Chen W, Tsiouris AJ, Wisnieff C, Spincemaille P, Prince MR, Wang Y. 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) Bauschke HH, Combettes PL. Convex analysis and monotone operator theory in Hilbert spaces. CMS books in mathematics,. p 1 online resource.

(3) Chambolle A, Pock T. A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision 2011;40(1):120-145.3. Wang Y, Liu T. Quantitative susceptibility mapping (QSM): Decoding MRI data for a tissue magnetic biomarker. Magn Reson Med 2015;73(1):82-101.

(4) Wang Y, Liu T. Quantitative susceptibility mapping (QSM): Decoding MRI data for a tissue magnetic biomarker. Magn Reson Med 2015;73(1):82-101.

(5) Liu T, Liu J, de Rochefort L, Spincemaille P, Khalidov I, Ledoux JR, Wang Y. Morphology enabled dipole inversion (MEDI) from a single-angle acquisition: comparison with COSMOS in human brain imaging. Magn Reson Med 2011;66(3):777-783.

(6) Trefethen LN, Bau D. Numerical linear algebra. Philadelphia: SIAM; 1997.

Figures

(a) At 3000 iterations, the relative errors were 8.93%, 8.84%, and 8.62% for GNCG (ε=10-6), GNCG (ε=10-8), PD, respectively. PD also computed the most accurate solution at 1000 iterations: the relative errors were 8.96%, 9.45%, and 8.67% for GNCG (ε=10-6), GNCG (ε=10-8), PD, respectively. The CPU running times of each algorithm per iteration were approximately 0.75sec. (b) Ground truth (GT), computed susceptibility maps (PD and GNCG (ε=10-6)) are shown in the first row. Difference images between GT and computed solutions from PD and GNCG are displayed in the bottom row. Close-up images correspond to the red bounding box in GT.

(a) Computed susceptibility maps from GNCG and PD on ex vivo data compared to COSMOS. (b) Accuracy assessment via linear regression and Bland-Altman plots between different algorithms and COSMOS. The CPU running times of each algorithm per iteration were 4.15sec and 4.05sec for GNCG and PD, respectively. The total number of iterations was 2000 and 2100 for GNCG was PD, respectively.

(a) Computed susceptibility maps from GNCG and PD on in vivo data compared to COSMOS. (b) Accuracy assessment via linear regression and Bland-Altman plots between different algorithms and COSMOS. For the subject 1 of size 230×230×46, the CPU running times of each algorithm per iteration were approximately 0.29sec. The total number of iterations was 1304 and 1000 for GNCG was PD, respectively.

Proc. Intl. Soc. Mag. Reson. Med. 25 (2017)
1961