Electrical Properties Tomography: A Methodological Review

Abstract

Electrical properties tomography (EPT) is an imaging method that uses a magnetic resonance (MR) system to non-invasively determine the spatial distribution of the conductivity and permittivity of the imaged object. This manuscript starts by providing clear definitions about the data required for, and acquired in, EPT, followed by comprehensively formulating the physical equations underlying a large number of analytical EPT techniques. This thorough mathematical overview of EPT harmonizes several EPT techniques in a single type of formulation and gives insight into how they act on the data and what their data requirements are. Furthermore, the review describes machine learning-based algorithms. Matlab code of several differential and iterative integral methods is available upon request.

1. Introduction

The electrical properties (EPs; conductivity $\sigma$ and permittivity $\epsilon$) of tissue have the potential to be used as biomarkers in many clinical applications. Tissue EPs depend on the tissue structure and composition. The conductivity varies largely as a function of fluid volumes and ionic concentrations, while the permittivity is largely influenced by the cellular membrane extent [1]. Cancer causes local changes of EPs relative to healthy tissues. The EPs of benign tissue compared to tumors are significantly different and have been reported to offer advantages in separating them from each other [2,3,4]. Similarly, the conductivity in cerebral ischemia is significantly decreased [5,6]. Conductivity measurements can therefore be helpful for better characterization of brain tumors [7,8], but they have also shown promising results for pelvic tumors [9], breast cancer [10] and ischemic stroke [11,12]. Knowledge of the EPs additionally allows for the calculation of the electromagnetic (EM) fields inside tissue. This makes them interesting for a wide range of clinical applications, such as electroencephalography (EEG) and electrocardiography (ECG) measurements to accurately localize internal electrical activities, deep brain stimulation to mitigate Parkinson’s disease symptoms, radio frequency (RF) ablation to remove arrhythmic genesis foci and RF hyperthermia for cancer treatment [13]. Additionally, they are critical to accurately determine the specific absorption rate (tissue heating) induced by EM waves [1].

Several EP map** approaches are explored to map the electrical properties of tissue in vivo. Electrical impedance tomography (EIT), for example, uses electrode mounting to detect currents injected into the sample [14]. This method is cost-effective and yields high temporal resolution, but poor spatial resolution due to the ill-posed nature of the inverse problem [13,15]. Magnetic induced tomography (MIT) applies an oscillating magnetic field to induce eddy currents in the object and detects the resulting magnetic fields outside the object [16]. However, it suffers from the same issues as EIT. Magnetic resonance electrical impedance tomography (MR-EIT) utilizes MRI to detect the magnetic field induced by the probing current [17]. This provides higher spatial resolution, but has a poor signal-to-noise ratio due to limitations on the amount of current injection [13,18,19]. Hall effect imaging (HEI) induces currents through surface electrodes and detects the emitted acoustic wave to reconstruct EPs [20]. This also has the potential to reach high resolution images, but all of the current injection based methods may suffer from shielding artifacts of non-conductivity tissue. Magneto-acoustic tomography with magnetic induction (MAT-MI) circumvents this shielding problem by inducing acoustic signals with time varying magnetic fields which are detected with ultrasound measurements [21]. However, methods that involve acoustic measurements are often limited to the surface of the object.

Electrical properties tomography (EPT) non-invasively images the conductivity and permittivity maps (simultaneously) in vivo from the radio frequency field signals obtained with MRI. The method does not require electrode mounting, does not induce additional external energy other than the inherent RF fields, and the RF fields can easily penetrate into most biological tissue. It uses a standard MRI system with regular RF coils. This concept was first introduced in 1991 by Haacke et al. [22] and first demonstrated in 2003 by Wen et al. [23]. The topic, however, only recently gained considerable interest by various research groups [1,13,19,24].

Several review papers discuss existing methods and review clinical applications [1,13,19,24]. These reviews, however, do not discuss the mathematical methodology in depth, which hampers the overview in terms of intrinsic assumptions. Here, however, a mathematical description of the acquired MR data and several differential and integral EPT approaches are described more thoroughly to give clear insights into the relations and differences between a large number of methods. This review thereby allows accurate comparisons between different methods and outlines their relative strengths and weaknesses. Extensions and generalizations are also mentioned. The EPT approaches are harmonized in terms of mathematical formulation, while maintaining as much as possible the structure of the original implementations to keep the transition from this manuscript to the references straightforward.

This review manuscript is organized as follows. First, general RF background information is presented together with a general formulation of an acquired MR image, which are necessary to understand some of the problems arising in EPT data acquisition. Second, some fundamental EPT equations are presented, from which the bulk of the analytical EPT approaches can be derived. The subsequent two sections discuss a large number of physical model-based EPT strategies, starting with methods that are based on transmit field data, followed by receive field-based methods. The review continues with a discussion about training-based EPT approaches. Finally, a general discussion is provided and EPT reconstruction examples are presented.

2. Phasor Representations for the RF Field

In EPT, knowledge about the RF field within the body is used to retrieve the dielectric properties (conductivity and permittivity) of tissue. This RF field is called the $B_1$ field and phasors are typically used to describe its behavior. What may lead to confusion is that there are actually two time conventions in common use to represent the RF field in terms of phasors. In particular, for a given time-domain RF field ${\mathit{B}}_1(\mathit{r},t)$ operating at a frequency $\omega >0$, phasors are introduced via the representation

$${\mathit{B}}_1(\mathit{r},t)=Re\left[{\widehat{\mathit{B}}}_1(\mathit{r},-\mathrm{i}\omega )exp(-\mathrm{i}\omega t)\right]$$

(1)

or alternatively the representation

$${\mathit{B}}_1(\mathit{r},t)=Re\left[{\widehat{\mathit{B}}}_1(\mathit{r},\mathrm{i}\omega )exp\left(\mathrm{i}\omega t\right)\right]$$

(2)

is used to describe the RF field. The vector ${\widehat{\mathit{B}}}_1(\mathit{r},-\mathrm{i}\omega )$ is the phasor of the RF field when the time factor $exp(-\mathrm{i}\omega t)$ is used, while ${\widehat{\mathit{B}}}_1(\mathit{r},\mathrm{i}\omega )$ is the phasor of the RF field in the situation where a time factor $exp\left(\mathrm{i}\omega t\right)$ is used.

For a given time convention, the phasor that corresponds to a given RF field is unique, and, since the time-domain RF field ${\mathit{B}}_1(\mathit{r},t)$ is real-valued, the phasors of the two representations are related by

$${\widehat{\mathit{B}}}_1^{\ast }(\mathit{r},-\mathrm{i}\omega )={\widehat{\mathit{B}}}_1(\mathit{r},\mathrm{i}\omega ),$$

(3)

where the asterisk denotes complex conjugation. In other words, for a given RF field, the phasors of the two representations are the complex conjugate of each other. We note that if Equation (2) is used to represent RF fields, the letter j is often used for the imaginary unit instead of the letter i. We adopt this notation here as well and write the representation of Equation (2) as

$${\mathit{B}}_1(\mathit{r},t)=Re\left[{\widehat{\mathit{B}}}_1(\mathit{r},\mathrm{j}\omega )exp\left(\mathrm{j}\omega t\right)\right].$$

(4)

Unless otherwise stated, we use the phasor representation of Equation (4) to describe the RF field.

Suppose now that we orient our reference frame such that the static background field ${\mathit{B}}_0$ is directed in the longitudinal z-direction (${\mathit{B}}_0=\pm B_0\left(\mathit{r}\right){\mathbf{i}}_z$, $B_0\left(\mathit{r}\right)>0$) and we have available a transverse RF field with x- and y-components only. The corresponding phasor of this RF field is given by

$${\widehat{\mathit{B}}}_1(\mathit{r},\mathrm{j}\omega )={\widehat{B}}_{1;x}(\mathit{r},\mathrm{j}\omega ){\mathbf{i}}_x+{\widehat{B}}_{1;y}(\mathit{r},\mathrm{j}\omega ){\mathbf{i}}_y,$$

(5)

which can be written as

$${\widehat{\mathit{B}}}_1(\mathit{r},\mathrm{j}\omega )={\widehat{\mathit{B}}}_1^{+}(\mathit{r},\mathrm{j}\omega )+{\widehat{\mathit{B}}}_1^{-}(\mathit{r},\mathrm{j}\omega )$$

(6)

with

$${\widehat{\mathit{B}}}_1^{+}(\mathit{r},\mathrm{j}\omega )={\widehat{B}}_1^{+}(\mathit{r},\mathrm{j}\omega )({\mathbf{i}}_x-\mathrm{j}{\mathbf{i}}_y)\mathrm{and}{\widehat{\mathit{B}}}_1^{-}(\mathit{r},\mathrm{j}\omega )={\left[{\widehat{B}}_1^{-}(\mathit{r},\mathrm{j}\omega )\right]}^{\ast }({\mathbf{i}}_x+\mathrm{j}{\mathbf{i}}_y),$$

(7)

where we have introduced the ${\widehat{B}}_1^{+}$ and ${\widehat{B}}_1^{-}$ fields defined as

$${\widehat{B}}_1^{+}(\mathit{r},\mathrm{j}\omega )=\frac{{\widehat{B}}_{1;x}(\mathit{r},\mathrm{j}\omega )+\mathrm{j}{\widehat{B}}_{1;y}(\mathit{r},\mathrm{j}\omega )}{2}$$

(8)

and

$${\widehat{B}}_1^{-}(\mathit{r},\mathrm{j}\omega )={\left[\frac{{\widehat{B}}_{1;x}(\mathit{r},\mathrm{j}\omega )-\mathrm{j}{\widehat{B}}_{1;y}(\mathit{r},\mathrm{j}\omega )}{2}\right]}^{\ast },$$

(9)

respectively. Substitution of Equation (6) into Equation (4) leads to the time-domain RF field decomposition

$${\mathit{B}}_1(\mathit{r},t)={\mathit{B}}_1^{+}(\mathit{r},t)+{\mathit{B}}_1^{-}(\mathit{r},t)$$

(10)

with

$${\mathit{B}}_1^{+}(\mathit{r},t)=Re\left[{\widehat{\mathit{B}}}_1^{+}(\mathit{r},\mathrm{j}\omega )exp\left(\mathrm{j}\omega t\right)\right]\mathrm{and}{\mathit{B}}_1^{-}(\mathit{r},t)=Re\left[{\widehat{\mathit{B}}}_1^{-}(\mathit{r},\mathrm{j}\omega )exp\left(\mathrm{j}\omega t\right)\right].$$

(11)

Finally, we decompose the scalar ${\widehat{B}}_1^{+}$ and ${\widehat{B}}_1^{-}$ fields into their real and imaginary parts as

$${\widehat{B}}_1^{\pm }=Re\left({\widehat{B}}_1^{\pm }\right)+\mathrm{j}Im\left({\widehat{B}}_1^{\pm }\right).$$

Using these decompositions in Equation (7) and substituting the results in the field expressions for ${\mathit{B}}_1^{+}(\mathit{r},t)$ and ${\mathit{B}}_1^{-}(\mathit{r},t)$ as given by Equation (11) gives

$${\mathit{B}}_1^{+}(\mathit{r},t)=Re\left({\widehat{B}}_1^{+}\right)\left[cos\left(\omega t\right){\mathbf{i}}_x+sin\left(\omega t\right){\mathbf{i}}_y\right]+Im\left({\widehat{B}}_1^{+}\right)\left[-sin\left(\omega t\right){\mathbf{i}}_x+cos\left(\omega t\right){\mathbf{i}}_y\right]$$

(12)

and

$${\mathit{B}}_1^{-}(\mathit{r},t)=Re\left({\widehat{B}}_1^{-}\right)\left[cos\left(\omega t\right){\mathbf{i}}_x-sin\left(\omega t\right){\mathbf{i}}_y\right]+Im\left({\widehat{B}}_1^{-}\right)\left[sin\left(\omega t\right){\mathbf{i}}_x+cos\left(\omega t\right){\mathbf{i}}_y\right].$$

(13)

From the above expressions, we observe that the ${\mathit{B}}_1^{+}$ vector traces out a circle in the transverse $xy$-plane and the radius of this circle is given by $|{\mathit{B}}_1^{+}|=|{\widehat{B}}_1^{+}|$. The ${\mathit{B}}_1^{-}$ vector also traces out a circle but this circle has a radius $|{\mathit{B}}_1^{-}|=|{\widehat{B}}_1^{-}|$ and is traversed in the opposite direction compared with the direction in which the circle of the ${\mathit{B}}_1^{+}$ field is traversed. Since in both cases the fields trace out a circle in the transverse plane, the ${\mathit{B}}_1^{+}$ and ${\mathit{B}}_1^{-}$ fields are called circularly polarized.

The direction of rotation depends upon the direction of the static background field. In particular, assume that our reference frame is such that the static ${\mathit{B}}_0$ field is in the negative z-direction: ${\mathit{B}}_0=-B_0\left(\mathit{r}\right){\mathbf{i}}_z$ with $B_0\left(\mathit{r}\right)>0$. From Equations (12) and (13), we observe that in this case the ${\mathit{B}}_1^{+}$ and ${\mathit{B}}_1^{-}$ fields rotate, respectively, in a left- and right-handed manner about the ${\mathit{B}}_0$ field. When the background field is directed in the positive ${\mathbf{i}}_z$-direction, the situation is reversed and the circularly polarized fields ${\mathit{B}}_1^{+}$ and ${\mathit{B}}_1^{-}$ rotate, respectively, in a right- and left-handed manner about the ${\mathit{B}}_0$ field.

To summarize, any transverse RF field can be decomposed into two circularly polarized fields, where one is polarized in a left-handed manner with respect to background field, while the other is polarized in a right-handed manner with respect to the background field. Explicitly, we have

$${\mathit{B}}_1(\mathit{r},t)={\mathit{B}}_1^{\mathrm{lh}}(\mathit{r},t)+{\mathit{B}}_1^{\mathrm{rh}}(\mathit{r},t),$$

(14)

where ${\mathit{B}}_1^{\mathrm{lh}}(\mathit{r},t)$ and ${\mathit{B}}_1^{\mathrm{rh}}(\mathit{r},t)$ rotate in a left- and right-handed manner about the ${\mathit{B}}_0$ field, respectively. In the case that this background field is in the negative ${\mathbf{i}}_z$-direction, we have

$${\mathit{B}}_1^{\mathrm{lh}}(\mathit{r},t)={\mathit{B}}_1^{+}(\mathit{r},t)\mathrm{and}{\mathit{B}}_1^{\mathrm{rh}}(\mathit{r},t)={\mathit{B}}_1^{-}(\mathit{r},t),$$

while if the background field is in the positive ${\mathbf{i}}_z$-direction we have

$${\mathit{B}}_1^{\mathrm{lh}}(\mathit{r},t)={\mathit{B}}_1^{-}(\mathit{r},t)\mathrm{and}{\mathit{B}}_1^{\mathrm{rh}}(\mathit{r},t)={\mathit{B}}_1^{+}(\mathit{r},t).$$

2.1. Transmit and Receive Fields

As is well known, a circularly polarized RF field that operates at the Larmor frequency and that also rotates in a left-handed manner about the ${\mathit{B}}_0$ field influences the orientation of the magnetization, which ultimately leads to measurable MR signals. During transmission then, the left-handed circularly polarized part of the RF field, ${\mathit{B}}_1^{\mathrm{lh}}$, is of interest. In the case that the background field is in the negative ${\mathbf{i}}_z$-direction, it is the scalar ${\widehat{B}}_1^{+}$ field that determines ${\mathit{B}}_1^{\mathrm{lh}}$, while, if the background field is in the positive ${\mathbf{i}}_z$-direction, it is the scalar field ${\widehat{B}}_1^{-}$ that determines ${\mathit{B}}_1^{\mathrm{lh}}$.

Now, in the MRI literature, the left-handed circularly polarized RF field is always described in terms of a scalar ${\widehat{B}}_1^{+}$ field, which seems to contradict the above observation that this field is described by ${\widehat{B}}_1^{-}$ in the case that the static background field is in the positive ${\mathbf{i}}_z$-direction. It is important, however, to realize that the scalar ${\widehat{B}}_1^{+}$ and ${\widehat{B}}_1^{-}$ fields are defined in terms of phasors that correspond to a particular time factor that is used to represent the RF field. Moreover, the phasors of the same RF field that correspond to the two different time factors are the complex-conjugate of each other (Equation (3)). Consequently, we have

