Reconstruction of dual-frequency conductivity by optimization of phase map in MREIT and MREPT

Abstract

Background

The spectroscopic conductivity distribution of tissue can help to explain physiological and pathological status. Dual frequency conductivity imaging by combining Magnetic Resonance Electrical Property Tomography (MREPT) and Magnetic Resonance Electrical Impedance Tomography (MREIT) has been recently proposed. MREIT can provide internal conductivity distributions at low frequency (below 1 kHz) induced by an external injecting current. While MREPT can provide conductivity at the Larmor frequency related to the strength of the magnetic field. Despite this potential to describe the membrane properties using spectral information, MREPT and MREIT techniques currently suffer from weak signals and noise amplification as they both reply on differentiation of measured phase data.

Methods

We proposed a method to optimize the measured phase signal by finding weighting factors according to the echo signal for MREPT and MREIT using the ICNE (Injected current nonlinear encoding) multi-echo pulse sequence. Our target weights are chosen to minimize the measured noise. The noise standard deviations were precisely analyzed for the optimally weighted magnetic flux density and the phase term of the positive-rotating magnetic field. To enhance the quality of dual-frequency conductivity images, we applied the denoising method based on the reaction-diffusion equation with the estimated noise standard deviations. A real experiment was performed with a hollow cylindrical object made of thin insulating film with holes to control the apparent conductivity using ion mobility and an agarose gel cylinder wrapped in an insulating film without holes to show different spectroscopic conductivities.

Results

The ability to image different conductivity characteristics in MREPT and MREIT from a single MR scan was shown by including the two objects with different spectroscopic conductivities. Using the six echo signals, we computed the optimized weighting factors for each echo. The qualities of conductivity images for MREPT and MREIT were improved by optimization of the phase map. The proposed method effectively reduced the random noise artifacts for both MREIT and MREPT.

Conclusion

We enhanced the dual conductivity images using the optimally weighted magnetic flux density and the phase term of positive-rotating magnetic field based on the analysis of the noise standard deviations and applying the optimization and denoising methods.

Background

The conductivity spectra of biological tissues can provide diagnostic medical information from the estimation of physiological and pathological conditions of in-vivo and ex-vivo tissue. However it is difficult to produce high resolution conductivity images inside the human body [1, 2]. The conventional conductivity imaging methods have limited spatial resolution and sensitivity inherited from the ill-posed nature of the problem [3]. In order to achieve sensitive conductivity images with high resolution, electric impedance imaging techniques based on a magnetic resonance imaging (MRI) have been propsed. These include magnetic resonance electrical impedance tomography (MREIT) and magnetic resonance electrical property tomography (MREPT) which are under active investigation [4–17]. Both methods use the internal magnetic field information obtained from the phase data of MRI scanner to reconstruct the internal conductivity image.

MREPT can provide the electrical conductivity information at the Larmor frequency by measuring the phase of positive rotating field due to the applied B ₁ field. It does not require application of external current and therefore directly recovers the conductivity distribution by taking the derivative of the measured phase signals twice. MREIT needs a pair of electrodes to inject current into an imaging object during MRI scan for measuring a magnetic flux density induced by the external injecting current. The MREIT technique used only the z-component of magnetic flux density, B _z , of B=(B _x ,B _y ,B _z ) to reconstruct the cross-sectional apparent conductivity image at a lower frequency range (below 1 kHz) [18–22]. The phase difference approach with an interleaved encoding scheme was adopted to cancel the systematic artifacts accumulated in phase signals and also reduce the random noise artifacts. Recently, a simultaneous conductivity imaging technique using a combination of MREPT and MREIT was proposed to provide the dual-frequency conductivities of tissue from a common MR scan [23]. MREPT using B ₁-mapping technique visualizes the conductivity and permittivity distributions at the Larmor frequency and MREIT recovers the apparent conductivity distribution when injecting low-frequency external current through the attached electrodes. Since the biological tissues show the frequency dependent conductivity property [1, 2], the simultaneous dual-frequency conductivity imaging using a single MR scan is beneficial to provide distinct electrical features of tissues quantitatively.

Both conductivity imaging techniques use the phase signals of measured MR data. MREIT and MREPT commonly suffer from weak signals and noise amplification by using the derivative of the measured signal. The noise level of B _z in MREIT is inversely proportional to the signal-to-noise ratio (Υ^j) of the MR magnitude image and the current injection pulse width. Because of the small amount of injection current and the poor quality of measured B _z change due to the injected current, it is difficult to perform in vivo human experiments using a conventional MR pulse sequence. To enhance the magnetic flux density due to the injected current in MREIT, the injected current nonlinear encoding (ICNE) method was introduced. This extended the duration of injecting current until the end of a readout gradient, and improved the signal by using a multi-echo train MR pulse sequence [24–26]. MREPT also suffers from low sensitivity due to the inherently poor signal to noise ratio and noise sensitive characteristics as it also needs the derivative of measured data, is very sensitive to the measured noise.

In this paper, we adopt a multiple spin echo MREIT pulse sequence based on the ICNE scheme to measure multiple phase data for MREIT and MREPT images. The acquired multiple phase data can be decomposed to the phase term reflecting the magnetic flux density signal induced by the injected current and the other phase term of positive rotating field due to the applied B ₁ field. We analyze the noise level of the two decomposed phase terms and minimize the measured random noise artifacts by applying the optimal combination of multiple phase terms. Also, we apply a denoising technique to the optimized magnetic flux density for MREIT and to the phase signal for MREPT in order to improve the quality of the reconstructed conductivity images. We prepared a conductivity phantom consisted of two different kinds of anomalies to show the difference of MREIT and MREPT. A phantom experiment is conducted to validate that the proposed method is able to improve the qualities of reconstructed dual-frequency conductivity images compared to the results of using conventional MREIT and MREPT reconstruction algorithms.

