## Introduction

Magnetic resonance imaging (MRI)-based magnetic resonance electrical property tomography (MREPT) (1, 2) is a strategy for noninvasively reconstructing electrical properties (EPs) (permittivity and electrical conductivity) of the human body without using additional hardware to a magnetic resonance (MR) scanner. It has been intensively developed in recent years for reconstructing the EPs of human tissues based on the measurable data (radiofrequency [RF] field distribution called B_{1}-map) (3) from an MRI scanner because the B_{1}-field carries the information of EP distribution in the human body.

MREPT is an important technology in various areas in medicine and biology. The constructed EP map can show the anatomical structure of the human body by the high contrast of EPs between different tissues (4). Cancerous tissues, including those at the early stages, are highly distinguishable from healthy ones (5, 6). Therefore, the EP map provides a good tool for early cancer detection. It can also provide guidance to design body-centric communications (7) and electromagnetics-based therapies (8). EP maps can also help in understanding the activities of cells, in particular with exposure to either static or dynamic electromagnetic waves (9). Moreover, an EP map is crucial for accurate calculation of specific absorption rate, which is a main parameter for assessing the risk of RF and microwave radiation and that of ultrahigh-field MRI scanning (10, 11).

MREPT was originally proposed by Haacke et al. (12) and first implemented by Wen (13). The Helmholtz equation for homogeneous material was proposed for calculating both conductivity and permittivity with positively circularly polarized component of the magnetic field (B1+), which corresponds to the RF transmit field. Another MREPT method was developed and systematically studied by Katscher et al. (14) and Voigt et al. (15). The integral forms of Faraday's law and Ampère's law were applied to formulate the governing equation based on B1+. However, these methods are based on an assumption of local homogeneity for simplification. This assumption leads to artifacts in the region where EPs vary either quickly or abruptly (3, 16, 17). The reconstruction errors are rigorously analyzed by Seo et al. (18). A recent work was focused on removing the local homogeneity assumption (19).

One full generalization of these methods was performed by Sodickson et al., and the technique was named Local Maxwell Tomography (LMT) (20, 21). The governing equation of this method was derived from Faraday's law and Ampère's law. The unmeasurable magnetic field components are eliminated by using Gauss's law. The method is physically feasible and free of assumptions. Magnetic fields measured using an MRI scanner with multiple channels are used to solve the unknown. Although the assumption of local homogeneity is removed in LMT, higher-order partial differentiation over the magnetic field is needed in the calculation, which significantly increases this method's sensitivity to noise.

Another MREPT method (22) was proposed on the basis of the fact that the magnetic field B1+ near the center of a birdcage or a transverse electromagnetic coil is homogeneous along its central axis (23, 24), and it is assumed that the contribution of the along-axis component of the magnetic field is negligible when the region of interest is close to the center of the coil. A first-order PDE including a gradient term of the unknown EP was formulated via properly manipulating Maxwell equations, and the gradient term describes the inhomogeneity of the object under investigation. This method is more general than the MREPT proposed by Katscher et al. (14) and has a simpler form than LMT.

Hafalir et al. (25) then named this method as convection–reaction MREPT (cr-MREPT), as its governing equation is similar to the convection–reaction equation (26). The first-order PDE was then solved by Hafalir et al. (25) for 2-dimensional reconstruction using a strong form with finite element method (FEM) (27). An obvious artifact was observed in the region where the gradient of |B1+| is close to zero. Although 2 methods, the double-excitation and the constraint MREPT, were further proposed by Hafalir et al. (25) to improve the results, artifacts reappear when noise exists. Meanwhile, cr-MREPT does not work with every phantom. Even with the phantoms used by the authors, the reconstructed EPs deviate seriously from the true values in a large region close to the central artifact. High sensitivity of noise was also observed.

In parallel with the work of Hafalir et al., Ammari et al. (19) also attempted to solve this PDE via an adjoint-based iterative reconstruction algorithm. The real and imaginary parts of the equation are separated to produce a coupled equation pair solved by the iterative scheme. This method can construct a blur image of the EPs without the aforementioned artifacts, but it is not applicable in the region where the gradient of |B1+| field is close to zero, and the produced linear system is ill-conditioned. The method also converges slower than the FEM-based method (19).

The solution of the first-order PDE shows persistent artifacts where ∇|B1+| approaches zero. Other than that, global spurious oscillations exist in the solution of this first-order PDE. Therefore, directly solving it with numerical methods is difficult. Here, the existing challenges for solving the PDE are discussed in detail. The existing artifacts in the numerical results are pointed out and classified. An effective viscosity-type regularization (28) solution is applied to mitigate the aforementioned problems. A Laplacian term is introduced into the first-order PDE for stabilizing the entire system, which can be solved by either the finite difference (FD) method or FEM. The governing PDE is in the first order, so its solution will be a linear function that depends on the space coordinates, and it can present sharp variations in space. Because the PDE is not a singularly disturbed equation (29, 26, 30), the introduced second-order Laplacian term will modify the PDE to obtain a second-order solution, which is an approximation of its exact solution and presents a smoother property in space.

The derivation of the formulation is briefly presented for the convenience of discussion; behaviors of the equation are analyzed in detail before introducing the viscosity-type regularization. Numerical experiments are used to facilitate the analysis of the existing problems and to show the effectiveness of the proposed approach. Finite difference method is applied for discretization. The proposed method is applied for reconstructing the EPs of different phantoms. Accuracy, stability, and noise tolerance of the method are illustrated with numerical results.

### Theory and Methodology

The governing PDE for cr-MREPT are given as follows:

This equation can be solved for reconstructing EPs of a 3-dimensional object (see online supplemental Appendix for the relating derivation of this equation). If γ does not vary severely along the *z* direction, ∂γ∂z becomes negligible. Equation (1) is then simplified as follows for a single section located at the center of birdcage coil:

Lxy and Mxy are vectors L and M without the *z*-component, respectively.

∇xy=∂∂xx∧+∂∂yy∧. Notice that ∇2H1+ still needs to be evaluated in a 3-dimensional form. Magnetic field *H*_{z} is not measurable using MRI, but its contribution is negligible compared with H1+, as *H*_{z} is generally close to zero in the center of the birdcage coil. Therefore equation (2) could be further simplified as follows:

For homogeneous regions, ∇xyγ=0; therefore, equation (3) can be simplified as γ=iωμH1+/∇2H1+, which is the governing equation proposed by Haacke et al. and Wen (12, 13), named, in this content, standard MREPT (stdMREPT).

Therefore, the original inverse problem equivalently becomes a problem of solving the first-order PDE (3). It has the form of a stationary convection–reaction equation, which is a mathematical model for fluid dynamics (31). Thus, the MREPT based on equation (3) is called cr-MREPT by Hafalir et al. (25).

