EPTlib: An Open-Source Extensible Collection of Electric Properties Tomography Techniques

Electric properties tomography (EPT) is a novel magnetic resonance imaging–based method to estimate non-invasively the distribution of the electric properties in the human body. In this paper, EPTlib, an open-source extensible C++ library collecting ready-to-use algorithms for electric properties tomography, is presented. Currently, EPTlib implements three techniques, named Helmholtz-EPT, convection-reaction-EPT and gradient-EPT, whose derivation and implementation is deeply discussed. Moreover, the configuration files needed by the terminal application included in EPTlib to apply the implemented techniques are outlined. The three techniques are applied to a couple of model problems in order to highlight their main features and the effects of the tunable parameters.


Introduction
Traditional magnetic resonance imaging (MRI) produces qualitative images that are weighted according to certain physical and operative parameters. The interpretation of these images is based only on the relative contrast between pixels of the same image and the comparison of the results of acquisitions conducted with different scanners or at different times is not straightforward, because of the lack of a common reference. Quantitative MRI, namely an MRI process that produces images which are pixel-by-pixel measurements of a physical quantity, is arousing increasing interest, because it would overcome the mentioned intercomparison issue. Examples of MRI-based quantitative imaging are magnetic resonance fingerprinting [1], which provides an estimate of the distribution of the relaxation times T 1 and T 2 of the biological tissues, and magnetic resonance elastography [2,3], which evaluates the spatial distribution of the elastic shear modulus within the patient body.
Electric properties tomography (EPT) is a MRI-based quantitative imaging technique which reconstructs the distribution of the electric properties (EPs), conductivity σ and relative permittivity ε r , of the biological tissues at the Larmor frequency of the used MRI scanner [4][5][6][7]. The benefit of a reliable EPT is twofold: on the one hand, EPT can be used as a diagnostic tool, since the electric conductivity of the biological tissues acts as a physical biomarker for some pathologies [8,9]; on the other hand, the knowledge of the actual distribution of the EPs would allow planning electromagnetic-based therapies, like oncological hyperthermia, specifically for each patient [10,11]. A plethora of different techniques implementing EPT has been proposed in the literature [7]. Despite there exist some peculiar approaches to EPT, like the one based on the water content estimated by T 1 -weighted images [12][13][14], almost all the proposed techniques map the EPs by elaborating the distribution of the rotating component B + 1 of the radiofrequency magnetic field B 1 generated by the MRI scanner [7]. The magnitude |B + 1 |, namely the transmit (Tx) sensitivity,

Software Architecture
The software solutions adopted in the development of the library EPTlib and its associated terminal application EPTlib_app are detailed in this section.