$$\begin{array}{ccc}\hfill {\widehat{B}}_1^{+}(\mathit{r},\mathrm{j}\omega )& =\frac{{\widehat{B}}_{1;x}(\mathit{r},\mathrm{j}\omega )+\mathrm{j}{\widehat{B}}_{1;y}(\mathit{r},\mathrm{j}\omega )}{2}\hfill & \\ & =\frac{{\widehat{B}}_{1;x}(\mathit{r},\mathrm{i}\omega )+\mathrm{i}{\widehat{B}}_{1;y}(\mathit{r},\mathrm{i}\omega )}{2}\hfill & \hfill (\mathrm{writing}\mathrm{i}\mathrm{instead}\mathrm{of}\mathrm{j})\\ & =\frac{{\widehat{B}}_{1;x}^{\ast }(\mathit{r},-\mathrm{i}\omega )+\mathrm{i}{\widehat{B}}_{1;y}^{\ast }(\mathit{r},-\mathrm{i}\omega )}{2}\hfill & \hfill (\mathrm{using}\mathrm{Equation}\left(3\right))\\ & ={\left[\frac{{\widehat{B}}_{1;x}(\mathit{r},-\mathrm{i}\omega )-\mathrm{i}{\widehat{B}}_{1;y}(\mathit{r},-\mathrm{i}\omega )}{2}\right]}^{\ast }={\widehat{B}}_1^{-}(\mathit{r},-\mathrm{i}\omega )\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill {\widehat{B}}_1^{-}(\mathit{r},\mathrm{j}\omega )& ={\left[\frac{{\widehat{B}}_{1;x}(\mathit{r},\mathrm{j}\omega )-\mathrm{j}{\widehat{B}}_{1;y}(\mathit{r},\mathrm{j}\omega )}{2}\right]}^{\ast }\hfill & \\ & =\frac{{\widehat{B}}_{1;x}^{\ast }(\mathit{r},\mathrm{j}\omega )+\mathrm{j}{\widehat{B}}_{1;y}^{\ast }(\mathit{r},\mathrm{j}\omega )}{2}\hfill & \\ & =\frac{{\widehat{B}}_{1;x}(\mathit{r},-\mathrm{j}\omega )+\mathrm{j}{\widehat{B}}_{1;y}^{\ast }(\mathit{r},-\mathrm{j}\omega )}{2}\hfill & \hfill (\mathrm{using}\mathrm{Equation}\left(3\right))\\ & =\frac{{\widehat{B}}_{1;x}(\mathit{r},-\mathrm{i}\omega )+\mathrm{i}{\widehat{B}}_{1;y}(\mathit{r},\mathrm{i}\omega )}{2}={\widehat{B}}_1^{+}(\mathit{r},-\mathrm{i}\omega )\hfill & \hfill (\mathrm{writing}\mathrm{i}\mathrm{instead}\mathrm{of}\mathrm{j}).\end{array}$$

In other words, the ${\widehat{B}}_1^{+}$ field always describes a left-handed circularly polarized field provided that the phasors of Equation (4) are used for ${\mathit{B}}_0$ defined in the negative z-direction, while the phasors of Equation (1) have to be used for ${\mathit{B}}_0$ defined in the positive z-direction. Since transmitting a left-handed circularly polarized field operating at the Larmor frequency enables us to manipulate the magnetization, the ${\widehat{B}}_1^{+}$ field is often referred to as the transmit field. Similarly, received signals can be expressed in terms of the right-handed circularly polarized field ${\mathit{B}}_1^{\mathrm{rh}}$, which is completely described by the ${\widehat{B}}_1^{-}$ field if the phasors of Equation (4) are used in its definition for ${\mathit{B}}_0$ defined in the negative ${\mathbf{i}}_z$-direction and the phasors of Equation (1) are used if the background field is in the positive ${\mathbf{i}}_z$-direction. For this reason, the ${\widehat{B}}_1^{-}$ field is often referred to as the receive field.

2.2. MR Imaging

The transmit field can be written in polar form as ${\widehat{B}}_1^{+}=\left|{\widehat{B}}_1^{+}\right|exp\left(\mathrm{j}{\widehat{\phi }}^{+}\right)$, where $\left|{\widehat{B}}_1^{+}\right|$ is the amplitude or magnitude of the transmit field and ${\widehat{\phi }}^{+}\in (-\pi ,\pi ]$ its phase. Similarly, the receive field can be written in polar form as ${\widehat{B}}_1^{-}=\left|{\widehat{B}}_1^{-}\right|exp\left(\mathrm{j}{\widehat{\phi }}^{-}\right)$, with $\left|{\widehat{B}}_1^{-}\right|$ its amplitude and ${\widehat{\phi }}^{-}\in (-\pi ,\pi ]$ its phase. Note that, to define a phase that is unique, we have restricted the transmit phase ${\widehat{\phi }}^{+}$ and the receive phase ${\widehat{\phi }}^{-}$ to the principle branch $(-\pi ,\pi ]$. Spatial information is encoded into the signal using magnetic field gradients, applied after the ${\widehat{B}}_1^{+}$ field has tipped the magnetization into the transverse plane. Due to the interaction with the body, the transmit field has a spatial dependence, denoted ${\widehat{B}}_1^{+}\left(\mathit{r}\right)$. The polar decomposition is used to express the acquired spatially dependent MR image as [13,19,24,25]

$$I\left(\mathit{r}\right)={\varrho }_0\left(\mathit{r}\right)sin\left(\gamma \tau \left|{\widehat{B}}_1^{+}\left(\mathit{r}\right)\right|\right)exp\left[\mathrm{j}{\widehat{\phi }}^{+}\left(\mathit{r}\right)\right]{\widehat{B}}_1^{-;\ast }\left(\mathit{r}\right),$$

(15)

with ${\varrho }_0$ the proton density, $\gamma$ the gyromagnetic ratio and $\tau$ the RF pulse duration and where ${\widehat{B}}_1^{-;\ast }$ is the complex conjugate of ${\widehat{B}}_1^{-}$. In this simplified expression for the acquired MR image, system dependent factors and contrast terms that underlie an MR image, such as $T_1$ and $T_2$ relaxation, are ignored. Of the transmit and receive fields, only the magnitude of the transmit field shows a non-linear impact on the MR image. This non-linear relation allows for the direct measurement of the transmit magnitude by combining images from different scans such that confounding factors cancel. However, the acquired phase is always the superposition of the phases of ${\widehat{B}}_1^{+}$ and ${\widehat{B}}_1^{-;\ast }$, called the transceive phase, which can not be disentangled from measurements and are therefore difficult to determine exactly. It has been observed that at 1.5 and 3 T the transmit phase closely resembles the phase of ${\widehat{B}}_1^{-;\ast }$ (see also the example given in Figure 1), and in those cases the transmit phase is therefore typically estimated as half the transceive phase: this is termed the transceive phase assumption [23,26]. Similarly, the (magnitude of the) receive field is weighted by the proton density, which is also difficult to disentangle. If the proton-density is not negligible, the proton-density or magnitude of the receive field can be extracted from their product term based on symmetry patterns of the transmit and receive fields in the case of a symmetrical object and imaging setup [27,28]. Additionally, the proton-density could be removed via suitable modeling based on image segmentation [27]. However, knowledge of the transmit phase, receive phase or receive magnitude individually is not always necessary, but could also potentially be determined through EPT.

2.3. Transmit and Receive Fields in Terms of Measurable Quantities

Consider a multi-element RF antenna with P transmit and Q receive channels. The transmit field from channel p (${\widehat{B}}_{1p}^{+}$) measured at receive channel q can then be written in measurable (known) and unknown terms as

$$\begin{array}{cc}\hfill {\widehat{B}}_{1p}^{+}& =\left|{\widehat{B}}_{1p}^{+}\right|exp\left[\mathrm{j}\left({\widehat{\phi }}_p^{+}-{\widehat{\phi }}_q^{-}\right)\right]exp\left(\mathrm{j}{\widehat{\phi }}_q^{-}\right)\hfill \\ & =\left|{\widehat{B}}_{1p}^{+}\right|exp\left(\mathrm{j}{\widehat{\phi }}_{pq}^{\pm }\right)exp\left(\mathrm{j}{\widehat{\phi }}_q^{-}\right)\hfill \\ & ={\widehat{B}}_{1p}^{+;\mathrm{TRX}(p,q)}exp\left(\mathrm{j}{\widehat{\phi }}_q^{-}\right)\hfill \end{array}$$

(16)

with ${\widehat{\phi }}_{pq}^{\pm }={\widehat{\phi }}_p^{+}-{\widehat{\phi }}_q^{-}$ the transceive phase (note that ${\widehat{\phi }}^{-}$ is sometimes defined as the argument of ${\widehat{B}}_1^{-;\ast }$, such that the transceive phase is given by ${\widehat{\phi }}^{\pm }={\widehat{\phi }}^{+}+{\widehat{\phi }}^{-}$), ${\widehat{\phi }}_p^{+}$ the absolute transmit phase of transmit channel p and ${\widehat{\phi }}_q^{-}$ the absolute receive phase of receive channel q. The term ${\widehat{B}}_{1p}^{+;\mathrm{TRX}(p,q)}=\left|{\widehat{B}}_{1p}^{+}\right|exp\left(\mathrm{j}{\widehat{\phi }}_{pq}^{\pm }\right)$ is the measurable term, while $exp\left(\mathrm{j}{\widehat{\phi }}_q^{-}\right)$ is the unknown term. Note that this formulation is applicable for RF coils in general with $P=Q=1$, while the subsequent two formulations in terms of relative transmit phases are only applicable for multi-element RF arrays with $P>1$ and/or $Q>1$. The transmit field from channel p can be written in terms of relative phase distributions as

$$\begin{array}{cc}\hfill {\widehat{B}}_{1p}^{+}& =\left|{\widehat{B}}_{1p}^{+}\right|exp\left[\mathrm{j}\left({\widehat{\phi }}_p^{+}-{\widehat{\phi }}_r^{+}\right)\right]exp\left(\mathrm{j}{\widehat{\phi }}_r^{+}\right)\hfill \\ & =\left|{\widehat{B}}_{1p}^{+}\right|exp\left(\mathrm{j}{\widehat{\phi }}_{pr}^{+}\right)exp\left(\mathrm{j}{\widehat{\phi }}_r^{+}\right)\hfill \\ & ={\widehat{B}}_{1p}^{+;\mathrm{rel}(p,r)}exp\left(\mathrm{j}{\widehat{\phi }}_r^{+}\right)\hfill \end{array}$$

(17)

with ${\widehat{\phi }}_{pr}^{+}={\widehat{\phi }}_p^{+}-{\widehat{\phi }}_r^{+}$ the transmit phase of channel p relative to the reference transmit phase ${\widehat{\phi }}_r^{+}$ of channel r. ${\widehat{B}}_{1p}^{+;\mathrm{rel}(p,r)}=\left|{\widehat{B}}_{1p}^{+}\right|exp\left(\mathrm{j}{\widehat{\phi }}_{pr}^{+}\right)$ is the measurable term, while $exp\left(\mathrm{j}{\widehat{\phi }}_r^{+}\right)$ is the unknown term. Additionally, the receive phase can also be written in a similar formulation. The measurable term is, however, weighted by the proton density. For the conjugate of the receive field, we have

$$\begin{array}{cc}\hfill {\widehat{B}}_{1q}^{-;\ast }& =\left|{\varrho }_0{\widehat{B}}_{1q}^{-}\right|exp\left[\mathrm{j}\left(-{\widehat{\phi }}_q^{-}+{\widehat{\phi }}_r^{-}\right)\right]\left|{\varrho }_0^{-1}\right|exp\left(-\mathrm{j}{\widehat{\phi }}_r^{-}\right)\hfill \\ & =\left|{\varrho }_0{\widehat{B}}_{1q}^{-}\right|exp\left(\mathrm{j}{\widehat{\phi }}_{rq}^{-}\right)\left|{\varrho }_0^{-1}\right|exp\left(-\mathrm{j}{\widehat{\phi }}_r^{-}\right)\hfill \end{array}$$

(18)

with ${\widehat{\phi }}_{rq}^{-}={\widehat{\phi }}_r^{-}-{\widehat{\phi }}_q^{-}$ the receive phase of reference channel r relative to channel q. $\left|{\varrho }_0{\widehat{B}}_{1q}^{-}\right|exp\left(\mathrm{j}{\widehat{\phi }}_{rq}^{-}\right)$ is the measurable term, while $\left|{\varrho }_0^{-1}\right|exp\left(-\mathrm{j}{\widehat{\phi }}_r^{-}\right)$ is the unknown term.

In summary, we use the phasor representation of Equation (4) to describe the RF field and orient the ${\mathit{B}}_0$ field in the negative ${\mathbf{i}}_z$-direction such that ${\widehat{B}}_1^{+}$ enables the manipulation of magnetization. Furthermore, we describe the transmit and receive fields with Equations (8) and (9), such that ${\widehat{B}}_1^{+}$ and ${\widehat{B}}_1^{-}$ described by the phasor representation of Equation (4) correspond to the ${\widehat{B}}_1^{-}$ and ${\widehat{B}}_1^{+}$ fields described by the phasor representation of Equation (1), respectively. Additionally, the MR image can be described by Equation (15), which shows that the transmit and receive field phases, as well as the proton density and the receive field magnitude, are entangled and therefore not directly available from MR acquisitions. Instead of making assumptions about the acquirable data to obtain absolute transmit or receive field maps, the transmit and receive fields can also be expressed in terms of known (directly derived from measurements) and unknown terms, as depicted in Equations (16)–(18).

3. Fundamental EPT Equations

Physical model-based EPT approaches all rely on a few fundamental equations from which their central equations are derived. To derive and understand the approaches, knowledge about the Maxwell’s equations, the Helmholtz equation and the scattering field formalism is required. These are summarized below.

3.1. First-Order Differential Equations: Maxwell’s Equations

Maxwell’s equations for time-harmonic fields are given by

$$\begin{array}{cc}\hfill -\mathbf{\nabla }\times \widehat{\mathit{H}}\left(\mathit{r}\right)+\eta \left(\mathit{r}\right)\widehat{\mathit{E}}\left(\mathit{r}\right)& =-{\widehat{\mathit{J}}}^{\mathrm{ext}}\left(\mathit{r}\right),\mathrm{and}\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \mathbf{\nabla }\times \widehat{\mathit{E}}\left(\mathit{r}\right)+\zeta \left(\mathit{r}\right)\widehat{\mathit{H}}\left(\mathit{r}\right)& =\mathbf{0},\hfill \end{array}$$

(20)

with $\eta \left(\mathit{r}\right)=\sigma \left(\mathit{r}\right)+\mathrm{j}\omega \epsilon \left(\mathit{r}\right)$ and $\zeta \left(\mathit{r}\right)=\mathrm{j}\omega \mu \left(\mathit{r}\right)$, which are, respectively, the per-unit-length admittance and impedance of the medium. Here, $\sigma$, $\epsilon$, $\mu$ and $\omega$ are the conductivity, permittivity, permeability and angular (RF) frequency, respectively. Additionally, in the MR setting, ${\widehat{\mathit{J}}}^{\mathrm{ext}}$ is an external current density distribution present on the MR coil that generates the EM fields. Since these sources are located outside the body and since the permeability of biological tissue is assumed to be constant and equal to that of vacuum, the RF field inside the body satisfies the Maxwell equations

$$\begin{array}{cc}\hfill -\mathbf{\nabla }\times \widehat{\mathit{B}}\left(\mathit{r}\right)+{\mu }_0\eta \left(\mathit{r}\right)\widehat{\mathit{E}}\left(\mathit{r}\right)& =\mathbf{0},\mathrm{and}\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill \mathbf{\nabla }\times \widehat{\mathit{E}}\left(\mathit{r}\right)+\mathrm{j}\omega \widehat{\mathit{B}}\left(\mathit{r}\right)& =\mathbf{0},\hfill \end{array}$$

(22)

with $\widehat{\mathit{B}}={\mu }_0\widehat{\mathit{H}}$. Furthermore, introducing the vectors

$${\mathbf{i}}^{+}=\frac{1}{2}({\mathbf{i}}_x+\mathrm{j}{\mathbf{i}}_y)\mathrm{and}{\mathbf{i}}^{-}=\frac{1}{2}({\mathbf{i}}_x-\mathrm{j}{\mathbf{i}}_y)$$

we have for the transmit and receive fields the expressions ${\widehat{B}}_1^{+}={\mathbf{i}}^{+}·\widehat{\mathit{B}}$ and ${\widehat{B}}_1^{-;\ast }={\mathbf{i}}^{-}·\widehat{\mathit{B}}$. Similarly, we define ${\widehat{E}}_1^{+}={\mathbf{i}}^{+}·\widehat{\mathit{E}}$ and ${\widehat{E}}_1^{-;\ast }={\mathbf{i}}^{-}·\widehat{\mathit{E}}$ and introduce the differentiation operators (Wirtinger derivatives)

$${\partial }^{+}={\mathbf{i}}^{+}·\mathbf{\nabla }=\frac{1}{2}({\partial }_x+\mathrm{j}{\partial }_y)\mathrm{and}{\partial }^{-}={\mathbf{i}}^{-}·\mathbf{\nabla }=\frac{1}{2}({\partial }_x-\mathrm{j}{\partial }_y).$$

(23)

Taking the inner product of ${\mathbf{i}}^{+}$ and the second Maxwell equation now gives an explicit expression for the transmit field, while taking the inner product of ${\mathbf{i}}^{-}$, and this second Maxwell equation gives an explicit expression for the receive field. Explicitly, we have

$${\widehat{B}}_1^{+}=\frac{1}{\omega }\left({\partial }^{+}{E}_z-{\partial }_z{\widehat{E}}_1^{+}\right)\mathrm{and}{\widehat{B}}_1^{-;\ast }=-\frac{1}{\omega }\left({\partial }^{-}{E}_z-{\partial }_z{\widehat{E}}_1^{-;\ast }\right).$$

(24)

