A mathematical and experimental analysis is proposed to investigate the effect of discretization of the gradient and divergence operators in Morphology Enabled Dipole Inversion (MEDI) for Quantitative Susceptibility Mapping (QSM). We show that the use of central differences leads mathematical consistencies that are absent using forward and backward differences. Experimentally, this leads to a suppression of checkerboard patterns and improves QSM reconstruction.

Let $$$\Omega \subset \mathbb{R}^3$$$ be a bounded open set; and let $$$\chi:\Omega \to \mathbb{R}$$$ be a susceptibility map. Then, a continuous formulation of MEDI for QSM is given by

$$\underset{\chi \in {\cal V}}{\operatorname{minimize}} \,\, \frac{\lambda}{2}\int_\Omega |w(d*\chi - b)|^2 + \int_\Omega g|\nabla \chi|,$$

where $$$w$$$ is a SNR weighting function, $$$d$$$ is the unit dipole kernel, $$$b$$$ is an extracted local field from a measured total field by background removal, and $$$g$$$ is a morphologically determined weighting function. A proper solution space $$$\cal V$$$ allows $$$\chi$$$ to have the set of discontinuities which models edges in susceptibility. Indeed, the energy functional is defined when the number of voxels in $$$\Omega$$$ is uncountably many. Then, the corresponding Euler$$$-$$$Lagrange equation is given by

$$\frac{\lambda}{2}\frac{\partial}{\partial \chi}\int_\Omega |w(d*\chi - b)|^2 - \operatorname{div}\left(g \frac{\nabla \chi}{|\nabla \chi|}\right)=0,$$

which is a parabolic partial differential equation. Note that the matrix formulation in MEDI is a discrete version of this continuous formulation where the discretization schemes are needed for the differential operators. Therefore, the accuracy of a computed solution$$$-$$$susceptibility map$$$-$$$substantially relies on how these operators are discretized. The lagged diffusivity fixed point iteration [3] is used in MEDI with numerical regularization where $$$|\nabla \chi|$$$ is replaced by $$$|\nabla \chi|_\epsilon := \sqrt{|\nabla \chi|_2^2 + \epsilon}$$$, $$$\epsilon >0$$$.

Next, we show that using the forward difference for $$$\nabla$$$ and backward difference for $$$\operatorname{div}$$$ [4] leads to consistent approximation of $$$\operatorname{div}(\nabla)$$$. The resulting approximation error at $$$(x,y,z) \in \Omega$$$ is $$\left| \operatorname{div}(\nabla \chi) - \left( \frac{\chi(x+h_x, y, z) - 2\chi(x, y, z) + \chi(x-h_x, y, z)}{4h_x^2} + \frac{\chi(x, y+h_y, z) - 2\chi(x, y, z) + \chi(x, y-h_y, z)}{4h_y^2} \right. \right. \nonumber\\\left. \left. + \frac{\chi(x, y, z+h_z) - 2\chi(x, y, z) + \chi(x, y, z-h_z)}{4h_z^2} \right) \right| = O(h_x^2) + O(h_y^2) + O(h_z^2) \to 0 \quad \text{as} \quad (h_x, h_y, h_z) \to 0.$$ Clearly, the above discrete approximation is consistent as it converges to $$$\operatorname{div}(\nabla \chi)$$$.

Let us take a look at its corresponding scheme implemented in MEDI [5] which implements central difference for both gradient and divergence. The approximation error can be computed by setting $$$h_x =2\delta_x$$$, $$$h_y =2\delta_y$$$, and $$$h_z =2\delta_z$$$ leading to: $$\left| \operatorname{div}(\nabla \chi) - \frac{1}{4}\left( \frac{\chi(x+2\delta_x, y, z) - 2\chi(x, y, z) + \chi(x-2\delta_x, y, z)}{4\delta_x^2} + \frac{\chi(x, y+2\delta_y, z) - 2\chi(x, y, z) + \chi(x, y-2\delta_y, z)}{4\delta_y^2} \right. \right. \nonumber\\\left. \left. + \frac{\chi(x, y, z+2\delta_z) - 2\chi(x, y, z) + \chi(x, y, z-2\delta_z)}{4\delta_z^2} \right) \right| = O(4\delta_x^2) + O(4\delta_x^2) + O(4\delta_x^2) \to 0 \quad \text{as} \quad (\delta_x, \delta_y, \delta_z) \to 0.$$ This shows the asymptotic inconsistency as it converges to $$$4\operatorname{div}(\nabla \chi)$$$ so that a numerical implementation based on this scheme may lead to artifacts.

The numerical phantom and one human in vivo data set from the MEDI Toolbox [1] were reconstructed using both difference schemes. A value of lambda corresponding to strong regularization was selected to emphasize the effect of the choice of difference scheme. An additional comparison was performed using rat brain GRE data acquired at 7T (Bruker, Billerica, MA).As can be seen in Figures 1, 2 and 3, susceptibility maps reconstructed using an inconsistent discretization suffers from checkerboard artifacts, which are suppressed when using the consistent discretization of the divergence and gradient operators.We have shown that the consistent method for the differential operators drastically reduces artifacts. On the other hand, such an inconsistent scheme seems to have been derived by applying the central difference method to the operators based on a Cartesian grid. While a similar issue was pointed out in [6], a different and more complicated approximation for the divergence was used. This work also presents a simple mathematical justification for the specific discretizations used.