Introduction

Electrical properties tomography (EPT) reconstructs admittivity, namely conductivity and permittivity, which are useful in oncology and dosimetry, from B1 field data¹. Conventional EPT methods can be divided into two types: partial-differential-equation-based (PDE-based) methods^2,3,4 and integral-equation-based (IE-based) methods^5,6. In this paper, we refer to these methods as local and global methods, respectively. Local methods are sensitive to noise due to numerical calculation of the Laplacian of the measured data, while global methods solve a nonlinear equation, possibly resulting in local minima. We previously proposed a direct and global method solving a linear IE of the electric field⁷. However, this method requires all the three components of the magnetic field $$$(H_{x},H_{y},H_{z})$$$, which are inaccessible in practical situations. In this paper, we propose a novel EPT method solving a linear IE for admittivity using only $$$H^{+}=(H_{x}+\mathrm{i}H_{y})/2$$$. We also propose a phase-based EPT method, in which conductivity is reconstructed from the measured B1 phase⁸, based on the same formulation. The proposed methods retrieve admittivity without calculating the Laplacian of the measured data, nor iteratively updating the estimate.

Method

According to Hafalir et al.², the following convection-reaction (cr) equation for impedivity $$$λ$$$ holds: $$λ{\Delta}H^{+}+{\nabla}λ\cdot\begin{pmatrix}2{\partial}H^{+}\\2\mathrm{i}{\partial}H^{+}\\\partial_{z}H^{+}\end{pmatrix}-\mathrm{i}ω_{0}μ_{0}H^{+}=0,\tag{1}$$ where $$$ω_{0}$$$ is Larmor's frequency, $$$μ_{0}$$$ is permeability of free space, and $$$\partial\equiv(\partial_{x}-\mathrm{i}\partial_{y})/2$$$. By introducing a new vector field $$$λ\nabla^{\mathrm{c}}H^{+}$$$, where $$$\nabla^{\mathrm{c}}\equiv\nabla-\mathrm{i}\nabla\times\hat{\boldsymbol{z}}$$$, Eq. (1) can be rewritten in a simple divergence-constraint form as follows: $$\nabla\cdot(λ\nabla^{\mathrm{c}}H^{+})=\mathrm{i}ω_{0}μ_{0}H^{+}.\tag{2}$$ Equation (2) is equivalent to the PDE derived in [4] to construct the modified finite difference (mFD) method.The cr/mFD method solves Eq. (1)/(2) using the finite element/difference method. In this paper, we derive a linear IE based on Helmholtz's identity by exploiting the fact that Eq. (2) is written in the divergence-constraint form.
Helmholtz's identity⁹ states that a smooth vector field $$$\boldsymbol{f}$$$ defined in a bounded region of interest $$$Ω$$$ can be represented by its divergence and curl as follows: $$\boldsymbol{f}(\boldsymbol{r})=-\oint_{{\partial}Ω}{\nabla}G(\boldsymbol{r}^{\prime};\boldsymbol{r})(\boldsymbol{n}\cdot\boldsymbol{f})\mathrm{d}S^{\prime}+\int_{Ω}{\nabla}G(\boldsymbol{r}^{\prime};\boldsymbol{r})(\nabla\cdot\boldsymbol{f})\mathrm{d}V^{\prime}+\oint_{{\partial}Ω}{\nabla}G(\boldsymbol{r}^{\prime};\boldsymbol{r})(\boldsymbol{n}\times\boldsymbol{f})\mathrm{d}S^{\prime}-\int_{Ω}{\nabla}G(\boldsymbol{r}^{\prime};\boldsymbol{r})\times(\nabla\times\boldsymbol{f})\mathrm{d}V^{\prime},\tag{3}$$ where $$$G(\boldsymbol{r}^{\prime};\boldsymbol{r})=1/(4π\lvert\boldsymbol{r}^{\prime}-\boldsymbol{r}\rvert)$$$. By applying integration by parts to the last term of Eq. (3), we obtain $$\boldsymbol{f}(\boldsymbol{r})=-\oint_{{\partial}Ω}{\nabla}G(\boldsymbol{r}^{\prime};\boldsymbol{r})(\boldsymbol{n}\cdot\boldsymbol{f})\mathrm{d}S^{\prime}+\int_{Ω}{\nabla}G(\boldsymbol{r}^{\prime};\boldsymbol{r})(\nabla\cdot\boldsymbol{f})\mathrm{d}V^{\prime}+\int_{Ω}(\boldsymbol{f}\times\nabla)\times{\nabla}G(\boldsymbol{r}^{\prime};\boldsymbol{r})\mathrm{d}V^{\prime},\tag{4}$$ which represents a vector field using its divergence and itself instead of its curl. Here, substituting Eq. (2) into Eq. (4) yields $$λ\nabla^{\mathrm{c}}H^{+}=-\oint_{\partial Ω}{\nabla}G(\boldsymbol{r}^{\prime};\boldsymbol{r})(\boldsymbol{n}{\cdot}λ\nabla^{\mathrm{c}}H^{+})\mathrm{d}S^{\prime}+\int_{Ω}{\nabla}G(\boldsymbol{r}^{\prime};\boldsymbol{r})\mathrm{i}ω_{0} μ_{0}H^{+}\mathrm{d}V^{\prime}+\int_{Ω}(λ\nabla^{\mathrm{c}}H^{+}\times\nabla)\times{\nabla}G(\boldsymbol{r}^{\prime};\boldsymbol{r})\mathrm{d}V^{\prime}.\tag{5}$$ Given the EP values on $$${\partial}Ω$$$, the surface integrals can be calculated. Equation (5) is a linear IE for impedivity including only first-order derivatives of $$$H^{+}$$$, which can be calculated using the Savitzky-Golay filter¹⁰, and can be efficiently solved by the FFT-CG method¹¹. In contrast to the previous method^7,12 in which we derived an IE for $$$\boldsymbol{E}=λ\nabla\times\boldsymbol{H}$$$ which requires all components of $$$\boldsymbol{H}$$$, the proposed method considers $$$λ\nabla^{\mathrm{c}}H^{+}$$$ including only $$$H^{+}$$$ measured by the B1 mapping.
A similar formulation can be derived for phase-based EPT. According to Gurler et al.¹³, the following cr equation for resistivity, $$$ρ$$$, holds: $$ρ{\Delta}ϕ+{\nabla}ρ\cdot{\nabla}ϕ-2ω_{0}μ_{0}=0,\tag{6}$$ where $$$ϕ$$$ is the measured B1 phase. By introducing a vector field $$$ρ{\nabla}ϕ$$$, Eq. (6) can also be rewritten in divergence-constraint form as follows: $$\nabla\cdot(ρ{\nabla}ϕ)=2ω_{0}μ_{0}.\tag{7}$$ Substituting Eq. (7) into Eq. (2) yields $$ρ{\nabla}ϕ=-\oint_{{\partial}Ω}{\nabla}G(\boldsymbol{r}^{\prime};\boldsymbol{r})(\boldsymbol{n}{\cdot}ρ{\nabla}ϕ)\mathrm{d}S^{\prime}+\int_{Ω}{\nabla}G(\boldsymbol{r}^{\prime};\boldsymbol{r})2ω_{0}μ_{0}\mathrm{d}V^{\prime}+\int_{Ω}(ρ{\nabla}ϕ\times\nabla)\times{\nabla}G(\boldsymbol{r}^{\prime};\boldsymbol{r})\mathrm{d}V^{\prime}.\tag{8}$$ Equation (8) is a linear IE for resistivity including only first-order derivatives of $$$ϕ$$$, and can also be solved by the FFT-CG method.
We tested the proposed method by a numerical simulation and a phantom experiment and compared the obtained results with those of the cr methods. We constructed a sphere model for the simulation and a cylindrical phantom for the experiment as shown in Fig. 1. The simulation was performed by COMSOL Multiphysics 5.5 (COMSOL Inc.) and $$$H^{+}$$$ was obtained with a resolution of $$$1.4\times1.4\times5~\mathrm{mm^{3}}$$$. The experiment was conducted by the 3 T MR scanner Magnetom Prisma (Siemens) at the University of Tokyo and $$$H^{+}$$$ was obtained with a resolution of $$$1.4\times1.4\times10~\mathrm{mm^{3}}$$$. For the B1 mapping, we used the double angle method¹⁴.