These relations tell us that the ${\widehat{B}}_1^{+}$ and ${\widehat{B}}_1^{-;\ast }$ fields result from a difference between transverse variations of the longitudinal electric field (as determined by the Wirtinger derivatives) and longitudinal variations of the transverse ${\widehat{E}}_1^{+}$ and ${\widehat{E}}_1^{-;\ast }$ fields. These equations are used as a starting point in the EPT method discussed in Section 4.7. For completeness, we mention that, if a similar procedure is followed for the first Maxwell equation, we obtain

$${\widehat{E}}_1^{+}=\frac{1}{\mathrm{j}{\mu }_0\eta }\left({\partial }^{+}{B}_z-{\partial }_z{\widehat{B}}_1^{+}\right)\mathrm{and}{\widehat{E}}_1^{-;\ast }=-\frac{1}{\mathrm{j}{\mu }_0\eta }\left({\partial }^{-}{B}_z-{\partial }_z{\widehat{B}}_1^{-;\ast }\right).$$

(25)

3.2. Second-Order Differential Equation: The Generalized Helmholtz Equation

Since the objective is to obtain the dielectric tissue parameters from magnetic field data, a second option is not to consider the electric field at all and to eliminate this field from the source-free first-order Maxwell system as given by Equations (21) and (22). To this end, we take the curl of Equation (21) and obtain

$$-\mathbf{\nabla }\times \mathbf{\nabla }\times \widehat{\mathit{B}}+{\mu }_0\mathbf{\nabla }\times \left(\eta \widehat{\mathit{E}}\right)=\mathbf{0}.$$

(26)

Since

$$\begin{array}{c}\hfill \mathbf{\nabla }\times \mathbf{\nabla }\times \widehat{\mathit{B}}=\mathbf{\nabla }\mathbf{\nabla }·\widehat{\mathit{B}}-{\nabla }^2\widehat{\mathit{B}}\end{array}$$

and

$$\begin{array}{c}\hfill \mathbf{\nabla }\times \left(\eta \widehat{\mathit{E}}\right)=\mathbf{\nabla }\eta \times \widehat{\mathit{E}}+\eta \mathbf{\nabla }\times \widehat{\mathit{E}}=\frac{1}{{\mu }_0}\frac{\mathbf{\nabla }\eta }{\eta }\times (\mathbf{\nabla }\times \widehat{\mathit{B}})-\mathrm{j}\omega \eta \widehat{\mathit{B}},\end{array}$$

Equation (26) can be written as

$$-\mathbf{\nabla }\mathbf{\nabla }·\widehat{\mathit{B}}+{\nabla }^2\widehat{\mathit{B}}+\frac{\mathbf{\nabla }\eta }{\eta }\times (\mathbf{\nabla }\times \widehat{\mathit{B}})-\eta \zeta \widehat{\mathit{B}}=\mathbf{0}.$$

Finally, taking the divergence of the second Maxwell equation (Equation (22)), we get $\mathbf{\nabla }·\widehat{\mathit{B}}=0$ and substituting this result in the above equation, we obtain the generalized Helmholtz equation

$${\nabla }^2\widehat{\mathit{B}}+\frac{\mathbf{\nabla }\eta }{\eta }\times \left(\mathbf{\nabla }\times \widehat{\mathit{B}}\right)+k^2\widehat{\mathit{B}}=\mathbf{0},$$

(27)

where $k={(-\eta \zeta )}^{1/2}={({\omega }^2{\mu }_0\epsilon -\mathrm{j}\omega {\mu }_0\sigma )}^{1/2}$ is the complex wave number with $\mathrm{Im}\left(k\right)\le 0$. Note that, for homogeneous media, $\eta$ is constant and the second term on the left-hand side vanishes. In this case, we have the Helmholtz equation

$${\nabla }^2\widehat{\mathit{B}}+k^2\widehat{\mathit{B}}=\mathbf{0}.$$

(28)

Taking the inner product of the vector ${\mathbf{i}}^{+}$ and Equation (28) gives the Helmholtz equation for the ${\widehat{B}}_1^{+}$ field

$${\nabla }^2{\widehat{B}}_1^{+}+k^2{\widehat{B}}_1^{+}=0,$$

(29)

which serves as a starting point for the EPT methods discussed in Section 4.1, Section 4.2, Section 4.3, Section 5.1 and Section 5.2. It is important to realize that the above Helmholtz equation is valid for homogeneous media ($\eta$ is constant) only. For general inhomogeneous media ($\eta$ is not constant), we have the generalized Helmholtz equation (Equation (27)). Dotting this equation with the vector ${\mathbf{i}}^{+}$, we end up with the generalized Helmholtz equation for the ${\widehat{B}}_1^{+}$ field given by

$${\nabla }^2{\widehat{B}}_1^{+}+{\mathbf{i}}^{+}·\left[\frac{\mathbf{\nabla }\eta }{\eta }\times \left(\mathbf{\nabla }\times \widehat{\mathit{B}}\right)\right]+k^2{\widehat{B}}_1^{+}=0.$$

(30)

This equation serves as a starting point for the EPT methods discussed in Section 4.4, Section 4.5 and Section 4.6, but with the second term on the left-hand side rewritten in terms of ${\widehat{B}}_1^{+}$ and ${\widehat{B}}_z$. Specifically, in the EPT methods of Section 4.4 and Section 4.5, the generalized Helmholtz equation is rewritten in terms of the gradient of ${\widehat{B}}_1^{+}$ and ${\widehat{B}}_z$, while, in the methods of Section 4.6, the generalized Helmholtz equation is written as a convection–reaction equation.

3.2.1. The Gradient-Type Generalized Helmholtz Equation

Let us first consider rewriting the generalized Helmholtz equation (Equation (30)) in terms of the gradient of the ${\widehat{B}}_1^{+}$ field and ${\widehat{B}}_z$, the z-component of the magnetic field. As a first step, we introduce the vector $\mathit{g}={\eta }^{-1}\mathbf{\nabla }\eta$ and write $g^{+}={\mathbf{i}}^{+}·\mathit{g}$. The second term on the left-hand side in Equation (30) can now be written as

$$\begin{array}{cc}\hfill {\mathbf{i}}^{+}·\left[\frac{\mathbf{\nabla }\eta }{\eta }\times \left(\mathbf{\nabla }\times \widehat{\mathit{B}}\right)\right]& ={\mathbf{i}}^{+}·\left[\mathit{g}\times \left(\mathbf{\nabla }\times \widehat{\mathit{B}}\right)\right]=\mathit{g}·\left({\partial }^{+}\widehat{\mathit{B}}-\mathbf{\nabla }{\widehat{B}}_1^{+}\right)\hfill \\ \hfill & =g_x({\partial }^{+}{\widehat{B}}_x-{\partial }_x{\widehat{B}}_1^{+})+g_y({\partial }^{+}{\widehat{B}}_y-{\partial }_y{\widehat{B}}_1^{+})+g_z({\partial }^{+}{\widehat{B}}_z-{\partial }_z{\widehat{B}}_1^{+}).\hfill \end{array}$$

Since

$${\partial }^{+}{\widehat{B}}_x-{\partial }_x{\widehat{B}}_1^{+}=-\mathrm{j}({\partial }^{+}{\widehat{B}}_y-{\partial }_y{\widehat{B}}_1^{+})=\frac{1}{2}\mathrm{j}({\partial }_y{\widehat{B}}_x-{\partial }_x{\widehat{B}}_y),$$

this can be written as

$${\mathbf{i}}^{+}·\left[\frac{\mathbf{\nabla }\eta }{\eta }\times \left(\mathbf{\nabla }\times \widehat{\mathit{B}}\right)\right]=(\mathrm{j}{g}_x-g_y)\frac{1}{2}({\partial }_y{\widehat{B}}_x-{\partial }_x{\widehat{B}}_y)+g_z({\partial }^{+}{\widehat{B}}_z-{\partial }_z{\widehat{B}}_1^{+}).$$

(31)

Furthermore, using $\mathbf{\nabla }·\widehat{\mathit{B}}=0$, we have

$$\frac{1}{2}({\partial }_y{\widehat{B}}_x-{\partial }_x{\widehat{B}}_y)=\mathrm{j}{\partial }_x{\widehat{B}}_1^{+}+{\partial }_y{\widehat{B}}_1^{+}+\frac{1}{2}\mathrm{j}{\partial }_z{\widehat{B}}_z$$

(32)

and substituting this relation in Equation (31) leads to

$${\mathbf{i}}^{+}·\left[\frac{\mathbf{\nabla }\eta }{\eta }\times \left(\mathbf{\nabla }\times \widehat{\mathit{B}}\right)\right]=-{\mathit{f}}^{+}·\mathbf{\nabla }{\widehat{B}}_1^{+}-{\mathit{h}}^{+}·\mathbf{\nabla }{\widehat{B}}_z$$

with

$${\mathit{f}}^{+}=4g^{+}{\mathbf{i}}^{-}+g_z{\mathbf{i}}_z\mathrm{and}{\mathit{h}}^{+}=-g_z{\mathbf{i}}^{+}+g^{+}{\mathbf{i}}_z.$$

With this result, we end up with the gradient-type generalized Helmholtz equation

$${\nabla }^2{\widehat{B}}_1^{+}-{\mathit{f}}^{+}·\mathbf{\nabla }{\widehat{B}}_1^{+}-{\mathit{h}}^{+}·\mathbf{\nabla }{\widehat{B}}_z+k^2{\widehat{B}}_1^{+}=0.$$

(33)

3.2.2. The Generalized Helmholtz Equation as a Convection–Reaction Equation

To arrive at the convection–reaction form of the generalized Helmholtz equation as used in EPT, we return to Equation (32) and rewrite this equation as

$$\frac{1}{2}({\partial }_y{\widehat{B}}_x-{\partial }_x{\widehat{B}}_y)=\mathrm{j}{\partial }_x{\widehat{B}}_1^{+}+{\partial }_y{\widehat{B}}_1^{+}+\frac{1}{2}\mathrm{j}{\partial }_z{\widehat{B}}_z=\mathrm{j}\left(2{\partial }^{-}{\widehat{B}}_1^{+}+\frac{1}{2}{\partial }_z{\widehat{B}}_z\right).$$

(34)

Substitution of this result in Equation (31) gives

$${\mathbf{i}}^{+}·\left[\frac{\mathbf{\nabla }\eta }{\eta }\times \left(\mathbf{\nabla }\times \widehat{\mathit{B}}\right)\right]=-{\mathbf{\beta }}^{+}·\mathit{g},$$

where

$${\mathbf{\beta }}^{+}=(2{\partial }^{-}{\widehat{B}}_1^{+}+\frac{1}{2}{\partial }_z{\widehat{B}}_z){\mathbf{i}}_x+\mathrm{j}(2{\partial }^{-}{\widehat{B}}_1^{+}+\frac{1}{2}{\partial }_z{\widehat{B}}_z){\mathbf{i}}_y+({\partial }_z{\widehat{B}}_1^{+}-{\partial }^{+}{\widehat{B}}_z){\mathbf{i}}_z.$$

(35)

The generalized Helmholtz equation now becomes

$${\nabla }^2{\widehat{B}}_1^{+}-{\mathbf{\beta }}^{+}·\mathit{g}+k^2{\widehat{B}}_1^{+}=0.$$

Dividing this equation by $\eta$ and using the definition of vector $\mathit{g}$, we arrive at our final form

$$u{\nabla }^2{\widehat{B}}_1^{+}+{\mathbf{\beta }}^{+}·\mathbf{\nabla }u-\zeta {\widehat{B}}_1^{+}=0$$

(36)

with $u={\eta }^{-1}$. Equation (36) is the generalized Helmholtz equation in convection–reaction form, where $u{\nabla }^2{\widehat{B}}_1^{+}-\zeta {\widehat{B}}_1^{+}$ is the reaction component and ${\mathbf{\beta }}^{+}$ is the convective field. Observe that the components of the convective field are directly related to the dielectric medium parameters and the z-component of the electric field strength via (cf. Equations (24) and (25))

$$2{\beta }_x=\mathrm{j}{\mu }_0\eta {\widehat{E}}_z,2{\beta }_y=\mathrm{j}2{\beta }_x=-{\mu }_0\eta {\widehat{E}}_z,\mathrm{and}{\beta }_z=-\mathrm{j}{\mu }_0\eta {\widehat{E}}^{+}.$$

(37)

These relations are used as a starting point in the EPT methods discussed in Section 4.7 and Section 4.8.

3.2.3. Helmholtz Equations for the Receive Field

For completeness, we mention that a similar procedure can be carried out for the ${\widehat{B}}_1^{-}$ field. In particular, taking the inner product of the vector ${\mathbf{i}}^{-}$ and Equation (27), we end up with

$${\nabla }^2{\widehat{B}}_1^{-;\ast }+{\mathbf{i}}^{-}·\left[\frac{\mathbf{\nabla }\eta }{\eta }\times \left(\mathbf{\nabla }\times \widehat{\mathit{B}}\right)\right]+k^2{\widehat{B}}_1^{-;\ast }=0,$$

(38)

which is the generalized Helmholtz equation for ${\widehat{B}}_1^{-;\ast }$. This equation can also be written in terms of gradients of the ${\widehat{B}}_1^{-;\ast }$ field and ${\widehat{B}}_z$ as

$${\nabla }^2{\widehat{B}}_1^{-;\ast }-{\mathit{f}}^{-}·\mathbf{\nabla }{\widehat{B}}_1^{-;\ast }-{\mathit{h}}^{-}·\mathbf{\nabla }{\widehat{B}}_z+k^2{\widehat{B}}_1^{-;\ast }=0,$$

(39)

with

$${\mathit{f}}^{-}=4g^{-}{\mathbf{i}}^{+}+g_z{\mathbf{i}}_z\mathrm{and}{\mathit{h}}^{-}=-g_z{\mathbf{i}}^{-}+g^{-}{\mathbf{i}}_z,$$

where we introduce $g^{-}={\mathbf{i}}^{-}·\mathit{g}$, or as a convection–reaction equation as

$$u{\nabla }^2{\widehat{B}}_1^{-;\ast }+{\mathbf{\beta }}^{-}·\mathbf{\nabla }u-\zeta {\widehat{B}}_1^{-;\ast }=0,$$

(40)

with

$${\mathbf{\beta }}^{-}=(2{\partial }^{+}{\widehat{B}}_1^{-;\ast }+\frac{1}{2}{\partial }_z{\widehat{B}}_z){\mathbf{i}}_x-\mathrm{j}(2{\partial }^{+}{\widehat{B}}_1^{-;\ast }+\frac{1}{2}{\partial }_z{\widehat{B}}_z){\mathbf{i}}_y+({\partial }_z{\widehat{B}}_1^{-;\ast }-{\partial }^{-}{\widehat{B}}_z){\mathbf{i}}_z.$$

For the vectorial Helmholtz equation of Equation (28), dotting with the vector ${\mathbf{i}}^{-}$ gives

$${\nabla }^2{\widehat{B}}_1^{-;\ast }+k^2{\widehat{B}}_1^{-;\ast }=0,$$

(41)

which is the Helmholtz equation for the ${\widehat{B}}_1^{-;\ast }$ field in the case of homogeneous media.

3.3. Volume Integral Equations

The fundamental integral equations are obtained through a scattering formalism by exploiting the linearity of Maxwell’s equations. Specifically, the total electromagnetic field in the presence of an object in an MR coil is denoted by $\left\{\widehat{\mathit{E}},\widehat{\mathit{H}}\right\}$, and this field is written as the sum of an incident and scattered field as

$$\left\{\widehat{\mathit{E}},\widehat{\mathit{H}}\right\}=\left\{{\widehat{\mathit{E}}}^{\mathrm{inc}},{\widehat{\mathit{H}}}^{\mathrm{inc}}\right\}+\left\{{\widehat{\mathit{E}}}^{\mathrm{sca}},{\widehat{\mathit{H}}}^{\mathrm{sca}}\right\},$$

(42)

where the incident field is defined as the field that is present in an empty (air-filled) RF coil. This incident field is generated by an external current density distribution ${\widehat{\mathit{J}}}^{\mathrm{ext}}$ representing the MR coil that occupies the bounded source domain $\mathbb{S}$. The governing equations for the incident field are

$${\widehat{\mathit{E}}}^{\mathrm{inc}}=\left(k_0^2+\mathbf{\nabla }\mathbf{\nabla }·\right){\widehat{\mathit{A}}}^{\mathrm{ext}}\mathrm{and}{\widehat{\mathit{H}}}^{\mathrm{inc}}={\eta }_0\mathbf{\nabla }\times {\widehat{\mathit{A}}}^{\mathrm{ext}},$$

(43)

with $k_0=\omega /c_0$ the wave number of the background medium, $c_0$ the electromagnetic wave speed of free space and ${\eta }_0=\mathrm{j}\omega {\epsilon }_0$ the admittance of the background medium. In the above field expressions, the vector potential ${\widehat{\mathit{A}}}^{\mathrm{ext}}$ is given by

$${\widehat{\mathit{A}}}^{\mathrm{ext}}\left(\mathit{r}\right)={\eta }_0^{-1}{\int }_{{\mathit{r}}^{\prime }\in \mathbb{S}}\widehat{G}\left(\mathit{r}-{\mathit{r}}^{\prime }\right){\widehat{\mathit{J}}}^{\mathrm{ext}}\left({\mathit{r}}^{\prime }\right)\mathrm{d}V,$$