The coefficients of PDE (3) are complex, making the PDE exhibit different mathematical and numerical behaviors from the stationary convection–reaction equations that have real coefficients. In particular, both equations show global spurious oscillations in their solutions (26). For convection–reaction equations (real coefficients), these global spurious oscillations are caused by the mesh density for discretization, and the results converge better when the mesh size decreases. However, this fact does not hold for equation (3), where the coefficients are complex. Instead, the spurious oscillation in equation (3) is caused by a fundamental mathematical drawback of itself; this is further discussed as follows.

Equation (3) can be reformed as follows:

with

F(r)=Lxy and

f(r,γ(r))=γ∇2H1+−iωμH1+. According to the Cauchy–Kowalevski theorem (

32), a necessary condition for equation (

4) to have a smooth solution is that

f(r,γ(r)) is analytic in region Ω (

19). However,

∇2H1+ is not smooth in the region where electrical parameters vary rapidly (as in

Figure 3H for a 2-layer cylindrical phantom with difference

ϵr and σ in different layers). Discontinuity of

∇2H1+ exists at the boundary where electrical parameters vary rapidly. Therefore, a smooth solution of σ and

ϵr cannot be obtained with equation (

3), and a better strategy is needed for good numerical stability, accuracy, and efficiency.

##### Figure 3.

σ (A) and ϵr (B) of phantom shown in Figure 2 are reconstructed by applying convection–reaction magnetic resonance electrical property tomography (MREPT) (cr-MREPT) based on equation (3). σ (C) and ϵr (D) reconstructed using equation (2). Real (E) and imaginary (F) parts of H1+, absolute value of Lx (G) and ∇2H1+ (H), are also presented.

Here, a viscosity-type regularization method (28) is applied to equation (3) to eliminate the global spurious oscillation in its solution. The governing equation becomes the following equation:

where

ρ∇x,y2γ is the additional Laplacian term, and ρ is a constant that is always <1 but >0. This approach does not lead to a singularly disturbed problem (

29,

26,

30). The additional Laplacian term plays a regularization role. This stabilized MREPT method is named as stabMREPT. There is no good way for deciding the optimum value of ρ yet, but it can be chosen as several times or a fraction of the maximum of

|Lx|, which is the

*x*-component of

Lxy. There are 2 ways to add the Laplacian term. One is adding it to the whole calculation domain, as shown in equation (

5), and the other is adding it to the region where it is necessary (

33) to produce more accurate results, but the optimal coefficients of the added Laplacian term need to be determined in a complex way, rendering it complicated. Here, only equation (

5) is implemented for showing the idea. Central difference formula of the FD method is used for discretization.

### Discretization Strategy

A square domain denoted as Ω and bounded by ∂Ω is discretized by small squares with Nx+2 and Ny+2 grid points along the *x* and *y* directions, respectively, giving (Nx+2)×(Ny+2) data points in space. Δx and Δy are the step lengths along *x* and *y* directions, respectively. Use *i* and *j* to index the mesh points along *x* and *y* with i:={i|i∈Z, 0≤i≤Nx+1} and *j* := {j|j∈Z, 0≤j≤Ny+1}, respectively. The values of H1+ and γ at the point (*i*, *j*) are represented as H{i,j} and γ{i,j}, respectively.

Allowing Lxy=[Lx, Ly] with Lx=−∂H1+/∂x+i∂H1+/∂y and Ly=−i∂H1+/∂x−∂H1+/∂y. The values of Lx and Ly at point (*i*, *j*) are represented with Lx{i,j} and Ly{i,j}, respectively. Using first-order central difference formula, Lx{i,j} and Ly{i,j} are given as follows:

Here, the imaginary unit *i* can be easily distinguished from the index *i* on the basis of its location. Representing the data of ∇2H1+ at the point (*i*, *j*) as d2H{i,j}, the second-order central difference formula is used to discretize ∇2H1+, which gives the following equation:

*k* is the index for sections along the *z*-direction. Two extra sections at −Δz and Δz away from the central section are used to calculate the partial differentiation of H1+ over *z*.

Another way for calculating Lx{i,j},Ly{i,j}, and d2H{i,j} is by using the Savitzky–Golay (SG) filter (34), which has been widely used to smooth noisy data. Here, H1+ is approximated with the following equation:

where (

*x*,

*y*) is the corresponding coordinate of

H1+.

*a*_{p},

*p* = 1, . ., 6 are polynomial coefficients, which need to be decided with the measured data of

H1+. As an example, at mesh grid point (

*i*,

*j*),

H{i,j} has 8 neighbors, as shown in

Figure 1. These neighbors or

H{i,j} itself could be used for interpolation with Stanley et al. (

9), which will produce an overfitting linear system about the unknown coefficients

*a*_{p},

*p* = 1, …, 6. A least-square method could then be used for solving the unknown.

##### Figure 1.

Neighbors of H{i,j} in a mesh grid.

##### Figure 4.

σ (A) and ϵr (B) reconstructed via least-square minimization with λ=0.01max{A} for the phantom shown in Figure 2. Central difference is used for ∇2H1+.

The derivatives over H1+ at the point (*i*, *j*) are then calculated using the following equations:

Here, (*x*_{i}, *y*_{j}) is the coordinate of the grid point (*i*, *j*). Notice that the Laplacian term over H1+ includes the second-order partial differentiation of H1+ over *z*, which is calculated with a central difference formula, with the data at sections located at −Δz and Δz. With equations (6)–(8), or (10)–(12), equation (5) can be easily discretized with the first- and second-order central difference formula, as follows:

A linear system is produced via assembling (13) for all points indexed by (*i*, *j*). The linear system is denoted as Ag=b. The sparsity of the tridiagonal matrix A allows efficient calculation with Gaussian elimination. The Dirichlet boundary condition is assumed here. Vector **g** represents the unknowns γ{i,j}. The relative permittivity and conductivity at (*i*, *j*) can be calculated based on the real and imaginary parts of γ{i,j}.

### Numerical Experiments

In the calculation, H1+ as input is simulated with FEM-based software COMSOL (35). A 16-leg high-pass shielded birdcage coil with a radius of 14.5 cm and a leg-length of 24 cm (36) is built to work at 127.74 MHz (3 T MRI system) to produce a homogeneous H1+ within the coil, and it is placed with its axis along the *z*-axis. The coil is driven by 2 ports with 500 V for each. Both ports are geometrically 90° apart with a 90° phase difference.

