Effects of Simulated Error-Sources on Different 3-D CSI-EPT Strategies

—Three-dimensional contrast source inversion-electrical properties tomography (3-D CSI-EPT) is an iterative reconstruction method that estimates the electrical properties of tissue from transmit ﬁeld magnetic resonance data. However, in order to bring 3-D CSI-EPT into practice for complex tissue structures and to understand the origin and effect of errors, insight in the sensitivities of reconstruction accuracy to the major error-sources is necessary. In this paper, different strategies for implementing 3-D CSI-EPT, including their iterative structure, are presented, of which the regularized implementation shows the most potential to be used in practice. Moreover, the inﬂuence of initialization, noise, stopping criteria, incident ﬁelds, B1-maps, transceive phase and domain truncation are discussed. We show that of all these different error-sources, initialization, accurate coil models and domain truncation have the most dramatic effect on electrical properties reconstructions in practice.


I. INTRODUCTION
K NOWLEDGE about tissue electrical properties (EPs), consisting of the conductivity (σ) and the permittivity (ε), is essential for determining accurate specific absorption rate levels in MRI and for patient-specific electromagnetic modeling that is performed prior to hyperthermia therapy treatment [1], [2].Furthermore, these parameters have the potential to be used as biomarkers in clinical applications, since it is possible to differentiate between ischemic and hemorrhagic strokes [3] and between benign and malignant breast lesions [4] based on the change in EPs.
The retrieval of EPs of biological tissue from measurements of magnetic fields generated by radiofrequency coils in a magnetic resonance imaging (MRI) scanner is called electrical properties tomography (EPT) and is a topic that has gained increasing attention in the last few years (see, for example, the three recent reviews [3], [5], [6]).Most EPT methods are based on a modified Helmholtz equation derived from Maxwell's equations in differential form, which shows the direct relationship between the EPs and the magnetic field strength (H), given by with k 2 = ω 2 μ(ε − jσ/ω) and η = σ + jωε, where ω denotes the angular frequency and μ the permeability.By assuming smooth transitions between tissue parameters, the first term on the right-hand side can be omitted and the conventional Helmholtz equation can be obtained.This allows for the computation of the EPs in a direct manner [7].The insight that the conductivity is primarily influenced by the phase, while the permittivity is mainly determined by the magnitude of the radiofrequency transmit field, resulted in phase-and magnitudeonly Helmholtz-based approaches [8]- [10].This additional simplification implies linearity of the phase-based equation which supersedes the requirement of the transceive phase assumption (the estimation of the transmit phase as half the transceive phase), and, since the magnitude map can be omitted, significantly speeds up the scan time which is invaluable for practical conductivity mapping implementations [11].However, since the method assumes spatial homogeneity of the underlying EPs, severe errors occur at boundaries between different tissue types.The additional homogeneity assumptions in the magnitude or phase of the transmit field degrade the accuracy in EP mapping even further.In addition, the second order differential operator makes them extremely sensitive to noise [12].More advanced reconstruction methods such as local Maxwell tomography [13], [14], convection-reaction EPT [15] and gradient-based EPT [16] solve Eq. ( 1) by deriving a set of linear equations, which enables improved reconstructions at tissue boundaries.The boundary problem is also solved in first-order induced-current EPT [17] and by making use of only first order differential operators and an integral formulation it reduces noise sensitivities.Cauchy-based EPT from [18] shows many similarities, but the EPs are reconstructed in a direct and explicit fashion.In global Maxwell tomography and contrast-source inversion EPT (CSI-EPT) integral formulations are used which avoid derivatives of the transmit MR data altogether.These methods were initially implemented in 2-D [19], [20], but have been extended to 3-D [21]- [24], which circumvents any assumption This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/ of longitudinal homogeneity of the object or electromagnetic fields, which is not applicable in general [24].
Three-dimensional CSI-EPT does not have to solve any timeconsuming forward problems, in contrast to global Maxwell tomography, and has been shown in simulations to be able to accurately reconstruct tissue transitions in three-dimensional settings [23].However, practical 3-D CSI-EPT implementations can produce several artifacts that are difficult to classify from measurement data only.In order to improve practical implementations, insights about the origin and effect of different error-sources are necessary.
In this manuscript, different strategies for implementing 3-D CSI-EPT are simulated: error-sources due to initialization, noise, iterative stopping criteria, incident fields, B 1 -mapping, transceive phase assumption and domain size are individually examined.The effect that regularization can have on the convergence of the cost functional and on overfitting to noise is also shown.Additionally, an overview is presented that shows where in the iterative process the computational expensive steps occur, and the time requirement of the most computationally expensive operators is discussed.