Results and Discussion

Figure 2 shows the reconstructed conductivity and permittivity maps of the sphere model when no noise was added. The cr method gives the results with spurious oscillations arising from the unstable nature of the cr equation¹⁵. By introducing viscosity-type regularization¹⁵, the cr method yields stable but blurred results. In contrast, the proposed method yields accurate and stable results even without additional regularization. The proposed method has the least L2 error as shown in the table. Figure 3 shows the reconstruction results when 1% Gaussian noise was added. In the conventional method, strong regularization is required to stabilize the reconstruction, resulting in over-smoothed estimates. The proposed method yields stable results while preserving the edge. Figure 4 shows the conductivity maps reconstructed by phase-based EPT. Similar to EPT, the proposed method yields a more stable result than the conventional method while preserving the edge. Figure 5 shows the reconstructed conductivity maps of the cylindrical phantom. The reference map obtained by the conductivity meter (Hanna Instruments, HI8733) is shown in Fig 5(a). The proposed EPT method successfully reconstructed conductivity and a similar map is obtained by the proposed phase-based EPT method using only the B1 phase.

Conclusion

We proposed novel EPT and phase-based EPT methods solving linear IEs for impedivity and resistivity derived from Helmholtz's identity. The proposed method could reconstruct EPs directly without second-order derivatives of the measured data. Numerical simulations and phantom experiments showed that the proposed method could retrieve EPs more robustly than the conventional method solving a linear PDE.

Acknowledgements

This work was financially supported by Canon Medical Systems Corporation.