Introduction

Global Maxwell Tomography (GMT) is a recently introduced volume-integral-equation-based method for non-invasive estimation of tissue electrical properties (EP), namely electric permittivity and conductivity, from magnetic resonance (MR) measurements [1,2]. GMT is posed as a volumetric inverse problem with respect to the EP. The cost function of GMT is minimized using L-BFGS-B, or L-BFGS with box constraints [3], a quasi-Newton algorithm that forms approximations of the Hessian. Initial numerical experiments with a human-head phantom showed an average error in the estimated permittivity and conductivity of 3.82% and 2.33%, respectively, after 1000 iterations, starting from a homogeneous initial guess. In this work, we present a new implementation of GMT that relies on the Projected Newton algorithm [4], with the goal of improving the accuracy of the EP reconstructions, while reducing the number of iterations.

Theory

We modified the original cost function of GMT [1,2] by removing a square-root. Furthermore, in the proposed implementation we operate over the pointwise fields which means that the magnetic and electric fields are multiplied by the inverse of the Grammian $$$\mathbf{\Delta}$$$:
$$\mathbf{h}_t=\mathbf{h}_i+\mathbf{\Delta}^{-1}\mathbf{K}\mathbf{J}\text{, and}$$
$$\mathbf{A}(\mathbf{\epsilon})\mathbf{J}=\left(\mathbf{I}-\mathbf{\Delta}^{-1}\mathbf{P}(\tfrac{\chi}{\epsilon})\mathbf{N}\right)\mathbf{J}=c_e\mathbf{P}(\tfrac{\chi}{\epsilon})\mathbf{e}_i$$
We solve the following optimization problem:
$$\min_{\epsilon\in\mathbb{C}^N}f(\epsilon)+\alpha{f_r(\epsilon)}\ \text{s.t.}\ \rm{Re}(\epsilon)\geq{1}\ \&\ \rm{Im}(\epsilon)\leq{0}\text{,}$$
where $$$f(\epsilon)$$$ is the data mismatch cost function and $$$f_r(\epsilon)$$$ is the Match Regularization cost function.
The data mismatch cost function operates over estimated and measured (with the hat) $$$\mathbf{b}_1^+$$$:
$$f\!\left(\mathbf{\epsilon},\mathbf{\overline{\epsilon}}\right)=\sum_k\sum_n\|\mathbf{\delta}_{kn}\|_2^2=\sum_k\sum_n\|\mathbf{b}_k\odot\overline{\mathbf{b}_n}-\hat{\mathbf{b}}_k\odot\overline{\hat{\mathbf{b}}_n}\|_2^2$$
where $$$\odot$$$ is the Hadamard product. The two summations in $$$n$$$ and $$$k$$$ iterate over all of the unique field maps of a multiple-channel transmit array. The gradient of the cost function in the CR-calculus sense is given by
$$\left(\frac{\partial f}{\partial\mathbf{\epsilon}}\right)^\ast=\left[\sum_k\mathbf{Q}^{\rm{T}}\!\left(\overline{\mathbf{\psi}_k}\odot\gamma_k\right)\right]\oslash\overline{\mathbf{\epsilon}}^2\text{,}$$
where
$$\mathbf{\psi}_k=\mathbf{\Delta}^{-1}\mathbf{N}\mathbf{J}_k+c_e\mathbf{e}_k^i$$
$$\mathbf{\gamma}_k=\mathbf{\Omega}^{\ast}\left(\sum_n\mathbf{b}_n\odot\mathbf{\delta}_{kn}\right)$$
$$\mathbf{\Omega}=\mu_0\mathbf{F}\mathbf{\Delta}^{-1}\mathbf{K}\mathbf{A}^{-1}\text{.}$$
In order to apply Newton to Global Maxwell Tomography, we require the second derivatives in the CR-calculus sense:
$$\frac{\partial}{\partial{\epsilon}}\frac{\partial f}{\partial\mathbf{\epsilon}}^\ast=\rm{D}^{-2}(\overline{\mathbf{\epsilon}})\mathbf{Q}^{\rm{T}}\sum_k\rm{D}(\overline{\mathbf{\psi}_k})\mathbf{\Omega}^{\ast}\!\left(\sum_n\rm{D}(\mathbf{\delta}_{kn})\mathbf{\Omega}\rm{D}(\mathbf{\psi}_n)+\rm{D}\left(|\mathbf{b}_n|^2\right)\mathbf{\Omega}\rm{D}(\mathbf{\psi}_k)\right)\mathbf{Q}\rm{D}^{-2}(\mathbf{\epsilon})$$
$$ \begin{split}\frac{\partial}{\partial\overline{\mathbf{\epsilon}}}\frac{\partial{f}}{\partial\mathbf{\epsilon}}^\ast= &-2\rm{D}^{-3}(\overline{\mathbf{\epsilon}})\rm{D}\!\left(\sum_k\mathbf{Q}^{\rm{T}}\!\left(\overline{\mathbf{\psi}_k}\odot\gamma_k\right)\right)\\&+\rm{D}^{-2}(\overline{\mathbf{\epsilon}})\sum_k\mathbf{Q}^{\rm{T}}\rm{D}(\overline{\mathbf{\psi}_k})\mathbf{A}^{-\ast}\!\left[\mathbf{N}^{\ast}\rm{D}(\mathbf{\gamma}_k)+\mu_0\mathbf{K}^{\ast}\mathbf{F}^{\ast}\rm{D}(\mathbf{b}_k)\sum_n\rm{D}(\mathbf{b}_n)\overline{\mathbf{\Omega}}\rm{D}(\overline{\mathbf{\psi}_n})\right]\mathbf{Q}\rm{D}^{-2}(\overline{\mathbf{\epsilon}})\\&+\rm{D}^{-2}(\overline{\mathbf{\epsilon}})\sum_k\mathbf{Q}^{\rm{T}}\rm{D}(\mathbf{\gamma}_k)\mathbf{\Omega}^{\rm{T}}\!\left[\sum_n\rm{D}(\overline{\mathbf{\delta}_{kn}})\overline{\mathbf{\Omega}}\rm{D}(\overline{\mathbf{\psi}_n})+\rm{D}\left(|\mathbf{b}_n|^2\right)\overline{\mathbf{\Omega}}\rm{D}(\overline{\mathbf{\psi}_k})\right]\mathbf{Q}\rm{D}^{-2}(\overline{\mathbf{\epsilon}})\end{split}
$$
When it comes to the Hessian matrix-vector product, we treat the Hessian as a real-valued matrix $$$\mathbf{H}\in\mathbb{R}^{2N\times2N}$$$. In order to switch between real- and complex-valued vectors, we define the operator $$$\mathbf{R}\colon\mathbb{C}^{N}\to\mathbb{R}^{2N}$$$ that concatenates the real and imaginary parts of the input and its inverse $$$\mathbf{R}^{-1}\colon\mathbb{R}^{2N}\to\mathbb{C}^{N}$$$. The Hessian matrix-vector product $$$\mathbf{H}\mathbf{v}$$$ is given by
$$\mathbf{H}(\epsilon)\mathbf{v}=\mathbf{R}\!\left( \frac{\partial}{\partial{\epsilon}}\frac{\partial f}{\partial\mathbf{\epsilon}}^\ast z +\frac{\partial}{\overline{\partial{\epsilon}}}\frac{\partial f}{\partial\mathbf{\epsilon}}^\ast\overline{z}\right)\text{,}$$
where $$$z=\mathbf{R}^{-1}\!\left(v\right)$$$.