where $\widehat{G}\left(\mathit{r}\right)$ is the Green’s function of the background medium given by

$$\widehat{G}\left(\mathit{r}\right)=\frac{exp\left(-\mathrm{j}{k}_0\left|\mathit{r}\right|\right)}{4\pi \left|\mathit{r}\right|},\mathrm{for}\mathit{r}\ne \mathbf{0}.$$

When there is an object present, scattered fields will be generated due to an induced scattering current density distribution ${\widehat{\mathit{J}}}^{\mathrm{sca}}$ having the object domain $\mathbb{D}$ as its support. The scattered fields are given by

$${\widehat{\mathit{E}}}^{\mathrm{sca}}=\left(k_0^2+\mathbf{\nabla }\mathbf{\nabla }·\right){\widehat{\mathit{A}}}^{\mathrm{sca}}\mathrm{and}{\widehat{\mathit{H}}}^{\mathrm{sca}}={\eta }_0\mathbf{\nabla }\times {\widehat{\mathit{A}}}^{\mathrm{sca}},$$

(44)

where the vector potential ${\widehat{\mathit{A}}}^{\mathrm{sca}}$ is given by

$${\widehat{\mathit{A}}}^{\mathrm{sca}}\left(\mathit{r}\right)={\eta }_0^{-1}{\int }_{{\mathit{r}}^{\prime }\in \mathbb{D}}\widehat{G}\left(\mathit{r}-{\mathit{r}}^{\prime }\right){\widehat{\mathit{J}}}^{\mathrm{sca}}\left({\mathit{r}}^{\prime }\right)\mathrm{d}V$$

with ${\widehat{\mathit{J}}}^{\mathrm{sca}}=\left(\eta -{\eta }_0\right)\widehat{\mathit{E}}$ the scattered current density distribution. Note that the scattering current density and consequently the scattered field vanish if the object is absent ($\eta ={\eta }_0$) and the total electromagnetic field is equal to the incident field. Finally, we mention that the ${\widehat{B}}_1^{+}$ field can be obtained from the vector potential as

$${\widehat{B}}_1^{+}={\widehat{B}}_1^{+;\mathrm{inc}}+{\widehat{B}}_1^{+;\mathrm{sca}}\mathrm{with}{\widehat{B}}_1^{+;\mathrm{inc}}=\frac{\omega }{c_0^2}\tilde{\mathbf{\nabla }}·{\widehat{\mathit{A}}}^{\mathrm{ext}}\mathrm{and}{\widehat{B}}_1^{+;\mathrm{sca}}=\frac{\omega }{c_0^2}\tilde{\mathbf{\nabla }}·{\widehat{\mathit{A}}}^{\mathrm{sca}},$$

(45)

where $\tilde{\mathbf{\nabla }}={\mathbf{i}}_z{\partial }^{+}-{\mathbf{i}}^{+}{\partial }_z$. The dielectric tissue parameters only influence ${\widehat{B}}_1^{+;\mathrm{sca}}$, that is, the effects of the medium parameters on the ${\widehat{B}}_1^{+}$ field have been separated from the excited incident ${\widehat{B}}_1^{+}$ field. These relations are used as starting point for the EPT methods discussed in Section 4.9, Section 4.10 and Section 4.11.

4. EPT Methods Requiring Transmit Field Map**

This section discusses analytical EPT approaches based on transmit field map**s. The section starts with direct local differential methods and roughly transitions to end with forward global integral methods. More specifically, the EPT methods discussed in this section are

• Section 4.1: Helmholtz-based EPT (H-EPT)
• Section 4.2: Simplified H-EPT (SH-EPT)
• Section 4.2.1: Poisson-based Conductivity Map** (P-CM)
• Section 4.3: Local Maxwell tomography (LMT)
• Section 4.4: Modified dual-excitation EPT (MDE-EPT)
• Section 4.5: Gradient-based EPT (G-EPT)
• Section 4.6: Convection–reaction EPT (CR-EPT)
• Section 4.6.1: Phase-only convection–reaction conductivity map** (PCR-CM)
• Section 4.7: Transverse EPT (T-EPT)
• Section 4.8: First-order induced-current EPT (foIC-EPT)
• Section 4.9: Variational Born iterative method EPT (VBIM-EPT)
• Section 4.10: Global Maxwell tomography (GMT)
• Section 4.11: Contrast source inversion EPT (CSI-EPT)

Other methods that do not require transmit field map** are discussed in Section 5. Machine-learning approaches are discussed in Section 6.

4.1. Helmholtz-Based EPT

Helmholtz-based EPT (H-EPT) assumes a homogeneous medium ($\nabla \eta =\mathbf{0}$) and is based on the Helmholtz equation (Equation (29)) [22,23,30]. Explicitly, assuming that the ${\widehat{B}}_1^{+}$ field is known, the tissue parameters are determined from

$$\frac{{\nabla }^2{\widehat{B}}_1^{+}}{{\widehat{B}}_1^{+}}=-k^2$$

(46)

and the definition of the wave number as

$$\sigma =\frac{1}{\omega {\mu }_0}Im\left(\frac{{\nabla }^2{\widehat{B}}_1^{+}}{{\widehat{B}}_1^{+}}\right)\mathrm{and}\epsilon =\frac{-1}{{\omega }^2{\mu }_0}Re\left(\frac{{\nabla }^2{\widehat{B}}_1^{+}}{{\widehat{B}}_1^{+}}\right).$$

(47)

This explicit method is extremely simple, easy to implement and fast to compute. However, the homogeneity assumption results in errors at tissue boundaries; the second-order derivative that acts on the data makes the method sensitive to noise [31,32,33]; and the method requires knowledge of the absolute transmit phase which is not directly available. To mitigate noise effects, filtered Laplacians with increased kernel size can be used, however, this leads to a severe numerical boundary error propagation [32,34]. The second-order differential has been reduced to first-order derivatives in an alternative formulation based on Gauss’ integral theorem, but image segmentation is required to implement this method [35,36]. Since the absolute transmit phase is in practice unavailable, it is typically estimated with the transceive phase assumption. However, since the Laplacian of a variable is the divergence of the gradient of the variable, this assumption can be prevented with multiple acquisitions from a multi-element array. This system namely allows for the determination of the gradient of the transmit phase of a reference transmit channel ($\mathbf{\nabla }{\widehat{\phi }}_r^{+}$) from relative transmit phases [37] or the gradient of the receive phase of the receive channel ($\mathbf{\nabla }{\widehat{\phi }}_q^{-}$) from transceive phase measurements [27].

4.2. Simplified H-EPT

Simplified H-EPT (SH-EPT) derives the conductivity and permittivity independently from the phase and magnitude of the ${\widehat{B}}_1^{+}$, respectively [8,35]. Starting point is again the Helmholtz equation (Equation (29)) for the ${\widehat{B}}_1^{+}$ field. However, here the polar decomposition of ${\widehat{B}}_1^{+}$ is substituted, which gives

$$\frac{{\nabla }^2\left|{\widehat{B}}_1^{+}\right|}{\left|{\widehat{B}}_1^{+}\right|}-{|\mathbf{\nabla }{\widehat{\phi }}^{+}|}^2+\mathrm{j}\left(2\frac{\mathbf{\nabla }\left|{\widehat{B}}_1^{+}\right|}{\left|{\widehat{B}}_1^{+}\right|}·\mathbf{\nabla }{\widehat{\phi }}^{+}+{\nabla }^2{\widehat{\phi }}^{+}\right)=-k^2$$

(48)

and, equating the real and imaginary parts in the above equation, we obtain

$$\sigma =\frac{1}{\omega {\mu }_0}\left(2\frac{\mathbf{\nabla }|{\widehat{B}}_1^{+}|·\mathbf{\nabla }{\widehat{\phi }}^{+}}{\left|{\widehat{B}}_1^{+}\right|}+{\nabla }^2{\widehat{\phi }}^{+}\right)\mathrm{and}\epsilon =\frac{-1}{{\omega }^2{\mu }_0}\left(\frac{{\nabla }^2\left|{\widehat{B}}_1^{+}\right|}{\left|{\widehat{B}}_1^{+}\right|}-{|\mathbf{\nabla }{\widehat{\phi }}^{+}|}^2\right).$$

(49)

Finally, assuming that ${\nabla }^2{\widehat{\phi }}^{+}>>2\frac{\mathbf{\nabla }\left|{\widehat{B}}_1^{+}\right|·\mathbf{\nabla }{\widehat{\phi }}^{+}}{\left|{\widehat{B}}_1^{+}\right|}$, we obtain

$$\sigma =\frac{1}{\omega {\mu }_0}{\nabla }^2{\widehat{\phi }}^{+}.$$

(50)

Note that, if the Helmholtz equation accurately describes the behavior of the ${\widehat{B}}_1^{+}$ field and if the above approximation holds, then only the phase of the ${\widehat{B}}_1^{+}$ field is required to determine the conductivity. We remark that, if we write the ${\widehat{B}}_1^{-}$ field in polar form as well and follow similar steps as for the ${\widehat{B}}_1^{+}$ field, we obtain from the Helmholtz Equation (41)

$$\sigma =-\frac{1}{\omega {\mu }_0}{\nabla }^2{\widehat{\phi }}^{-},$$

where ${\widehat{\phi }}^{-}$ is the phase of the ${\widehat{B}}_1^{-}$ field. Consequently, if the transceive phase ${\widehat{\phi }}^{\pm }={\widehat{\phi }}^{+}-{\widehat{\phi }}^{-}$ is available, we have

$$\sigma =\frac{1}{2\omega {\mu }_0}{\nabla }^2{\widehat{\phi }}^{\pm }.$$

Similarly, the assumption $\frac{{\nabla }^2\left|{\widehat{B}}_1^{+}\right|}{\left|{\widehat{B}}_1^{+}\right|}>>{|\mathbf{\nabla }{\widehat{\phi }}^{+}|}^2$ results in

$$\epsilon =\frac{-1}{{\omega }^2{\mu }_0}\frac{{\nabla }^2\left|{\widehat{B}}_1^{+}\right|}{\left|{\widehat{B}}_1^{+}\right|}.$$

(51)

Clearly, in this case only the magnitude of the ${\widehat{B}}_1^{+}$ field is required to determine the permittivity.

This EPT method is similar to H-EPT, but allows for conductivity or permittivity map** without requiring the availability of both the magnitude and phase of the ${\widehat{B}}_1^{+}$ field if the corresponding additional assumptions hold. The validity of these assumptions need to be investigated further. If only one of the EP maps is required, this approach enables, for example, shorter acquisition times or an increase in signal-to-noise ratio (SNR) of the transmit field map.

4.2.1. Poisson-Based Conductivity Map**

Poisson-based conductivity map** (P-CM) considers Equation (50) as a Poisson equation for the phase. More precisely, in P-CM, we consider the Poisson equation

$${\nabla }^2{\widehat{\varphi }}^{+}=\omega {\mu }_0\sigma$$

on ${\mathbb{R}}^3$ and observe that the right-handed side has the object domain $\mathbb{D}$ as its support. Requiring that $\widehat{\phi }$ decays sufficiently fast at infinity ($|\widehat{\phi }|$ decreases as $1/|\mathit{r}|$ uniformly in $\mathit{r}/|\mathit{r}|$ as $|\mathit{r}|\to \infty$), we have

$${\widehat{\varphi }}^{+}\left(\mathit{r}\right)=-\omega {\mu }_0{\int }_{{\mathit{r}}^{\prime }\in \mathbb{D}}{\widehat{G}}_{\mathrm{P}}\left(\mathit{r}-{\mathit{r}}^{\prime }\right)\sigma \left({\mathit{r}}^{\prime }\right)\mathrm{d}V$$

(52)

where ${\widehat{G}}_{\mathrm{P}}\left(\mathit{r}\right)$ is the 3D static Green’s function given by

$${\widehat{G}}_{\mathrm{P}}\left(\mathit{r}\right)=\frac{1}{4\pi |\mathit{r}|},\mathrm{for}\mathit{r}\ne \mathbf{0}.$$

In P-CM, we assume that the phase of the transmit field ${\widehat{\phi }}^{+}$ is known within the object, let $\mathit{r}\in \mathbb{D}$ in Equation (52), set ${\widehat{\varphi }}^{+}\left(\mathit{r}\right)={\widehat{\phi }}^{+}\left(\mathit{r}\right)$ for $\mathit{r}\in \mathbb{D}$ and retrieve a conductivity profile by minimizing Equation (52) in a least-squares sense.

P-CM is an integral formulation of the methods described in [38,39]. Its global integral approach has an inherent noise suppression effect which makes this method more robust to noise than local differentiation methods. Additionally, the minimization process allows for the inclusion of regularization as well. However, the method has an increased computational complexity compared to differential Helmholtz-based EPT approaches.

4.3. Local Maxwell Tomography

The simplified form of Local Maxwell Tomography (LMT) assumes the availability of a multi-element array and substitutes the polar decomposition of the ${\widehat{B}}_1^{+}$ field as presented in Equation (16) into the Helmholtz equation (Equation (29)) [40]. We then obtain

$$\frac{{\nabla }^2{\widehat{B}}_{1p}^{+;\mathrm{TRX}(p,q)}}{{\widehat{B}}_{1p}^{+;\mathrm{TRX}(p,q)}}=2\left(\mathbf{\nabla }{\widehat{\phi }}_{pq}^{\pm }-\mathrm{j}\frac{\mathbf{\nabla }\left|{\widehat{B}}_{1p}^{+}\right|}{\left|{\widehat{B}}_{1p}^{+}\right|}\right)·\mathbf{\nabla }{\widehat{\phi }}_q^{-}+\left({\left|\mathbf{\nabla }{\widehat{\phi }}_q^{-}\right|}^2-\mathrm{j}{\nabla }^2{\widehat{\phi }}_q^{-}\right)-k^2.$$

(53)

This local equation is assumed to hold inside the object domain $\mathbb{D}$ and can be written as

$${\mathsf{a}}^T\left(\mathit{r}\right)\mathsf{x}\left(\mathit{r}\right)=b\left(\mathit{r}\right)\mathrm{with}\mathit{r}\in \mathbb{D},$$

(54)

with

$$\mathsf{a}\left(\mathit{r}\right)=\left[\begin{array}{c}2{\partial }_x{\widehat{\phi }}_{pq}^{\pm }\left(\mathit{r}\right)-2\mathrm{j}\frac{{\partial }_x\left|{\widehat{B}}_{1p}^{+}\left(\mathit{r}\right)\right|}{\left|{\widehat{B}}_{1p}^{+}\left(\mathit{r}\right)\right|}\\ 2{\partial }_y{\widehat{\phi }}_{pq}^{\pm }\left(\mathit{r}\right)-2\mathrm{j}\frac{{\partial }_y\left|{\widehat{B}}_{1p}^{+}\left(\mathit{r}\right)\right|}{\left|{\widehat{B}}_{1p}^{+}\left(\mathit{r}\right)\right|}\\ 2{\partial }_z{\widehat{\phi }}_{pq}^{\pm }\left(\mathit{r}\right)-2\mathrm{j}\frac{{\partial }_z\left|{\widehat{B}}_{1p}^{+}\left(\mathit{r}\right)\right|}{\left|{\widehat{B}}_{1p}^{+}\left(\mathit{r}\right)\right|}\\ 1\\ 1\end{array}\right],\mathsf{x}\left(\mathit{r}\right)=\left[\begin{array}{c}{\partial }_x{\widehat{\phi }}_q^{-}\left(\mathit{r}\right)\\ {\partial }_y{\widehat{\phi }}_q^{-}\left(\mathit{r}\right)\\ {\partial }_z{\widehat{\phi }}_q^{-}\left(\mathit{r}\right)\\ {\left|\mathbf{\nabla }{\widehat{\phi }}_q^{-}\left(\mathit{r}\right)\right|}^2-\mathrm{j}{\nabla }^2{\widehat{\phi }}_q^{-}\left(\mathit{r}\right)\\ -k^2\left(\mathit{r}\right)\end{array}\right],$$

and

$$b\left(\mathit{r}\right)=\frac{{\nabla }^2{\widehat{B}}_{1p}^{+;\mathrm{TRX}(p,q)}\left(\mathit{r}\right)}{{\widehat{B}}_{1p}^{+;\mathrm{TRX}(p,q)}\left(\mathit{r}\right)},$$

Requiring that Equation (54) holds at N different locations with position vectors ${\mathit{r}}_n\in \mathbb{D}$, $n=1,2,\dots ,N$ (e.g., with N the total number of pixels/voxels and ${\mathit{r}}_n$ the position vector of the center of the nth pixel/voxel), we obtain the set of equations ${\mathsf{a}}^T\left({\mathit{r}}_n\right)\mathsf{x}\left({\mathit{r}}_n\right)=b\left({\mathit{r}}_n\right)$ for $n=1,2,\dots ,N$, which can be written as an underdetermined system $\mathsf{A}\mathsf{x}=\mathsf{b}$, where $\mathsf{A}$ is an N-by-$5N$ matrix given by

$$\mathsf{A}=\left[\begin{array}{cccc}{\mathsf{a}}^T\left({\mathit{r}}_1\right)& & & \\ & {\mathsf{a}}^T\left({\mathit{r}}_2\right)& & \\ & & \ddots & \\ & & & {\mathsf{a}}^T\left({\mathit{r}}_N\right)\end{array}\right]$$