II. THEORY
In this section, the integral representations for the EM field and the CSI-EPT method are briefly summarized.The background is assumed to consists of free space, and the spatially-varying conductivity and permittivity of objects are isotropic at the Larmor frequency ω.The permeability μ is assumed to be tissueindependent and equal to that of free space.Throughout this manuscript, the time factor convention exp (+jωt) is adopted.

A. Fundamental Electromagnetic Field Representations
Let the magnetic resonance coil occupy a bounded source domain S and let an object occupy a disjoint object domain D. An externally applied current density distribution J ext on the coil generates an electromagnetic field {E, H} which is defined by the superposition of the incident and scattered fields.
The incident electromagnetic field is the field that is present inside the coil in the absence of an object and is given by and where A ext is the vector potential given by with G the scalar Green's function of the homogeneous background medium defined as with x = (x, y, z) the position vector, η 0 = jωε 0 the per-unitlength admittance of the background medium and k 0 = ω/c 0 the wavenumber of the background medium, with c 0 the electromagnetic wave speed in vacuum.Note that this Green's function assumes that there are no other current sources in infinite space other than the currents on the coil.The scattered field due to the presence of an object is given by and where A sca is the vector potential given by with w the contrast source, defined as w(x) = χ(x)E(x), with χ(x) = η(x)/η 0 − 1 the contrast function describing the object and containing the information about the EPs [25]- [27].The contrast source is the scattering current density distribution normalized by the free space per-unit-length admittance, including both the conduction and displacement current density distributions, that results due to the presence of the object.The scattering current density and consequently the scattered field would vanish if the object is absent (η = η 0 ).The contrast function is the additional per-unit-length admittance imposed on top of the free space per-unit-length admittance, normalized by free space per-unit-length admittance.Once the contrast source is established due to the object, then it can be viewed as a source situated in free space.