Methods

We used the Magnetic Resonance Integral Equation (MARIE) simulation tool [5] for the forward problem of GMT, employing the piecewise linear basis functions for the volumetric currents [6]. Additionally, as in the original GMT publication [1], we neglected the coupling between the head and the incident fields, which is expected to have a negligible effect.
We used Projected Newton [4] as the optimizer when estimating the electrical properties. Projected Newton is a modification of Newton's Method that allows for bound constraints on the decision variables. In addition, when solving for the Newton step $$$\mathbf{H}\mathbf{\Delta}\mathbf{\epsilon}=-\mathbf{g}(\mathbf{\epsilon})$$$, we used the minres algorithm for a loose tolerance of $$$10^{-1}$$$ and a maximum of 7 iterations because the Hessian is not positive semi-definite.
We performed numerical experiments on the same four-compartment tissue-mimicking phantom (4,096 voxels at a 6 mm isotropic resolution) and realistic "Billie" head model [7] (25,575 voxels at a 5 mm isotropic resolution) used in previous work [1], with the same parameters for the regularizer (Figure 1). We used $$$8$$$ idealized excitations [1] as incident fields, with a peak SNR of 50 and 200 for the four-compartment and head phantom, respectively. The operating frequency and the electrical properties were set for $$$B_0 = 7$$$Tesla.

Results

Figure 2 compares the initial and final EP estimates with the ground truth, for the four-compartment phantom. Starting from a homogeneous initial guess of EP, GMT converged in 74 iterations. The mean errors were 0.27% and 0.48% for $$$\epsilon_{\rm{R}}$$$ and $$$\sigma$$$, respectively. The maximum errors were 1.05% and 1.72% for $$$\epsilon_{\rm{R}}$$$ and $$$\sigma$$$, respectively.
For the head phantom, starting again from a worst–case homogeneous guess, GMT converged after 103 iterations. The mean errors were 4.5% and 3% for $$$\epsilon_{\rm{R}}$$$ and $$$\sigma$$$, respectively. The maximum errors were 98.6% and 97.7% for $$$\epsilon_{\rm{R}}$$$ and $$$\sigma$$$, respectively. These peak errors were due to one-off voxels in regions of very high contrast, as suggested by Figure 3.

Discussion & Conclusion

We demonstrated a new implementation of GMT that uses a Projected Newton method as the optimizer. Similar results to the original GMT, based on L-BFGS-B method, were obtained in considerably fewer iterations. For the 4-compartment phantom, we saw a reduction in iterations by roughly a factor of two. For the realistic head phantom, the Projected Newton approach converged in 103 iterations, yielding results with comparable mean errors to 1000 L-BFGS-B iterations. The peak errors were considerably lower with Projected Newton: 98.6% vs. 262% for $$$\mathbf{\epsilon}_{\rm{R}}$$$. Additionally, the errors in our new approach were due to one-off voxels, whereas previously the large errors in estimated EP were artifacts that did not correspond to any tissue values. Therefore, for complex geometries, Projected Newton not only drastically reduces the number of iterations, but also yields more accurate and artifact-free solutions. Future work on GMT will try to address one-off errors in regions of extremely high contrast using more sophisticated regularization and Machine Learning methods. In conclusion, this work shows that the application of Projected Newton further cements GMT as one of the most accurate MR-based EP estimation algorithms implemented to date.

Acknowledgements

This work was supported by NIH R01 EB024536 and by NSF 1453675. It was performed under the rubric of the Center for Advanced Imaging Innovation and Research (CAI$$$^2$$$R, www.cai2r.net), an NIBIB Biomedical Technology Resource Center (NIH P41 EB017183).