(55)

and

$$\mathsf{x}={\left[\begin{array}{c}{\mathsf{x}}^T\left({\mathit{r}}_1\right),{\mathsf{x}}^T\left({\mathit{r}}_2\right),\dots ,{\mathsf{x}}^T\left({\mathit{r}}_N\right)\end{array}\right]}^T\mathrm{and}\mathsf{b}={\left[\begin{array}{c}b\left({\mathit{r}}_1\right),b\left({\mathit{r}}_2\right),\dots b\left({\mathit{r}}_N\right)\end{array}\right]}^T.$$

(56)

Since there are five unknowns associated with each point of interest ${\mathit{r}}_n$, at least five linearly independent transmit field measurements are carried out, producing the set of equations ${\mathsf{A}}_i\mathsf{x}={\mathsf{b}}_i$, $i=1,2,\dots ,I$, where $I\ge 5$ is the total number of transmit field measurements. The total set of field equations can now be written as

$$\left[\begin{array}{c}{\mathsf{A}}_1\\ {\mathsf{A}}_2\\ ⋮\\ {\mathsf{A}}_I\end{array}\right]\mathsf{x}=\left[\begin{array}{c}{\mathsf{b}}_1\\ {\mathsf{b}}_2\\ ⋮\\ {\mathsf{b}}_I\end{array}\right]$$

(57)

and this square ($I=5$) or overdetermined ($I>5$) system is solved in the least-squares sense to obtain vector $\mathsf{x}$. Finally, the EPs at location ${\mathit{r}}_n$ can be obtained by equating the fifth entry in $\mathsf{x}\left({\mathit{r}}_n\right)$ to $-k^2\left({\mathit{r}}_n\right)$.

This method requires no knowledge of the unavailable absolute transmit phase. However, since there are multiple unknowns for each point of interest, several independent transmit field measurements are required, which are typically only available on 7 T MRI systems, to derive a unique solution. The amount of required transmit fields can be reduced by extending the method with receive field measurements. The same procedure as described above can be carried out for the receive field in terms of measurable quantities, as presented in Equation (18), and, under the assumption of homogeneous proton density, a similar equation may be derived [40]. When more field maps are available, this last assumption can be prevented and the gradient and Laplacian of the proton density can also be determined [40]. LMT can be further generalized to also take the spatial variations of the tissue EPs into account, such that it becomes free from object and field assumptions. This, however, comes at the cost of increasing the number of unknowns and therefore requiring a larger amount of transmit and/or receive field maps [41].

4.4. Modified Dual-Excitation EPT

Modified dual-excitation EPT (MDE-EPT) uses Equation (33) as a starting point and it assumes that the gradient of the z-component of the magnetic flux density vanishes [28]. We then obtain

$${\nabla }^2{\widehat{B}}_1^{+}=\mathit{f}·\mathbf{\nabla }{\widehat{B}}_1^{+}-k^2{\widehat{B}}_1^{+}=4{\partial }^{-}{\widehat{B}}_1^{+}{g}^{+}+g_z{\partial }_z{\widehat{B}}_1^{+}-k^2{\widehat{B}}_1^{+}.$$

(58)

This local equation is assumed to hold inside the object domain $\mathbb{D}$ and can be written as

$${\mathsf{a}}^T\left(\mathit{r}\right)\mathsf{x}\left(\mathit{r}\right)=b\left(\mathit{r}\right),\mathit{r}\in \mathbb{D},$$

(59)

where $\mathsf{a}\left(\mathit{r}\right)$ and $\mathsf{x}\left(\mathit{r}\right)$ are 3-by-1 vectors given by

$$\mathsf{a}\left(\mathit{r}\right)=\left[\begin{array}{c}{\partial }^{-}{\widehat{B}}_1^{+}\left(\mathit{r}\right)\\ {\partial }_z{\widehat{B}}_1^{+}\left(\mathit{r}\right)\\ {\widehat{B}}_1^{+}\left(\mathit{r}\right)\end{array}\right]\mathrm{and}\mathsf{x}\left(\mathit{r}\right)=\left[\begin{array}{c}4g^{+}\left(\mathit{r}\right)\\ g_z\left(\mathit{r}\right)\\ -k^2\left(\mathit{r}\right)\end{array}\right]$$

and $b\left(\mathit{r}\right)={\nabla }^2{\widehat{B}}_1^{+}\left(\mathit{r}\right)$. Requiring that Equation (59) holds at N different locations with position vectors ${\mathit{r}}_n\in \mathbb{D}$, $n=1,2,\dots ,N$, leading to a system of equations $\mathsf{A}\mathsf{x}=\mathsf{b}$, where the N-by-$3N$ matrix $\mathsf{A}$, the $3N$-by-1 vector $\mathsf{x}$ and the N-by-1 vector $\mathsf{b}$ are of similar form as in LMT (cf. Equations (55) and (56)). Since there are three unknowns ($g^{+}\left({\mathit{r}}_n\right)$, $g_z\left({\mathit{r}}_n\right)$, and $-k^2\left({\mathit{r}}_n\right)$) associated with each point of interest ${\mathit{r}}_n$, at least three linearly independent transmit field measurements are carried out, producing the set of equations ${\mathsf{A}}_i\mathsf{x}={\mathsf{b}}_i$, $i=1,2,\dots ,I$, where $I\ge 3$ is the total number of transmit field measurements. The total set of field equations can now again be written as Equation (57) and this square ($I=3$) or overdetermined ($I>3$) system is solved in the least-squares sense to obtain vector $\mathsf{x}$. Finally, the EPs at location ${\mathit{r}}_n$ can be obtained by equating the third entry in $\mathsf{x}\left({\mathit{r}}_n\right)$ to $-k^2\left({\mathit{r}}_n\right)$.

This approach does not require homogeneity of the object, which allows for improved tissue boundary reconstructions. However, since there are three unknowns, at least three independent transmit fields are required. Additionally, the method is restricted to regions with spatially invariant z-component of the magnetic field. Note that the original form, dual-excitation EPT, assumed knowledge of the unavailable x- and y-components of the magnetic fields [42] and therefore required only two linearly independent excitations/measurements to determine the EP maps. MDE-EPT, however, can be extended by including Equation (39), again under the assumption of vanishing gradient of the z-component of the magnetic field, into the system of equations. The name dual-excitation is then again justified, in the sense that two independent excitations result in four equations if both the transmit and receive fields are acquired [43].

4.5. Gradient-Based EPT

Gradient-based EPT (G-EPT) continues from Equation (58) and writes it in terms of absolute and relative transmit phases with respect to a reference element, as presented in Equation (17) [44,45,46,47,48]. We then obtain

$$\begin{array}{cc}\hfill {\nabla }^2{\widehat{B}}_{1p}^{+;\mathrm{rel}(p,r)}& =-2\mathrm{j}\mathbf{\nabla }{\widehat{B}}_{1p}^{+;\mathrm{rel}(p,r)}·\mathbf{\nabla }{\widehat{\phi }}_r^{+}-{\widehat{B}}_{1p}^{+;\mathrm{rel}(p,r)}\left(-|\mathbf{\nabla }{\widehat{\phi }}_r^{+}{|}^2+\mathrm{j}{\nabla }^2{\widehat{\phi }}_r^{+}\right)\hfill \\ \hfill & +4\left({\partial }^{-}{\widehat{B}}_{1p}^{+;\mathrm{rel}(p,r)}+\mathrm{j}{\widehat{B}}_{1p}^{+;\mathrm{rel}(p,r)}{\partial }^{-}{\widehat{\phi }}_r^{+}\right)g^{+}\hfill \\ \hfill & +\left({\partial }_z{\widehat{B}}_{1p}^{+;\mathrm{rel}(p,r)}+\mathrm{j}{\widehat{B}}_{1p}^{+;\mathrm{rel}(p,r)}{\partial }_z{\widehat{\phi }}_r^{+}\right)g_z-k^2{\widehat{B}}_{1p}^{+;\mathrm{rel}(p,r)},\hfill \end{array}$$

(60)

with ${\widehat{\phi }}_r^{+}$ the unknown absolute transmit phase of reference channel r. First, the gradient $g^{+}$ is determined. Similar to LMT and MDE-EPT, the above equation is written in the form

$${\mathsf{a}}^T\left(\mathit{r}\right)\mathsf{x}\left(\mathit{r}\right)=b\left(\mathit{r}\right)\mathrm{with}\mathit{r}\in \mathbb{D},$$

(61)

where

$$\mathsf{a}\left(\mathit{r}\right)=\left[\begin{array}{c}{\widehat{B}}_{1p}^{+;\mathrm{rel}(p,r)}\left(\mathit{r}\right)\\ {\partial }_x{\widehat{B}}_{1p}^{+;\mathrm{rel}(p,r)}\left(\mathit{r}\right)\\ {\partial }_y{\widehat{B}}_{1p}^{+;\mathrm{rel}(p,r)}\left(\mathit{r}\right)\\ {\partial }_z{\widehat{B}}_{1p}^{+;\mathrm{rel}(p,r)}\left(\mathit{r}\right)\end{array}\right],$$

$$\mathsf{x}\left(\mathit{r}\right)=\left[\begin{array}{c}|\mathbf{\nabla }{\widehat{\phi }}_r^{+}{\left(\mathit{r}\right)|}^2-\mathrm{j}{\nabla }^2{\widehat{\phi }}_r^{+}\left(\mathit{r}\right)+4\mathrm{j}{\partial }^{-}{\widehat{\phi }}_r^{+}\left(\mathit{r}\right)g^{+}\left(\mathit{r}\right)+\mathrm{j}{\partial }_z{\widehat{\phi }}_r^{+}\left(\mathit{r}\right)g_z\left(\mathit{r}\right)-k^2\left(\mathit{r}\right)\\ -2\mathrm{j}{\partial }_x{\widehat{\phi }}_r^{+}\left(\mathit{r}\right)+4g^{+}\left(\mathit{r}\right)\\ -2\mathrm{j}{\partial }_y{\widehat{\phi }}_r^{+}\left(\mathit{r}\right)-4\mathrm{j}{g}^{+}\left(\mathit{r}\right)\\ -2\mathrm{j}{\partial }_z{\widehat{\phi }}_r^{+}\left(\mathit{r}\right)+g_z\left(\mathit{r}\right)\end{array}\right],$$

and $b\left(\mathit{r}\right)={\nabla }^2{\widehat{B}}_{1p}^{+;\mathrm{rel}(p,r)}\left(\mathit{r}\right)$. Equation (61) is required to hold at N different locations of interest with position vectors ${\mathit{r}}_n\in \mathbb{D}$, $n=1,2,\dots ,N$, leading to a system of equations $\mathsf{A}\mathsf{x}=\mathsf{b}$, where the N-by-$4N$ matrix $\mathsf{A}$, the $4N$-by-1 vector $\mathsf{x}$ and the N-by-1 vector $\mathsf{b}$ are of a similar form as in LMT and MDE-EPT (cf. Equations (55) and (56)).

Since there are four unknowns associated with each point of interest (the elements of vector $\mathsf{x}\left({\mathit{r}}_n\right)$), at least four linearly independent transmit field measurements are carried out and these produce the set of equations ${\mathsf{A}}_i\mathsf{x}={\mathsf{b}}_i$, $i=1,2,\dots ,I$, with $I\ge 4$ the total number of transmit field measurements. The total set of equations can now again be written as in Equation (57) and this square ($I=4$) or overdetermined ($I>4$) system is solved in the least-squares sense to obtain vector $\mathsf{x}$. From this vector, $g^{+}\left({\mathit{r}}_n\right)$ can be determined from the second or third entry of $\mathsf{x}\left({\mathit{r}}_n\right)$.

Second, the gradient is integrated using the definition $g^{+}={\partial }^{+}ln\left(\eta \right)$ and an additional least-squares minimization process, where seed points (point belonging to a subdomain of the object domain where the EPs are known) are used to obtain absolute EP maps.

The additional integration step in G-EPT acts as a low-pass filter and makes the approach relatively robust to noise. Additionally, using the relative phase has the benefit that influences of receive field, chemical shift, magnetic susceptibility and eddy currents on the phase are mitigated. The method, however, requires multiple transmit elements as well as knowledge of seed points to derive absolute EP maps. The seed points can be derived by surrounding the object with a gel with known EPs, (dubbed boundary informed G-EPT) [48]. Additionally, since the transverse gradients of the absolute phase of the reference channel can also be derived from $\mathsf{x}\left(2\right)$ and $\mathsf{x}\left(3\right)$ in the first step of G-EPT, the seed points can be selected in an automated fashion by using the Helmholtz-based EPT approach in homogeneous regions (dubbed automated G-EPT) [47].

4.6. Convection–Reaction EPT

Convection–reaction EPT (CR-EPT) [49,50,51] assumes that the ${\widehat{B}}_1^{+}$ field is known and solves the generalized Helmholtz equation in convection–reaction form (Equation (36)), for convenience repeated here,

$$u{\nabla }^2{\widehat{B}}_1^{+}+{\mathbf{\beta }}^{+}·\mathbf{\nabla }u-\zeta {\widehat{B}}_1^{+}=0,$$

(62)

in a least-squares sense for the inverse admittance parameter u under the assumption of invariance of the z-component of the magnetic flux density in the convective field (Equation (35)), and derives the tissue parameters as

$$\sigma =Re\left(\frac{1}{u}\right),\mathrm{and}\epsilon =\frac{1}{\omega }Im\left(\frac{1}{u}\right).$$

(63)

This method is again not restricted to regions with homogeneous tissue structures, does not require seed points and does not require a multi-element array. However, the absolute transmit phase is again required, which is not directly available from measurements but can be accurately estimated for many cases as half the transceive phase. Furthermore, the method is restricted to regions with spatially invariant z-component of the magnetic field. Additionally, the method suffers from a reconstruction artifact in the region with low convective field.

4.6.1. Phase-Only Convection–Reaction Conductivity Map**

Phase-only convection–reaction conductivity map** (PCR-CM) [52] simplifies the generalized Helmholtz equation in convection–reaction form (Equation (36)) by dividing it by ${\widehat{B}}_1^{+}$ and assuming $\mathbf{\nabla }\left|{\widehat{B}}_1^{+}\right|=\mathbf{0}$ and $\mathbf{\nabla }{\widehat{B}}_z=\mathbf{0}$, which gives

$$u\left(-{\left|\mathbf{\nabla }{\widehat{\phi }}^{+}\right|}^2+\mathrm{j}{\nabla }^2{\widehat{\phi }}^{+}\right)+{\tilde{\mathbf{\beta }}}^{+}·\mathbf{\nabla }u-\zeta =0,$$

(64)

with

$${\tilde{\mathbf{\beta }}}^{+}=2\mathrm{j}{\partial }^{-}{\widehat{\phi }}^{+}{\mathbf{i}}_x-2{\partial }^{-}{\widehat{\phi }}^{+}{\mathbf{i}}_y+\mathrm{j}{\partial }_z{\widehat{\phi }}^{+}{\mathbf{i}}_z.$$

The same procedure can be performed for the generalized Helmholtz equation in convection–reaction form in terms of receive fields (Equation (40)), which yields

$$u\left(-{\left|\mathbf{\nabla }{\widehat{\phi }}^{-}\right|}^2-\mathrm{j}{\nabla }^2{\widehat{\phi }}^{-}\right)+{\tilde{\mathbf{\beta }}}^{-}·\mathbf{\nabla }u-\zeta =0,$$

(65)

with

$${\tilde{\mathbf{\beta }}}^{-}=-2\mathrm{j}{\partial }^{+}{\widehat{\phi }}^{-}{\mathbf{i}}_x-2{\partial }^{+}{\widehat{\phi }}^{-}{\mathbf{i}}_y-\mathrm{j}{\partial }_z{\widehat{\phi }}^{-}{\mathbf{i}}_z.$$

The addition of these two equations gives

$$u\left(-{\left|\mathbf{\nabla }{\widehat{\phi }}^{\pm }\right|}^2-2\mathbf{\nabla }{\widehat{\phi }}^{+}·\mathbf{\nabla }{\widehat{\phi }}^{-}+\mathrm{j}{\nabla }^2{\widehat{\phi }}^{\pm }\right)+{\tilde{\mathbf{\beta }}}^{\pm }·\mathbf{\nabla }u-2\zeta =0,$$

(66)

with

$${\tilde{\mathbf{\beta }}}^{\pm }=\left[{\partial }_y\left({\widehat{\phi }}^{+}+{\widehat{\phi }}^{-}\right)+\mathrm{j}{\partial }_x{\widehat{\phi }}^{\pm }\right]{\mathbf{i}}_x-\left[{\partial }_x\left({\widehat{\phi }}^{+}+{\widehat{\phi }}^{-}\right)+\mathrm{j}{\partial }_y{\widehat{\phi }}^{\pm }\right]{\mathbf{i}}_y+\mathrm{j}{\partial }_z{\widehat{\phi }}^{\pm }.$$

In the case that $Im\left({\tilde{\mathbf{\beta }}}^{\pm }\right)·Re\left(\mathbf{\nabla }u\right)>>Re\left({\tilde{\mathbf{\beta }}}^{\pm }\right)·Im\left(\mathbf{\nabla }u\right)$, the imaginary part of this equation can be written as a convection–reaction equation in terms of the resistivity $\rho ={\sigma }^{-1}$ as

