# Magnetic Resonance-Electrical Properties Tomography by Directly Solving Maxwell's Curl Equations
## 2 Theory and Method
### Optimization  
EPs were estimated through an iterative optimization procedure. The optimization variables (EPs) are \(\mathbf{x}_{\text{EP}}=[\epsilon_{\text{tissue}}\ \mathbf{\sigma}_{\text{tissue}}]\), that is, the permittivity and conductivity values of each voxel of the imaging subject. Using the developed FDTD solver described in Section 2.1, the tissue dielectric properties variable \(\mathbf{x}_{\text{EP}}\) were mapped to the \(B_{1}^{+}\) field as:  
\[B_{1}^{+}=\mathcal{T}_{\text{FDTD}}(\mathbf{x}_{\text{EP}}), \tag{1}\]  
where \(\mathcal{T}_{\text{FDTD}}\) denotes the FDTD process.  
Figure 2: The finite-difference-time-domain (FDTD) configuration. Here, a 12-leg birdcage coil was modelled and loaded with a subject (phantom or head); a perfect electric conductor (PEC) was placed around the birdcage coil to simulate the RF shield, and perfectly matched layers (PML) boundary conditions were placed at the top and bottom of the computation domain. With the graphics processing unit (GPU) computing, the FDTD calculation of the electromagnetic (EM) fields only took about 18 s for the given case involving 1.4 million Yee-cells.  
Then, a cost function which describes the residual between the mapped \(B_{1}^{+}\) fields and the measured \(B_{1}^{+}\) field can be built as:  
\[C=\left\|B_{1}^{+}\right|-\left|B_{1;\text{meas}}^{+}\right|\right\|_{2}, \tag{2}\]  
Therefore, an unconstrained optimization problem which aims to minimize the cost function can be formed as:  
\[\min_{x}\left\|\left|\mathcal{F}_{\text{FDTD}}(\mathbf{x}_{\text{EP}})\right|- \left|B_{1;\text{meas}}^{+}\right|\right\|_{2}, \tag{3}\]  
The problem in Equation (3) can be efficiently solved by using optimization methods such as quasi-Newton or trust-region techniques. The EPs of the imaging subject can be obtained after the optimum \(\mathbf{x}_{\text{EP}}\) is found from Equation (3).  
In this work, the genetic algorithm (GA, standard MATLAB function) was used to find the optimal solution of \(\mathbf{x}_{\text{EP}}\). In the optimization procedure, an initial EP distribution denoted as \(\mathbf{x}_{\text{O}}\) was set based on the LHA. Inside the uniform region, the EP distribution was directly derived from the LHA-based result. On the tissue boundaries, average EPs of the neighboring tissues were used. This operation avoided the boundary artifacts in the LHA-based results. The variable range in the GA was set as \(\mathbf{x}_{\text{O}}\pm 0.3\overline{\mathbf{x}}_{\text{O}}\), where \(\overline{\mathbf{x}}_{\text{O}}\) is the mean value of \(\mathbf{x}_{\text{O}}\). Alternatively, the initial EPs of each tissue can be set based on the data from Reference [32], and the variable range was set as \(\mathbf{x}_{\text{O}}\pm 0.3\mathbf{x}_{\text{O}}\), where \(\mathbf{x}_{\text{O}}\) are the reported values in Reference [33]. The GA-based optimization was implemented to satisfy one of the following stop criteria: (1) the fitness, \(C\), was less than the expected fitness (least-square residual error), \(C_{\text{desired}}\), and (2) the number of iterations exceeds certain pre-set limits, \(N_{\text{desired}}\). Once condition (1) or (2) was satisfied, the optimization process could be terminated. The parameter setup of \(C_{\text{desired}}\) and \(N_{\text{desired}}\) is problem-dependent. In this study, it was observed that the variation of residual error was negligible after 50 iterations, and when it reached a value of around \(10^{-3}\). Therefore, we chose \(C_{\text{desired}}\) and \(N_{\text{desired}}\) as \(10^{-3}\) and 100, respectively.  
The iterative optimization procedure is implemented in MATLAB (R2015a, 64 bits, The MathWorks, Inc., Natick, MA, USA).