The phantom of 2 coaxial cylinders used for simulation is shown in Figure 2A. The cylinders have the same height h=24cm. The radii of the exterior and interior ones are r1=5cm and r2=2.5cm, respectively. The exterior cylinder is filled by a material with ϵr1=75 and σ1=0.5 S/m, whereas for the interior one, ϵr2=50 and σ2=1 S/m.

##### Figure 2.

The phantom of 2 cylinders with the 3D view (A) and its cross section view (B). *h* = 24 cm, *r*_{1} = 5 cm, and *r*_{2} = 2.5 cm. Dashed-line square defines region of interest. Unit of sigma is S/m.

##### Figure 5.

σ (A) and ϵr (B) of the phantom in Figure 2 are reconstructed by cr-MREPT with *N*_{x} = *N*_{y} = 199.

##### Figure 6.

Reconstructed results of σ (A) and ϵr (B) for the 2-cylinder phantom shown in Figure 2 with stabilized MREPT (stabMREPT). δ=0.37, *N*_{x} = *N*_{y} = 199.

The phantom is placed in the birdcage coil, with its axis and center coinciding with those of the coil. A Dell Precision workstation (Round Rock, Texas) with Intel Xeon CPU with 2 processors (4 cores for each) and 48 GB RAM is used for simulation. The phantom is meshed with tetrahedral elements. The maximum mesh size for −0.5cm<z<0.5cm is 1 mm. The mesh size is then increased gradually inside out. This strategy guarantees the accuracy of the calculated H1+ at the sections of interest and minimizes the number of elements to reduce the required memory. H1+ of the sections at *z* = 0, *z* = −0.25 cm, and z=0.25cm is extracted as inputs for the reconstruction.

The domain Ω is defined as the square inscribed to the cross section of the exterior cylinder, indicated by the dashed lines in Figure 2B. It is then meshed with Nx=Ny=99. The data of H1+ are directly extracted from COMOSL.

The mesh used here will lead to a space resolution of 0.75 mm. Because simulation data are used here, not much concerns are needed. But, in practice, this resolution will need a long MR scan and produce low signal-to-noise ratio (SNR). A coarse mesh can be used first to perform the MR scan. With the obtained H1+, one can easily apply the cyclic regularized Savitzky–Golay filter (37) to reconstruct H1+ of a finer mesh with the required resolution. This will also help in reducing noise.

Relative permittivity ϵr and conductivity δ of phantom shown in Figure 2 are first reconstructed by equation (3), the cr-MREPT method. Differentiations of H1+ are calculated with the central difference formula. As shown in Figure 3, A and B, significant artifacts are observed around the center of the interior cylinder and in the regions where σ and ϵr are discontinuous. Furthermore, global spurious oscillations are also clearly seen. For checking the effect of neglecting the contribution of *H*_{z}, equation (2) is also applied for the reconstruction. Results are shown in Figure 3, C and D, and not much improvement is seen. The global spurious oscillations and the artifacts still exist.

Hafalir et al. (25) report that the central artifacts are caused by the fact that Lx and Ly are close to 0. Figure 3, E and F shows the real and imaginary parts of the simulated H1+ at the central section (*z* = 0), denoted as RH1+ and IH1+, achieve to their extrema in the vicinity of the center, as shown in the figures, indicating that ∂H1+/∂x and ∂H1+/∂y are close to zero. This results in Lx{i,j}≈0 and Ly{i,j}≈0 for the points in the center region. The absolute values of Lx are shown in Figure 3G (|Lx|=|Ly|).

A minimization method is also applied to mitigate the global spurious oscillation when solving equation (3). The original linear system produced by equation (3) is written as Ag=b, which is then transformed into a least-square minimization problem min||b−Ag||2+λ||g||2 with a penalty term λ||g||2. The damping parameter λ>0. In this way, the solution for **g** is regularized. It is then solved with the iterative algorithm LSQR, which is based on the Golub–Kahan bidiagonalization process (38). The condition number of A here is about 10^{3}, so (ATA+λI)g=ATb,A being not symmetric, is not suitable for the calculation (it is equivalent to the applied minimization problem), because the 2-norm condition number of ATA is exactly the square of the one of A.

This minimization is conducted with λ=0.01max||A||2, which is a large value for λ. As shown in Figure 4, the results converge better to the true values, but an obvious “4-blade” artifact is observed. These global spurious oscillations remain.

Moreover, this global spurious oscillation is independent of grid sizes Δx and Δy. Figure 5, A and B show the results calculated with Nx=Ny=199, which is a much finer mesh in the region Ω. Global spurious oscillation still exists. Many artifacts appear with very large values. The reconstructed values of |ϵr| and |σ| achieve to 500 and 10.

Therefore, neither directly the solution of equations (2) and (3) nor that of the least-square minimization method produces accurate results. This is expected by the fundamental drawback existing in the governing PDE, which is explained in section Theory and Methodology. However, reconstruction with the proposed stabMREPT based on the new second-order PDE (5) shows much better performance, which is shown as follows.

The same phantom (shown in Figure 2) is used for the proposed stabMREPT. The region Ω is discretized with the same mesh density (Nx=Ny=199). The parameter ρ=0.37, which is 7 times of the maximum value of |Lx|. Figure 6 shows the reconstructed conductivities and permittivities. The normalized root-mean-square error (NRMSE) is defined as %error=100||(etrue−erecon)||2/||etrue||2. The reconstructed results of σ and ϵr are very close to the true values. Comparison of the reconstructed σ and ϵr with their true values at the line *y* = 0 yields the plots shown in Figure 7. A very good match is observed. Relative errors are all <1% except the boundary where the electrical parameters are discontinuous. Moreover, including *H*_{z} into equation (5) produces the same results.

##### Figure 7.

Comparison between the true values and the reconstructed ones of σ (A) and ϵr (B) along *y* = 0. stabMREPT is applied for the 2-cylinder phantom. The regularization parameter δ=0.37, *N*_{x} = *N*_{y} = 199.

##### Figure 8.

The map of normalized root-mean-square error (NRMSE) for σ reconstructed with stdMREPT (A) and stabMREPT (B). The 2-cylinder phantom is used. δ=0.37, *N*_{x} = *N*_{y} = 199.

##### Figure 9.

Simplified head model with 5 different regions characterized with electrical parameters of scalp (ϵr=62,σ=0.54 S/m, gray), skull (ϵr=21,σ=0.12 S/m, deep blue), gray matter (ϵr=73,σ=0.59 S/m, yellow), white matter (ϵr=52,σ=0.34 S/m, light blue), and cerebrospinal fluid (ϵr=84,σ=2.14 S/m, red).