$$\rho {\nabla }^2{\widehat{\phi }}^{\pm }+\mathbf{\nabla }{\widehat{\phi }}^{\pm }·\mathbf{\nabla }\rho -2\omega {\mu }_0=0,$$

(67)

which can be solved in a least-squares sense. This equation is in the form of a convection–diffusion–reaction equation with zero diffusion term. A diffusion term would act as a low-pass filter and increases numerical stability of the approach. To suppress spurious oscillations, an artificial diffusion term $c{\nabla }^2\rho$ is typically added to the fundamental equation, where c is an empirically determined constant diffusion coefficient. The conductivity $\sigma$ is finally retrieved as the inverse of the resistivity $\rho$.

This method can be seen as a generalized version of phase-only Helmholtz-based EPT implementation as discussed in Section 4.2, which allows for large spatial variations of the tissue conductivity. However, the method has an increased computational complexity and the required assumptions do not hold for high field strengths.

4.7. Transverse EPT

Transverse EPT (T-EPT) [53] assumes that the RF field has a so-called E-polarized field structure within a certain transverse plane, by which we mean that longitudinal variations of the transverse electric field and the longitudinal variation of the magnetic field essentially vanish within this plane (${\partial }_z{\widehat{E}}_x={\partial }_z{\widehat{E}}_y=0$, and ${\partial }_z{\widehat{B}}_z=0$ for $z=\mathrm{constant}$). Usually, the plane $z=\mathrm{constant}$ is taken to be the midplane of a birdcage coil, since it has been observed that the RF field has an approximate E-polarized field structure within this midplane [54]. Note that, for two-dimensional configurations with no spatial variations in the z-direction and a z-directed external electric current source, the E-polarized field structure is exact.

Taking the E-polarized field assumption into account, it follows from Maxwell’s equations that (cf. Equations (24) and (37))

$$\frac{4}{\mathrm{j}{\mu }_0}{\partial }^{-}{\widehat{B}}_1^{+}=\eta {\widehat{E}}_z\mathrm{and}{\widehat{B}}_1^{+}=\frac{1}{\omega }{\partial }^{+}{\widehat{E}}_z,\mathrm{within}\mathrm{the}\mathrm{plane}z=\mathrm{constant}.$$

(68)

In T-EPT, these two equations are combined into a single normalized functional given by

$$F\left({\widehat{E}}_z,\eta \right)=\frac{{∥\frac{4}{\mathrm{j}{\mu }_0}{\partial }^{-}{\widehat{B}}_1^{+}-\eta {\widehat{E}}_z∥}_{\mathbb{D}}^2}{{∥\frac{4}{\mathrm{j}{\mu }_0}{\partial }^{-}{\widehat{B}}_1^{+}∥}_{\mathbb{D}}^2}+\frac{{∥{\widehat{B}}_1^{+}-\frac{1}{\omega }{\partial }^{+}{\widehat{E}}_z∥}_{\mathbb{D}}^2}{{∥{\widehat{B}}_1^{+}∥}_{\mathbb{D}}^2},$$

(69)

where ${\parallel ·\parallel }_{\mathbb{D}}$ is an $L^2$-norm defined on $\mathbb{D}$. Ths functional is iteratively minimized in an alternate fashion using conjugate-gradient-type update formulas for ${\widehat{E}}_z$, followed by a direct update of $\eta$. This two-step update procedure is repeated until convergence or a maximum number of iterations is reached. Finally, the conductivity and permittivity reconstructions follow from the reconstructed admittance $\eta$ as

$$\sigma =Re\left(\eta \right)\mathrm{and}\epsilon =\frac{1}{\omega }Im\left(\eta \right).$$

(70)

In the remainder of this manuscript, these multi-step inversion methods are summarized in a listing (see Listing 1 for the update process of T-EPT).

This method has no second-order but only first-order derivatives that act on the measurement data, increasing noise robustness. Additionally, the method computes the z-component of the electric field strength, which can be helpful in SAR computations. However, the method is restricted to regions where the RF field is approximately E-polarized, such as in the midplane of a birdcage RF coil.

Listing 1. Transverse EPT (T-EPT).

Given initial guesses ${\tilde{\eta }}^{\left[0\right]}$ for the admittance and ${\tilde{E}}_z^{\left[0\right]}$ for the electric field strength For $n=1,2,\dots$ (a) Fix the admittance ${\tilde{\eta }}^{[n-1]}$ and update the electric field strength according to the update formula $${\tilde{E}}_z^{\left[n\right]}={\tilde{E}}_z^{[n-1]}+{\alpha }^{\left[n\right]}{v}^{\left[n\right]}$$ where $\alpha$ are the update coefficients and v the Polak-Ribière update directions [55]. (b) Update the admittance according to $${\tilde{\eta }}^{\left[n\right]}=\frac{4}{\mathrm{j}{\mu }_0}\frac{{\tilde{E}}_z^{\left[n\right]\ast }{\partial }^{-}{\widehat{B}}_1^{+}}{{\left|{\tilde{E}}_z^{\left[n\right]}\right|}^2}.$$ (c) Stop if objective function is smaller than user specified tolerance level, or if maximum number of iterations has been reached. End

4.8. First-Order Induced-Current EPT

First-order induced-current EPT (foIC-EPT) considers E-polarized RF fields and thus assumes that the electric field strength is mainly directed in the longitudinal z-direction [56]. For E-polarized fields, we have (cf. Equations (37) and (44))

$$\eta {\widehat{E}}_z=\frac{4}{\mathrm{j}{\mu }_0}{\partial }^{-}{\widehat{B}}_1^{+}\mathrm{and}{\widehat{E}}_z^{\mathrm{sca}}=-\zeta {\mathcal{G}}_{\mathrm{A}}\left\{\eta {\widehat{E}}_z\right\}-k_0^2{\mathcal{G}}_{\mathrm{A}}\left\{{\widehat{E}}_z\right\},$$

(71)

where we introduce the vector potential operator

$${\mathcal{G}}_{\mathrm{A}}\left\{\mathit{x}\right\}\left(\mathit{r}\right)={\eta }_0^{-1}{\int }_{{\mathit{r}}^{\prime }\in \mathbb{D}}\widehat{G}\left(\mathit{r}-{\mathit{r}}^{\prime }\right)\mathit{x}\left({\mathit{r}}^{\prime }\right)\mathrm{d}V.$$

Combining these two equations, together with the linearity of Maxwell’s equations (see Equation (42)), gives

$${\widehat{E}}_z^{\mathrm{inc}}-4\omega {\mathcal{G}}_{\mathrm{A}}\left\{{\partial }^{-}{\widehat{B}}_1^{+}\right\}={\widehat{E}}_z+k_0^2{\mathcal{G}}_{\mathrm{A}}\left\{{\widehat{E}}_z\right\},$$

(72)

which can be solved for the z-component of the electric field strength, if its incident component is known. Finally, the conductivity and permittivity can be derived via

$$\sigma =Re\left(\frac{4}{\mathrm{j}{\mu }_0}\frac{{\widehat{E}}_z^{\ast }{\partial }^{-}{\widehat{B}}_1^{+}}{|{\widehat{E}}_z{|}^2}\right)\mathrm{and}\epsilon =\frac{1}{\omega }Im\left(\frac{4}{\mathrm{j}{\mu }_0}\frac{{\widehat{E}}_z^{\ast }{\partial }^{-}{\widehat{B}}_1^{+}}{|{\widehat{E}}_z{|}^2}\right).$$

(73)

This approach has only first-order derivatives that act on the measured transmit field data, and an integral formulation for the electric field strength determination, making the method robust to noise. However, the method is again restricted to a region with an E-polarized field structure. Additionally, the method requires knowledge of the incident field, which cannot be measured directly. Incident fields are typically estimated from a simulation setup or from a reference scan of a phantom with known EPs. Note that the formulation is presented as a three-dimensional problem, but, in the case of an E-polarized field structure, it can be simplified to a two-dimensional setting. Cauchy-based EPT shares a lot of similarities with foIC-EPT, but the electric field strength is derived via a Cauchy integral which allows for the computation of the EPs in a direct manner through complex analysis [57,58,59,60].

4.9. Variational Born Iterative Method EPT

The variational Born iterative method EPT (VBIM-EPT) is a volumetric integral method that iteratively updates the tissue parameters based on improved estimations of the transmit field by solving forward and inverse problems [61,62]. Given knowledge of the incident fields and an initial estimation for the contrast function, the electric field strength is derived from (cf. Equations (42) and (44))

$${\widehat{\mathit{E}}}^{\mathrm{inc}}\left(\mathit{r}\right)=\widehat{\mathit{E}}\left(\mathit{r}\right)-\left(k_0^2+\mathbf{\nabla }\mathbf{\nabla }·\right){\int }_{{\mathit{r}}^{\prime }\in \mathbb{D}}\widehat{G}\left(\mathit{r}-{\mathit{r}}^{\prime }\right)\chi \left({\mathit{r}}^{\prime }\right)\widehat{\mathit{E}}\left({\mathit{r}}^{\prime }\right)\mathrm{d}V,$$

(74)

with $\chi =\left(\eta -{\eta }_0\right){\eta }_0^{-1}$. Equation (74) is a forward problem. Based on the derived electric field strength and the estimated contrast function, an estimate of the scattered part of the transmit field is computed as (cf. Equation (45))

$${\tilde{B}}_1^{+;\mathrm{sca}}\left(\mathit{r}\right)=\frac{\omega }{c_0^2}\tilde{\mathbf{\nabla }}·{\int }_{{\mathit{r}}^{\prime }\in \mathbb{D}}\widehat{G}\left(\mathit{r}-{\mathit{r}}^{\prime }\right)\chi \left({\mathit{r}}^{\prime }\right)\widehat{\mathit{E}}\left({\mathit{r}}^{\prime }\right)\mathrm{d}V,$$

(75)

and the residual $\delta {\widehat{B}}_1^{+;\mathrm{sca}}$ is derived according to $\delta {\widehat{B}}_1^{+;\mathrm{sca}}={\widehat{B}}_1^{+}-{\widehat{B}}_1^{+;\mathrm{inc}}-{\tilde{B}}_1^{+;\mathrm{sca}}$. The residual in the contrast function is then determined by solving

$$\delta {\widehat{B}}_1^{+;\mathrm{sca}}\left(\mathit{r}\right)=\frac{\omega }{c_0^2}\tilde{\mathbf{\nabla }}·{\int }_{{\mathit{r}}^{\prime }\in \mathbb{D}}\widehat{G}\left(\mathit{r}-{\mathit{r}}^{\prime }\right)\delta \chi \left({\mathit{r}}^{\prime }\right)\widehat{\mathit{E}}\left({\mathit{r}}^{\prime }\right)\mathrm{d}V$$

(76)

for $\delta \chi$, which is an inverse problem. The contrast function is then updated as ${\chi }^{[n+1]}={\chi }^{\left[n\right]}+\delta \chi$. Based on this new estimation of the contrast function, the procedure is repeated until a convergence criterion has been reached (see Listing 2). Finally, the conductivity and permittivity maps are derived via

$$\sigma =-\omega {\epsilon }_{0}Im\left(\chi \right),\mathrm{and}\epsilon ={\epsilon }_{0}Re\left(\chi +1\right)$$

(77)

This method does not apply any derivatives on the measured transmit field. Instead, it makes use of an integral formulation, making the method noise robust. However, the method requires knowledge of the incident fields, and solving the forward and inverse problems iteratively is computationally prohibitively expensive.

Listing 2. Variational Born Iterative Method-EPT (VBIM-EPT).

Given initial guesses ${\tilde{\chi }}^{\left[0\right]}$ for the contrast function For $n=1,2,\dots$ (a) Fix the contrast function to ${\tilde{\chi }}^{[n-1]}$ and determine the electric field strength ${\tilde{\mathit{E}}}^{\left[n\right]}$ by solving Equation (74) for $\widehat{\mathit{E}}$ (solve the forward problem). (b) Knowing contrast function ${\tilde{\chi }}^{[n-1]}$ and corresponding electric field strength ${\tilde{\mathit{E}}}^{\left[n\right]}$, compute the scattered magnetic flux density ${\tilde{B}}_1^{+;\mathrm{sca};\left[n\right]}$ according to Equation (75). (c) Compute the residual $\delta {\widehat{B}}_1^{+;\mathrm{sca};\left[n\right]}$ according to $$\delta {\widehat{B}}_1^{+;\mathrm{sca};\left[n\right]}={\widehat{B}}_1^{+}-{\widehat{B}}_1^{+;\mathrm{inc}}-{\tilde{B}}_1^{+;\mathrm{sca};\left[n\right]}.$$ (d) Fix the data residual to $\delta {\widehat{B}}_1^{+;\mathrm{sca};\left[n\right]}$ and the electric field strength to ${\tilde{\mathit{E}}}^{\left[n\right]}$ and determine the contrast residual $\delta {\chi }^{\left[n\right]}$ by solving Equation (76) for $\delta \chi$ (solve the inverse problem). (e) Update the contrast function according to the update formula $${\tilde{\chi }}^{\left[n\right]}={\tilde{\chi }}^{[n-1]}+\delta {\chi }^{\left[n\right]}.$$ (f) Stop if $\delta {\widehat{B}}_1^{+;\mathrm{sca};\left[n\right]}$ is smaller than user specified tolerance level, or if maximum number of iterations has been reached. End

4.10. Global Maxwell Tomography

Global Maxwell tomography (GMT) is a volumetric integral method that iteratively updates the tissue parameters based on improved estimations of the transmit field by solving a forward problem and minimizing an objective function [63,64,65]. GMT makes use of the identity [66]

$$\left(k_0^2+\mathbf{\nabla }\mathbf{\nabla }·\right){\int }_{{\mathit{r}}^{\prime }\in \mathbb{D}}\widehat{G}\left(\mathit{r}-{\mathit{r}}^{\prime }\right)\widehat{\mathit{w}}\left({\mathit{r}}^{\prime }\right)\mathrm{d}V=\mathbf{\nabla }\times \mathbf{\nabla }\times {\int }_{{\mathit{r}}^{\prime }\in \mathbb{D}}\widehat{G}\left(\mathit{r}-{\mathit{r}}^{\prime }\right)\widehat{\mathit{w}}\left({\mathit{r}}^{\prime }\right)\mathrm{d}V-\widehat{\mathit{w}},\mathit{r}\in \mathbb{D},$$

and transforms the electric field integral representation of Equation (74) into a current density volume integral representation, given by

$${\eta }_0\chi \left(\mathit{r}\right){\widehat{\mathit{E}}}^{\mathrm{inc}}\left(\mathit{r}\right)=\left(1+\chi \left(\mathit{r}\right)\right){\widehat{\mathit{J}}}^{\mathrm{sca}}\left(\mathit{r}\right)-\chi \left(\mathit{r}\right)\mathbf{\nabla }\times \mathbf{\nabla }\times {\eta }_0^{-1}{\int }_{{\mathit{r}}^{\prime }\in \mathbb{D}}\widehat{G}\left(\mathit{r}-{\mathit{r}}^{\prime }\right){\widehat{\mathit{J}}}^{\mathrm{sca}}\left(\mathit{r}\right)\mathrm{d}V.$$

(78)

Given knowledge of the incident field and an initial contrast function, this equation is solved for the scattered current density distribution (a forward problem), which is used to estimate the scattered component of the transmit field ${\tilde{B}}_1^{+;\mathrm{sca}}$ via Equation (45). An objective function is introduced

$$F\left(\chi \right)=\frac{{∥{\widehat{B}}_1^{+;\mathrm{sca}}-{\tilde{B}}_1^{+;\mathrm{sca}}\left({\widehat{\mathit{J}}}^{\mathrm{sca}}\left(\chi \right)\right)∥}_{\mathbb{D}}^2}{{∥{\widehat{B}}_1^{+;\mathrm{sca}}∥}_{\mathbb{D}}^2}.$$

(79)

Based on its gradient with respect to $\chi$ the contrast function is updated. This process is iterated until a convergence criterion has been reached (see Listing 3). Finally, the EPs are derived from the contrast function via Equation (77).

This method is similar to VBIM-EPT, but removes the inverse problem in every iteration. Even though a computational expensive inverse problem is removed, the method remains computationally expensive since the gradient updates typically require a large amount of iterations. The presented formulation still requires knowledge of the absolute transmit phase; however, this has been addressed by reformulating the objective function to only consider the magnitude of the transmit field or, in the case of a multi-element transmit system, by reformulating it in terms of magnitude and relative phases [64,65].

Listing 3. Global Maxwell Tomography (GMT).

Given initial guess ${\tilde{\chi }}^{\left[0\right]}$ for the contrast function For $n=1,2,\dots$ (a) Fix the contrast function to ${\tilde{\chi }}^{[n-1]}$ and determine the scattered current density distribution ${\tilde{J}}^{\mathrm{sca};\left[n\right]}$ by solving Equation (78) for ${\widehat{J}}^{\mathrm{sca}}$ (solve the forward problem). (b) Fix the scattered current density distribution to ${\tilde{\mathit{J}}}^{\mathrm{sca};\left[n\right]}$ and update the contrast function according to the update formula $${\tilde{\chi }}^{\left[n\right]}={\tilde{\chi }}^{[n-1]}+{\beta }^{\left[n\right]}{d}^{\left[n\right]}.$$ where $\beta$ are the update coefficients and d the Polak-Ribière update directions [55] (c) Stop if objective function of Equation (79) is smaller than user specified tolerance level, or if maximum number of iterations has been reached. End