Methods

Governing equation

We denote an imaging object as Ω. The admittivity in Ω is κ=σ(ω)+i ω ε(ω), where σ(ω) and ε(ω) are the conductivity and permittivity, respectively, at the angular frequency ω. Time-harmonic Maxwell’s equations relate the electric field E and magnetic flux field H:

$$\nabla \times \mathbf{E}\left(\mathbf{r}\right)=-\mathrm{i\omega }{\mu }_0\mathbf{H}\left(\mathbf{r}\right)\text{and}\nabla \times \mathbf{H}\left(\mathbf{r}\right)=\kappa \left(\mathbf{r}\right)\mathbf{E}\left(\mathbf{r}\right),\mathbf{r}\in \Omega ,$$

(1)

where the current J and the electric field intensity E satisfy the relation J=κ E by Ohm’s law.

In MREIT, we can assume that the injecting current is sufficiently low frequency to meet i ω ε(ω)≈0 and neglect i ω μ ₀≈0, Maxwell’s equation (1) satisfies

$$\begin{array}{ll}-{\nabla }^2\mathbf{H}\left(\mathbf{r}\right)& =\nabla \times \nabla \times \mathbf{H}\left(\mathbf{r}\right)=\nabla \times \left({\sigma }_L\right(\mathbf{r}\left)\mathbf{E}\right(\mathbf{r}\left)\right)\\ =\nabla {\sigma }_L\left(\mathbf{r}\right)\times \mathbf{E}\left(\mathbf{r}\right)=\frac{\nabla {\sigma }_L\left(\mathbf{r}\right)}{{\sigma }_L\left(\mathbf{r}\right)}\times (\nabla \times \mathbf{H}(\mathbf{r}\left)\right),\mathbf{r}\in \mathrm{\Omega .}\end{array}$$

(2)

where σ _L =σ(ω) denotes the low-frequency electrical conductivity responding to the externally injecting current. Since the electric field intensity E=−∇u is a gradient form under negligible angular frequency ω, we have