To further check the performance of stabMREPT compared with stdMREPT at both the boundary and at the homogeneous regions (no abrupt change in EPs), the NRMSE is calculated here for the reconstructed σ of the 2-cylinder phantom. The parameters for the reconstruction are the same as those used for Figure 6. The map of NRMSE is shown in Figure 8. A reconstruction error of stdMREPT at the boundary of interior cylinder is >100%, whereas a large error is also observed in the region far from the boundary. A reconstruction error at the boundary of the cylinder is also observed in the results obtained with stabMREPT, which is about 45%, and this error decreases gradually along the radius. But the reconstruction error at the region far away from the boundary is much smaller than the one at stdMREPT, where the largest error is <3%. Therefore, stabMREPT performs better than stdMREPT in both the boundary and homogeneous regions.

To further test the proposed approach, a simplified human head phantom as shown in Figure 9 is simulated. The model has 5 different regions indicated by different colors. Each region relates to a tissue type of the human brain (4). These tissues are, from external to internal, scalp (ϵr=62,σ=0.54S/m), skull (ϵr=21,σ=0.12S/m), gray matter (ϵr=73,σ=0.59S/m), white matter (ϵr=52,σ=0.34S/m), and cerebrospinal fluid (ϵr=84,σ=2.14S/m) at 127.74 MHz (39). For the convenience of calculation, the region indicated with gray color is set to a 12- × 14-cm (length × width) square). The 2-dimensional section shown in Figure 9 is drawn and then extruded in COMSOL to produce a 1-cm-thick model. The thickness is primarily limited by the workstation's RAM. The curved edges and the central red region are meshed with an element size smaller than 1 mm, and other regions are meshed with an element size smaller than 1.5 mm.

The region Ω defined by −6 cm < *x* < 6 cm and −7 cm < *y* < 7 cm is meshed with *N*_{x} = 150 and *N*_{y} = 175, respectively. The parameter ρ=4.1e−3, which is 0.7 times the maximum value of |Lx|. The calculation is finished within 1 minute. Results are given in Figure 10. True values are shown in Figure 10, A and B. The reconstructed σ with stabMREPT (Figure 10C) matches quite well to true values with NRMSE <5%, although boundary artifacts are observed in the results of ϵr (Figure 10D). These results are better than those computed with cr-MREPT (Figure 10, E and F) and stdMREPT (Figure 10, G and H). Figure 11 shows a comparison of the variation of σ and ϵr along the line of *x* = 0 for stabMREPT and stdMREPT. A strong oscillation exists in the stdMREPT results.

##### Figure 10.

True σ and ϵr of the head model are given in (A) and (B). Results are shown in stabMREPT (C) and (D). Results of cr-MREPT are given in (E) and (F). Results of stdMREPT are given in (G) and (H). The first and second rows are results of σ and ϵr, respectively.

##### Figure 11.

The variation of conductivity σ (A) and relative permittivity ϵ_{r} (B) along *x* = 0 for different methods.

In practical measurement, the data of H1+ are always polluted by noise at different levels. The noise tolerance of stabMREPT is tested using the simplified human brain model. The SNR is defined as SNR=20log10|H1+|/|N|, where H1+ is the simulated data. *N* is the added complex white Gaussian noise generated with a given noise level in decibel. SNR is set to be 60 dB and ρ=4.1×10−3 in equation (5). The SG filter is used for calculating the derivatives of H1+ to smooth noisy data. Figure 12, A and B shows the results reconstructed by cr-MREPT. The accuracy of the results obtained is low. Results of stdMREPT are given in Figure 12, C and D. The results are severely degraded by the noise. Very large values and global oscillations are observed in the results. On the contrary, the noise has no strong effect on the results of σ reconstructed with stabMREPT, although the results for ϵr are slightly affected, as shown in Figure 12, E and F. Moreover, the reconstructed results are severely affected when the noise level is increased to be ≥50 dB. However, the effect of noise on stabMREPT becomes rapidly insignificant with an increase in SNR. Examples are given in Figure 12, G and H for stabMREPT with SNR=65dB, a considerable improvement is observed. In contrast, performance of cr-MREPT and stdMREPT did not improve when SNR is increased from 60 to 65 dB. The effect of noise on stabMREPT is negligible when SNR is smaller than 70 dB.

##### Figure 12.

The noise tolerance of cr-MREPT (A and B), standard MREPT (stdMREPT) (C and D), and stabMREPT (E and F). The first and the second rows are reconstructed results of σ and ϵ_{r}, respectively. From left to right, the first 3 columns are reconstructed results of cr-MREPT, stdMREPT, and stabMREPT with 60 dB noise. The last row shows the results of stabMREPT with 65 dB noise (G and H). The Savitzky–Golay (SG) filter is used to calculate partial differentiations over H1+ for smoothing noisy data.

Research Articles

Download PDF (4.21 MB)

TOMOGRAPHY, March 2017, Volume 3, Issue 1:50-59

DOI: 10.18383/j.tom.2016.00283

## An MR-Based Viscosity-Type Regularization Method for Electrical Property Tomography

Changyou Li

^{1}, Wenwei Yu^{2}, Shao Ying Huang^{3}^{1}School of Electronics and Information, Northwestern Polytechnical University, China;^{2}Center for Frontier Medical Engineering, Chiba University, Japan; and^{3}Bio-Medical Group, Engineering Product Development, Singapore University of Technology and Design, Singapore## Abstract

Here, a method based on viscosity-type regularization is proposed for magnetic resonance electrical property tomography (MREPT) to mitigate persistent artifacts when it is used to reconstruct a map of electrical properties based on data from a magnetic resonance imaging scanner. The challenges for solving the corresponding partial differential equation (PDE) are discussed in detail. The existing artifacts in the numerical results are pointed out and classified. The methods in the literature for MREPT are mainly based on an assumption of local homogeneity, which makes the approach simple but leads to artifacts in the transition region where electrical properties vary rapidly. Recent work has focused on eliminating the assumption of local homogeneity, and one of the solutions is convection–reaction MREPT that is based on a first-order PDE. Numerical solutions of the PDE have persistent artifacts in certain regions and global spurious oscillations. Here, a method based on viscosity-type regularization is proposed to effectively mitigate the aforementioned problems. Finite difference method is used for discretizing the governing PDE. Numerical experiments are presented to analyze the problem in detail. Electrical properties of different phantoms are successfully retrieved. The efficiency, accuracy, and noise tolerance of the proposed method are illustrated with numerical results.

## Introduction

