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 B1-map) (3) from an MRI scanner because the B1-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 (), 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 . 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 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 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 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 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, becomes negligible. Equation (1) is then simplified as follows for a single section located at the center of birdcage coil:
and are vectors and without the z-component, respectively.
. Notice that still needs to be evaluated in a 3-dimensional form. Magnetic field Hz is not measurable using MRI, but its contribution is negligible compared with , as Hz is generally close to zero in the center of the birdcage coil. Therefore equation (2) could be further simplified as follows:
For homogeneous regions, ; therefore, equation (3) can be simplified as , 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:32), a necessary condition for equation (4) to have a smooth solution is that is analytic in region Ω (19). However, is not smooth in the region where electrical parameters vary rapidly (as in Figure 3H for a 2-layer cylindrical phantom with difference and σ in different layers). Discontinuity of exists at the boundary where electrical parameters vary rapidly. Therefore, a smooth solution of σ and cannot be obtained with equation (3), and a better strategy is needed for good numerical stability, accuracy, and efficiency. 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 , which is the x-component of . 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.
A square domain denoted as Ω and bounded by is discretized by small squares with and grid points along the x and y directions, respectively, giving data points in space. and are the step lengths along x and y directions, respectively. Use i and j to index the mesh points along x and y with and j := , respectively. The values of and γ at the point (i, j) are represented as and , respectively.
Allowing with and . The values of and at point (i, j) are represented with and , respectively. Using first-order central difference formula, and 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 at the point (i, j) as , the second-order central difference formula is used to discretize , which gives the following equation:
k is the index for sections along the z-direction. Two extra sections at and away from the central section are used to calculate the partial differentiation of over z.
Another way for calculating , and is by using the Savitzky–Golay (SG) filter (34), which has been widely used to smooth noisy data. Here, is approximated with the following equation:Figure 1. These neighbors or itself could be used for interpolation with Stanley et al. (9), which will produce an overfitting linear system about the unknown coefficients ap, p = 1, …, 6. A least-square method could then be used for solving the unknown.
The derivatives over at the point (i, j) are then calculated using the following equations:
Here, (xi, yj) is the coordinate of the grid point (i, j). Notice that the Laplacian term over includes the second-order partial differentiation of over z, which is calculated with a central difference formula, with the data at sections located at and . 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 . The sparsity of the tridiagonal matrix allows efficient calculation with Gaussian elimination. The Dirichlet boundary condition is assumed here. Vector g represents the unknowns . The relative permittivity and conductivity at (i, j) can be calculated based on the real and imaginary parts of .
In the calculation, 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 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 . The radii of the exterior and interior ones are and , respectively. The exterior cylinder is filled by a material with and S/m, whereas for the interior one, and S/m.
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 is 1 mm. The mesh size is then increased gradually inside out. This strategy guarantees the accuracy of the calculated at the sections of interest and minimizes the number of elements to reduce the required memory. of the sections at z = 0, z = −0.25 cm, and 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 . The data of 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 , one can easily apply the cyclic regularized Savitzky–Golay filter (37) to reconstruct of a finer mesh with the required resolution. This will also help in reducing noise.
Relative permittivity and conductivity δ of phantom shown in Figure 2 are first reconstructed by equation (3), the cr-MREPT method. Differentiations of 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 are discontinuous. Furthermore, global spurious oscillations are also clearly seen. For checking the effect of neglecting the contribution of Hz, 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 and are close to 0. Figure 3, E and F shows the real and imaginary parts of the simulated at the central section (z = 0), denoted as and , achieve to their extrema in the vicinity of the center, as shown in the figures, indicating that and are close to zero. This results in and for the points in the center region. The absolute values of are shown in Figure 3G ().
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 , which is then transformed into a least-square minimization problem with a penalty term . The damping parameter . 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 here is about 103, so being not symmetric, is not suitable for the calculation (it is equivalent to the applied minimization problem), because the 2-norm condition number of is exactly the square of the one of .
This minimization is conducted with , 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 and . Figure 5, A and B show the results calculated with , 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 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 (). The parameter , which is 7 times of the maximum value of . Figure 6 shows the reconstructed conductivities and permittivities. The normalized root-mean-square error (NRMSE) is defined as . The reconstructed results of σ and are very close to the true values. Comparison of the reconstructed σ and 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 Hz into equation (5) produces the same results.
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 (), skull (), gray matter (), white matter (), and cerebrospinal fluid () 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 Nx = 150 and Ny = 175, respectively. The parameter , which is 0.7 times the maximum value of . 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 , although boundary artifacts are observed in the results of (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 along the line of x = 0 for stabMREPT and stdMREPT. A strong oscillation exists in the stdMREPT results.
In practical measurement, the data of 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 , where is the simulated data. N is the added complex white Gaussian noise generated with a given noise level in decibel. is set to be 60 dB and in equation (5). The SG filter is used for calculating the derivatives of 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 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 , 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.
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 , which is the x-component of . 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.
Although 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 reconstructing from the measured noisy data, which will further reduce the influence of noise and increase the noise tolerance of the method.