B. Contrast Source Inversion-Electrical Properties Tomography
In EPT the goal is to reconstruct the contrast function χ from B + 1 = μ 0 2 (H x + jH y ) field measurements.Several EPT approaches have operators that act directly on this field data to retrieve the EPs [3], [5], while in CSI-EPT operators effectively operate on the scattered B + 1 field component.The scattered component is denoted by and this field is taken as a starting point.In CSI-EPT, the incident field is assumed to be a known quantity which can be obtained from simulations based on a known coil geometry, or derived from phantom experiments with known electrical properties.The scattered field of Eq. ( 9) is then easily derived by subtracting the incident field from the measured B + 1 field.1) Naive CSI-EPT: CSI-EPT makes use of the fact that Eq. ( 9) can also be computed from knowledge of the contrast source w, via Eq.(7), and searches for a contrast source w such that the mismatch between the measured field and the predicted fields is minimized.This mismatch is represented by the data residual, given by In this equation, G B is the data operator that computes B +;sca 1 from estimations of the contrast source via with ∇ = − 1 2 (i x + ji y )∂ z + 1 2 (∂ x + j∂ y )i z .This mismatch can be minimized by iteratively minimizing the data functional [28] where || • || denotes the norm defined on D, which is introduced in Appendix I.
Once the optimal contrast source is obtained, the contrast function is obtained by minimizing the discrepancy in Maxwell's equations.This is realized by introducing the object residual, given by where G E is the object operator that computes E sca , given by This object residual is minimized when where the index k ranges over the set {x, y, z}.Hereafter, this method is referred to as N-CSI.
2) Traditional CSI-EPT: Traditional CSI-EPT ensures that Maxwell's equations are fulfilled at every iteration by enforcing an object functional, given by which is combined with the data functional into a total functional as which gives rise to a non-linear optimization problem in which the contrast source and contrast function are updated in a two-step manner at each iteration [23].The update step of the contrast source is presented in Appendix II-A and two different contrast function update steps, i.e. the traditionally used direct method as well as a conjugate gradient method, are presented in Appendix II-B.Throughout this manuscript, the minimization of this total functional F T is denoted by T-CSI.Specifically, T-CSI-dir denotes the application of the direct update step for the contrast function and T-CSI-cg denotes the use of the conjugate gradient update step.
3) Regularized CSI-EPT: When the data is perturbed by noise, the approach has the tendency to overfit.To reduce this effect, a multiplicative total variation regularization parameter is added to T-CSI.Since the regularizing objective function is multiplied by the sum of the data and objective functional instead of added to these functionals, this implementation has the advantage that no artificial regularization parameter is required, which otherwise would require extensive tuning and vary from problem to problem.It has been shown to be effective for 2-D CSI-EPT implementations [20], [29].The regularization functional F R is given by TV , (18) where the total variation factor is defined as [30] with the steering parameter chosen as [31] (δ Since the total variation factor is independent of the contrast source, the regularization does not change the update for w, but it does adjust the update for χ.The altered conjugate gradient updates for χ are presented in Appendix II-C.In this manuscript, this third regularized implementation of CSI-EPT is denoted by R-CSI.
The different cost functionals represent the error that still remains in the data and object equation.CSI-EPT continues the iterative process for example as long as the cost functional is larger than a predefined value, or a predefined maximum number of iterations has not yet been reached.The latter is used throughout this manuscript.
The three methods and their steps as described in Appendix II are presented as a flowchart in Fig. 1.After initialization, the iterative process of updating the contrast source (and the electric field and contrast function) is started, until a stopping criterion has been reached.The shades of red indicate the computational cost.Darker colors indicate a higher computational complexity.

III. SIMULATION SETUP
This study was performed in a 7 T neuroimaging configuration, where a head-sized birdcage transmit coil driven in quadrature at 300 MHz was simulated using XFdtd (XF7.5, Remcom State College, PA, USA).The birdcage coil model has 16 rungs, a length of 195 mm and a diameter of 300 mm, and is surrounded by an RF shield which is 220 mm in length and 360 mm in diameter.A resonant model of the coil was numerically tuned using 7.1 pF capacitors, and loaded by either the male body model Duke or the female body model Ella from the virtual family dataset with their corpus callosum centered in the coil [32].An ideal version of the coil was also simulated, where the capacitors were replaced by 50 Ohm voltage sources with 1 V amplitudes driven such that a circular polarization was obtained [33].The resulting current density on the coil structure was subsequently used to derive incident and total fields required for CSI-EPT using Eqs.( 2)-( 5), for each of the simulated configurations.The coils were also driven in anti-quadrature mode with an opposite phase shift between the sources, and the same steps were performed to obtain the receive fields that were used to obtain the transceive phase.All custom codes were implemented using MATLAB 2019b (The MathWorks, Inc., Natick, Massachusetts, United States).The The shades of red show an estimation of the most computational expensive steps.Darker red means a higher computational complexity.Note that the computation of E is also a computational expensive step in N-CSI.However, it only has to be computed once.total fields are computed for head models derived from the Duke body model.These head models occupy domains that range from 105×83×43 up to 105×83×73 voxels, each with dimensions of 2.5×2.5×3.5 mm 3 .
IV. STRATEGIES Fig. 2 shows the conductivity and relative permittivity maps of the true model and those from different reconstruction strategies, evaluated on the three-dimensional Duke head model with a domain size of 105 × 83 × 43 voxels, and initialized with a homogeneous mask containing the average values of the EPs (σ = 0.58 S/m, ε r = 43) followed by a forward computation to retrieve the initialization of the electric field strength and contrast source.This initialization has been used throughout this manuscript, unless stated otherwise.The depicted transversal and coronal slices are taken through the center of the head coil.Figs.2(b,g) show the reconstructions that are derived from N-CSI after 10 000 iterations.The approach converges to a reconstruction that is an extremely poor representation of the true EP model.The reconstructions with T-CSI-dir after 10 000 iterations, presented in Figs.2(c,h), give an improved reconstruction, but the quality is still poor.By updating the contrast function not via the direct method, but with the conjugate gradient method as discussed in Appendix II-B, a significant improvement in the reconstructions is achieved, see Figs. 2(d,i).Figs.2(e,j) show that R-CSI results in a similar reconstruction result as T-CSI-cg, as can be expected in the absence of noise.In these two approaches the low E-field issue results in a major reconstruction artifact [20], [23].Its location is for clarification indicated with a black arrow in the R-CSI permittivity map.Line profiles through the center of the object for the T-CSI-cg and R-CSI maps, as well as line profiles through the low E-field artifact for the R-CSI map are presented in Fig. S2 in the Supplementary Information.The mean and standard deviation of the three main tissue regions (white matter, gray matter and cerebrospinal fluid), as well as the relative residual error (RRE; see Appendix III-A) of the whole volume from the reconstruction results of Fig. 2 are presented in Table S2 V. ERROR-SOURCES Due to the fact that T-CSI-cg and R-CSI show an overall superior EP map reconstruction, we focus on these methods in our error-source analysis.

A. Initialization State
One of the key error-sources in CSI-EPT is the initialization, i.e. the starting condition for the iterative process.Proper initialization can prevent CSI-EPT from reaching a local minimum and reduce the computation time that is necessary to reach the tolerance level of the cost function.Reconstruction results from different initializations are presented in Section S1 in the Supplementary Information.The data show that initializing the reconstruction technique with, for example, backprojection can result in a stagnation of the iterative process, where a local minimum is reached that corresponds to inaccurate EPs.Certain local minima can be prevented by incorporating a priori information into the initialization, such as the average expected values for the EPs.The determination of a good initialization is, however, not straightforward.

B. Noise
A second important error-source in CSI-EPT is noise in the B + 1 map.Fig. 3 shows the T-CSI-cg and R-CSI reconstructions for noisy data (Gaussian noise added to the real and imaginary part of B + 1 ) with a signal-to-noise ratio (SNR; see Appendix III-B) of 75 and 50, to assess their noise sensitivity.The presence of noise degrades the reconstruction of smaller tissue structures, but the main structures are properly reconstructed for both SNR levels after 1000 iterations (see Figs. 3(a-d,g-j)).Note that the noise seems to reduce the low E-field artifact in the regularized reconstruction.
The degradation due to the presence of noise is also depicted in Fig. 4, which shows the relative residual error as a quantification of the results from the different strategies after different numbers  (i,j) in the Supplementary Information, which show similar results to the noiseless case shown in Figs.S2(e,f).

C. Stopping Criterion
Another parameter that influences CSI-EPT reconstructions is the stopping criterion.Fig. 4(b) additionally shows that, in the case of noise, an increase in RRE can occur with more iterations.The results of T-CSI-cg and R-CSI at 10 000 iterations of this noisy case are presented in Figs.3(e-f,k-l) and compared to the previously discussed 1000 iterations with the same noise level (Figs.3(c-d,i-j)).The T-CSI-cg reconstruction becomes significantly more noisy with more iterations, while the same effect is not observed for the R-CSI reconstruction.The differences with the noiseless case indicate that overfitting to noise is substantial in T-CSI-cg, while this is considerably less in R-CSI.

D. B1-Mapping
A measured B 1 -map (B + 1 magnitude map) can be perturbed by several inaccuracies, such as global offsets, low flip angle bias and other error-propagation mechanisms.The effect of an incorrect global magnitude scaling is shown in Fig. 5    1 magnitude map can result in inaccuracies in the CSI reconstruction.The heading represent the scaling factor that is applied on B + 1 , from 20% underestimation (a,f) up to 20% overestimation (e,j).Results are for R-CSI after 10 000 iterations.Conductivity (a-e) and relative permittivity (f-j).
32 and 17 for the scaling of 0.8, 0.9, 1.1 and 1.2, respectively.Large under-and overestimations show blurring in the reconstructed images, and the loss of detailed tissue structures.With a 20% under-or overestimation, global tissue structures remains visible, and with a 10% under-or overestimation relatively small tissue structures remain distinguishable.This shows that R-CSI reconstructions are fairly insensitive to small under-and overestimations.Spatially non-uniform noise, which would be more realistic as the data quality in areas of low excitation flip are known to suffer from increased B 1 mapping inaccuracies, resulted in very similar reconstruction quality and behavior (data not shown).

E. Coil Loading
The incident field is influenced by coil-subject interactions, called the loading effect, which is not explicitly accessible.Fig. 6.Incident fields and TPA effects.R-CSI reconstructions on noiseless data after 10 000 iterations when different incident fields are used (a-d,g-j) and when the TPA is applied to the simulated transceive phase (e-f,k-l).Using the incident field from the tuned coil model loaded with Ella (a,g), loaded with a sphere with average expected EP values (b,h), unloaded (c,i), and from the ideal coil loaded with the correct Duke model (d,j).Using the estimated transmit phase through the TPA in the case of the tuned coil setting (e,k) and in the case of the ideal coil setting (f,l).
Different incident fields are therefore simulated to assess the sensitivity of the R-CSI reconstruction to coil loading variations.While providing the total field data obtained in the Duke body model, the incident fields were generated using either the Ella body model, a homogeneous sphere (radius = 75 mm, σ = 0.58 S/m, ε r = 43), an empty coil, as well as the ideal coil model to determine the resulting error in the reconstructed EPs.The incident fields, with the exception of that of the empty coil, are scaled such that the mean of their corresponding absolute transmit field in the center region (3×3×3 voxels) matches the mean absolute transmit field of the tuned coil loaded with the Duke model in the same region.The reconstruction results are shown in Figs.6(a-d,g-j).It shows that the effects of slight inter-subject loading variations are not severe (Figs.6(a,g); RRE is 0.47 and 0.37 for the conductivity and permittivity, respectively).Loading the coil with the sphere results in blurring and the RRE increases to 0.52 and 0.42 for the conductivity and permittivity, respectively (Figs. 6(b,h)).When the incident field from an empty (unloaded) tuned coil is used, most of the detailed anatomical structure is lost, and the reconstructed images are very blurred (Figs.6(c,i)).Additionally, the exact coil model is typically not available.As an illustration, the reconstruction when the incident field from an ideal version of the coil is used is shown in Figs.6(d,j).The reconstruction results deviate severely from the ground truth, failing to capture any structural detail.

F. Transceive Phase
In EPT, the transceive phase, the additive combination of the transmit and the receive phase, is typically acquired.Figs.6(e,k) show the results when the transceive phase assumption (TPA), estimating the transmit phase as half the transceive phase, is applied on the 7 T tuned coil example.A similar reconstruction quality as using the true transmit field phase can be observed in the transversal conductivity slice, while several artifacts are observed in the coronal slice.In the permittivity map more substantial artifacts are observed, but the general outline of the tissue structure remains intact.Line plots are provided in Figs.S2(k,l) in the Supplementary Information.Figs.6(f,l) show the results when the TPA is applied on the ideal coil setup as introduced in Section III.Compared to the tuned coil model case, a lower reconstruction quality in the center transversal conductivity slice is observed.However, the coronal slice, as well as the permittivity results show a higher accuracy in the outer regions.The mean and standard deviation in the different segmentation regions are of comparable quality (see Table S2), but the RRE of the conductivity and permittivity maps are 0.53 and 0.40 for the tuned and 0.50 and 0.33 for the ideal setting, respectively, indicating a generally higher accuracy in the ideal coil setting.

G. Domain Truncation
In EM scattering formulations such as those underlying CSI-EPT, one typically assumes that the object is surrounded by air, i.e. that the object is fully enclosed within the object domain.In practise, however, B + 1 data is acquired within a region that is smaller than the entire subject.To evaluate the influence of outof-volume scattering on the reconstruction quality, we truncated the reconstructed domain to different sizes.The body model with a size as shown in Fig. 7(a) results in a |B + 1 | field as shown in Fig. 7(b).If the reconstructed domain is smaller, i.e. truncated to e.g.3.5, 7 or 10.5 cm, the |B + 1 | fields from forward simulations correspond to those presented in Figs.7(c-e), respectively.If the reconstruction domain is chosen as large as the object, the assumption of a vanishing object is satisfied, and an accurate reconstruction is achieved, as shown by the absolute contrast function in Fig. 7(g).However, if the assumption is violated, by reducing the reconstruction domain by e.g.3.5, 7 or 10.5 cm, reconstructions as shown in Figs.7(h-j) are achieved.Quantification values can be found in Table S2 and an additional line plot corresponding to Fig. 7(j) can be found in Figs.S2(m,n) in the Supplementary Information.The mean and standard deviation of the segmented regions as well as the line plot indicate a relatively good (but smoothed) reconstruction around the center region of the coil.The large RRE values, that refer to the general quality of the whole domain, indicate significant errors in other regions.Fig. 7 shows that these substantial errors occur mostly around the boundary region where the assumption of a vanishing object is violated.Note that in the previous results the |B +  1 | fields were taken equal to the one shown in Fig. 7(e), corresponding to the domain size 105 × 83 × 41 voxels, and thus neglecting the effects from the region below the brain.

VI. DISCUSSION AND CONCLUSION
In order to get more insight in the origin and effect of some of the most important error-sources on 3-D CSI-EPT that occur in practice, different update approaches for the contrast source as well as for the contrast function are compared, and the effects of (suboptimal) initialization, noise, stopping criteria, incident fields, B 1 -mapping, transceive phase and domain size on the reconstructions are addressed.Furthermore, the analytical implementation of different CSI-EPT strategies are summarized, and an overview is given that shows where the computationally expensive steps occur in the iterative process.
N-CSI applies the conjugate gradient approach to a linear problem, which gives fast results but is prone to reaching a local minimum.T-CSI-dir applies a two-step update approach to a non-linear problem which shows improved reconstructions: however, it is still prone to local minima and the direct update approach of the contrast function does not necessarily reduce errors [29].Improvements can be obtained by updating the contrast function more gradually in the iterative process as done in T-CSI-cg.In R-CSI multiplicative total variation regularization is employed and only minor differences are observed in the noiseless case compared to T-CSI-cg.Note that the RRE of the permittivity is generally lower compared to that of the conductivity, which is probably due to the B + 1 being mainly influenced by displacements currents and not conduction currents at 7 T.
When noise is present in the data, R-CSI outperforms the other strategies.Applying multiplicative regularization to CSI-EPT makes the method significantly more noise robust, while the others already are noise robust compared to non-integral equation-based EPT approaches [3], [5].Note that the stopping criterion is also a form of regularization.Throughout this manuscript the stopping criterion is selected to be a maximum amount of iterations: however, a predefined tolerance level for the cost functional, or step size could be used as well.R-CSI reduces the tendency of overfitting to noise at a high number of iterations and therefore does not require fine-tuning of the stopping criterion.
Spatially non-uniform noise distributions can occur due to B 1 mapping inaccuracies in areas of low excitation flip angle.These noise distributions resulted in similar reconstruction quality and behavior.Additionally, global scaling offsets in B 1 -mapping can occur due to unknown RF losses along the RF chain.Reconstructions with different offsets show that R-CSI is not very sensitive to small under-or overestimations of the B + 1 map.Offsets might be improved by calibrating the map with a reference phantom and/or by using directional couplers at the input of the coil.Note that due to the linearity of Maxwell's equations, the effect of inaccurate magnitude scaling can be reflected as an incorrectly scaled incident magnetic field strength, which therefore has similar effects on the reconstruction.CSI-EPT requires knowledge of the incident electric and magnetic fields.Here the incident fields are derived from the currents that run along the (loaded) coil and the shield in XFdtd simulations.In this way, the loading effect is taken into account, and the scattered fields created by the currents on the shield are interpreted as known incident fields such that the homogeneous Green's function is applicable.In practice, an accurate model of the coil setup is typically not available.In experimental settings it is possible to approximate components of the incident fields via measurements of phantoms with known EP values [34].However, reconstruction errors remain due to inaccuracies in the incident field because of coil loading variations, for example.Simulation results of the brain indicate that as long as a comparable reference object is used for determination of the incident fields, the remaining loading differences do not introduce severe errors in the reconstruction results.Since in the head model setting the incident field from the tuned coil loaded with a homogeneous spherical phantom did not show severe reconstruction differences compared to using the true incident field, we expect that an accurate body model would even suffice in regions with larger inter-subject variations, such as the abdominal region.
In EPT, the absolute transmit phase map cannot be acquired from a single-element coil, such as an ordinary birdcage coil.In practice the transceive phase, the superposition of the transmit and receive phase, is measured.For cylindrical objects and at low frequencies, the transmit phase from the coil in transmit mode is almost identical in structure to the receive phase of the coil in receive mode.In those cases, the transceive phase assumption can be used to estimate the transmit phase.The TPA results at 7 T show a smaller RRE in the ideal coil setting than in the tuned coil setting.This is as expected, since the transmit and receive phase are more symmetrical in the ideal setup.Our TPA results show a fairly good agreement with the ground truth, even at the relatively high field strength of 7 T.The results are expected to improve with lower field strengths, which are more often used in practice.Additional correction strategies can be applied, such as iterative updates of the receive phase, which has been shown to work for two-dimensional implementations of CSI-EPT [35].Some CSI-EPT implementations circumvent the transceive phase aspect by reformulating CSI-EPT such that it uses amplitude data only [36], [37].Other possible solutions can be sought by using transceiver arrays and formulating the problem in terms of a relative phase, as employed in gradientbased EPT [16], [38].
CSI-EPT implementations assume that the reconstruction domain fully embeds the object that contributes to the scattered B + 1 field to be able to make efficiently use of the FFT [20], [23].However, in practice, the region of acquisition is smaller than the full coverage of the coil in order to prevent excessive acquisition times.This results in a truncation of the reconstructed domain, which introduces errors due to the violation of the vanishing object assumption.These errors are most pronounced at the region close to the truncation, therefore the reconstructed domain should in practice always be larger than the domain of interest.In practice one should attempt to include all regions that contribute significantly to the B + 1 data in the acquisition and reconstruction domain.An increased domain leads to an increase in acquisition time as well as an increase in reconstruction time.However, due to the linearity of Maxwell's equations, solutions can be sought by projecting contributions of the scattered fields generated by an estimated (fictional) object outside the reconstruction domain onto the incident fields.
From the T-CSI-cg and R-CSI reconstructions a low electric field artifact can be observed, which typically occurs around the center of a volume RF coil [15], [23], [39].Reconstructions at those regions might be improved by incorporating complementary antenna settings or through active or passive shimming [11], [39].Note that the slices shown are through the center of the birdcage coil, which show the largest sensitivity to the low electric field issue, and other slices might therefore show better results.
Fig. 1 shows the differences between the various updating schemes of the discussed strategies and gives insight into the computation costs of each step.The total computation time for 1000 iterations for the object domain with 105×83×43 voxels are approximately 16, 37 and 44 minutes for N-CSI, T-CSI and R-CSI, respectively.More details about the computationally cost for different time consuming operators are presented in Section S2 in the Supplementary Information.Note that a significant computational speed-up can be achieved by making use of a GPU [40].
The fact that R-CSI outperforms the other EPT strategies discussed in this paper, without a substantial increase in computation time, gives this approach great potential for practical applications.To implement the method in practice, the initialization, coil model and domain truncation should be carefully taken care of, while noise level, stopping criterion, B 1 -mapping, coil loading and transceive phase problem play a more subordinate role in CSI-EPT.

APPENDIX
The steps in the contrast source inversion method are well described in the literature [30], [41], [42], but discussed once more in order to give a complete overview in the context of 3-D CSI-EPT.

APPENDIX I DEFINITION OF THE INNER PRODUCT, ADJOINT OPERATORS AND NORM
The inner product on D of two vector functions u and v is defined as where the overbar denotes complex conjugation.This inner product induces a norm given by ||u|| = u, u 1/2 .The adjoints of the operators G B (see Eq. ( 11)) and G E (see Eq. ( 14)) with respect to the above inner product are denoted by G * B and G * E , respectively, and are given by and x ∈D G(x − x)r(x ) dV. (23) Note that the gradient operator in Eq. ( 22) acts on the integral, and not on ρ as described in previous implementations [28], [42], which saves two computationally expensive convolutions.

A. Update of Contrast Source
Keeping the contrast source fixed, the contrast source is updated as w [n] = w [n−1] + α [n] v [n] , (24) where the step length α is given by α with η B = (||B +;sca 1 || 2 ) −1 and η [n] E = (||χ [n−1] E inc || 2 ) −1 .The update direction v are taken to be the Polak-Ribière update directions, given by v [n] = g [n]  w + g [n] w , g [n] w − g Finally, the gradient of the functional with respect to w is given by with χ the complex conjugate of the contrast function.Note that this is the contrast source update that is used in T-CSI and R-CSI.The updates for N-CSI are similar, and are implemented by setting η E to zero in the equations, saving a considerable computational cost.

B. Update of Contrast Function
The contrast function can be determined by minimizing Eq. ( 16).By assuming that the denominator of Eq. ( 16) is independent of χ, the new contrast function can be found by simply minimizing the nominator of the object functional, leading to where the index k ranges over the set {x, y, z}.However, due to the dependency of χ in the denominator of Eq. ( 16), the updating scheme might not reduce errors and it might be better to update the contrast function using the conjugate-gradient update formula [29] χ [n] = χ [n−1] + β [n] d [n] , (29) where the step length β is given by β [n] = −(aC − Ac) + (aC − Ac) 2 − 4(aB − Ab)(bC − Bc) 2(aB − Ab) , (30) with and the update directions for χ are given by χ , g [n] χ − g Finally, the gradient of the functional with respect to χ is given by which is the gradient when η E is again assumed to be independent of χ.

C. Update of Contrast Function in Case of Regularization
In the case of the multiplicative total variation factor, the step length β is given by the real solution of the three roots of the third order polynomial aB + bA + 2(aC + bB + cA)β + 3(bC + cB)β 2 + 4cCβ 3 , (34) with a = F B w [n] + F E w [n] , χ [n−1] , b = 2η [n] E χ [n−1] E [n] − w [n] , d [n] , c = η

Fig. 1 .
Fig. 1.Flowchart of the steps that are taken in the different CSI-EPT methods: naive CSI-EPT (N-CSI), traditional CSI-EPT (T-CSI) and regularized CSI-EPT (R-CSI).The shades of red show an estimation of the most computational expensive steps.Darker red means a higher computational complexity.Note that the computation of E is also a computational expensive step in N-CSI.However, it only has to be computed once.

Fig. 3 .
Fig. 3. Noise sensitivity.The reconstructions of T-CSI-cg and R-CSI after 1000 iterations when noise with an SNR of 75 is implemented (a,g and b,h, respectively), and the reconstructions at different iteration numbers when an SNR of 50 is implemented (c-f and i-l).Conductivity (a-f) and relative permittivity (g-l).

Fig. 4 .
Fig. 4. Relative residual error for the conductivity and permittivity at different iterations for different different CSI-EPT strategies.Noiseless case (a) and with uniform Gaussian noise with an SNR of 50 (b).

Fig. 5 .
Fig. 5. Scaling dependency.Incorrect scaling of the B +1 magnitude map can result in inaccuracies in the CSI reconstruction.The heading represent the scaling factor that is applied on B + 1 , from 20% underestimation (a,f) up to 20% overestimation (e,j).Results are for R-CSI after 10 000 iterations.Conductivity (a-e) and relative permittivity (f-j).

Fig. 7 .
Fig. 7. Domain truncation.Model of the tissue structure of the head (a) and the |B + 1 | field for different sizes of the head model (b-e).The ground truth absolute contrast function (f) and the corresponding |χ| reconstructions when the field of (b) is truncated to different sizes (g-j).The data is noiseless and reconstruction results are after 10 000 iterations with R-CSI.The horizontal lines indicate the position of the end ring of the coil.
and the corresponding quantification values are presented in Table S2 in the Supplementary Information.The data have an SNR of 12, 26,