Magnetic resonance imaging (MRI)-based magnetic resonance electrical property tomography (MREPT) (1, 2) is a strategy for noninvasively reconstructing electrical properties (EPs) (permittivity and electrical conductivity) of the human body without using additional hardware to a magnetic resonance (MR) scanner. It has been intensively developed in recent years for reconstructing the EPs of human tissues based on the measurable data (radiofrequency [RF] field distribution called B

_{1}-map) (3) from an MRI scanner because the B_{1}-field carries the information of EP distribution in the human body.MREPT is an important technology in various areas in medicine and biology. The constructed EP map can show the anatomical structure of the human body by the high contrast of EPs between different tissues (4). Cancerous tissues, including those at the early stages, are highly distinguishable from healthy ones (5, 6). Therefore, the EP map provides a good tool for early cancer detection. It can also provide guidance to design body-centric communications (7) and electromagnetics-based therapies (8). EP maps can also help in understanding the activities of cells, in particular with exposure to either static or dynamic electromagnetic waves (9). Moreover, an EP map is crucial for accurate calculation of specific absorption rate, which is a main parameter for assessing the risk of RF and microwave radiation and that of ultrahigh-field MRI scanning (10, 11).

MREPT was originally proposed by Haacke et al. (12) and first implemented by Wen (13). The Helmholtz equation for homogeneous material was proposed for calculating both conductivity and permittivity with positively circularly polarized component of the magnetic field (B 1 + ), which corresponds to the RF transmit field. Another MREPT method was developed and systematically studied by Katscher et al. (14) and Voigt et al. (15). The integral forms of Faraday's law and Ampère's law were applied to formulate the governing equation based on B 1 + . However, these methods are based on an assumption of local homogeneity for simplification. This assumption leads to artifacts in the region where EPs vary either quickly or abruptly (3, 16, 17). The reconstruction errors are rigorously analyzed by Seo et al. (18). A recent work was focused on removing the local homogeneity assumption (19).

One full generalization of these methods was performed by Sodickson et al., and the technique was named Local Maxwell Tomography (LMT) (20, 21). The governing equation of this method was derived from Faraday's law and Ampère's law. The unmeasurable magnetic field components are eliminated by using Gauss's law. The method is physically feasible and free of assumptions. Magnetic fields measured using an MRI scanner with multiple channels are used to solve the unknown. Although the assumption of local homogeneity is removed in LMT, higher-order partial differentiation over the magnetic field is needed in the calculation, which significantly increases this method's sensitivity to noise.

Another MREPT method (22) was proposed on the basis of the fact that the magnetic fieldB 1 + near the center of a birdcage or a transverse electromagnetic coil is homogeneous along its central axis (23, 24), and it is assumed that the contribution of the along-axis component of the magnetic field is negligible when the region of interest is close to the center of the coil. A first-order PDE including a gradient term of the unknown EP was formulated via properly manipulating Maxwell equations, and the gradient term describes the inhomogeneity of the object under investigation. This method is more general than the MREPT proposed by Katscher et al. (14) and has a simpler form than LMT.

Hafalir et al. (25) then named this method as convection–reaction MREPT (cr-MREPT), as its governing equation is similar to the convection–reaction equation (26). The first-order PDE was then solved by Hafalir et al. (25) for 2-dimensional reconstruction using a strong form with finite element method (FEM) (27). An obvious artifact was observed in the region where the gradient of| B 1 + | is close to zero. Although 2 methods, the double-excitation and the constraint MREPT, were further proposed by Hafalir et al. (25) to improve the results, artifacts reappear when noise exists. Meanwhile, cr-MREPT does not work with every phantom. Even with the phantoms used by the authors, the reconstructed EPs deviate seriously from the true values in a large region close to the central artifact. High sensitivity of noise was also observed.

In parallel with the work of Hafalir et al., Ammari et al. (19) also attempted to solve this PDE via an adjoint-based iterative reconstruction algorithm. The real and imaginary parts of the equation are separated to produce a coupled equation pair solved by the iterative scheme. This method can construct a blur image of the EPs without the aforementioned artifacts, but it is not applicable in the region where the gradient of| B 1 + | field is close to zero, and the produced linear system is ill-conditioned. The method also converges slower than the FEM-based method (19).

The solution of the first-order PDE shows persistent artifacts where∇ | B 1 + | approaches zero. Other than that, global spurious oscillations exist in the solution of this first-order PDE. Therefore, directly solving it with numerical methods is difficult. Here, the existing challenges for solving the PDE are discussed in detail. The existing artifacts in the numerical results are pointed out and classified. An effective viscosity-type regularization (28) solution is applied to mitigate the aforementioned problems. A Laplacian term is introduced into the first-order PDE for stabilizing the entire system, which can be solved by either the finite difference (FD) method or FEM. The governing PDE is in the first order, so its solution will be a linear function that depends on the space coordinates, and it can present sharp variations in space. Because the PDE is not a singularly disturbed equation (29, 26, 30), the introduced second-order Laplacian term will modify the PDE to obtain a second-order solution, which is an approximation of its exact solution and presents a smoother property in space.

The derivation of the formulation is briefly presented for the convenience of discussion; behaviors of the equation are analyzed in detail before introducing the viscosity-type regularization. Numerical experiments are used to facilitate the analysis of the existing problems and to show the effectiveness of the proposed approach. Finite difference method is applied for discretization. The proposed method is applied for reconstructing the EPs of different phantoms. Accuracy, stability, and noise tolerance of the method are illustrated with numerical results.

## Theory and Methodology

The governing PDE for cr-MREPT are given as follows:

## (1)

This equation can be solved for reconstructing EPs of a 3-dimensional object (see online supplemental Appendix for the relating derivation of this equation). If γ does not vary severely along the∂ γ ∂ z becomes negligible. Equation (1) is then simplified as follows for a single section located at the center of birdcage coil:

zdirection,## (2)

z-component, respectively.H_{z}is not measurable using MRI, but its contribution is negligible compared withH_{z}is generally close to zero in the center of the birdcage coil. Therefore equation (2) could be further simplified as follows:## (3)

For homogeneous regions,∇ x y γ = 0 ; therefore, equation (3) can be simplified as γ = i ω μ H 1 + / ∇ 2 H 1 + , which is the governing equation proposed by Haacke et al. and Wen (12, 13), named, in this content, standard MREPT (stdMREPT).

Therefore, the original inverse problem equivalently becomes a problem of solving the first-order PDE (3). It has the form of a stationary convection–reaction equation, which is a mathematical model for fluid dynamics (31). Thus, the MREPT based on equation (3) is called cr-MREPT by Hafalir et al. (25).

