Inversion without Explicit Jacobian Calculations in Electrical Impedance Tomography

Electrical impedance tomography (EIT) is the inverse problem of finding the internal conductivity distribution of a medium given boundary electrical measurements performed via an electrode array onto its surface. Conventional inversion schemes adopt Tikhonov regularized Newton-type methods. Following a transport back-transport approach, we develop in this work an adjoint approach which allows reducing computational burden in enabling inversion without explicit Jacobian calculation. Forward and back-projection operators are defined from potential gradients, along with their efficient implementation. These derivations allow the transparent use of inversion algorithms. We first check the implementation of operators. We then evaluate how reconstructions perform on simulated noisy data using a preconditioned conjugate gradient. We eventually practice our inversion framework on experimental data acquired in vitro from a saline phantom.


Introduction
Biological tissues and fluids exhibit specific electrical properties, conductivity and permittivity , which may be used to infer their internal structure. Medical Electrical Impedance Tomography (EIT) is a soft-field, non-invasive imaging modality, which aims at reconstructing the distribution of internal electrical properties from boundary electrical measurements [1]. It works by injecting known low intensity alternating electrical currents through the body via an array of surface electrodes. The distribution of admittivity , at frequency , dictates internal current flows. Resulting voltages are measured between several pairs of electrodes.
EIT uses low frequency electrical currents from 1Hz to 10MHz. In this frequency range, current flows cannot be confined to a plane by a set of external electrodes; all measurements may be affected by a change in conductivity anywhere in the domain, not only those in the ray path. Hence, electrical imaging requires three dimensional modelling and needs adequate computational processing [2].
Using electrical measurements, reconstruction of admittivity distribution involves a two-step mathematical framework. The forward problem consists in building a numerical model able to predict expected measurements at the electrodes, given a distribution of admittivity. The inverse problem then estimates internal electrical properties through an optimization function which adjusts admittivity distribution in order to minimise the distance between measurements and calculated values.
Because the EIT inverse problem is non-linear and ill-posed, it cannot be easily solved. The standard reconstruction approach uses a family of regularized Newton-type methods [2] where admittivity distribution is represented using a piecewise constant model over voxels. Forward problem is linearized at some initial admittivity estimate and a Jacobian matrix is constructed. This matrix, also called sensitivity matrix, relates the admittivity of each internal voxel to the boundary measurements. Then, admittivity reconstruction involves the inversion of the ill-conditioned Jacobian matrix through generalized Tikhonov regularization. This requires inverting large linear systems of size equal to the number of voxels in the mesh times the number of measurements.
In order to optimise computational efficiency by avoiding explicit calculations of Jacobian, we propose here a transport back-transport method [3] that defines an adjoint reconstruction framework for EIT. Matrices involved in transport and back-transport operations present a reduced size by a factor of the number of electrodes, typically 16 or more in modern EIT systems. This approach opens new opportunities to use a wide range of inversion algorithms for EIT. It also offers a way to reduce complexity of EIT calculations.
In this paper, we focus only on conductivity reconstructions. The derivation of the forward and back-projection operators is presented. The reconstruction framework is then evaluated using the case of a 2D preconditioned conjugate gradient (PCG) inversion scheme. Conductivity reconstructions are performed against experimental data acquired in vitro on a saline phantom.

Sensitivity calculations
Optimization based methods are used in EIT to infer internal electrical property distribution. They require determining the derivative of the voltage measurements with respect to a conductivity parameter. The Jacobian matrix, or sensitivity matrix, contains all these partial derivatives.
The classic estimation of the Jacobian is based on the perturbation approach where two configurations are considered ( Figure 1). The first configuration indexed by is the actual measurement configuration, in which the source is used. The second configuration indexed by is a virtual measurement configuration in which source and detector have been interchanged.

Figure 1. Actual (left) and virtual (right) measurement configurations for source and detector
Suppose the domain under study has been meshed into a partition of elements. The perturbation of the voltage , measured in the situation where the source is and detector , due to a localized perturbation of conductivity in the element , for a fixed injected current of amplitude is [4]: in which ( ) ( ) are the potential gradients, respectively in actual and virtual configurations. With electrodes, when rotating detector, we usually count E measurements per injection configuration: considering all measurements, the Jacobian matrix has dimension . Jacobian coefficients are calculated considering a piecewise linear potential, which leads to constant potential gradients per element. Hence, equation (1) reduces to: where | | is the volume of element in 3D. Potential gradients over each element of the mesh can be computed from standard finite element derivation of EIT physics [5].

Adjoint reconstruction framework
In this work, we follow a transport back-transport approach [3] which uses adjoint fields for reconstruction. Starting from an initial guess of conductivity, the framework calculates the difference between the computed and the actual measurements by solving a direct transport problem, and then transports these residuals back into the medium by solving a corresponding adjoint transport problem. In order to avoid explicit calculations of the Jacobian and its adjoint, the Jacobian matrix may be factorized under the form of equation (3) with matrices , given in equation (4), containing potential gradients relative to coordinate , respectively for source injection configurations and virtual detector injection configurations.
This factorization leads to express the direct and adjoint transport operators without explicit calculations of Jacobians but only of matrices and , possibly equal.

Direct transport: forward operator
The direct transport problem predicts expected measurements given a map of conductivity. It performs the operation . From the perturbation formula of equation (2) under the assumption of piecewise linear potential, the operation can be written in the formalism which uses gradient matrices, with :

Adjoint transport: back-projection operator
The adjoint transport problem transports back to the image the residuals between computed and measured values. It performs the adjoint operation ( ). Using gradient matrices, it follows:

Implementation of inversion
The adjoint reconstruction framework depicted in the previous section allows using numerous inversion algorithms. In this work, standard PCG algorithm is used. CG features intrinsic regularizing properties and offers a good compromise between robustness, convergence, computation and image quality. Left preconditioning is used to weight measurements by the noise covariance matrix. Numerical simulations are conducted using standard Finite Element Method in EIT with the Complete Electrode Model [1]. Custom Matlab routines adapted from the EIDORS library [6] allow the calculations of potential gradients along each mesh element and the implementations of forward and back-projection operators. Experimental measurements are performed on a saline phantom of 4cm diameter featuring 14 equallyspaced copper electrodes ( Figure 2) with a custom-built EIT system [7]. Measurements are carried out using a four-electrode method with an injected current of 100µApp at 98kHz. Inhomogeneities are created by the presence of metallic cylinders of either 1cm or 5mm diameter, i.e. ⁄ or ⁄ of phantom diameter. A 5mm diameter hollow insulating cylinder is also used.

Reconstruction approach
Numerical derivations depicted in previous sections are first validated on simple test cases in 2D before considering image reconstruction of simulated and experimental measurements. A 2D mesh consisting of elements will be used throughout the rest of this paper.

Implementation validation
Two reference configurations, similar to the ones in Figure 2, are used to validate both forward and back-projection implementations. Background conductivity is 1S/m. Inclusions either centred or decentred exhibit a conductivity contrast of 0.1. Using the EIDORS library [6], the Jacobian is calculated along with the residuals, obtained by inversion of the FEM admittance matrix. Proposed operators can be checked by comparing their results with the application of or .

Simulated noisy measurements
Difference reconstructions are conducted on simulated numerical measurements, perturbed with a 5% additive Gaussian noise on every measurement. These reconstructions allow studying the reconstruction framework behavior in presence of noise. The mesh used for the simulation of measurements (18,647 elements) is different from the mesh used for reconstructions (6,713 elements), in order to avoid committing an inverse crime [2].

Experimental measurements
Difference reconstructions are performed given sets of data acquired from the experimental device. Data collection is performed 10 times in order to evaluate the mean value of measurements and their variance. The latter is incorporated in the left preconditioner.

Validation of reconstruction framework implementation
Comparison between calculated gradients and theoretical ones, for a given potential distribution ( ) ( ) , shows a mean value of 1 and a standard deviation of 0.05. This indicates that the implementation of elemental gradient calculations works fine. Both projection and back-projection ratios exhibit no discrepancies (mean value of 1 and corresponding standard deviation inferior to ), which validates the right implementation of corresponding operators.

2D simulation results
Difference reconstructions are performed against simulated measurements (Figure 3). be recovered, it is a classical behaviour of EIT [1]. In presence of noisy data, the PCG algorithm is able to detect the inclusions along with their conductivity contrast. But, spatial resolution has been lost in the reconstruction process, which is also a known feature of EIT [1]. Future developments of the adjoint framework presented in this work might consider inversion algorithms that favour sparse solutions.
The back-projection operator readily offers determination of sensitivity maps for any desired source-detector configuration (Figure 4). These maps present both positive sensitivity areas (yellow, red) and negative ones (blue). The latter explains the lack of spatial resolution put forward in previous reconstruction in the presence of noise.

Phantom experiments
Difference reconstructions from experimental data are now performed. Experiments are first led with conductive inclusion ( ⁄ of phantom diameter) featuring a contrast of with the saline solution background ( Figure 5). The PCG algorithm and subsequent forward and back-projection operators perform well in experimental conditions, with an estimated maximum noise of 5% in the corresponding measurements. Reconstruction drift in position seems to be handled quite properly.
Experiments are then conducted with two conductive inclusions ( ⁄ and ⁄ of phantom diameter) or with a hollow insulating cylinder ( ⁄ of phantom diameter) ( Figure 6).  The framework behaves properly on both conductivity distributions. Recovering an inclusion whose dimensions equal ⁄ of phantom diameter almost matches the spatial resolution limit of 14 electrode EIT.
Sharp contrasts are not recovered in this work since it requires using non-linear inversion algorithms.

Discussion
In this paper, we presented an inversion scheme which does not require explicit calculations of Jacobian. Instead, we practice a transport back-transport framework using an adjoint method and define two operators: forward and back-projection. Both can be efficiently determined from elemental gradient matrices detailed in previous sections. Such an approach allows using a wide range of inversion algorithms in the reconstruction process of EIT.
The validity of the implementation we propose was first checked. We then ran 2D difference reconstructions on simulated numerical data in which a 5% additive Gaussian noise was added. Classical behaviours of EIT reconstruction were observed, namely the loss of spatial resolution with noisy measurements. Eventually, we tested the method on experimental data acquired in vitro on a saline phantom.
In this scope, the framework allows seamlessly to determine estimated conductivity profiles and measurement configuration sensitivity maps. With the 14 electrode phantom used in this work, gains in matrices size used in calculations involve a reduction of 14 in one dimension, i.e. manipulated matrices were of size 6,713 14 instead of standard 6,713 14².
During the inversion process, the PCG algorithm easily allows incorporating prior information about the solution. In this work, measurement variance was incorporated in the left preconditioner. Consequent images prove to be qualitatively superior to the ones obtained using conventional regularization schemes, i.e. generalized Tikhonov. This aspect deserves more investigation.
Further studies will translate these developments in 3D and work on recovering more complicated conductivity distributions. Future developments of the adjoint framework presented in this work might consider reconstruction processes that favour sparse solutions and use non-linear inversion algorithms.