4.11. Contrast Source Inversion EPT

Contrast-Source Inversion EPT (CSI-EPT) formulates the inversion problem as a purely optimization problem in which a single functional is iteratively minimized [67,68]. CSI-EPT combines the multiplication of the contrast function and the electric field strength into a single variable, the so-called contrast source $\widehat{\mathit{w}}=\chi \widehat{\mathit{E}}$. The scattered electric field strength is then given by (cf. Equation (44))

$${\widehat{\mathit{E}}}^{\mathrm{sca}}\left(\mathit{r}\right)=\left(k_0^2+\mathbf{\nabla }\mathbf{\nabla }·\right){\int }_{{\mathit{r}}^{\prime }\in \mathbb{D}}\widehat{G}\left(\mathit{r}-{\mathit{r}}^{\prime }\right)\widehat{\mathit{w}}\left({\mathit{r}}^{\prime }\right){\mathcal{G}}_{\mathrm{E}}\left\{\widehat{\mathit{w}}\right\}\left(\mathit{r}\right),$$

and the scattered transmit field operator is then given by (cf. Equation (45))

$${\widehat{B}}_1^{+;\mathrm{sca}}\left(\mathit{r}\right)=\frac{\omega }{c_0^2}\tilde{\mathbf{\nabla }}·{\int }_{{\mathit{r}}^{\prime }\in \mathbb{D}}\widehat{G}\left(\mathit{r}-{\mathit{r}}^{\prime }\right)\widehat{\mathit{w}}\left({\mathit{r}}^{\prime }\right){\mathcal{G}}_{\mathrm{B}}\left\{\widehat{\mathit{w}}\right\}\left(\mathit{r}\right),$$

which are used to set up an objective functional (cf. Equation (74))

$$F\left(\widehat{\mathit{w}},\chi \right)=\frac{{∥\chi {\widehat{\mathit{E}}}^{\mathrm{inc}}-\widehat{\mathit{w}}+\chi {\mathcal{G}}_{\mathrm{E}}\left\{\widehat{\mathit{w}}\right\}∥}_{\mathbb{D}}^2}{{∥\chi {\widehat{\mathit{E}}}^{\mathrm{inc}}∥}_{\mathbb{D}}^2}+\frac{{∥{\widehat{B}}_1^{+;\mathrm{sca}}-{\mathcal{G}}_{\mathrm{B}}\left\{\widehat{\mathit{w}}\right\}∥}_{\mathbb{D}}^2}{{∥{\widehat{B}}_1^{+;\mathrm{sca}}∥}_{\mathbb{D}}^2},$$

(80)

which is minimized in a two-step “fix-one-minimize-for-the-other” update process. First, the contrast function is fixed and the contrast source is updated from the gradient of the cost function with respect to $\widehat{\mathit{w}}$. Once the contrast source is updated, the electric field strength is calculated as

$$\widehat{\mathit{E}}={\widehat{\mathit{E}}}^{\mathrm{inc}}+{\mathcal{G}}_{\mathrm{E}}\left\{\widehat{\mathit{w}}\right\},$$

and the contrast function is updated by solving the least-squares problem ${∥\tilde{\chi }\tilde{\mathit{E}}-\tilde{\mathit{w}}∥}_{\mathbb{D}}^2$, which gives

$$\chi =\frac{\widehat{\mathit{w}}·{\widehat{\mathit{E}}}^{\ast }}{{\left|\widehat{\mathit{E}}\right|}^2},$$

or by fixing the contrast source and updating the contrast function from the gradient of the cost functional with respect to $\chi$. This two-step update procedure is iteratively repeated until a stop** criterion has been reached (see Listing 4). Finally, the tissue parameters are derived from the contrast function via Equation (77).

This approach only applies fast forward computations and does not have to solve any forward problem as in VBIM-EPT or GMT. However, since the approach typically still requires a lot of iterations, the method remains time consuming compared to direct methods such as H-EPT. Additionally, the transmit phase remains required, which can only be accurately approximated in specific cases. Naive [69,70,71,72], two-dimensional [73,74,75,76,77,78], magnitude-based [79] and segmented [80,81] implementations have been proposed to improve for example convergence or applicability.

Listing 4. Contrast Source Inversion-EPT (CSI-EPT).

Given initial guesses ${\tilde{\chi }}^{\left[0\right]}$ and ${\tilde{\mathit{w}}}^{\left[0\right]}$ for the contrast function and contrast source, respectively For $n=1,2,\dots$ (a) Fix the contrast function to ${\tilde{\chi }}^{[n-1]}$ and update the contrast source according to the update formula [67] $${\tilde{\mathit{w}}}^{\left[n\right]}={\tilde{\mathit{w}}}^{[n-1]}+{\alpha }^{\left[n\right]}{\mathbf{v}}^{\left[n\right]}.$$ (b) Compute the corresponding electric field strength ${\widehat{\mathit{E}}}^{\left[n\right]}$ according to $${\tilde{\mathit{E}}}^{\left[n\right]}={\widehat{\mathit{E}}}^{\mathrm{inc}}+{\mathcal{G}}_{\mathrm{E}}\left\{{\tilde{\mathit{w}}}^{\left[n\right]}\right\}.$$ (c) Compute the contrast function according to $$\tilde{\chi }=\frac{\tilde{\mathit{w}}·{\tilde{\mathit{E}}}^{\ast }}{{\left|\tilde{\mathit{E}}\right|}^2},$$ or fix the contrast source to ${\tilde{\mathit{w}}}^{\left[n\right]}$ and update the contrast function according to the update formula [67] $${\tilde{\chi }}^{\left[n\right]}={\tilde{\chi }}^{[n-1]}+{\beta }^{\left[n\right]}{d}^{\left[n\right]}.$$ (d) Stop if objective function of Equation (80) is smaller than user specified tolerance level, or if maximum number of iterations has been reached. End

5. EPT Methods Not Requiring Transmit Field Map**

The previously discussed EPT approaches can be extended with receive fields. However, there are also methods that do not require transmit fields. The followin EPT approaches are discussed in this section:

• Section 5.1: Single-acquisition EPT (SA-EPT)
• Section 5.2: Image-based EPT (I-EPT)

5.1. Single-Acquisition EPT

Single-acquisition EPT (SA-EPT) [82] rewrites the Helmholtz equation for the receive field (Equation (41)) as

$$\begin{array}{cc}\hfill \frac{{\nabla }^2{\widehat{B}}_1^{-;\ast }}{{\widehat{B}}_1^{-;\ast }}& =\frac{\mathbf{\nabla }·\left({\widehat{B}}_1^{-;\ast }\frac{\mathbf{\nabla }{\widehat{B}}_1^{-;\ast }}{{\widehat{B}}_1^{-;\ast }}\right)}{{\widehat{B}}_1^{-;\ast }},\hfill \\ \hfill & =\frac{\mathbf{\nabla }{\widehat{B}}_1^{-;\ast }}{{\widehat{B}}_1^{-;\ast }}·\frac{\mathbf{\nabla }{\widehat{B}}_1^{-;\ast }}{{\widehat{B}}_1^{-;\ast }}+\mathbf{\nabla }·\left(\frac{\mathbf{\nabla }{\widehat{B}}_1^{-;\ast }}{{\widehat{B}}_1^{-;\ast }}\right)=-k^2,\hfill \end{array}$$

(81)

which shows that knowledge of the rate of change of the field is sufficient to derive the EPs. By introducing the receive field of channel q relative to reference channel r as

$${\widehat{B}}_{1qr}^{-;\ast }=\frac{{\widehat{B}}_{1q}^{-;\ast }}{{\widehat{B}}_{1r}^{-;\ast }},$$

(82)

we obtain through the Laplacian of the relative receive field [82]

$$\begin{array}{cc}\hfill {\nabla }^2{\widehat{B}}_{1qr}^{-;\ast }& =\frac{{\widehat{B}}_{1q}^{-;\ast }}{{\widehat{B}}_{1r}^{-;\ast }}\left(\frac{{\nabla }^2{\widehat{B}}_{1q}^{-;\ast }}{{\widehat{B}}_{1q}^{-;\ast }}-\frac{{\nabla }^2{\widehat{B}}_{1r}^{-;\ast }}{{\widehat{B}}_{1r}^{-;\ast }}\right)-2\mathbf{\nabla }{\widehat{B}}_{1qr}^{-;\ast }·\frac{\mathbf{\nabla }{\widehat{B}}_{1r}^{-;\ast }}{{\widehat{B}}_{1r}^{-;\ast }},\hfill \\ \hfill & =-2\mathbf{\nabla }{\widehat{B}}_{1qr}^{-;\ast }·\frac{\mathbf{\nabla }{\widehat{B}}_{1r}^{-;\ast }}{{\widehat{B}}_{1r}^{-;\ast }},\hfill \end{array}$$

(83)

since each element measures the same EPs. This local equation is assumed to hold inside the object domain $\mathbb{D}$ and can be written as

$${\mathsf{a}}^T\left(\mathit{r}\right)\mathsf{x}\left(\mathit{r}\right)=b\left(\mathit{r}\right),\mathit{r}\in \mathbb{D},$$

(84)

where $\mathsf{a}\left(\mathit{r}\right)$ and $\mathsf{x}\left(\mathit{r}\right)$ are 3-by-1 vectors given by

$$\mathsf{a}\left(\mathit{r}\right)=\left[\begin{array}{c}-2{\partial }_x{\widehat{B}}_{1qr}^{-;\ast }\left(\mathit{r}\right)\\ -2{\partial }_y{\widehat{B}}_{1qr}^{-;\ast }\left(\mathit{r}\right)\\ -2{\partial }_z{\widehat{B}}_{1qr}^{-;\ast }\left(\mathit{r}\right)\end{array}\right]\mathrm{and}\mathsf{x}\left(\mathit{r}\right)=\left[\begin{array}{c}\frac{{\partial }_x{\widehat{B}}_{1r}^{-;\ast }}{{\widehat{B}}_{1r}^{-;\ast }}\left(\mathit{r}\right)\\ \frac{{\partial }_x{\widehat{B}}_{1r}^{-;\ast }}{{\widehat{B}}_{1r}^{-;\ast }}\left(\mathit{r}\right)\\ \frac{{\partial }_x{\widehat{B}}_{1r}^{-;\ast }}{{\widehat{B}}_{1r}^{-;\ast }}\left(\mathit{r}\right)\end{array}\right]$$

and $b\left(\mathit{r}\right)={\nabla }^2{\widehat{B}}_{1qr}^{-;\ast }\left(\mathit{r}\right)$. Equation (84) is required to hold at N different locations of interest with position vectors ${\mathit{r}}_n\in \mathbb{D}$, $n=1,2,\dots ,N$, leading to a system of equations $\mathsf{A}\mathsf{x}=\mathsf{b}$, where the N-by-$3N$ matrix $\mathsf{A}$, the $3N$-by-1 vector $\mathsf{x}$ and the N-by-1 vector $\mathsf{b}$ are of a similar form as in LMT, MDE-EPT and G-EPT (cf. Equations (55) and (56)). Once the $\frac{\mathbf{\nabla }{\widehat{B}}_1^{-;\ast }}{{\widehat{B}}_1^{-;\ast }}$ term is obtained, the tissue parameters can be determined from Equation (81).

This method does not require absolute transmit field data, but relies only on relative receive fields which are directly available. This results in the elimination of specific artifacts, since common terms for the different elements can be eliminated. These receive fields can be derived from a single acquisition; however, this requires a multi-element array with a minimum of four receive elements, since at least three linearly independent relative receive fields are required to determine a unique solution. Additionally, since the gradient term is derived in a minimization process based on second-order derivatives and an additional divergence is applied on the gradient term, the method relies on third-order derivatives giving strong dependence on the SNR of the images.

5.2. Image-Based EPT

Image-based EPT (I-EPT) uses the acquired MR image directly for the reconstruction of the tissue parameters. For a low flip angle $\alpha =\gamma \tau \left|{\widehat{B}}_1^{+}\right|$, we have $sin\alpha \approx \alpha$, thus the MR image for any low-flip-angle sequence is essentially given by (cf. Equation (15))

$$I={\varrho }_0\gamma \tau {\widehat{B}}_1^{+}{\widehat{B}}_1^{-;\ast }.$$

(85)

In I-EPT, information from this image is used, instead of estimated transmit or receive fields. The relevant equation applied to this image is derived by multiplying Equation (29) with ${\widehat{B}}_1^{-;\ast }$ and Equation (41) with ${\widehat{B}}_1^{+}$ and adding them together, which gives

$${\nabla }^2\left({\widehat{B}}_1^{+}{\widehat{B}}_1^{-;\ast }\right)+2k^2{\widehat{B}}_1^{+}{\widehat{B}}_1^{-;\ast }-2\mathbf{\nabla }{\widehat{B}}_1^{+}·\mathbf{\nabla }{\widehat{B}}_1^{-;\ast }=0.$$

(86)

By taking ${\widehat{B}}_1^{+}{\widehat{B}}_1^{-;\ast }={\left(\sqrt{{\widehat{B}}_1^{+}{\widehat{B}}_1^{-;\ast }}\right)}^2$ and defining $a=\sqrt{{\widehat{B}}_1^{+}{\widehat{B}}_1^{-;\ast }}$ and $b=\sqrt{\frac{{\widehat{B}}_1^{-;\ast }}{{\widehat{B}}_1^{+}}}$, we obtain

$${\nabla }^2{a}^2+2k^2{a}^2-2\mathbf{\nabla }\frac{a}{b}·\mathbf{\nabla }\left(ab\right)=0,$$

and by using the product rule of the scalar Laplacian (${\nabla }^2{a}^2=2a{\nabla }^{2}a+2\mathbf{\nabla }a·\mathbf{\nabla }a$) and dividing by $2a^2$, we find

$$\frac{{\nabla }^{2}a}{a}+k^2+\underset{(\ast )}{\underbrace{\frac{1}{a^2}\left(\mathbf{\nabla }a·\mathbf{\nabla }a-\mathbf{\nabla }\frac{a}{b}·\mathbf{\nabla }\left(ab\right)\right)}}=0.$$

(87)

The underbraced term denoted by $(\ast )$ is an error term that can be simplified to $\frac{1}{4}{\left|\mathbf{\nabla }ln\frac{{\widehat{B}}_1^{-;\ast }}{{\widehat{B}}_1^{+}}\right|}^2$ [83] and can be neglected when ${\widehat{B}}_1^{+}$ and ${\widehat{B}}_1^{-;\ast }$ are similar to each other. This results in the following Helmholtz equation

$$\frac{{\nabla }^2\sqrt{{\widehat{B}}_1^{+}{\widehat{B}}_1^{-;\ast }}}{\sqrt{{\widehat{B}}_1^{+}{\widehat{B}}_1^{-;\ast }}}=-k^2.$$

(88)

Equation (88) remains valid when the ${\widehat{B}}_1^{+}{\widehat{B}}_1^{-;\ast }$ term is multiplied by a constant, since it drops out of the equation. The variables in front of the ${\widehat{B}}_1^{+}{\widehat{B}}_1^{-;\ast }$ term in Equation (88) are relatively constant throughout space in regions where the Helmholtz equation applies, and the image as described in Equation (85) can therefore be applied in Equation (88).

This method does not require the acquisition of transmit and receive fields, which results in reduced scan time and an increase in SNR, since the image SNR is greater than that of transmit or receive field maps. Additionally, in this formulation, the errors resulting from ${\widehat{B}}_1^{+}$ and ${\widehat{B}}_1^{-}$ differences are reduced to a first-order effect with respect to the difference compared to the conventional H-EPT method. A zero echo-time sequence has been proposed due to its immunity to eddy current and static magnetic field ($B_0$) inhomogeneity-induced phase changes, as well as its speed and SNR efficiency [83]. The method has also been proposed with a fast spin echo sequence together with a $T_2$ relaxation pattern between echoes to increase noise robustness [84]. A generalized image-based EPT form which includes the gradient of the EPs has also been proposed [85].

6. Data Driven Deep Learning Approaches for Solving Inverse Problems

Solving inverse problem by data-driven deep learning approaches is an emerging field with recent examples from the fields of gravitation [86], radioastronomy [87], medical imaging reconstruction [88] and electromagnetics [89,90]. The basic advantage of these data-driven approaches is that it allows insertion of more tailored a priori information about the specific inverse problem under study. A key aspect is the learning-based approach where during a training phase deep neural networks learn to perform a specific task in the inverse process by feeding them with many ground truth examples. After training the neural network, the inference is extremely fast, sometimes only a few seconds. This constitutes a clear computational gain over more conventional iterative approaches to solve inverse problems. Other advantages include a higher level of tolerance to noise on the input data and a higher flexibility on the required input data as neural networks can act as learned surrogate models.

6.1. Convolutional Neural Networks