Library
EPTlib is an open-source C++ library collecting EPT techniques that includes a terminal application, EPTlib_app, for a ready-to-use adoption of the implemented techniques. The sources and the releases of EPTlib are shared through a GitHub repository, accessible from the website of the project https://eptlib.github.io. Currently, EPTlib implements three EPT techniques, but, by taking advantage of the object-oriented programming paradigm of C++, it can be easily extended with other implementations in the future.
The library has some external dependency, as depicted in Figure 1. In particular, EPTlib depends privately on the interface library Eigen 3 (https://eigen.tuxfamily.org), for linear algebra computations, and publicly on the library HDF 5 (https://www.hdfgroup.org/ solutions/hdf5), for the input and output of datasets, and on the libraries Dynamic Bitset and Operators of the collection Boost (https://www.boost.org), for handling binary maps. Optionally, EPTlib uses the framework GoogleTest (https://google.github.io/googletest) for unit testing. The structure of the library is summarized in Figure 2, where the collaboration diagram of the class implementing Helmholtz-EPT is reported. An identical diagram holds also for the classes implementing convection-reaction-EPT and gradient-EPT. The content of the entire library is protected by the namespace eptlib, that will be omitted when referring to the classes from now on. All the EPT techniques implemented in EPTlib derive from the common abstract class EPTInterface, which provides the elementary methods to handle input and output streams of data and the signature of the method that launches the EPT technique. The abstract class EPTInterface is expected to provide also the post-processing procedures to refine the EPs maps, which can be adopted independently of the used EPT technique. Currently, a median filter is implemented through the class MedianFilter, which reshapes the filter window according to the anatomical geometry deduced from a reference image [9]. The number of input variables can change between different EPT techniques; however, the Tx sensitivity |B + 1 | and the TRx phase ϕ ± (expressed in radians) are common inputs and, as such, are defined directly in EPTInterface. The number of Tx and receive (Rx) channels can change as well from a technique to another, opening the possibility to exploit multiple Tx sensitivities and TRx phases in input. Similarly, although some techniques could provide additional output variables, the electrical conductivity σ and the relative permittivity ε r are the common output of EPT and so are defined in EPTInterface. Both the input and output common variables are stored as three-dimensional real-valued images through the specialized class template Image<double>, whose import and export as three-dimensional datasets in HDF5 files is handled by library routines.
All the EPT techniques currently implemented in EPTlib are differential techniques that require the computation of the derivatives of the input data. This computation is based on a noise-robust second-order finite difference Savitzky-Golay filter [36] implemented in the class FDSavitzkyGolayFilter. The filter approximates within a moving window the quantity to derive with a second-order polynomial function, which is then derived analytically. In case of a wrapped phase map ϕ, the derivative is computed observing that and in order to avoid the appearance of singularities around the discontinuities of 2π. Both the class MedianFilter and the class FDSavitzkyGolayFilter define the filter window by means of an instance of the class Shape, which specializes the features of the Dynamic Bitset of the collection Boost.

Configuration File 1 Common settings for EPTlib_app
The configuration file has a number of settings common to all the EPT techniques and a parameter section that depends on the selected technique. The common settings, collected in the example Configuration File 1, are described in the following. The first three lines of Configuration File 1 are for basic setups: • title and description provide a brief description of the case study.
• method is an integer value that selects the EPT technique to use according to the following list: 0 for Helmholtz-EPT, 1 for convection-reaction-EPT and 2 for gradient-EPT.
The section mesh describes geometrically the input three-dimensional images: • size is the number of voxels of the input data in the three dimensions. • step is the length in meters of the voxels in the three dimensions.
The section input defines the input data: • frequency is the Larmor frequency in hertz of the input data. • tx-channels and rx-channels are the number of Tx and Rx channels in the input data, respectively. They must assume values in certain ranges depending on the adopted EPT technique. • tx-sensitivity and trx-phase are the addresses of the Tx sensitivity |B + 1 | and the TRx phase ϕ ± datasets in suitable HDF5 files. The address is provided as a string in which the HDF5 file address is separated from the internal dataset identifier by a colon. Moreover, in case of multiple Tx channels, the wildcard character '>' can be used to distinguish the data from each channel. It will be replaced by an increasing integer number from zero up to the number of Tx channels minus one. Similarly, the wildcard character '<' is used for multiple Rx channels. Thus, the string "in.h5:/trx_phase><" denotes the datasets "/trx_phase00", "/trx_phase10", "/trx_phase01", "/trx_phase11", "/trx_phase02" and so on in the file "in.h5".
• wrapped-phase is true if the input TRx phase maps are wrapped. This setting is optional and false by default.
The optional section input.wildcard can be used to change the wildcard characters for the input address definition: • tx-character and rx-character are the above-mentioned wildcard character. By default '>' and '<', respectively. • start-from is the integer number from which the wildcard character replacement starts. By default 0. • step is the integer step in the wildcard character replacement. By default it is 1.
Lastly, the section output defines the output addresses: • electric-conductivity and relative-permittivity are the addresses of the datasets in HDF5 files in which the electric conductivity σ and the relative permittivity ε r estimated by the EPT technique will be stored, respectively. The address is provided as a string in which the HDF5 file address is separated from the internal dataset identifier by a colon. Thus, the string "out.h5:/sigma" denotes the dataset "/sigma" in the file "out.h5".
For an up-to-date description of the configuration file settings, check the website of EPTlib, https://eptlib.github.io.

EPT Techniques
Each EPT technique implemented in EPTlib is accurately described in this section, highlighting the underlying mathematical assumptions and the operative conditions under which it can be applied. The specific settings of the configuration file for EPTlib_app are reported as well.

Helmholtz-EPT
Assuming a constant magnetic permeability equal to the one of vacuum µ 0 [37], the time-harmonic RF magnetic flux density B 1 generated by the MRI scanner in the biological tissues satisfies where ω = 2π f is the angular frequency, being f the Larmor frequency, andε = ε r ε 0 − iσ/ω is the complex permittivity, being ε 0 the electric permittivity of vacuum. This equation can be written equivalently as whose left hand side is deduced from the vector calculus identity ∇ × (∇ × B 1 ) = ∇(∇ · B 1 ) −∇ 2 B 1 and from the Gauss's law for magnetism ∇ · B 1 = 0. In the regions where the EPs are almost homogeneous, namely ∇ε −1 ∼ = 0, the latter coincides with the vector Helmholtz equation.
Since the positively rotating component B + 1 of the RF coil is defined as [38] B + 1 = the Cartesian components of the vector Helmholtz equation can be combined linearly to obtain the Helmholtz equation for B + 1 , which inverted algebraically with respect toε leads to the fundamental equation of Helmholtz-EPT [17], Clearly, where the EPs does not satisfy the assumption of homogeneity, like at the tissue boundaries, this equation introduces systematic errors.
By denoting B + 1 = |B + 1 |e iϕ + , an explicit equation for each EP can be deduced. Precisely, and where the Tx phase ϕ + is estimated as half of the measurable TRx phase ϕ ± , when the same RF coil is used both in transmission and in reception with a polarization switch [5].
It is worth noting that if the Rx sensitivity |B − 1 | is almost homogeneous, like in the case of a volume coil used in reception or of an optimized combination of Rx channels [40], the electric conductivity can be estimated by the analogous relation where the complex Rx sensitivity The combination of (10) and (11) leads to a phase-based approximation that uses directly the measurable TRx phase ϕ ± = ϕ + + ϕ − , with no need of estimating ϕ + , Helmholtz-EPT is implemented in EPTlib through the class EPTHelmholtz. Since EPTlib applies a Savitzky-Golay filter [36] to estimate the derivatives of the input maps, the terminal application EPTlib_app requires additional settings in the configuration file to determine the filter window. They belong to the optional section parameter.savitzky-golay and are shown in the example Configuration File 2: • size is the number of voxels along each semi-axis of the filter window (cf. Figure 3).
shape is an integer value that selects the shape of the filter window: 0 for a cross, 1 for an ellipsoid and 2 for a cuboid (cf. Figure 3). By default 0.
Moreover, to apply Helmholtz-EPT it is necessary to select in the configuration file the method 0 and to set tx-channels and rx-channels in the input section equal to 1. If both the input addresses tx-sensitivity and trx-phase are defined in the input section, then the complete Helmholtz-EPT (6) is applied. Otherwise, the magnitude-based approximation (9) or the phase-based approximation (13) is applied, depending on the provided input address. 0 (Cross) 1 (Ellipsoid) 2 (Cuboid)

Convection-Reaction-EPT
Convection-reaction-EPT overcomes the homogeneity assumption of Helmholtz-EPT by elaborating Equation (4) with no assumptions on the distribution of the EPs. The linearly combination of the x and y components of the vector Equation (4) according to the Tx sensitivity (5) leads to the following equation [19], where the variable γ =ε −1 has been introduced and the transverse derivatives of the longitudinal component of B 1 (i.e., ∂B 1,z /∂x and ∂B 1,z /∂y) have been assumed to be negligible. This assumption holds, for instance, around the mid-plane of a birdcage coil, in which case also the longitudinal derivative of B 1,z would be negligible, and so By substituting (15) in (14), the following equation, in which only the measurable component B + 1 of the magnetic flux density appears, is obtained [19], The conservative form of the latter equation is [41] being β + = ∇B + 1 − i∇ × B + 1ẑ the convective term. The latter equation is treated as a complex-valued partial differential equation (PDE) with respect to γ and complemented with Dirichlet boundary conditions, namely the value of γ at the boundary of the computational domain is established a priori. Since, in case of a birdcage coil, Equation (17) holds only around the mid-plane of the coil, where ∇B 1,z ∼ = 0, the determination of the boundary conditions on the top and bottom slices of the region of interest could be tricky. This issue could be avoided by solving the PDE in a single slice with the two-dimensional approximation obtained assuming the longitudinal invariance of the EPs (i.e., ∂γ/∂z ∼ = 0), Similarly to Helmholtz-EPT, also convection-reaction-EPT admits a phase-based approximation [23]. It can be deduced by observing that if the Tx sensitivity |B + 1 | is almost homogeneous, then Equation (17) is approximated by being β + ϕ = ∇ϕ + − i∇ × (ϕ +ẑ ). A similar equation can be deduced for the Rx sensitivity under the assumption of a Rx magnetic flux density whose longitudinal component has negligible derivatives and of an almost homogeneous Rx sensitivity |B − 1 |. It reads . The sum of Equation (19) and (20) leads to where the TRx phase ϕ ± = ϕ + + ϕ − appears. Finally, assuming that σ/ω ε r ε 0 , as it actually is for many biological tissues at 64 MHz (the Larmor frequency of 1.5 T MRI scanners) [18], it is possible to approximate the unknown as γ ∼ = iω/σ. Thus, the imaginary part of (21) reads [23,42] ∇ · ρ∇ϕ ± = 2ωµ 0 , where the resistivity ρ = σ −1 has been introduced. Similarly to Equation (18), also (22) can be approximated in two-dimensions assuming that ∂ρ/∂z ∼ = 0. The PDEs on which convection-reaction-EPT is based are of difficult solution because of the lack of a diffusion term. The introduction of an artificial diffusion term in the equations acts like a viscosity-type regularization of the problem [43] both for the complex-valued Equation (17) and its two-dimensional approximation (18) and for the real-valued phase-based approximation (22). The regularized PDE for the complete convection-reaction-EPT is where λ denotes the artificial diffusion coefficient. An analogous correction is introduced in the other equations.
In EPTlib, the class EPTConvReact manages the convection-reaction-EPT technique. The PDEs are solved numerically according to the finite difference method. The complexvalued PDEs are approximated with second-order centred finite differences, whereas the real-valued PDEs derived by the phase-based approximation are discretized with stable first-order upwind finite differences. In all the cases, the artificial diffusion term is discretized with second-order centred finite differences. The obtained linear system is solved using the biconjugate gradient stabilized method implemented in Eigen 3. volume-tomography = false 13: imaging-slice = 3 14: artificial-diffusion = true 15: artificial-diffusion-coefficient = 0.1

Configuration File 3 Convection-reaction-EPT specific settings for EPTlib_app
The first-and second-order derivatives of the Tx sensitivity and the TRx phase are estimated by EPTlib applying a Savitzky-Golay filter [36]. Thus, EPTlib_app requires the optional section parameter.savitzky-golay in the configuration file, as described above for Helmholtz-EPT. In addition, it requires the settings reported in the example Configuration File 3 for determining the Dirichlet boundary conditions and selecting the proper PDE. Precisely, the optional section parameter.dirichlet set the Dirichlet boundary conditions as constant values in the entire boundary: • electric-conductivity is the electric conductivity in siemens per meter forced at the domain boundary. By default 0. It is worth noting that this condition is used also by the phase-based approximated convection-reaction-EPT (22), in which case the assigned value must be different from zero in order to allow the correct definition of the resistivity ρ. • relative-permittivity is the relative permittivity forced at the domain boundary. By default 1.
The optional section parameter selects the two-dimensional approximation (18) and the artificial diffusion regularization (23): • volume-tomography is false if the two-dimensional approximation is selected. By default false. • imaging-slice is the index of the slice selected for the two-dimensional tomography. This setting is ignored if volume-tomography is true. The mid-plane of the input data is selected by default. If the number of slices is even, then the lower of the two planes in the middle is selected by default. • artificial-diffusion is true if the artificial diffusion regularization is adopted. By default false. • artificial-diffusion-coefficient is the value of the artificial diffusion coefficient λ. This setting is ignored if artificial-diffusion is false. By default 0.
Moreover, to apply convection-reaction-EPT it is necessary to select in the configuration file the method 1 and to set tx-channels and rx-channels in the input section equal to 1. If both the input addresses tx-sensitivity and trx-phase are defined in the input section, then the complete convection-reaction-EPT (17) is applied. Otherwise, if only trx-phase is defined, the phase-based approximation (22) is applied.

Gradient-EPT
Gradient-EPT is developed starting from the same equation of convection-reaction-EPT (17), obtained without the assumption of homogeneity of the EPs used by Helmholtz-EPT. However, the computation is significantly different between the two techniques. In particular, gradient-EPT assumes that a parallel transmission (pTx) system with multiple Tx channels is used. In order to distinguish the quantities associated to each Tx channel, they are denoted in the following by the subscript i with reference to the i-th Tx channel.
Being based on Equation (17), gradient-EPT works under the hypothesis that the gradient of the longitudinal component of B 1,i is negligible, for any channel i. This is verified around the mid-plane of longitudinally symmetric sources, which are usually adopted for volume coils also in pTx systems.
By noting that 1 Equation (16) can be conveniently arranged in the equivalent form [20] ∇B + 1,i · (g + , −ig where g = ∇ log(ε) and g + = g x + ig y have been introduced. In a pTx system, the positively rotating component B + 1,i cannot be estimated easily, since the TRx phase does not approximates, in general, the Tx phase. However, assuming that the same Rx coil, or combination of Rx coils, has been used to acquire the signal of each Tx channel, the measured TRx phases ϕ ± i = ϕ + i + ϕ − can be combined to estimate the relative Tx phases ϕ +,r i with respect to a reference Tx phase ϕ + 0 , Thus, the positively rotating component of the i-th channel can be written as being B +,r 1,i = |B + 1,i |e iϕ +,r i the measurable part and ϕ + 0 the unknown reference phase. Substituting Equation (27) in (25) and cancelling out the common factor e iϕ + 0 , the following equation is obtained [20], −i2∇B +,r 1,i · ∇ϕ + 0 + ∇B +,r 1,i · (g + , −ig + , g z ) + B +,r 1,i θ + = ∇ 2 B +,r 1,i , where θ + = −ω 2 µ 0ε + |∇ϕ + 0 | 2 − i∇ 2 ϕ + 0 + i∇ϕ + 0 · (g + , −ig + , g z ) combines the EPs with the other unknown terms. A system of linear equations in the algebraically independent unknowns ∇ϕ + 0 , g + , g z and θ + is obtained in each voxel by writing Equation (28) for each Tx channel i and solved in the least-squares sense. In order to reduce the presence of numerical artefacts, a linear system is solved using each Tx channel as the reference one; the results are then averaged weighted by the Tx sensitivity magnitude of the corresponding reference channel [20]. An estimate of the EPs can be deduced by combining the unknown terms.
The solution of (28) is the local step of gradient-EPT, which can be followed by a noise-filtering global step, consisting in the solution in the computational domain Ω of the minimization problem [20] where g +,0 and g z,0 denotes the solutions of the local step in each voxel, whilst g + (ε) and g z (ε) are the derivatives of the guess distributionε. Considering u = log(ε) as the unknown, problem (29) Since it involves only the derivatives of u, this equation has infinitely many solutions amongst which the correct one must be selected by means of a constraint or a regularization term. In the former case, a number of seed points in which the EPs are constrained to assume certain values is introduced [20]. Numerically, this operation is analogous to setting Dirichlet boundary conditions, although in internal isolated points. In the latter case, a regularization term is added to the minimization problem, which becomes [44] miñ ε g + (ε) − g +,0 with λ the regularization coefficient,ε 0 the solution of the local step and Ω 0 ⊂ Ω the sub-domain where the regularization is applied. The Euler equation associated to the regularized minimization problem (31) for a generic test function v is being the omitted parts equal to (30). Since the quality ofε 0 is higher in the homogeneous regions, the sub-domain Ω 0 is defined as the union of the voxels where the EPs are almost homogeneous. The region is identified based on the distribution of the gradient magnitude G = |g + | 2 + |g z | 2 . Precisely, given a tolerance k, Ω 0 is the union of all and only the voxels where G < k max(G). Other definitions of Ω 0 can be found in the literature, like for example definitions based on its extension with respect to the entire computational domain Ω [45].
Because of the large number of small linear systems to be solved in the local step and of the large linear system in the global step, gradient-EPT is computationally the most expensive of the presented techniques. The computational burden can be reduced with a two-dimensional approximation of the equations, which can be obtained under the assumption of the longitudinal invariance of the EPs (i.e., g z ∼ = 0) and of the reference phase (i.e., ∂ϕ + 0 /∂z ∼ = 0). In that case, the equation for the local step (28) becomes −i2∇ xy B +,r 1,i · ∇ xy ϕ + 0 + ∇ xy B +,r 1,i · (g + , −ig + ) + B +,r 1,i θ + xy = ∇ 2 B +,r 1,i , with θ + xy = −ω 2 µ 0ε + |∇ xy ϕ + 0 | 2 − i∇ 2 xy ϕ + 0 + i∇ xy ϕ + 0 · (g + , −ig + ), whereas the global step keeps the same minimization problem (29) without the term relative to g z .

Configuration File 4 Gradient-EPT specific settings for
The example Configuration File 4 collects the additional settings of EPTlib_app specific for gradient-EPT. As already described for Helmholtz-EPT, the optional section parameter.savitzky-golay select the window of the Savitzky-Golay filter [36] adopted for the estimation of the first-and second-order derivatives of the Tx sensitivities and the relative phases in the local step. The optional section parameter.seed-point sets a seed point constraint in the global step (29): • use-seed-point is true if the seed point constraint is selected. By default false. • coordinates is the list of the integer coordinates of the seed points in the mesh of voxels (index 0 is used for the first voxel). This setting is mandatory if use-seed-point is true, otherwise it is ignored. • electric-conductivity is the list of the electric conductivity in siemens per meter forced in each seed point. This setting is mandatory if use-seed-point is true, otherwise it is ignored. • relative-permittivity is the list of the relative permittivity forced in each seed point. This setting is mandatory if use-seed-point is true, otherwise it is ignored.
The lists of coordinates, electric-conductivity and relative-permittivity must have as many elements as the number of selected seed points. The optional section parameter.regularization sets the regularization term in the global step (31) if the seed point constraint is not used (this section is ignored if use-seed-point in section parameter.seed-point is true): • regularization-coefficient is the value of the regularization coefficient λ. By default 1. • gradient-tolerance is the value of the tolerance k for determining the sub-domain Ω 0 . By default 0. • output-mask is the address of the dataset in HDF5 file in which the characteristic function of Ω 0 will be stored. The address is provided as a string in which the HDF5 file address is separated from the internal dataset identifier by a colon. Thus, the string "out.h5:/mask" denotes the dataset "/mask" in the file "out.h5". If this setting is omitted, the data will not be stored.
The optional section parameter selects the two-dimensional approximation (33) and if the execution must end after the local step (28) or proceed with the global-step (29): • volume-tomography is false if the two-dimensional approximation is selected. By default false. • imaging-slice is the index of the slice selected for the two-dimensional tomography. This setting is ignored if volume-tomography is true. The mid-plane of the input data is selected by default. • full-run is false if the execution must stop after the local step. By default true.
Moreover, to apply gradient-EPT it is necessary to select in the configuration file the method 2 and to set rx-channels in the input section equal to 1. Since in the local step (28) there are nine real-valued unknowns and each Tx channel leads to a complex-valued equation, a pTx system with at least five Tx channels is needed by gradient-EPT. Thus, the setting tx-channels in the input section must be at least equal to 5. Table 1 summarizes the input requirements and the output capabilities of the EPT techniques implemented in EPTlib, distinguishing between possible phase-based and magnitude-based approximations. Moreover, it is worth recalling the main hypotheses on which the three techniques are based:

1.
Helmholtz-EPT assumes that the EPs are locally homogeneous, leading to large systematic errors at the tissue boundaries.

2.
The complete versions of both Helmholtz-EPT and convection-reaction-EPT estimate the Tx phase ϕ + as half of the measurable TRx phase ϕ ± . This is a reasonable estimation when the same RF coil is used both in transmission and in reception with a polarization switch [5], or when special care is taken in the phase acquisition, like for example adopting a CLEAR reconstruction [46].

3.
Both convection-reaction-EPT and gradient-EPT assume that the longitudinal component of B 1 has a negligible gradient. This assumption limits the applicability of these techniques around the mid-plane of longitudinally symmetric sources, which are usually adopted for volume coils.

4.
Gradient-EPT requires acquisitions performed with a pTx system with at least five Tx channels. Currently, this kind of equipment is mainly employed in the research of ultra high field scanners.

Examples of Application
The implemented Helmholtz-EPT and convection-reaction-EPT techniques were applied on virtual data obtained from the simulation of a 1.5 T low-pass 16-leg birdcage body-coil (radius = 36 cm, height = 45 cm). The body-coil was used both in transmission and in reception with polarization switch. The gradient-EPT technique, instead, was applied on data from the simulation of a 7 T 8-channel head-coil made of eight loops arranged elliptically [47] (height = 16 cm, semi-axes = 11.2 cm and 13.2 cm). The same head-coil was used in transmission and reception and the raw data acquired by the eight channels were combined according to a circular polarization. Both the coil models were used to perform head imaging of the anatomical human model Duke [48], as illustrated in Figure 4. In order to test the EPT techniques in realistic conditions, the simulated sensitivities were sampled on an isotropic Cartesian mesh with sides of 2 mm and were corrupted with additive spatially-uncorrelated Gaussian noise, with a uniform standard deviation and an average signal-to-noise ratio (SNR) of 200 and 100. For the 1.5 T case, the simulated distributions of Tx sensitivity |B + 1 | and TRx phase ϕ ± are shown in Figure 5, as well as the noisy samples. The electric conductivity σ recovered in the mid-plane by applying Helmholtz-EPT on the virtual data generated by the 1.5 T body-coil is reported in Figure 6, where also the actual distribution of σ is shown for comparison. The results were collected for different SNR in the input data and for windows of the Savitzky-Golay filter with various shapes and sizes. The noiseless case (Figure 6b) is elaborated with a minimal filter-the cross of size [1,1,1]-corresponding to the traditional second-order centred finite difference scheme. The large errors introduced at tissue boundaries (including those developed longitudinally and invisible in the two-dimensional map) were due to the assumption of homogeneity of the EPs, on which Helmholtz-EPT was based. The presence of noise in the case with SNR 200 (Figure 6c) forced us to adopt a wider window for the Savitzky-Golay filter. Precisely, a cuboid of size [2,2,2] was used. This led to a blurred reconstruction, which managed the input noise at the cost of a lower quality. With an amount of noise such that the input SNR was 100, the size of the filter window was further increased to [3,3,3]. In this case, the effect of different window shapes was analyzed by comparing the cross (Figure 6d), the ellipsoid ( Figure 6e) and the cuboid (Figure 6f). The cuboid, being the shape that relied on the largest number of voxels, had the strongest damping effect on the input noise, but, at the same time, reduced the effective image resolution by blurring many details.
The considered problem satisfied with large accuracy the assumption of homogeneity of the Tx and the Rx sensitivities. On the one hand, this allowed us to obtain almost identical reconstructions as those reported in Figure 6 if the phase-based approximated Helmholtz-EPT technique was adopted instead of the complete one. On the other hand, the estimation of the relative permittivity ε r was significantly affected by the Laplacian of the Tx sensitivity (7), whose small variations were easily suppressed by the noise. As a consequence, the distributions of ε r recovered by applying Helmholtz-EPT with the adopted filter windows on the noisy input data were extremely noisy and completely unreadable.

Convection-Reaction-EPT
The convection-reaction-EPT technique applied on the complex input data of the simulated 1.5 T body-coil led to extremely ill-conditioned linear systems unless an overregularization was introduced, making the method unfeasible even in the noiseless case. The results reported in Figure 7 are obtained by the phase-based approximation of convection-reaction-EPT with artificial diffusion regularization applied on a three-dimensional portion of data. Precisely, 21 slices (10 above and 10 below the mid-plane) were considered and found to be enough to reduce the influence of the boundary conditions on the top and bottom slices on the reconstruction in the mid-plane. Figure 7 collects the distribution of σ recovered for different SNR in the input data and different values of the artificial diffusion coefficient λ. All the results were obtained using the minimal Savitzky-Golay filter, namely the cross of size [1,1,1]. The reference distribution of σ was the one reported in Figure 6a. The reconstructions obtained operating on noiseless data show that a negligible loss of details occurred when the coefficient λ was increased from 0.001 up to 0.003. In particular, a small reduction in the value of the cerebrospinal fluid (CSF) conductivity was observed, with no sensible consequences on the effective image resolution. Moreover, the distributions of σ obtained by convection-reaction-EPT with noiseless input (cf. Figure 7) were significantly improved with respect to the distribution estimated by Helmholtz-EPT (cf. Figure 6b), especially because of the absence of large errors at the tissue boundaries.
When the input data were noisy, the conditioning of the problem at the basis of convection-reaction-EPT worsened, making its resolution more time-consuming and affecting the quality of the result. However, the artificial diffusion term restored the conditioning with a regularization effect that became more relevant when the coefficient λ increased, as can be observed in Figure 7 for the cases with SNR 200 and, more evidently, with SNR 100. It is worth noting that, differently from the increased windows of the Savitzky-Golay filter used for Helmholtz-EPT, the artificial diffusion term managed the noise propagation without affecting significantly the image resolution. The adoption of a large window for the Savitzky-Golay filter in conjunction with the artificial diffusion regularization could help when dealing with very noisy input data.

7 T 8-Channel Head-Coil
The virtual data generated by simulating the 7 T 8-channel head-coil were used to test the gradient-EPT technique. A three-dimensional portion of data of 11 slices (five above and five below the mid-plane) was used. The obtained distributions of the electric conductivity σ are collected in Figure 8, where the actual distribution of σ at 300 MHz (the Larmor frequency for 7 T MRI) is reported as well for comparison. The results were collected for different SNR in the input data and were obtained by solving the global step of gradient-EPT with the seed point constraint or with the regularization term. In the former case, a cross of size [2,2,2] was used as window for the Savitzky-Golay filter for any input SNR, whereas in the latter a cuboid of size [2,2,2] was adopted to obtain reasonable outcomes. Seven seed points distributed in the white matter and located in mid-plane were used for the seed point constrained reconstructions and are pictured as red dots on the recovered maps in Figure 8. The actual EPs of the white matter adopted in the simulations were assigned to the seed points. Their adoption led to an extremely detailed reconstruction of the distribution of σ, even when noisy input data with SNR 100 were used. Regarding the reconstructions obtained with the regularization term, a regularization coefficient λ = 10,000 was used and the regularization mask was determined with a tolerance k = 0.02 for the noiseless input and k = 0.03 for the noisy input. The resulting masks and the recovered distributions of σ, both reported in Figure 8, were very robust with respect to the input noise, although the adoption of a larger window for the Savitzky-Golay filter reduced the effective resolution of the maps with respect to the one obtained with the seed point constraint.
The results obtained by the two-dimensional approximation of gradient-EPT are collected in Figure 9. Besides reducing the computational burden, the two-dimensional approximation simplified the solution of the global step with the regularization term, which in this case was performed, as for the reconstructions obtained with the seed point constraint, using a cross window of size [2,2,2] for the Savitzky-Golay filter. This made the regularized reconstruction more detailed than in the three-dimensional case, although the effect of noise was more relevant. The same seed points used in the three-dimensional case were used also in the two-dimensional approximation with seed point constraint, which led to a reconstruction robust with respect to noise, despite the systematic errors due to the two-dimensional approximation. Figure 9. Electric conductivity distributions recovered by the two-dimensional approximation of gradient-EPT, to be compared with the actual distribution at 300 MHz reported in Figure 8. All the distributions are recovered using a cross window of size [2,2,2] for the Savitzky-Golay filter. The seed points selected for the constraint are marked by red dots. The distributions recovered with regularization used a regularization coefficient λ = 1 000. The masks adopted for regularization were obtained with a tolerance k = 0.02 in the noiseless case, k = 0.03 with SNR 200 and k = 0.05 with SNR 100. The resulting masks are reported. Results obtained with different amount of noise (noiseless, SNR = 200 and SNR = 100) were collected.

Conclusions
In this paper, EPTlib, a novel open-source C++ library collecting EPT techniques, has been presented. EPTlib provides a common platform for the implementation and diffusion of EPT techniques at the state of the art, and makes them available and easily accessible to the scientists interested in their applications. The collected techniques are accessible in a ready-to-use form in the terminal application included in EPTlib, which is deeply described in this paper. Updates on the library and its terminal application can be found on the library website https://eptlib.github.io, where examples of usage and datasets are made available as well.
Currently, EPTlib implements three EPT techniques, whose performances when applied on a couple of model problems have been reported to highlight their features and the role of the parameters that can be tuned in the terminal application. Helmholtz-EPT is affected by large errors at the tissue boundaries and, in presence of noise, it requires the adoption of large windows for the Savitzky-Golay filter, reducing the effective resolution of the reconstructed maps. Convection-reaction-EPT, in its phase-based approximated version, provides smoother results than Helmholtz-EPT and reaches higher effective resolutions by handling the noisy input with an artificial diffusion regularization term. Finally, gradient-EPT significantly improves the results of the other techniques by relying on multiple input data. The parameters of each technique should be tuned depending on the noise distribution of the input data. In general, a higher SNR leads to a better effective resolution of the conductivity map, since it allows the adoption of smaller filter windows or lower regularization coefficients. The SNR of the input data could be improved, for instance, by averaging many acquisitions, adopting proper sequences or relying on high field scanners. Despite the differences between the techniques, all the reconstructions are performed in acceptable computational times, since in a laptop PC equipped with an Intel Core i7 at 1.80 GHz they took at most 30 seconds to solve the model problems.
In the future, EPTlib could become a reference platform for EPT development, helping in its widespread adoption and growing with the implementation of new techniques and the improvement of the ones already implemented.

Funding:
The results here presented have been developed in the framework of the EMPIR Project 18HLT05 QUIERO. This project 18HLT05 QUIERO has received funding from the EMPIR programme co-financed by the Participating States and from the European Union's Horizon 2020 research and innovation programme.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
Publicly available open-source software were used in this study. This software can befound here: https://eptlib.github.io.