Correction of electrode modelling errors in multi-frequency EIT imaging

The differentiation of haemorrhagic from ischaemic stroke using electrical impedance tomography (EIT) requires measurements at multiple frequencies, since the general lack of healthy measurements on the same patient excludes time-difference imaging methods. It has previously been shown that the inaccurate modelling of electrodes constitutes one of the largest sources of image artefacts in non-linear multi-frequency EIT applications. To address this issue, we augmented the conductivity Jacobian matrix with a Jacobian matrix with respect to electrode movement. Using this new algorithm, simulated ischaemic and haemorrhagic strokes in a realistic head model were reconstructed for varying degrees of electrode position errors. The simultaneous recovery of conductivity spectra and electrode positions removed most artefacts caused by inaccurately modelled electrodes. Reconstructions were stable for electrode position errors of up to 1.5 mm standard deviation along both surface dimensions. We conclude that this method can be used for electrode model correction in multi-frequency EIT.


Background
Multifrequency electrical impedance tomography (MFEIT) is a method for imaging biological tissues with frequency-dependent conductivity. While time-difference (TD) EIT is used to image conductivity changes between a baseline voltage measurement and a data measurement, MFEIT differentiates tissues based on their conductivity spectra at different modulation frequencies of the applied current. Therefore MFEIT does not require a baseline measurement and can be used as a diagnostic tool for multiple applications , Malich et al 2003, including stroke type differentiation (Holder 1992, Romsauerova et al 2006, Packham et al 2012.
While cerebral haemorrhages require surgery, ischaemic strokes could be treated with a clot dissolving drug if diagnosed within 3 h of onset. The current diagnostic procedure is to take a computed tomography (CT) scan, and results in only 2.5-6% of the 80% of ischaemic strokes to be treated in time (Power 2004, Saver et al 2013. While EIT cannot compete with CT in terms of image quality, the small size and low cost of an EIT system make it feasible to equip ambulances for early ischaemia diagnosis and thrombolytic treatment. A novel non-linear method for performing MFEIT using spectral constraints was recently proposed by Malone et al (2014a), because the complicated geometry of the head excludes linear multi-frequency reconstruction algorithms, such as weighted frequency difference (Jun et al 2009). In the new MFEIT algorithm, the inverse problem is reformulated to express the conductivity at each frequency in terms of tissue volume fractions and known tissue conductivity spectra. This separation of variables makes it possible to use measurements at different frequencies simultaneously, since the tissue fractions are frequency independent. An analysis of the influence of modelling errors on images reconstructed with this method has shown that inaccurately modelled electrode positions strongly affect the image quality, whereas wrongly modelled contact impedances were negligible and spectral errors were tolerable as long as tissue spectra did not overlap (Malone et al 2014b).
Methods for correction of electrode modelling inaccuracies are recently gaining interest in the EIT community (Soleimani et al 2006, Dardé et al 2012, and have already been applied to a realistic three dimensional head model for simultaneous TD recovery of conductivity changes and electrode movement (Jehl et al 2015b). In this paper, the first application of simultaneous MFEIT image reconstruction and electrode model correction is demonstrated on a realistic 3D head model with skull and scalp. The implementation is discussed and the performance is evaluated for different levels of electrode position errors.

Purpose
The purpose of this study is to evaluate the performance of simultaneous recovery of electrode positions and conductivity spectrum changes. Specifically, the following two questions will be addressed: (i) Does the simultaneous recovery of conductivity changes and electrode modelling errors remove image artefacts caused by inaccurately modelled electrode positions? (ii) At which magnitude of electrode position errors does the proposed algorithm begin to fail?
To answer these questions, multi-frequency boundary voltages were simulated on a fine 5 million element head mesh with different levels of electrode position errors. Reconstructions were made with and without the proposed addition of electrode modelling corrections and the resulting images were compared. It was found that the proposed algorithm could stably recover simulated strokes in the presence of electrode modelling errors of up to 1.5 mm standard deviation.

Tissue fraction reconstruction with electrode position correction
The general structure of the used multi-frequency reconstruction algorithm was identical to the one published by Malone et al (2014a). To be able to use conductivity measurements at different frequencies ω = … i W , 1, , i simultaneously in a single image reconstruction, the conductivity spectrum of one element in the mesh was described as a linear combination of the known spectra of present tissues = … t j T , 1, , j . Instead of reconstructing conductivities directly, the prior knowledge of the conductivity spectrum ε ij of all tissues was therefore used to assign fractions f nj of these tissues to all finite elements = … n N 1, , , such that . The modified Jacobian matrix at each frequency was obtained from the non-linear EIT forward map A using the chain rule where = … r R 1, , are the lines in the protocol, i.e. the combinations of different current injection and voltage measurement electrode pairs d and m. The Jacobian matrix ( ) σ J i relating voltage changes to changes in conductivity σ i at each frequency ω i , was computed using the adjoint fields method (derived e.g. in the appendix of Polydorides and Lionheart (2002)), giving one entry for element n as is the electric potential emerging when the drive current I d is applied to the electrodes and ( ) ∈ Ω u H m 1 the electric potential when a unit current is applied to the two measurement electrodes.
In order to correct for wrongly modelled electrode positions, this 'traditional' Jacobian matrix was augmented by an electrode boundary Jacobian ∈ × R EBJ R E , relating electrode boundary changes (in our case electrode movement) to voltage changes. Given a continuous vector field v on the boundary of electrode = … e E 1, , , one entry of the EBJ can be computed similarly to the conductivity Jacobian (Dardé et al 2012, Jehl et where z e is the contact impedance, U e d and U e m the drive and measurement electrode potentials and ∂ n E the outward normal of the electrode boundary, tangential to the head surface. The vector field v describes the studied change in the electrode boundary, e.g. to compute the EBJ with respect to movement along one direction, the vector field was chosen to point homogeneously in this direction. The Jacobian matrices were combined and the unknown electrode position errors p were appended to the vector of the tissue fractions to be recovered and p consisted of two variables per electrode that described electrode movements along both surface directions x s and y s , The regularised objective function to be minimised was then with ( ) ω v i being the boundary voltages measured at frequency ω i and the regularisation term The regularisation matrix D comprised one Laplacian matrix per recovered tissue and one identity matrix for the electrode movement variables. All components were scaled according to the expected standard deviation of the corresponding variable changes = std 0.01 The minimisation of the objective function (5) was performed by alternating steps of gradient projection) and damped Gauss-Newton algorithms. The gradient projection (Nocedal and Wright 1999) step was used to quickly move to the neighbourhood of the minimum, while considering the constraints on the fractions. This was done by computing the step sizes along the gradient, at which one of the fraction values reaches a constraint, i.e. 0 or 1. The change of the objective function value along each of the resulting intervals was approximated quadratically using Taylor series. Once one Taylor approximation found a minimum on an interval, this so-called Cauchy point x c was chosen, otherwise the gradient projection algorithm continued on the next interval until a minimum was found.
The subsequent Gauss-Newton step was only applied to the components that did not reach a constraint during the gradient projection. The search direction was calculated by solving using a generalised minimal residual algorithm in order to avoid the explicit calculation of the Hessian matrix H, which is the second derivative of ( ) Φ x (which was approximated by disregarding the second order derivative of the residual error). The step width along direction d was determined using the Brent line-search method (Brent 1973) and the resulting minimum x g was projected back to the fraction constraints. The point, x c or x g , that gave a smaller function value was chosen for the next iteration.
The reconstruction of the fractions was constrained to the closed interval [0, 1] and the constraint ∑ = = f 1 j T nj 1 was enforced by substituting = − ∑ = f f 1 j T j 1 2 . After two iterations of gradient projection and Gauss-Newton, the electrode positions had normally converged and were subsequently kept fixed for the remaining iterations. The number of iterations of this reconstruction algorithm was fixed to 10 for all image reconstructions. To avoid the inverse crime (Lionheart 2004) and speed up image reconstruction, all reconstruction were made on a coarse 180 thousand element mesh on which the skull and scalp were kept fixed and it was assumed the inside of the skull was occupied by either the brain or the stroke with the initial guess being the healthy brain. The regularisation parameter of τ = × − 8 10 10 was chosen empirically and was the same for all reconstructed images. After each iteration of gradient projection and Gauss-Newton minimisation, the regularisation parameter was halved in the case of ischaemias, and given that the spectral contrast was lower, divided by three for haemorrhages (Viklands andGulliksson 2001, Malone et al 2014b).

Data simulation
The simulation parameters in this study were identical to the ones used in the preliminary feasibility study (Malone et al 2014b). Boundary voltages were computed on a fine 5 million element mesh, which was created from a CT scan of a human head and included three homogeneous tissues: brain, skull and scalp. Computation of the boundary voltages was done with Peits (Jehl et al 2015a) on all 16 processors of a workstation with two 2.4 GHz Intel Xeon CPUs with eight cores and 20 MB cache each. It took less than 2 min to compute the required 31 forward solutions for each frequency. 32 electrodes of diameter 10 mm were placed on the surface of the model in the same positions used to acquire EEG measurements (Tidswell et al 2001). The electrodes were modelled using the complete electrode model (CEM) (Somersalo et al 1992), and the contact impedance was set to Ω ⋅| | E 1k for all electrodes, where | | E was the electrode area. The amplitude of the current was set to 140 μA and twelve frequencies between 5 Hz and 5 kHz were used ( figure 1(a)). 31 linearly independent current injection pairs were created by finding the maximum spanning tree of the electrode positions, thereby maximising the distance between injecting electrodes. Voltage measurements were made for each injection on all adjacent pairs not involved in delivering current. The total number of measurements acquired for each frequency was 869.
Strokes were simulated by changing the conductivities of all elements within a 1.5 cm radius of the stroke location. Simulated locations were set in a posterior (figure 1(c)) or lateral position in the head ( figure 1(b)), and stroke conductivities were set to the spectral values of ischaemic tissue (figure 1(a)) or to the conductivity of blood, 0.697 S m −1 , for haemorrhage (Horesh 2006). Both proportional and additive noise was added to all simulated voltages: where ( ) ς rand indicates a random number drawn from a Gaussian distribution with zero mean and standard deviation ς. The standard deviation of the proportional noise was ς = 0.02% p and the standard deviation of the additive noise was ς µ = 5 V a , which correspond to human measurements (Goren et al 2015).
Electrode positions can currently be measured to around 1 mm precision using photogrammetry (Qian and Sheng 2011). Other technologies, such as the commercial MicroScribe, laser 3D scanners, or electrode helmets, can achieve an even higher precision in electrode localisation. Simulated electrode position errors were therefore created by drawing two random numbers for each electrode from Gaussian distributions with zero mean and standard deviations 0.5 mm, 1 mm, 1.5 mm and 2 mm. According to these drawn random numbers, the electrodes were then moved along a two dimensional surface coordinate system Jehl et al (2015b). Consequently, the overall electrode position errors were the combined errors of the displacement along the two surface dimensions, and deviations of up to three times the standard deviation of the error were expected in the majority of cases. Previous methods for computing the electrode movement Jacobian include the computationally intensive differential approximations by moving nodes in the mesh (Soleimani et al 2006), which are limited to relatively coarse meshes, direct methods based on the mesh geometry (Gómez-Laberge and Adler 2008) and the approximation error approach (Nissinen et al 2011). However, the analytical formulation of the Fréchet derivative with respect to the electrode boundary presented by Dardé et al (2012) is the most flexible approach, since it can be used for different electrode characteristics independent of mesh refinement and can be implemented in a very fast and memory efficient way (Jehl et al 2015b). Another advantage of this implementation is the possibility to move electrodes without altering the mesh geometry, by changing the assignment of the surface facets to the electrodes instead. This allows for larger movement of electrodes, while maintaining a good mesh quality and refinement.

Image quantification
The quality of an image was objectively quantified in terms of the ability to distinguish the stroke from the brain. This was done by assessing the fraction f s corresponding to the tissue that made up the anomaly. First, the images reconstructed on the 180 thousand element mesh were averaged onto cubic voxels with 0.5 cm sides. The volume P corresponding to the reconstructed perturbation was identified as the largest connected cluster of voxels with values larger than 50% of the maximum of the image (Malone et al 2014a, Jehl et al 2015b. Three measures of image errors were defined: (i) Location error: ratio between the distance ∥( ) ∥ x y z , , P P P of the centre of mass of the reconstructed perturbation P from the actual position, and the average dimension of the head ( ) d d d mean , ,

Multi-frequency tissue fraction reconstructions
When electrode positions were modelled accurately in the reconstruction algorithm, image reconstructions with electrode modelling correction were slightly worse than without correction ( figure 2(a)). The average image error of reconstructions without electrode correction was 10% and with correction 17%; particularly the reconstruction of the lateral haemorrhage was not good (24%). The reason for the decreased image quality was the larger number of variables to recover, which slightly increased the ill-posedness of the inverse problem. Consequently, conductivity changes were sometimes explained by electrode movements, if this reduced the value of the objective function.
With the traditional fraction reconstruction MFEIT algorithm, already 0.5 mm of electrode position modelling errors made stroke detection impossible ( figure 3(a)). The average image error of the reconstructions without electrode modelling correction was 55% ( figure 3(b)).
Simultaneous recovery of stroke tissue fractions and electrode positions significantly improved image quality in the presence of electrode modelling errors ( figure 4(a)). The average image error of the reconstructed images with correction was 23% ( figure 4(b)), excluding the three outliers (numbers 5, 10 and 14) only 16%. Remarkably, all three bad reconstructions occurred for ischaemic strokes, suggesting that local conductivity spectrum changes caused by ischaemia were more difficult to separate from electrode position errors than for haemorrhage. Such image errors are the result of the combination of errors added to the simulated voltages and the combined uncertainty on all electrode positions, and are consequently difficult to characterise. The shape error of the reconstruction of a lateral ischaemia with electrode movement of 2 mm (number 13) is misleadingly small, because the recovered perturbation was a diagonal disc with very similar x-y-z dimensions than the simulated stroke.

Electrode placement correction
The 2-norm of the difference between recovered and simulated electrode position mismatch for both surface dimensions ( x d and dy) was computed as ( ) , for electrodes i. While the electrode position correction was more accurate for ischaemic strokes when electrodes were correctly modelled, for position mismatch of 1 mm standard deviation the  correction was better in the presence of haemorrhages (table 1). The accuracy of the electrode position recovery tended to correlate with the quality of the reconstructed image (figures 2(a) and 4(a)). This is intuitive and was already observed for time-difference electrode movement corrections (Jehl et al 2015b).
The reason for the worse 2-norm of ischaemic electrode position correction of 1 mm errors, was electrode 17 (entries 33 and 34 on x-axis of figure 5). Interestingly, after one iteration of the proposed algorithm, the electrode position recovery of this electrode was still accurate. Only in the second iteration, the electrode was moved several millimetres in both surface dimensions. Since electrode 17 was located laterally, 5.4 cm from the centre of the lateral ischaemia, this affected the reconstruction of the lateral ischaemia more than the reconstruction of the posterior ischaemia ( figure 4(a)).

Discussion
The simultaneous recovery of stroke tissue fractions and electrode positions significantly improved image quality in the presence of electrode modelling errors and had only small negative effects on the image quality when electrode were modelled accurately. For movements between 0.5 mm and 2 mm standard deviation along both surface dimensions, the average image error was 23% compared to 55% without electrode correction. While reconstructions of haemorrhagic strokes were even successful in the presence of 2 mm of electrode errors, ischaemia detection was less reliable from 1 mm onwards. The lower reliability could be caused by a worse differentiability of ischaemic changes to electrode movements, but this could so far not be confirmed. Ideally, several reconstructions would have been made for each electrode movement level in order to characterise the effect of electrode modelling errors over a number of samples. However, the computational expense of multiple repetitions was prohibitive, since reconstruction of a single image took around 6 h. For the same reason, only two stroke positions were Table 1. 2-norm of the difference in simulated and recovered electrode position errors when electrodes were accurately modelled in the reconstruction algorithm (first row) and when there was a position mismatch of 1 mm standard deviation along both surface dimensions (second row).

Ischaemia Haemorrhage
Lateral Posterior Lateral Posterior 0 mm ⋅ − 0.8 10 3 ⋅ − 1.6 10 3 ⋅ − 1.7 10 3 ⋅ − 3.7 10 3 1 mm ⋅ − 13.0 10 3 ⋅ − 11.2 10 3 ⋅ − 7.1 10 3 ⋅ − 5.8 10 3  studied. Nonetheless, the produced images clearly illustrate the advantage of simultaneous tissue fraction and electrode position recovery in the presence of electrode modelling errors. Many methods for computing the Jacobian matrix with respect to electrode movement could have been used on the multi-frequency data. The presented method for correction of inaccurate electrode modelling had the advantage over previous approaches, that the mesh geometry did not have to be altered. This allowed for the iterative recovery of movements of several millimetres on fine meshes, which is not possible with traditional differential approaches for electrode movement recovery (Soleimani et al 2006). Furthermore, the implementation based on the Fréchet derivative (Dardé et al 2012) used here has been shown to be fast and memory efficient (Jehl et al 2015b).

Conclusion
Simultaneous iterative electrode position correction with the fraction reconstruction method using spectral constraints was applied to a numerical head phantom with realistic conductivities. Realistic noise was added to the simulated voltages to investigate the robustness of the proposed method. The results show that (i) the simultaneous recovery of tissue volume fractions and electrode position errors removed most image artefacts caused by inaccurately modelled electrodes. (ii) while haemorrhagic strokes could be reconstructed with electrode position errors up to 2 mm standard deviation in both surface dimensions, the reconstruction of ischaemic strokes was less reliable from electrode movements of 1 mm onwards.
Further work is required to understand why ischaemic stroke reconstructions were less reliable with the proposed method, and to correct for it. Additionally, it has so far not been studied how stable non-linear multi-frequency reconstruction methods are in the presence of geometric modelling errors, such as skull shape. We plan to validate the presented results in tank experiments with 3D printed head shaped tanks and skull (Jehl et al 2015b) and recommend the presented algorithm for any planned human studies.