The coefficients of PDE (3) are complex, making the PDE exhibit different mathematical and numerical behaviors from the stationary convection–reaction equations that have real coefficients. In particular, both equations show global spurious oscillations in their solutions (26). For convection–reaction equations (real coefficients), these global spurious oscillations are caused by the mesh density for discretization, and the results converge better when the mesh size decreases. However, this fact does not hold for equation (3), where the coefficients are complex. Instead, the spurious oscillation in equation (3) is caused by a fundamental mathematical drawback of itself; this is further discussed as follows.

Equation (3) can be reformed as follows:

## (4)

## Figure 3.

σ (A) andϵ r (B) of phantom shown in Figure 2 are reconstructed by applying convection–reaction magnetic resonance electrical property tomography (MREPT) (cr-MREPT) based on equation (3). σ (C) and ϵ r (D) reconstructed using equation (2). Real (E) and imaginary (F) parts of H 1 + , absolute value of L x (G) and ∇ 2 H 1 + (H), are also presented.

Here, a viscosity-type regularization method (28) is applied to equation (3) to eliminate the global spurious oscillation in its solution. The governing equation becomes the following equation:

## (5)

x-component of## Discretization Strategy

A square domain denoted as Ω and bounded by∂ Ω is discretized by small squares with N x + 2 and N y + 2 grid points along the ( N x + 2 ) × ( N y + 2 ) data points in space. Δ x and Δ y are the step lengths along i : = { i | i ∈ Z , 0 ≤ i ≤ N x + 1 } and { j | j ∈ Z , 0 ≤ j ≤ N y + 1 } , respectively. The values of H 1 + and γ at the point (H { i , j } and γ { i , j } , respectively.

xandydirections, respectively, givingxandydirections, respectively. Useiandjto index the mesh points alongxandywithj:=i,j) are represented asAllowingL x y = [ L x , L y ] with L x = − ∂ H 1 + / ∂ x + i ∂ H 1 + / ∂ y and L y = − i ∂ H 1 + / ∂ x − ∂ H 1 + / ∂ y . The values of L x and L y at point (L x { i , j } and L y { i , j } , respectively. Using first-order central difference formula, L x { i , j } and L y { i , j } are given as follows:

i,j) are represented with## (6)

## (7)

Here, the imaginary unit∇ 2 H 1 + at the point (d 2 H { i , j } , the second-order central difference formula is used to discretize ∇ 2 H 1 + , which gives the following equation:

ican be easily distinguished from the indexion the basis of its location. Representing the data ofi,j) as## (8)

kis the index for sections along thez-direction. Two extra sections atz.Another way for calculatingL x { i , j } , L y { i , j } , and d 2 H { i , j } is by using the Savitzky–Golay (SG) filter (34), which has been widely used to smooth noisy data. Here, H 1 + is approximated with the following equation:

## (9)

x,y) is the corresponding coordinate ofa_{p},p= 1, . ., 6 are polynomial coefficients, which need to be decided with the measured data ofi,j),a_{p},p= 1, …, 6. A least-square method could then be used for solving the unknown.## Figure 1.

Neighbors ofH { i , j } in a mesh grid.

## Figure 4.

σ (A) andϵ r (B) reconstructed via least-square minimization with λ = 0.01 max { A } for the phantom shown in Figure 2. Central difference is used for ∇ 2 H 1 + .

The derivatives overH 1 + at the point (

i,j) are then calculated using the following equations:## (10)

## (11)

## (12)

Here, (H 1 + includes the second-order partial differentiation of H 1 + over − Δ z and Δ z . With equations (6)–(8), or (10)–(12), equation (5) can be easily discretized with the first- and second-order central difference formula, as follows:

x_{i},y_{j}) is the coordinate of the grid point (i,j). Notice that the Laplacian term overz, which is calculated with a central difference formula, with the data at sections located at## (13)

A linear system is produced via assembling (13) for all points indexed by (A g = b . The sparsity of the tridiagonal matrix A allows efficient calculation with Gaussian elimination. The Dirichlet boundary condition is assumed here. Vector γ { i , j } . The relative permittivity and conductivity at (γ { i , j } .

i,j). The linear system is denoted asgrepresents the unknownsi,j) can be calculated based on the real and imaginary parts of## Numerical Experiments

In the calculation,H 1 + as input is simulated with FEM-based software COMSOL (35). A 16-leg high-pass shielded birdcage coil with a radius of 14.5 cm and a leg-length of 24 cm (36) is built to work at 127.74 MHz (3 T MRI system) to produce a homogeneous H 1 + within the coil, and it is placed with its axis along the

z-axis. The coil is driven by 2 ports with 500 V for each. Both ports are geometrically 90° apart with a 90° phase difference.The phantom of 2 coaxial cylinders used for simulation is shown in Figure 2A. The cylinders have the same heighth = 24 cm . The radii of the exterior and interior ones are r 1 = 5 cm and r 2 = 2.5 cm , respectively. The exterior cylinder is filled by a material with ϵ r 1 = 75 and σ 1 = 0.5 S/m, whereas for the interior one, ϵ r 2 = 50 and σ 2 = 1 S/m.

## Figure 2.

The phantom of 2 cylinders with the 3D view (A) and its cross section view (B).

h= 24 cm,r_{1}= 5 cm, andr_{2}= 2.5 cm. Dashed-line square defines region of interest. Unit of sigma is S/m.## Figure 5.

σ (A) andϵ r (B) of the phantom in Figure 2 are reconstructed by cr-MREPT with

N_{x}=N_{y}= 199.## Figure 6.

Reconstructed results of σ (A) andϵ r (B) for the 2-cylinder phantom shown in Figure 2 with stabilized MREPT (stabMREPT). δ = 0.37 ,

N_{x}=N_{y}= 199.The phantom is placed in the birdcage coil, with its axis and center coinciding with those of the coil. A Dell Precision workstation (Round Rock, Texas) with Intel Xeon CPU with 2 processors (4 cores for each) and 48 GB RAM is used for simulation. The phantom is meshed with tetrahedral elements. The maximum mesh size for− 0.5 cm < z < 0.5 cm is 1 mm. The mesh size is then increased gradually inside out. This strategy guarantees the accuracy of the calculated H 1 + at the sections of interest and minimizes the number of elements to reduce the required memory. H 1 + of the sections at z = 0.25 cm is extracted as inputs for the reconstruction.

z= 0,z= −0.25 cm, andThe domain Ω is defined as the square inscribed to the cross section of the exterior cylinder, indicated by the dashed lines in Figure 2B. It is then meshed withN x = N y = 99 . The data of H 1 + are directly extracted from COMOSL.

The mesh used here will lead to a space resolution of 0.75 mm. Because simulation data are used here, not much concerns are needed. But, in practice, this resolution will need a long MR scan and produce low signal-to-noise ratio (SNR). A coarse mesh can be used first to perform the MR scan. With the obtainedH 1 + , one can easily apply the cyclic regularized Savitzky–Golay filter (37) to reconstruct H 1 + of a finer mesh with the required resolution. This will also help in reducing noise.