$$\left\{\begin{array}{cccc}\nabla ·({\sigma }_L\nabla u(\mathbf{r}\left)\right)& =& 0& \text{in}\Omega \\ -{\sigma }_L\nabla u\left(\mathbf{r}\right)·\mathbf{n}\left(\mathbf{r}\right)& =& g\left(\mathbf{r}\right)& \text{on}\mathrm{\partial \Omega }\end{array}\right.$$

(3)

where n is the outward normal vector on the surface ∂ Ω and g is the current density on the surface.

In MREPT, by assuming the locally homogeneity, the magnetic field H can be expressed as a simple form:

$$-{\nabla }^2\mathbf{H}\left(\mathbf{r}\right)=\nabla \kappa \left(\mathbf{r}\right)\times \mathbf{E}\left(\mathbf{r}\right)-\mathrm{i\omega }{\mu }_0\kappa \left(\mathbf{r}\right)\mathbf{H}\left(\mathbf{r}\right)=-\mathrm{i\omega }{\mu }_0\kappa \left(\mathbf{r}\right)\mathbf{H}\left(\mathbf{r}\right),\mathbf{r}\in \mathrm{\Omega .}$$

(4)

The phase terms retrieved from the measured phase are the combination of the phase terms of positive $H^{+}=\frac{H_x+{\mathit{\text{iH}}}_y}{2}$ and negative $H^{-}=\frac{H_x-{\mathit{\text{iH}}}_y}{2}$ rotating field. For a restricted situation such as the usage of transmit-receive coil, the phase term of positive rotating field can be retrieved from the measured phase. The transverse field of H can be decomposed into the positively rotating field $H^{+}=\frac{H_x+{\mathit{\text{iH}}}_y}{2}$ and the negatively rotating field $H^{-}=\frac{H_x-{\mathit{\text{iH}}}_y}{2}$. We denote φ⁺ and φ⁻ as the phase term of H⁺ and H⁻, respectively. When single channel transmit-receive coils are used, the transmit φ⁺ can be reasonably similar to φ⁻ [16]. The measurable positively rotating field $H^{+}=\frac{H_x+{\mathit{\text{iH}}}_y}{2}$ of H=(H _x ,H _y ,H _z ) also satisfies the relation under the locally homogeneity condition

$$-{\nabla }^2{H}^{+}\left(\mathbf{r}\right)=-\mathrm{i\omega }{\mu }_0\kappa \left(\mathbf{r}\right)H^{+}\left(\mathbf{r}\right),\mathbf{r}\in \mathrm{\Omega .}$$

(5)

By separating the real and imaginary parts of κ and assuming φ⁺≈φ⁻, the phase term of H⁺ can estimate the conductivity σ _H =σ(ω) which reflects the high-frequency property responding to the RF pulse [16]:

$${\sigma }_H\left(\mathbf{r}\right)\approx \frac{{\nabla }^2{\phi }^{+}\left(\mathbf{r}\right)}{{\mu }_0\omega },\mathbf{r}\in \mathrm{\Omega .}$$

(6)

ICNE multi-echo pulse sequence and multiple phase data

The ICNE multi-echo pulse sequence is adopted to acquire multiple echoes measured at the echo time $T_{E_j},j=1,\cdots ,N_E$, where N _E is the number of echoes per repetition time T _R . Figure 1 shows its schematic diagram where I⁺ and I⁻ are sequentially injected positive and negative currents synchronized with the applied 180° rephasing pulses with N _E times. The j-th complex image, ${\zeta }_j^{\pm }$, corresponding to the j-th readout gradient generated with the positive and negative injecting currents I⁺ and I⁻, can be expressed as

$${\zeta }_j^{\pm }\left(\mathbf{r}\right)=M_{\mathit{\text{xy}},j}^{\alpha }\left(\mathbf{r}\right)e^{i{\delta }_{\epsilon }}{e}^{i{\phi }^{+}}{e}^{\pm {(-1)}^j\mathrm{i\gamma }{T}_{c^j}{B}_{z^j}\left(\mathbf{r}\right)},j=1,\cdots ,N_E$$

(7)

Figure 1

Schematic diagram for ICNE multi-echo pulse sequence to acquire dual-frequency conductivity images. ‘RF’ shows the timing for excitation of 90° and 180° RF pulses. G _z , G _y and G _x are slice selection gradient, phase encoding gradient and readout gradient, respectively. ${\zeta }_j^{\pm }$ is the j-th echo signal. ‘ICNE’ shows the alternating injection currents in the form of pulses synchronized with the RF pulse.

where $M_{\mathit{\text{xy}},j}^{\alpha }\left(\mathbf{r}\right)=M_{\mathit{\text{xy}}}^{\alpha }\left(\mathbf{r}\right)e^{-T_{E_j}/T_2\left(\mathbf{r}\right)}$ is the j-th transverse magnetization at a flip angle α, δ _ε is the systematic phase artifact, γ is the gyromagnetic ratio of hydrogen, T ₂ is the transverse relaxation time, and $T_{c^j}$ is the j-th duration of the injecting current.

The j-th transverse magnetization $M_{\mathit{\text{xy}},j}^{\alpha }\left(\mathbf{r}\right)$ can be rewritten in detail as

$$M_{\mathit{\text{xy}},j}^{\alpha }\left(\mathbf{r}\right)=C_1{M}_0\left(\mathbf{r}\right)H^{-}\left(\mathbf{r}\right)sin\left(C_2\alpha \right|H^{+}\left(\mathbf{r}\right)\left|\right)e^{-\frac{T_{E_j}}{T_2\left(\mathbf{r}\right)}},j=1,\cdots ,N_E$$

(8)

where C ₁ and C ₂ are system-dependent constants. Assuming δ _ε ≈0, the phase signal can be separated as

$${\mathcal{P}}_j^{\pm }\left(\mathbf{r}\right):=\text{arg}\left({\zeta }_j^{\pm }\left(\mathbf{r}\right)\right)=2{\phi }^{+}\left(\mathbf{r}\right)\pm {(-1)}^j\gamma T_{c^j}{B}_{z^j}\left(\mathbf{r}\right),j=1,\cdots ,N_E.$$

(9)

By subtracting and adding the phase signals ${\mathcal{P}}_j^{\pm }$, we can obtain the magnetic flux density B _z induced by the externally injecting current at the low-frequency and the phase term of H⁺ at the Larmor frequency simultaneously.

$$\begin{array}{cccc}{\mathcal{P}}_j^{+}\left(\mathbf{r}\right)+{\mathcal{P}}_j^{-}\left(\mathbf{r}\right)& =& 4{\phi }_j^{+}\left(\mathbf{r}\right),& j=1,\cdots ,N_E\\ {\mathcal{P}}_j^{+}\left(\mathbf{r}\right)-{\mathcal{P}}_j^{-}\left(\mathbf{r}\right)& =& 2\gamma T_{c^j}{B}_{z^j}\left(\mathbf{r}\right),& j=1,\cdots ,N_E\end{array}$$

(10)

From the relation 10, we can compute the j-th phase signals for MREPT and the j-th z-component of magnetic flux density for MREIT:

$$\begin{array}{cccc}{\phi }_j^{+}& =& \frac{{\mathcal{P}}_j^{+}+{\mathcal{P}}_j^{-}}{4},& j=1,\cdots ,N_E\\ B_{z^j}& =& \frac{{\mathcal{P}}_j^{+}-{\mathcal{P}}_j^{-}}{2\gamma T_{c^j}},& j=1,\cdots ,N_E\end{array}$$

(11)

The noise standard deviation of the measured magnetic flux density $B_{z^j}$ in 11 is inversely proportional to the current injection time $T_c^j$ and the signal-to-noise ratio (Υ^j) of MR magnitude image following [27, 28]

$${\mathit{\text{sd}}}_{B_{z^j}}\left(\mathbf{r}\right)=\frac{1}{2\gamma T_c^j{{\rm Y}}^j\left(\mathbf{r}\right)},j=1,\cdots ,N_E.$$

(12)

A similar analysis provides the noise standard deviation of φ j+ as

$${\mathit{\text{sd}}}_{{\phi }_j^{+}}\left(\mathbf{r}\right)=\frac{1}{4{{\rm Y}}^j\left(\mathbf{r}\right)},j=1,\cdots ,N_E.$$

(13)

Optimal combination of multiple $B_{z^j}$

In order to enhance the conductivity image at the low frequency using MREIT, we have to improve the quality of the B _z signal. An alternative to reduce the noise of B _z is to find an optimal combination of multiple $B_{z^j},j=1,\cdots ,N_E,$ by determining the weighting factors which satisfied conditions of ξ _j (r)>0 and $\sum _{j=1}^{N_E}{\xi }_j\left(\mathbf{r}\right)=1$. When we assume that the measured k-space signals are contaminated with independent identically distributed (IID) complex Gaussian random noise as same amount of expected noise standard deviation of each $B_{z^j}$ in 12, the noise variance of the weighted $B_{z^j}$ satisfies the following relation

$$\text{Var}\left(B_z^{\xi }\right(\mathbf{r}\left)\right)=\sum _{j=1}^{N_E}{\xi }_j^2\left(\mathbf{r}\right)\text{Var}\left(B_{z^j}\right(\mathbf{r}\left)\right)\propto \sum _{j=1}^{N_E}\frac{{\xi }_j^2\left(\mathbf{r}\right)}{T_{c^j}^2\left|{\zeta }_j^{\pm }\right(\mathbf{r}){|}^2}.$$

(14)

To determine the optimized weighting factor ξ _j (r), the noise variance of the combined $B_z^{\xi }$ is required to minimize the following conditions:

$$\left\{\begin{array}{l}\underset{{\xi }_j,j=1,\cdots ,N_E}{min}{\mathcal{F}}_{B_z}({\xi }_1,\cdots ,{\xi }_{N_E})\left(\mathbf{r}\right)\\ \text{subject to}{\mathcal{G}}_{B_z}({\xi }_1,\cdots ,{\xi }_{N_E})\left(\mathbf{r}\right)=1,{\xi }_j\left(\mathbf{r}\right)>0\end{array}\right.$$

(15)

where ${\mathcal{F}}_{B_z}({\xi }_1,\cdots ,{\xi }_{N_E}):=\sum _{j=1}^{N_E}\frac{{\xi }_j^2\left(\mathbf{r}\right)}{T_{c^j}^2{\left|{\zeta }_j^{\pm }\right(\mathbf{r}\left)\right|}^2}$ and ${\mathcal{G}}_{B_z}({\xi }_1,\cdots ,{\xi }_{N_E}):=\sum _{j=1}^{N_E}{\xi }_j\left(\mathbf{r}\right)$. The method of Lagrange multipliers to solve the optimization problem in 15 introduces a new variable λ and requires the condition that the gradients of ${\mathcal{F}}_{B_z}$ and ${\mathcal{G}}_{B_z}$ are parallel. Using the method of Lagrange multipliers, the optimal weighting factor ξ _j (r) can be obtained as

$${\xi }_j\left(\mathbf{r}\right)=\frac{{\Psi }_j\left(\mathbf{r}\right)}{\sum _{j=1}^{N_E}{\Psi }_j\left(\mathbf{r}\right)}$$

(16)

where ${\Psi }_j\left(\mathbf{r}\right):=T_{c^j}^2\left|{\zeta }_j^{\pm }\right(\mathbf{r}){|}^2$.

The optimized magnetic flux density to produce a conductivity image in MREIT can be generated by the weighted average of multiple $B_{z^j}$ with the determined weighting factor of ξ _j (r) as $B_z^{\xi }\left(\mathbf{r}\right)=\sum _{j=1}^{N_E}{\xi }_j\left(\mathbf{r}\right)B_{z^j}\left(\mathbf{r}\right)$ using the ICNE multi-echo pulse sequence. The noise standard deviation of $B_z^{\xi }$ can be precisely described by substituting 16 into the noise variance of $B_z^{\xi }$ in 14:

$${\mathit{\text{sd}}}_{B_z^{\xi }}\left(\mathbf{r}\right)=\frac{{\mathcal{N}}_{M_{\mathit{\text{xy}}}}\left(\mathbf{r}\right)}{2\gamma M_{\mathit{\text{xy}}}^{\alpha }\left(\mathbf{r}\right)\sqrt{\sum _{j=1}^{N_E}{T}_{c^j}^2{e}^{-\frac{2T_{E_j}}{T_2\left(\mathbf{r}\right)}}}}.$$

(17)

where ${\mathcal{N}}_{M_{\mathit{\text{xy}}}}$ denotes the noise level of magnitude image.

Optimal combination of multiple φ j+

The noise variance of the weighted φ j+ satisfies the following relation similar as 14.

$$\text{Var}\left({\phi }^{+,\chi }\right)=\sum _{j=1}^{N_E}{\chi }_j^2\left(\mathbf{r}\right)\text{Var}\left({\phi }_j^{+}\right)\propto \sum _{j=1}^{N_E}\frac{{\chi }_j^2\left(\mathbf{r}\right)}{\left|{\zeta }_j^{\pm }\right(\mathbf{r}){|}^2}$$

(18)

The noise variance of the combined φ⁺ with weighting factors is required to minimize the following conditions,

$$\left\{\begin{array}{l}\underset{{\chi }_j,j=1,\cdots ,N_E}{min}{\mathcal{F}}_{{\phi }^{+}}({\chi }_1,\cdots ,{\chi }_{N_E})\left(\mathbf{r}\right)\\ \text{subject to}{\mathcal{G}}_{{\phi }^{+}}({\chi }_1,\cdots ,{\chi }_{N_E})\left(\mathbf{r}\right)=1,{\chi }_j\left(\mathbf{r}\right)>0\end{array}\right.$$

(19)

where ${\mathcal{F}}_{{\phi }^{+}}({\chi }_1,\cdots ,{\chi }_{N_E}):=\sum _{j=1}^{N_E}\frac{{\chi }_j^2\left(\mathbf{r}\right)}{\left|{\zeta }_j^{\pm }\right(\mathbf{r}){|}^2}$ and ${\mathcal{G}}_{{\phi }^{+}}({\chi }_1,\cdots ,{\chi }_{N_E}):=\sum _{j=1}^{N_E}{\chi }_j\left(\mathbf{r}\right)$. The optimal weighting factor, χ _j (r), corresponding to each echo phase can be determined by using the method of the Lagrange multipliers as

$${\chi }_j\left(\mathbf{r}\right)=\frac{{\Phi }_j\left(\mathbf{r}\right)}{\sum _{j=1}^{N_E}{\Phi }_j\left(\mathbf{r}\right)}$$

(20)

where ${\Phi }_j\left(\mathbf{r}\right):=\left|{\zeta }_j^{\pm }\right(\mathbf{r}){|}^2$.

The optimally weighted average of multiple phase signal φ j+ with the determined weighting factor of χ _j (r), ${\phi }^{+,\chi }\left(\mathbf{r}\right)=\sum _{j=1}^{N_E}{\chi }_j\left(\mathbf{r}\right){\phi }_j^{+}\left(\mathbf{r}\right)$, can produce an enhanced conductivity image in MREPT. The noise standard deviation of the optimized phase signal, φ^+,χ, is computed by substituting (20) into the noise variance of φ^+,χ in (18) to estimate pixel-by-pixel noise precisely:

$${\mathit{\text{sd}}}_{{\phi }^{+,\chi }}\left(\mathbf{r}\right)=\frac{{\mathcal{N}}_{M_{\mathit{\text{xy}}}}\left(\mathbf{r}\right)}{4M_{\mathit{\text{xy}}}^{\alpha }\left(\mathbf{r}\right)\sqrt{\sum _{j=1}^{N_E}{e}^{-\frac{2T_{E_j}}{T_2\left(\mathbf{r}\right)}}}}$$

(21)

Considering only random noise effects, the optimal weighting factor χ _j (r) reduces the noise level, depending on the echo number and T ₂-decay rate at each imaging pixel.

Denoising method using the estimated ${\mathit{\text{sd}}}_{B_z^{\xi }}$ and ${\mathit{\text{sd}}}_{{\phi }^{+,\chi }}$

A denoising technique is applied to gradients of the optimized magnetic flux density for MREIT image and the optimized phase signal for MREPT image since the reconstruction procedures only require the differentiated measured data. Also, it is advantageous to denoise $\nabla B_z^{\xi }$ and ∇φ^+,χ because the measured magnetic flux density of $B_z^{\xi }$ and the phase of H⁺ are inherently continuous without conventional edge information.

To denoise the optimized $\nabla B_z^{\xi }$ and ∇φ^+,χ, we use the following conventional reaction-diffusion equation [29, 30]:

$$\left\{\begin{array}{lll}\frac{\mathrm{\partial v}}{\mathrm{\partial t}}(\mathbf{r},t)& =& \nabla ·\left(\psi \left(\mathbf{r}\right)\tilde{\nabla }v(\mathbf{r},t)\right)-\beta \left(\mathbf{r}\right)\left(v\right(\mathbf{r},t)-f(\mathbf{r}\left)\right)\\ v(\mathbf{r},0)& =& f\left(\mathbf{r}\right)\hfill \end{array}\right.$$

(22)

where the parameter β is a fidelity term and the function ψ satisfies

$$\psi \left(s\right)>0,\underset{s\to \infty }{lim}\psi \left(s\right)=0.$$

(23)

Here, the initial state f denotes each component of the optimized magnetic flux density $\nabla B_z^{\xi }$ or the phase signal ∇φ^+,χ, which is to be denoised. The time dependent solution v(r,t) is the denoised image of the initial state image f.

To remove the random noise artifact while preserving the edge information, it is important to determine the diffusion function ψ and β. Using the estimated noise level of ${\mathit{\text{sd}}}_{{\phi }^{+,\chi }}$, the diffusion function ψ and the fidelity function β can be specified as

$$\psi \left(\mathbf{r}\right)\propto \frac{1}{|\nabla v(\mathbf{r},t\left)\right|+{\varsigma }_{\epsilon }}\text{and}\beta \left(\mathbf{r}\right)\propto {\mathit{\text{sd}}}_{{\phi }^{+,\chi }}\left(\mathbf{r}\right)$$

(24)

where ς _ε is a small parameter to guarantee the positive sign of |∇v(r,t)|+ς _ε . For the magnetic flux density $B_z^{\xi }$, the diffusion function ψ can be selected similarly.

The quality of magnetic flux density signal in MREIT depends on the magnitude intensity and width of injected current simultaneously, while the quality of phase signal in MREPT only depends on the magnitude intensity. Therefore, there are common effects due to the decay of the magnitude intensity in MREPT and MREIT processing. There are also different effects caused by externally injecting current in MREIT. These occur when we calculate each of the optimal weighting factors and the noise standard deviation in the optimal combination of multiple phase and magnetic flux density images. The proposed optimization method uses the decay rate of magnitude intensities to determine the weighting factors in MREIT and MREPT, which relate with the T ₂ values. However, the determined weighting factors in 16 and 20 only include the measured magnitude intensity at each echo time $T_{E_j},j=1,\cdots ,N_E$. The proposed method therefore does not need estimates T ₂ values to optimize the multiple phase signals.

Phantom design and experimental setup

In order to evaluate the proposed optimizing method using ICNE multi-echo MR pulse sequence, we designed a phantom which can show different characteristics of reconstructed conductivity images by MREIT and MREPT. Figure 2(a) and (b) show the phantom configuration at the middle slice of the object and the photo view of the phantom, respectively. We built a cylindrical phantom filled with a saline of 0.2 Sm⁻¹ conductivity (0.3 g/L NaCl and 1 g/L CuSO₄). The diameter and height of the phantom were 11 cm and 14 cm, respectively. We attached four carbon hydrogel electrodes (HUREV Co. Ltd, Korea) on the side of the acrylic container for injecting current.

Figure 2

Experimental set-up for a phantom. (a) The phantom configuration at the middle slice of the testing object, (b) a cylindrical phantom with saline solution of 0.2 Sm⁻¹ including a thin film object and an agarose anomaly wrapped by an insulating film, and (c) a magnitude image at the middle slice acquired at T _E =15 ms.

Two different objects were positioned inside phantom. The left one was a thin hollow cylindrical object with the diameter of 4 cm using an insulating thin film of 0.4 mm thickness. We punched four holes which had 2 mm diameter with equally spaced around the circumference. Because we filled the same saline of 0.2 Sm⁻¹ inside and outside of the hollow cylindrical object with holes, the apparent electrical conductivities inside and outside of object were determined by the movement of ions through the holes [31]. The right cylindrical object with the diameter of 3 cm was made of an agarose gel (1 g/L CuSO₄, 2.1 g/L NaCl, 15 g/L Agar) to generate different conductivity, spin density and T ₂-decay compared with the background saline. The conductivity of agarose gel was 1.10 Sm⁻¹. It was wrapped in a thin insulating film without holes.

A transversal injecting current of 10 mA was introduced into the phantom via a pair of recessed carbon hydrogel electrodes attached at the middle of the phantom. Figure 2(c) shows the magnitude image at the middle slice. We used the ICNE multi-echo MR pulse sequence by injecting current of alternating polarity synchronized with the multiple refocusing pulses as shown in Figure 1.

Imaging parameters in a 3T MRI scanner (Achieva, Philips) with birdcage transmit-receive (Tx/Rx) RF head coil were as follows: repetition time T _R =1200 ms, data acquisition time width T _s =3.584 ms, echo-spacing $T_{E_{\mathit{\text{sp}}}}=15$ ms, total number of echo N _E =6 and number of averaging = 4. The reconstructed image matrix was 128×128, with a FOV of 180×180 mm², 5 mm slice thickness. Total scanning time for both vertical and horizontal injection currents was 20 minutes with an interleaved phase encoding acquisition.

Results

Figure 3 shows the multiple MR magnitude images measured at the echo time $T_{E_j}=15\times j$ ms for j=1,⋯,6, where the right anomaly shows rapid decay of the transverse magnetization due to magnetic field non-uniformity and spin-spin transverse relaxation.

Figure 3

MR magnitude images measured at each echo time. Echo number marked on the upper left part of each image as (1), (2), ⋯, (6). Magnitude images using ICNE multi-echo pulse sequence acquired at the echo time $T_{E_j}=15\times j$ ms for j=1,⋯,6.

Figure 4 shows the phase images of φ j+ of H⁺ map and magnetic flux density $B_{z^j}$ induced by the vertically injected current using ICNE multi-echo pulse sequence at the echo time $T_{E_j}=15\times j$ ms, j=1,⋯,6. The images in Figure 4(a) and (b) were estimated by adding and subtracting the phase signals of ${\zeta }_j^{\pm }$, respectively. Since the electromagnetic wave at the Larmor frequency of 128 MHz could penetrate the thin insulating film and there was same saline inside and outside of the thin insulating wall of the left hollow cylinder, the measured phase values of φ j+ for the left object in MREPT could not produce conductivity difference between inside and outside. However, the measured $B_{z^j}$ images at the low frequency in MREIT showed the different value inside the left anomaly because MREIT measures the apparent magnetic flux density corresponding to the ion mobility and intrinsic conductivities of composite materials. Note that the slope of magnetic flux density was changed abruptly due to the high current density around the holes in the same direction of the current injection. Regarding the different conductivity agarose gel cylinder wrapped in the thin insulating film without holes, both images in Figure 4(a) and (b) could present distinctions between inside of the object and the background saline.

Figure 4

Phase images and magnetic flux density images of each echo. (a) Phase images φ j+ of H⁺ map in radian and (b) magnetic flux density images $B_{z^j}$ in ’nT’ unit using ICNE multi-echo pulse sequence at the echo time $T_{E_j}=15\times j$ ms for j=1,⋯,6.

Using the proposed optimization method with the estimated noise standard deviations ${\mathit{\text{sd}}}_{{\phi }_j^{+}}$ and ${\mathit{\text{sd}}}_{B_z^j},j=1,\cdots ,6$, we found the optimal weighting factors for the phase of H⁺ and the magnetic flux density by the injected current, respectively. Figure 5(a) and (b) show the determined weighting factors for dual-frequency conductivity images in MREPT and MREIT combined method, respectively. In MREPT, the intensity of weighting factors was monotonically decreased because the noise level of φ j+ only depended on T ₂-decay rate of ${\zeta }_j^{\pm },j=1,\cdots ,6$, whereas the weighting factors in MREIT show the different characteristics because they depended on the combination of T ₂-decay rate and the duration of injected current simultaneously. In Figure 5(b), the intensity of weighting factors in the background area was monotonically increased because the effect of time increasing for the injection current was greater than the magnitude attenuation effect of echo signals. However, the weighting factors in the right agarose gel anomaly region were different due to the combination of both effects. The weighting factor values at the second echo were higher than those of other echoes.

Figure 5

Images of weighting factors corresponding to each echo. (a) Images of weighting factors corresponding to each echoes for φ j+ and (b) images for magnetic flux density $B_{z^j}$ to produce an optimal φ^+,χ and $B_z^{\xi }$, respectively. The weighting factors were scaled from 0 to 1.

Figure 6 shows the reconstructed dual-frequency conductivity images using the multiple echoes φ j+ and $B_{z^j}^d$ for j=1,⋯,6 and d=1,2, where $B_{z^j}^1$ and $B_{z^j}^2$ were the measured magnetic flux density data by the vertically and horizontally injected currents, respectively. To reconstruct the high-frequency conductivity distribution in Figure 6(a), we used phase for each echo signal as

$${\sigma }_H^j\left(\mathbf{r}\right)=\frac{{\nabla }^2{\phi }_j^{+}\left(\mathbf{r}\right)}{{\mu }_0\omega },\mathbf{r}\in \Omega ,j=1,\cdots ,6.$$

(25)

Figure 6

Reconstructed conductivity images of each echo. (a) Reconstructed conductivity images using the multiple echoes φ j+ for j=1,⋯,6 at Larmor frequency and (b) conductivity images using $B_{z^j}^d$ for j=1,⋯,6 at the low frequency.

Here, we could not detect the hollow cylinder in the left side made of the thin insulating film with holes. The high-frequency conductivity inside and outside of the left object was estimated to be the same value. The quality of reconstructed conductivity image for each echo in Figure 6(a) was monotonically degraded corresponding to the weighting factors shown in Figure 5(a). Table 1 shows the values of weighting factor inside and outside of the right agarose object wrapped by insulating film.

Table 1 Weighting factor values for MREPT inside ( ${\mathcal{R}}_{\mathit{\text{in}}}$ ) and outside ( ${\mathcal{R}}_{\mathit{\text{out}}}$ ) the agarose anomaly

We estimated the projected current density J^P, which was the optimal divergence-free component from the measured z-component magnetic flux density in order to reconstruct the low-frequency conductivity images [32]. The projected current density consists of the background current density J⁰=−∇α and${\nabla }_{\mathit{\text{xy}}}^{\perp }\eta =\left(\frac{\mathrm{\partial \eta }}{\mathrm{\partial y}},-\frac{\mathrm{\partial \eta }}{\mathrm{\partial x}},0\right)$, where the potential η satisfies at the imaging slice Ω _t ⊂Ω

$$\begin{array}{cccc}{\nabla }_{\mathit{\text{xy}}}^2\eta & =& \frac{1}{{\mu }_0}{\nabla }^2{B}_z& \text{in}{\Omega }_t\\ \eta & =& 0& \text{on}\partial {\Omega }_t.\end{array}$$

(26)

where ∇ _xy denotes the two-dimensional gradient. The absolute conductivity distribution can be written as a matrix form using the two measured$B_z^d,d=1,2,$data

$$\left({\mathbf{A}}_{\mathbf{0}}+{\mathbf{A}}^{\ast }\right){\left(\frac{\mathrm{\partial \tau }}{\mathrm{\partial x}},\frac{\mathrm{\partial \tau }}{\mathrm{\partial y}}\right)}^T=\mathbf{b}\text{in}{\Omega }_t$$

(27)

where Ω _t is an imaging slice, τ= logσ, and

$${\mathbf{A}}_{\mathbf{0}}=\left(\begin{array}{cc}\frac{-\partial {\alpha }^1}{\mathrm{\partial y}}& \frac{\partial {\alpha }^1}{\mathrm{\partial x}}\\ \frac{-\partial {\alpha }^2}{\mathrm{\partial y}}& \frac{\partial {\alpha }^2}{\mathrm{\partial x}}\end{array}\right),{\mathbf{A}}^{\ast }=\left(\begin{array}{cc}\frac{\partial {\eta }^1}{\mathrm{\partial x}}& -\frac{\partial {\eta }^1}{\mathrm{\partial y}}\\ \frac{\partial {\eta }^2}{\mathrm{\partial x}}& -\frac{\partial {\eta }^2}{\mathrm{\partial y}}\end{array}\right),\mathbf{b}=\frac{1}{{\mu }_0}\left(\begin{array}{c}{\nabla }^2{B}_z^1\\ {\nabla }^2{B}_z^2\end{array}\right).$$

(28)

Figure 6(b) shows the reconstructed low-frequency conductivity images by solving 27 using the j-th measured two magnetic flux densities$B_{z^j}^d,d=1,2$obtained with the vertically and horizontally injected currents. The reconstructed conductivity values inside the hollow cylinder with holes in the left side reflected the externally injected current flowed through the four holes. Regarding the agarose gel object wrapped in the insulating film without holes, the apparent conductivity reconstructed by solving 27 shows the same anomaly as an insulator since the measured B _z data only reflected the current density distribution due to the externally injected current through the attached electrodes. The quality of conductivity image using the sixth echo data was poor inside the agarose gel anomaly region due to relatively short T ₂-decay relaxation time. It was consistent with the corresponding weighting factor in Figure 5(b). Table 2 shows the values of weighting factor inside and outside of the right agarose object wrapped by insulating film.

Table 2 Weighting factor values for MREIT inside ( ${\mathcal{R}}_{\mathit{in}}$ ) and outside ( ${\mathcal{R}}_{\mathit{out}}$ ) the agarose anomaly

Figure 7(a) and (e) show the normalized noise standard deviations of optimally weighted φ^+,χand$B_z^{\xi }$, respectively. In previous approaches without using multi-echo data with optimizing method, the conductivity images for MREPT and MREIT were reconstructed as shown in Figure 7(b) and (f). Due to the different T ₂-decay rate and the duration of injected current, the intensities of${\mathit{\text{sd}}}_{{\phi }^{+,\chi }}$and${\mathit{\text{sd}}}_{B_z^{\xi }}$represented different characteristics. The reconstructed dual-frequency conductivity images with the optimally weighted φ^+,χand$B_z^{d,\xi },d=1,2,$were displayed in Figure 7(c) and (g). Solving the reaction-diffusion equation 22 with the estimated noise standard deviations${\mathit{\text{sd}}}_{{\phi }^{+,\chi }}$and${\mathit{\text{sd}}}_{B_z^{\xi }}$, we denoised ∇φ^+,χand$\nabla B_z^{d,\xi }$; the iteration number was 400 and the regularization function of β was$\frac{0.05}{\mathcal{N}\left({\mathit{\text{sd}}}_{{\phi }^{+,\chi }}\right)}$for MREPT case, where$\mathcal{N}\left({\mathit{\text{sd}}}_{{\phi }^{+,\chi }}\right)$was the normalized noise standard deviation in Figure 7(a). The denoising parameters for MREIT were the same as those for MREPT. Figure 7(d) and (h) show the reconstructed dual-frequency conductivity images after denoising of ∇φ^+,χand$\nabla B_z^{d,\xi }$, respectively.

Figure 7

Optimized dual-frequency conductivity images with and without applying the optimization and denoising method. (a) and (e) Normalized noise standard deviation of the optimally weighted phase of H⁺and$B_z^{\xi }$, respectively. (b) and (f) Reconstructed conductivity images without using the optimization and denoising method. (c) and (d) Reconstructed conductivity distributions at Larmor frequency using the optimized φ^+,χwithout denoising and with the proposed denoising technique applying to ∇φ^+,χ, respectively. (g) and (h) Reconstructed low-frequency conductivity distributions without denoising and with the proposed denoising technique applying to$\nabla B_z^{d,\xi },d=1,2$, respectively.

Discussion

MREPT and MREIT used an internal magnetic flux density information from a phase imaging can provide conductivity images with high resolution compared to the conventional electrical impedance imaging techniques. Moreover, the interleaved measurement technique with alternating injected current polarity based on a multi-spin-echo pulse sequence, enables simultaneous dual conductivity images. From a single MR scan, we could obtain the low-frequency conductivity images by subtracting and the high-frequency conductivity images by adding the measured phase terms. Despite their advantages and usefulness, MREPT and MREIT techniques commonly suffer from noise amplifications in phase imaging. The proposed method used a ICNE multi-echo pulse sequence based on a spin echo pulse sequence in order to suppress the background field inhomogeneity. This also influenced to the reconstructed high-frequency conductivity using the MREPT technique. Even though the estimated weighting factors for MREIT and MREPT optimally reduce the random noise artifacts, non-uniform noise artifacts from various sources can severely deteriorate the reconstructed conductivity distributions. To be a clinical device used practically, it is important to develop noise reduction techniques taking into account the physical properties and experimental environments of MREIT and MREPT.

The designed phantom consists of two different kinds of anomalies to show the frequency dependent conductivity information. The stable apparent conductivity contrast of hollow cylindrical object made of thin insulating transparency film was only controlled by ion mobility through holes on the film wall. It excluded other effects caused by any ion concentration gradient. Therefore, MREIT images could only present the hollow cylindrical object with holes when comparing the phase terms${\mathcal{P}}^{+}$and${\mathcal{P}}^{-}$in 10. The apparent conductivity of the region inside cylinder was a non-zero conductivity value because the externally injected current entered the anomaly through the holes. In MREPT, the reconstructed high-frequency conductivity was same inside and outside of the anomaly since the measured phase of H⁺only reflected the material conductivities inside and outside of the cylindrical film. Another agarose gel cylinder wrapped in an insulating film without holes shows different conductivity characteristics in MREPT and MREIT. Due to the short T ₂-decay relaxation time of agarose gel, the noise levels of${\mathit{\text{sd}}}_{{\phi }^{+,\chi }}$and${\mathit{\text{sd}}}_{B_z^{\xi }}$were relatively severe compared to those in the other regions. The noise level of each${\mathit{\text{sd}}}_{{\phi }_j^{+}}$was strictly increasing with respect to the echo number depending on the T ₂-decay rate, whereas those of${\mathit{\text{sd}}}_{B_{z^j}}$had different characteristics depending on$T_{c^j}$and the T ₂-decay rate simultaneously.

MREIT provides the conductivity distribution at a low frequency, whereas MREPT produces the conductivity at the Larmor frequency of 128 MHz at 3 T. The cell membranes consisted of the phospholipid bilayer with embedded proteins behave as a capacitor or insulating film. It shows a complicated pattern of conductivity depending on the measured frequencies. For the biological tissues, the membranes restrict the flow of current of low frequency. MREIT images reflect the membrane properties quantitatively and MREIT has a potential to visualize the anisotropic conductivity tensor map. On the other hand, the membranes become transparent at high frequency and provide a relatively degraded sensitivity in the conductivity image. MREPT may provide different conductivity characteristics compared to MREIT results. Based on the previous studies related to the complex conductivity spectra measured at multiple frequencies within the range of 10 Hz to several kHz, spectroscopic complex conductivity distribution can contribute to explain the physiological and pathological status of internal tissues [33]. The dual-frequency conductivity imaging by combination of MREPT and MREIT is an advanced technique to show different information in spite of displaying two extreme cases at low and the Larmor frequencies. It is meaningful to develop a method to produce multi-frequency conductivity image spectra with high spatial resolution and sensitivity using MR scanner within the range of DC to the Larmor frequency since significantly distinguishable signal changes in biological tissues depending on the applied frequency range.

We have a plan to support the clinically applicable combined conductivity imaging method based on MR techniques. We may apply this method to detect the destroyed cell membrane or the modification of cell or tissue structure due to necrosis or apoptosis initiated by inflammatory response inside the body. Since RF ablation and cryoablation for cancer treatments are accompanied by the destruction of cell membranes, it will be a good method to diagnose ablated lesions based on physiological status of tissue to improve the safety and predict the local recurrence after ablation.

Conclusions

Cross-sectional conductivity imaging methods for high spatial resolution and sensitivity inside the human body has been actively investigated. MREIT and MREPT techniques show different electrical properties at low frequency (below 1 kHz) as MREIT uses an externally injecting current while MREPT uses the Larmor frequency of an MRI scanner. Recently, the dual-frequency conductivity imaging from a single MR scan simultaneously was proposed based on combination of MREPT and MREIT. Even though it produces conductivity spectral information, MREPT and MREIT have commonly suffered from weak signals and noise amplification since the procedures to reconstruct the conductivity in MREIT and MREPT need to differentiate the measured phase signals. We suggested the optimization method to find weighting factors according to echo signals for MREPT and MREIT using the ICNE multi-echo pulse sequence, which minimized the noise artifacts in the measured phase data. The noise standard deviations were precisely analyzed for the optimally weighted magnetic flux density and the phase term of positive-rotating magnetic field. We applied the denoising method based on the reaction-diffusion equation with the estimated noise standard deviations to enhance the quality of dual-frequency conductivity images. A real phantom experiment was performed to validate the proposed method using common measured data to reconstruct the dual-frequency conductivity distributions for MREPT and MREIT.