The most popular for image processing and reconstruction are convolutional neural networks (CNN) [91,92,93]. These neural networks employ convolutions in their architecture whose kernels are optimized during training for their given task. The training phase is a very computationally intensive process involving the backpropagation of the errors during the training (in fact, a large scale optimization process itself) over the various weights between the nodes (often >1 million) in the network. This process takes place on a GPU and is facilitated by mature software packages that are able to exploit the parallel nature of the GPU in an efficient manner.

To obtain an idea on what information a trained CNN triggers, it is insightful to process an input image through the trained network and generate the output of the various units in the layers as images (so-called feature maps) comprising of local and global predictive information to perform the task it was trained for. Although CNNs are also heavily used for classification problems such as image segmentation, for use in the EPT reconstruction, only the regression task is relevant. In a regression CNN, where the famous U-NET [93] is a prime example, an encoder and the image information is processed by various convolutional kernels and pooling operations and fed through activation layers which constitute non-linear elements and are essential for the ability to learn. During the downsampling path, spatial contextual information is learnt. In regression problems, this is followed by the decoder where convolution and upsampling takes place to recover the spatial information of the desired output matrix size. In fully convolutional networks, the architecture solely employs operations such as convolution, pooling, activation and upsampling. Avoiding fully connected layers makes the inference much faster as fewer weights are needed, and the network can work regardless of the original image size. Skip connections short-circuiting corresponding layers in the encoder and decoder such as in the popular U-NET [93] are essential to recover fine-grained spatial information lost in the pooling or downsampling layers.

6.2. Deep Learning for EPT Reconstruction: Single Feedforward Approaches

An important distinction can be made between deep learning inverse approaches where single feed forward networks are employed and hybrid approaches where deep neural networks are themselves embedded in the iterative optimization process of solving the inverse problem. In principle, both approaches can be applied to EPT. Using the feedforward approach, Mandija et al. [94] employing convolutional neural networks demonstrated that deep learning EPT (DL-EPT) can reconstruct more noise robust dielectric parameter maps than conventional Helmholtz based EPT. An essential element of this feed forward approach is that the network constitutes a surrogate EPT reconstruction model implicitly learnt from the training data and takes the measured complex transmit field information as the input. This learning-based approach creates more flexibility than state-of-the-art MR-EPT techniques, which require electromagnetic quantities dictated by electromagnetic first principles which are not accessible with MRI. As example, in DL-EPT, a feedforward network can be trained with MR accessible quantities (e.g., the transceive phase) only. Interestingly, Mandija et al. [94] demonstrated that, also for a deep learning approach, almost all predictive features to reconstruct electrical conductivity are contained in the transceive phase maps in accordance with our insight from electromagnetic principles underpinning conventional EPT.

6.3. Training Data and Generalization to Unseen Data

Essential of course is the availability of training data which can nowadays be easily generated by electromagnetic simulations including realistic RF coil models, phantoms and body models. In this way, a high degree of a priori knowledge, such as the specific MRI coil setup, can be introduced. The advantage of using this simulation-based approach to generate the training transmit field data with superimposed artificial noise, is that the ground truth is available. However, in silica training data might not reflect realistic experimental conditions. Therefore, approaches that employ training data reconstructed with more conventional EPT reconstruction schemes are also used [95]. In the work of Gavazzi et al. [96], a 3D patch neural network approach was used where the receptive field is more local (size of the 3D patch) forcing it to perform dielectric parameter estimation from more local ${\widehat{B}}_1^{+}$ magnitude and phase information. A further advantage is that it can work with varying matrix size of the input data. A key question is of course how a single feed forward neural network approach behaves when it is tested on unseen input data that was not directly included in the training data, e.g., pathologic tissue with different dielectric parameter values or in the presence of motion artifacts on the ${\widehat{B}}_1^{+}$ maps. An obvious mitigation is to augment the training data sufficiently (e.g., part of the transmit maps can be artificially corrupted with motion artifacts) to obtain more robust results.

6.4. Deep Learning EPT: Integrating Deep Learning into Iterative EPT Schemes

Another option to improve the generalization to unseen data is to retain the physics in the reconstruction framework while still benefiting from the advantages of deep learning. A first approach was published for EPT by Leijsen et al. [97], who demonstrated that initial estimates provided by deep learning led to better convergence for 3D CSI-EPT. This integration can be further improved. New hybrid approaches are now emerging in medical image construction where neural networks are embedded in conventional iterative reconstruction schemes [98,99,100,101,102]. The physics related to the reconstruction problem is still explicitly included by means of a physics-based forward operator (e.g., Fourier transform in case of MR image reconstruction). Experiences from the medical image reconstruction demonstrated that this leads to much better generalization to unseen data in the training [100]. The neural networks can be inserted in the iterative procedure for various tasks. For example, neural networks can be trained to learn regularization filters much more tailored to the specific application than applying standard regularization kernels [101]. Alternatively, the networks can be used to perform the update task, i.e., determining the update direction based on the data mismatch and the regularization term [98]. Employing such an approach, the convergence is often much faster as a priori information on the optimization landscape is learned in the training phase, enabling faster convergence. These hybrid approaches should also be possible to combine with iterative EPT schemes such as 3D CSI-EPT where the physics is included by a forward operator (e.g., Green’s function approach) linking a certain electrical property distribution to the measured data (${\widehat{B}}_1^{+}$ magnitude and transceive phase information). Such a methodology would be an ideal scenario as it would harvest the power of deep learning to accelerate reconstruction and include tailored a priori information from the learning phase, while still retaining to the physics-based modeling and the data consistency.

6.5. Outlook

It is clear that deep learning offers much benefit for EPT in terms of achieving higher quality reconstructions. The feedforward approaches using CNNs have demonstrated clear potential in terms of noise robustness, flexibility on inputs and computational speed. A key question is the generalization to data not encountered in the training. Interestingly, also in EPT’s sister field of Quantitative Susceptibility Map** deep learning improves the quality of susceptibility reconstruction over conventional methodology [103,104]. Promising results have been achieved here on the generalization issue. An attractive alternative might be the integration of learned networks into a conventional iterative EPT approach as occurs in medical image reconstruction. This still retains the physics, offering better generalization, while also being able to include more a priori information and providing faster reconstructions.

7. Discussion and Conclusions

This paper presents a mathematical analysis of a large number of different methods for EPT, each with their own relative strengths and weaknesses. By comparing the results from each approach, one can make a number of general statements, the most important of which are listed below.

7.1. Approach Description

EPT approaches can be sorted into several different categories. These categories can give some general insight into how the methods work, and what kind of restrictions exist. Here, we sort the methods into three categories.

• Differential methods or integral methods
• Local methods that reconstruct the EPs at a specific location by only taking the information from the direct neighbourhood into account, or global methods that take the whole imaging domain into account to reconstruct the EP maps as a whole
• Direct methods that act directly on the data to reconstruct the EPs, also called backward methods since they run ‘backwards’ from the measured field map to the underlying EPs, or forward methods that employ forward models or solve forward problems in the inversion scheme and act ‘indirectly’ on the data

For each of the transmit field-based methods discussed in the manuscript, the approach descriptions are assigned in

Table 1. Approach description for the EPT approaches as described in this manuscript.

	H-EPT	SH-EPT	LMT	MDE-EPT	G-EPT	CR-EPT	T-EPT	foIC-EPT	VBIM-EPT	GMT	CSI-EPT	SA-EPT	I-EPT
Differential (order)	✓(2)	✓(2)	✓(2)	✓(2)	✓(2)	✓(2)	✓(1)	✓(1)	✗	✗	✗	✓(3)	✓(2)
Integral	✗	✗	✗	✓	✗	✗	✗	✓	✓	✓	✓	✗	✗
Local	✓	✓	✓	✓	✓	✗	✓	✓	✗	✗	✗	✓	✓
Global	✗	✗	✗	✗	✓	✓	✓	✓	✓	✓	✓	✗	✗
Direct (backward)	✓	✓	✓	✓	✓	✓	✓	✓	✗	✗	✗	✓	✓
Forward (indirect)	✗	✗	✗	✗	✗	✗	✓	✓	✓	✓	✓	✗	✗

.

7.1.1. Differential vs. Integral

A general observation is that differential approaches have an inherent noise amplification, while integral approaches are more noise robust due to the inherent low-pass filtering properties of the relevant integrals. The higher is the order of the differentials acting on the data, the larger is the noise amplification. A comparison between a second-order differential approach and an integral approach is shown for simulated three-dimensional noisy data (SNR = 100, as defined in [94]) in Figure 2. It shows that integral methods are in essence more noise robust and that typical noise reduction implementations such as regularization do not overcome this disadvantage.

7.1.2. Local vs. Global

A commonality among many approaches is that there appears to be a trade-off between having an adverse noise effect in local methods due to higher-order derivatives (second order and up) acting on the data or having the bias effect of the EM field structure in global methods. A comparison between two-dimensional implementations of local and global approaches that assume knowledge of the complex transmit signal is shown for simulated two-dimensional noiseless data in Figure 3. The reconstructions of the local method H-EPT shows boundary errors due to assumed homogeneity of the underlying tissue, but the method is accurate in regions which have locally spatially invariant tissue properties. The reconstructions of the global methods (T-EPT, foIC-EPT and CSI-EPT) take the inhomogeneity of the EPs into account, but suffer from a bias related to the low electric field strength (low convective field). Additionally, global methods have the potential to reach local minima in their optimization process.

Global methods allow for the inclusion of regularization in the optimization problem which can be employed to correct for the bias, resolve local minima or improve noise robustness.

7.1.3. Direct vs. Forward

A general observation is that direct methods are relatively fast, while forward methods tend to be computational expensive and time consuming, especially those that require the results of forward and/or inverse problems iteratively. Forward methods have yet to be demonstrated to be clinically feasible. Forward methods, however, typically simultaneously reconstruct additional field maps, such as the electric field strength which would be useful for SAR computations.

7.2. Data Requirements

EPT approaches can also be categorized on which type of data they require, or what kind of assumptions about the data are required. For each of the transmit field-based methods discussed in the manuscript, data requirements are assigned in

Table 2. Data requirements for the EPT approaches as described in this manuscript. Note that data requirements can be influenced by extensions or generalizations of the methods.

	H-EPT	SH-EPT	LMT	MDE-EPT	G-EPT	CR-EPT	T-EPT	foIC-EPT	VBIM-EPT	GMT	CSI-EPT	SA-EPT	I-EPT
${\widehat{B}}_1$ term	${\widehat{B}}_1^{+}$	$\left|{\widehat{B}}_1^{+}\right|,{\widehat{\phi }}^{+}$	$\left|{\widehat{B}}_1^{+}\right|,{\widehat{\phi }}_{pq}^{\pm }$	${\widehat{B}}_1^{+}$	$\left|{\widehat{B}}_1^{+}\right|,{\widehat{\phi }}_{pr}^{+}$	${\widehat{B}}_1^{+}$	${\widehat{B}}_1^{+}$	${\widehat{B}}_1^{+}$	${\widehat{B}}_1^{+}$	${\widehat{B}}_1^{+}$	${\widehat{B}}_1^{+}$	$\left|{\widehat{B}}_1^{-}\right|,{\widehat{\phi }}_{pq}^{-}$	${\widehat{B}}_1^{+}{\widehat{B}}_1^{-;\ast }$
Multi-element array	✗	✗	✓	✓	✓	✗	✗	✗	✗	✗	✗	✓	✗
Incident fields	✗	✗	✗	✗	✗	✗	✗	✓	✓	✓	✓	✗	✗
Seed Points	✗	✗	✗	✗	✓	✗	✗	✗	✗	✗	✗	✗	✗
$\mathbf{\nabla }\eta =0$	✓	✓	✓	✗	✗	✗	✗	✗	✗	✗	✗	✓	✓
$\mathbf{\nabla }\left|{\widehat{B}}_1^{+}\right|=0$	✗	✓	✗	✗	✗	✗	✗	✗	✗	✗	✗	✗	✗
$\mathbf{\nabla }{\widehat{\phi }}^{+}=0$	✗	✓	✗	✗	✗	✗	✗	✗	✗	✗	✗	✗	✗
${\partial }_x{\widehat{B}}_z=0$	✗	✗	✗	✓	✓	✓	✗	✗	✗	✗	✗	✗	✗
${\partial }_y{\widehat{B}}_z=0$	✗	✗	✗	✓	✓	✓	✗	✗	✗	✗	✗	✗	✗
${\partial }_z{\widehat{B}}_z=0$	✗	✗	✗	✓	✓	✓	✓	✓	✗	✗	✗	✗	✗
${\partial }_z{\widehat{E}}_x=0$	✗	✗	✗	✗	✗	✗	✓	✓	✗	✗	✗	✗	✗
${\partial }_z{\widehat{E}}_y=0$	✗	✗	✗	✗	✗	✗	✓	✓	✗	✗	✗	✗	✗
${\widehat{B}}_1^{+}={\widehat{B}}_1^{-;\ast }$	✗	✗	✗	✗	✗	✗	✗	✗	✗	✗	✗	✗	✓

.

7.2.1. Measurable and Non-Measurable Data

For accurate reconstruction of the EPs, ideally measurements of all three components of the ${\mathit{B}}_1$ field would be possible. However, the z-component cannot be measured. Additionally, the x- and y-components cannot be acquired in a direct fashion, and determination would require the absolute transmit field as well as the absolute receive field. The absolute magnitude of the transmit field can be acquired in a direct fashion, while the absolute transmit phase and absolute receive field can not. Measured phases are always a superposition of the transmit and receive phase, and the magnitude of the receive field is always weighted by the proton-density. This data unavailability is one of the fundamental challenges that makes EPT complicated.

EPT approaches that assume availability of the complex transmit field typically estimate the absolute transmit phase by applying the transceive phase estimation (TPA), while methods that incorporate absolute receive fields go hand-in-hand with homogeneity or symmetry assumptions to eliminate the proton-density bias. To bypass these assumptions, solutions are often sought by reformulating the problem in terms of only directly-available field quantities. With regular RF coils, this comes down to formulating the inverse problem in only magnitude data or (in combination with) transceive phase data. However, with multi-element RF coil arrays, the acquisition of relative fields are possible, by dividing the complex signal measured in an element by the signal obtained in a particular reference element. This allows the derivation of the (gradient of the) transmit (and receive) phase, as well as EPT formulations based on the relative phase, instead of absolute or transceive phase, or formulations based on receive fields only. Additionally, this type of coil can eliminate specific artifacts since common terms for the different channels can be eliminated. However, the acquisition of multiple ${\widehat{B}}_1^{+}$ fields requires lengthy scans which can compromise patient comfort, throughput or SNR. Moreover, these multi-element RF coil elements are not yet widely available in clinical settings. About 50% of clinical MR scanners have a field strength of 1.5 T and have a body coil with a single transmit channel. About 45% of clinical MR scanners have a field strength of 3 T, of which the older ones have the same arrangement, and the newer one typically have two independent transmit channels that can produce different degrees of elliptically polarized RF fields. High field scanners, such as 7 T scanners, can have up to eight transmit channels with independent magnitude and phase control; however, there are only about 100 of these worldwide available.

Integral methods often require knowledge of the incident electric and magnetic field strength which are inaccessible with MRI. Typically, a reference scan from a phantom with known EPs or a simulation setup is used for estimation of the incident fields. However, since the incident fields are dependent on the loading of the coil, they are to a certain level patient-specific. Patient-specific coil–subject interactions remain unknown and a source of error.

7.2.2. Field/Object Structure

The inverse problem in EPT can be significantly simplified by assuming local homogeneity of the object. Methods that apply this local homogeneity assumption (LHA) are most often used in clinical studies due to their simplicity and ease of implementation. These methods however suffer from significant errors at tissue boundaries where the LHA is violated, making them impractical in regions with small tissue structures. For larger tissue structures, tissue segmentation can be used to improve boundary reconstructions.

The EPT problem can also be simplified by assuming an E-polarized field structure, i.e., assuming negligible (gradients of the) longitudinal component of the magnetic flux density, sometimes in combination with the assumption of vanishing (gradients of the) transverse components of the electric fields strength. This approximation is typically applied in the transverse midplane of a birdcage coil, where they are relatively small [13].

7.3. State Of Development

One might select a method based on different criteria, for example based on the SNR level, on the availability of multiple transmit or receive elements or incident fields, or whether or not the region of interest contains a homogeneous medium, a low electric field region or E-polarized fields. Due to the large number of EPT approaches with a large variety in requirements, assumptions and complexity, and since the EPT field is relatively new, most methods are not at the stage of clinical use yet (see

Table 3. Type of experiments performed for the EPT approaches as described in this manuscript. The table is constructed to the best of the authors knowledge. For references, see the main text.

	H-EPT	SH-EPT	LMT	MDE-EPT	G-EPT	CR-EPT	T-EPT	foIC-EPT	VBIM-EPT	GMT	CSI-EPT	SA-EPT	I-EPT
Simulation	✓	✓	✓	✓	✓	✓	✓	✓	✓	✓	✓	✓	✓
Phantom	✓	✓	✗	✗	✓	✓	✗	✓	✗	✓	✓	✓	✓
in vivo	✓	✓	✗	✗	✓	✗	✗	✗	✗	✗	✗	✗	✓
Clinical	✗	✓	✗	✗	✗	✗	✗	✗	✗	✗	✗	✗	✗

). To facilitate development, comparison and prototy** of EPT approaches, MATLAB code of the approaches presented in Figure 2 and Figure 3 is available upon request from the corresponding author.