Relative permittivityϵ r and conductivity δ of phantom shown in Figure 2 are first reconstructed by equation (3), the cr-MREPT method. Differentiations of H 1 + are calculated with the central difference formula. As shown in Figure 3, A and B, significant artifacts are observed around the center of the interior cylinder and in the regions where σ and ϵ r are discontinuous. Furthermore, global spurious oscillations are also clearly seen. For checking the effect of neglecting the contribution of

H_{z}, equation (2) is also applied for the reconstruction. Results are shown in Figure 3, C and D, and not much improvement is seen. The global spurious oscillations and the artifacts still exist.Hafalir et al. (25) report that the central artifacts are caused by the fact thatL x and L y are close to 0. Figure 3, E and F shows the real and imaginary parts of the simulated H 1 + at the central section (R H 1 + and I H 1 + , achieve to their extrema in the vicinity of the center, as shown in the figures, indicating that ∂ H 1 + / ∂ x and ∂ H 1 + / ∂ y are close to zero. This results in L x { i , j } ≈ 0 and L y { i , j } ≈ 0 for the points in the center region. The absolute values of L x are shown in Figure 3G (| L x | = | L y | ).

z= 0), denoted asA minimization method is also applied to mitigate the global spurious oscillation when solving equation (3). The original linear system produced by equation (3) is written asA g = b , which is then transformed into a least-square minimization problem min | | b − A g | | 2 + λ | | g | | 2 with a penalty term λ | | g | | 2 . The damping parameter λ > 0 . In this way, the solution for A here is about 10( A T A + λ I ) g = A T b , A being not symmetric, is not suitable for the calculation (it is equivalent to the applied minimization problem), because the 2-norm condition number of A T A is exactly the square of the one of A .

gis regularized. It is then solved with the iterative algorithm LSQR, which is based on the Golub–Kahan bidiagonalization process (38). The condition number of^{3}, soThis minimization is conducted withλ = 0.01 max | | A | | 2 , which is a large value for λ. As shown in Figure 4, the results converge better to the true values, but an obvious “4-blade” artifact is observed. These global spurious oscillations remain.

Moreover, this global spurious oscillation is independent of grid sizesΔ x and Δ y . Figure 5, A and B show the results calculated with N x = N y = 199 , which is a much finer mesh in the region Ω. Global spurious oscillation still exists. Many artifacts appear with very large values. The reconstructed values of | ϵ r | and | σ | achieve to 500 and 10.

Therefore, neither directly the solution of equations (2) and (3) nor that of the least-square minimization method produces accurate results. This is expected by the fundamental drawback existing in the governing PDE, which is explained in section Theory and Methodology. However, reconstruction with the proposed stabMREPT based on the new second-order PDE (5) shows much better performance, which is shown as follows.

The same phantom (shown in Figure 2) is used for the proposed stabMREPT. The region Ω is discretized with the same mesh density (N x = N y = 199 ). The parameter ρ = 0.37 , which is 7 times of the maximum value of | L x | . Figure 6 shows the reconstructed conductivities and permittivities. The normalized root-mean-square error (NRMSE) is defined as % error = 100 | | ( e true − e recon ) | | 2 / | | e true | | 2 . The reconstructed results of σ and ϵ r are very close to the true values. Comparison of the reconstructed σ and ϵ r with their true values at the line

y= 0 yields the plots shown in Figure 7. A very good match is observed. Relative errors are all <1% except the boundary where the electrical parameters are discontinuous. Moreover, includingH_{z}into equation (5) produces the same results.## Figure 7.

Comparison between the true values and the reconstructed ones of σ (A) andϵ r (B) along δ = 0.37 ,

y= 0. stabMREPT is applied for the 2-cylinder phantom. The regularization parameterN_{x}=N_{y}= 199.## Figure 8.

The map of normalized root-mean-square error (NRMSE) for σ reconstructed with stdMREPT (A) and stabMREPT (B). The 2-cylinder phantom is used.δ = 0.37 ,

N_{x}=N_{y}= 199.## Figure 9.

Simplified head model with 5 different regions characterized with electrical parameters of scalp (ϵ r = 62 , σ = 0.54 S/m, gray), skull (ϵ r = 21 , σ = 0.12 S/m, deep blue), gray matter (ϵ r = 73 , σ = 0.59 S/m, yellow), white matter (ϵ r = 52 , σ = 0.34 S/m, light blue), and cerebrospinal fluid (ϵ r = 84 , σ = 2.14 S/m, red).

To further check the performance of stabMREPT compared with stdMREPT at both the boundary and at the homogeneous regions (no abrupt change in EPs), the NRMSE is calculated here for the reconstructed σ of the 2-cylinder phantom. The parameters for the reconstruction are the same as those used for Figure 6. The map of NRMSE is shown in Figure 8. A reconstruction error of stdMREPT at the boundary of interior cylinder is >100%, whereas a large error is also observed in the region far from the boundary. A reconstruction error at the boundary of the cylinder is also observed in the results obtained with stabMREPT, which is about 45%, and this error decreases gradually along the radius. But the reconstruction error at the region far away from the boundary is much smaller than the one at stdMREPT, where the largest error is <3%. Therefore, stabMREPT performs better than stdMREPT in both the boundary and homogeneous regions.

To further test the proposed approach, a simplified human head phantom as shown in Figure 9 is simulated. The model has 5 different regions indicated by different colors. Each region relates to a tissue type of the human brain (4). These tissues are, from external to internal, scalp (ϵ r = 62 , σ = 0.54 S/m ), skull (ϵ r = 21 , σ = 0.12 S/m ), gray matter (ϵ r = 73 , σ = 0.59 S/m ), white matter (ϵ r = 52 , σ = 0.34 S/m ), and cerebrospinal fluid (ϵ r = 84 , σ = 2.14 S/m ) at 127.74 MHz (39). For the convenience of calculation, the region indicated with gray color is set to a 12- × 14-cm (length × width) square). The 2-dimensional section shown in Figure 9 is drawn and then extruded in COMSOL to produce a 1-cm-thick model. The thickness is primarily limited by the workstation's RAM. The curved edges and the central red region are meshed with an element size smaller than 1 mm, and other regions are meshed with an element size smaller than 1.5 mm.

The region Ω defined by −6 cm <ρ = 4.1 e − 3 , which is 0.7 times the maximum value of | L x | . The calculation is finished within 1 minute. Results are given in Figure 10. True values are shown in Figure 10, A and B. The reconstructed σ with stabMREPT (Figure 10C) matches quite well to true values with NRMSE < 5 % , although boundary artifacts are observed in the results of ϵ r (Figure 10D). These results are better than those computed with cr-MREPT (Figure 10, E and F) and stdMREPT (Figure 10, G and H). Figure 11 shows a comparison of the variation of σ and ϵ r along the line of

x< 6 cm and −7 cm <y< 7 cm is meshed withN_{x}= 150 andN_{y}= 175, respectively. The parameterx= 0 for stabMREPT and stdMREPT. A strong oscillation exists in the stdMREPT results.## Figure 10.

True σ andϵ r of the head model are given in (A) and (B). Results are shown in stabMREPT (C) and (D). Results of cr-MREPT are given in (E) and (F). Results of stdMREPT are given in (G) and (H). The first and second rows are results of σ and ϵ r , respectively.

## Figure 11.

The variation of conductivity σ (A) and relative permittivity ϵ

_{r}(B) alongx= 0 for different methods.In practical measurement, the data ofH 1 + are always polluted by noise at different levels. The noise tolerance of stabMREPT is tested using the simplified human brain model. The SNR is defined as SNR = 20 log 10 | H 1 + | / | N | , where H 1 + is the simulated data. SNR is set to be 60 dB and ρ = 4.1 × 10 − 3 in equation (5). The SG filter is used for calculating the derivatives of H 1 + to smooth noisy data. Figure 12, A and B shows the results reconstructed by cr-MREPT. The accuracy of the results obtained is low. Results of stdMREPT are given in Figure 12, C and D. The results are severely degraded by the noise. Very large values and global oscillations are observed in the results. On the contrary, the noise has no strong effect on the results of σ reconstructed with stabMREPT, although the results for ϵ r are slightly affected, as shown in Figure 12, E and F. Moreover, the reconstructed results are severely affected when the noise level is increased to be ≥50 dB. However, the effect of noise on stabMREPT becomes rapidly insignificant with an increase in SNR. Examples are given in Figure 12, G and H for stabMREPT with SNR = 65 dB , a considerable improvement is observed. In contrast, performance of cr-MREPT and stdMREPT did not improve when SNR is increased from 60 to 65 dB. The effect of noise on stabMREPT is negligible when SNR is smaller than 70 dB.

Nis the added complex white Gaussian noise generated with a given noise level in decibel.## Figure 12.

The noise tolerance of cr-MREPT (A and B), standard MREPT (stdMREPT) (C and D), and stabMREPT (E and F). The first and the second rows are reconstructed results of σ and ϵH 1 + for smoothing noisy data.

, respectively. From left to right, the first 3 columns are reconstructed results of cr-MREPT, stdMREPT, and stabMREPT with 60 dB noise. The last row shows the results of stabMREPT with 65 dB noise (G and H). The Savitzky–Golay (SG) filter is used to calculate partial differentiations over_{r}## Discussion and Conclusion

Here, a viscosity-type regularization method is used to effectively improve the performance of the formerly proposed cr-MREPT, which is based on a first-order complex coefficient PDE. The proposed method is then named as stabMREPT. The regularization is implemented on the basis of the fact that adding the Laplacian term does not lead to a singularly disturbed problem. The PDE is discretized using FD method.

The reconstruction method based on the same PDE was named as cr-MREPT by Hafalir et al. (25). However, several drawbacks exist in this method, although it is rarely reported in the literature. Here, these drawbacks are analyzed in detail. Different numerical simulations are given for the purpose of illustration; 2 PDE properties are coupled. The governing equation of cr-MREPT is equation (1), which is in the form of pure convection–reaction PDE. It is well known in the literature that it is tricky to numerically solve this kind of PDE (40). With the complex coefficients, this kind of PDE has a smooth solution with some conditions, as discussed in Theory and Methodology. Therefore, directly using the convection–reaction PDE for reconstruction can be difficult.

The viscosity-type regularization method is then suggested. As illustrated in the paper, this method can efficiently mitigate the persistent artifacts existing in the convection–reaction PDE. The regularization parameter ρ needs to be positive and smaller than 1. It can be chosen as several times or a fraction of the maximum of| L x | , which is the L x y . Blurring can be observed at the boundaries of different tissues, but the reconstruction accuracy can be considerably better. The accuracy of the proposed method is shown by numerical results obtained by using simple and complex phantoms. Good performances in several numerical experiments are observed by comparing them with other existing methods.

x-component ofAlthough the noise tolerance is not yet good for practical use, it is already considerably better than stdMREPT and cr-MREPT. The noise tolerance of this method was compared with that of other methods. Another algorithm based on FEM will be implemented in the next step to add the regularization term only when it is needed. A regularization-based SG method (37) can be used for reconstructingH 1 + from the measured noisy data, which will further reduce the influence of noise and increase the noise tolerance of the method.

## Supplemental Materials

## Supplemental Appendix:

http://dx.doi.org/10.18383/j.tom.2016.00283.sup.01

## Notes

[1] Abbreviations:

MREPT

Magnetic resonance electrical property tomography

PDE

partial differential equation

MRI

magnetic resonance imaging

EPs

electrical properties

MR

magnetic resonance

LMT

Local Maxwell Tomography

cr-MREPT

convection–reaction MREPT

FD

finite difference

FEM

finite element method

stdMREPT

standard MREPT

stabMREPT

stabilized MREPT

SNR

signal-to-noise ratio

NRMSE

normalized root-mean-square error

SG

Savitzky–Golay

cr-MREPT

convection–reaction MREPT

## References

## Journal Information

Journal ID (nlm-ta): tom

Journal ID (publisher-id): TOMOG

Title: Tomography

Subtitle: A Journal for Imaging Research

Abbreviated Title: Tomography

ISSN (print): 2379-1381

ISSN (electronic): 2379-139X

Publisher: Grapho Publications, LLC (Ann Abor, Michigan)

## Article Information

Copyright statement: © 2017 The Authors. Published by Grapho Publications, LLC

Copyright: 2017, Grapho Publications, LLC

License (open-access, http://creativecommons.org/licenses/by-nc-nd/4.0/):

This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Publication date (print): March 2017

Volume: 3

Issue: 1

Pages: 50-59

Publisher ID: TOM-00283-16

DOI: 10.18383/j.tom.2016.00283

## Supplemental Media

Supplemental Media:Appendix A: Derivation of the PDE for cr-MREPT.View this media larger in a new window

## PDF

Download the article PDF (4.21 MB)

Download the full issue PDF (141.48 MB)

## Mobile-ready Flipbook

View the full issue as a flipbook (Desktop and Mobile-ready)