Correcting electrode modelling errors in EIT on realistic 3D head models

Electrical impedance tomography (EIT) is a promising medical imaging technique which could aid differentiation of haemorrhagic from ischaemic stroke in an ambulance. One challenge in EIT is the ill-posed nature of the image reconstruction, i.e. that small measurement or modelling errors can result in large image artefacts. It is therefore important that reconstruction algorithms are improved with regard to stability to modelling errors. We identify that wrongly modelled electrode positions constitute one of the biggest sources of image artefacts in head EIT. Therefore, the use of the Fréchet derivative on the electrode boundaries in a realistic three-dimensional head model is investigated, in order to reconstruct electrode movements simultaneously to conductivity changes. We show a fast implementation and analyse the performance of electrode position reconstructions in time-difference and absolute imaging for simulated and experimental voltages. Reconstructing the electrode positions and conductivities simultaneously increased the image quality significantly in the presence of electrode movement.


Background
Electrical impedance tomography (EIT) is a three dimensional (3D) medical imaging technology which relies on a small, low cost device and is therefore of interest in many medical diagnostic fields. Compared to other medical imaging technologies such as MRI and CT, EIT has a low resolution due to the ill-posed nature of the image reconstruction problem. Small measurement noise and modelling errors introduce large artefacts into the images and can make detection of physiological changes impossible. Consequently, there is a lot of interest in finding ways to account for different sources of errors, such as varying electrode contact impedances, object boundary shape and electrode positions. Contact impedances were found to have particularly severe effects on image quality for current injection schemes where voltages are measured on injecting electrodes (Kolehmainen et al 1997), but less so when injecting current through electrode pairs while measuring on all other electrodes (Malone et al 2014).
The reconstruction of the boundary shape of the imaged object is of particular interest for thoracic imaging studies (e.g. lung ventilation, gastric emptying or heart cycles), because breathing changes the shape of the thorax significantly. Therefore, a simultaneous reconstruction of the conductivities and the boundary shape would improve the image quality (Boyle et al 2012, Dardé et al 2013. Given the difficulties in parametrising a realistic 3D surface, it may be more realistic to focus on mitigating geometrical errors in the electrode positions, of which there are relatively few parameters. Initial approaches to correcting electrode positions include differential approximations of the electrode displacement Jacobian (Soleimani et al 2006), direct methods based on the mesh geometry (Gómez-Laberge and Adler 2008) and the approximation error approach (Nissinen et al 2008). The most fundamental and straightforward way of computing a Jacobian matrix predicting voltage changes due to electrode boundary changes was shown by Dardé et al (2012). They derived an explicit formula for the Fréchet derivative with respect to the electrode boundary, which can be used to compute the Jacobian matrix with respect to electrode size, position and shape. Both, Soleimani et al (2006) and Dardé et al (2012), applied electrode movement corrections to cylindrical 3D saline phantoms with positive results, but the performance in anatomically realistic 3D problems has not been studied.
While absolute image reconstruction is theoretically possible in EIT, in realistic 3D applications it is prone to image artefacts caused by model inaccuracies and measurement errors (Kolehmainen et al 1997, Heikkinen et al 2002, Malone et al 2013. Consequently, the most commonly used time-difference method subtracts baseline measurements from the measurements of physiological changes, in an effort to cancel the effects of incorrect modelling and static instrumentation noise. The focus of this paper is on the imaging of head injury. Time-difference (TD) imaging can be used for long term monitoring of haemorrhagic transformations, where electrode movement regularly introduces artefacts (an introduction to the use of EIT for imaging haemorrhagic transformation can be found in Xu et al (2010)). For stroke type differentiation, time difference imaging is not possible and absolute or multi-frequency (MF) reconstruction methods have to be used (Jun et al 2009, Malone et al 2013. We chose to assess the quality of reconstructed images in TD and absolute (or static) imaging, to study the feasibility of including electrode model correction into existing MF methods.

Purpose
The purpose of this work was to develop a method for correction of errors in electrode modelling in EIT of the human head. The questions to be answered were: (1) which aspect of the electrodes-position, size, shape-most affects the boundary voltages and how accurately does the electrode boundary Jacobian matrix describe this aspects influence? (2) Does the simultaneous reconstruction of conductivities and electrode parameters remove image artefacts originating from imprecise electrode modelling in simulations and experiments? (3) Does the correction work in absolute imaging, where electrodes are modified iteratively?
These were addressed using a fast calculation of the Jacobian matrix with respect to electrode boundary perturbations, which was computed on a 4 million element mesh in less than a minute (described in section 2). Characterising the different aspects of the electrodes, it was found that the electrode position was the dominating variable and that it was accurately modelled by the electrode boundary Jacobian (section 3). The correction of electrode positions from time-difference data was therefore applied to an anatomically realistic 3D head model in simulations (section 4) and experiments on a 3D printed head shaped saline tank with skull (section 5). Finally, the correction for electrode movements was applied to absolute imaging (section 6).

Mathematical formulation
Computation of a Jacobian matrix in EIT typically requires the modelling of current flow in an object, the so-called forward problem. Forward solutions are commonly computed based on the complete electrode model (CEM), where the electric potential u solves the following problem on a body Ω with boundary Γ and electrodes E m m Here, σ is the conductivity distribution in Ω, z m the contact impedance of the electrodes and ν ∂ the outward unit normal to Γ. The CEM is proven to have a unique solution u H 1 ( ) ∈ Ω which depends continuously on I Somersalo et al 1992). The 'traditional' Jacobian matrix J σ which translates a change in conductivity to a change in measured voltages can then be computed with the efficient adjoint fields method (derived e.g. in the appendix of Polydorides and Lionheart (2002) Ω is the electric potential emerging when the drive current I d is applied to the electrodes and I u H a 1 ( ) ( ) ∈ Ω the electric potential when a unit current is applied to the two measurement electrodes. V da δ ∈ R is then the linearly approximated voltage change between the two measurement electrodes when the conductivity changes by L ( ) δσ ∈ Ω ∞ . The detailed derivation and proof of the electrode boundary Jacobian (EBJ), which describes voltage changes caused by electrode boundary movement, can be found in Dardé et al (2012); here we outline the parts relevant for the implementation. Electrode boundary changes can be characterised by C 1 vector fields on the electrode boundaries E ∂ as is the vector field defining the change in electrode boundary (e.g. change in electrode size, position or shape).
The measurement map including perturbed electrodes can then be considered as the operator and its Fréchet derivative with respect to the vector field v is fulfils where v τ is the component of v which is tangential to Γ, n E ∂ the outward normal of E ∂ which is tangential to Γ and u U , R . It is interesting to note that the computation of the EBJ resembles the computation of the traditional Jacobian matrix (with respect to the conductivity) using the adjoint field method (Polydorides and Lionheart 2002) in that the results of 'drive current' and 'measurement current' injections are used. Therefore, the calculation of the EBJ does not require any additional forward simulations to those performed for the calculation of the traditional Jacobian matrix.

Implementation
To describe the movement of an electrode along the surface of a mesh, a surface coordinate system is required. We created a surface coordinate system x y s s [ ] × as follows. Starting from the surface node with the largest y-coordinate (top of the head-(0, 0) in figure 1(a)), we moved along the surface in x and z direction and found surface nodes that are a defined distance d ε ± away from the previous point, where ε should be large enough that surface nodes are found. Distance d was chosen such that the surface coordinate points were far enough from each other to describe the surface curvature without significant influence of the mesh discretisation. If no more surface nodes could be found that either fulfilled these requirements or were above a certain y coordinate limit, we stopped. Once the x s and y s axis were found, the four quadrants were filled with points in the same fashion, by averaging the coordinates of all surface vertices within the area spanned by d ε ± (figure 1(a)). To describe directions of movement for each electrode we subsequently found the closest surface coordinate points to that electrode (figures 1(b) and (c)). A mapping back from surface coordinates to mesh nodes was not required. This approach had the advantage that the surface coordinate system could be stored together with the mesh and used for any electrode positions. The Jacobian matrix with respect to electrode radius was computed with a vector field v pointing from the electrode centre to the center of each electrode boundary edge.
For small meshes the electrode boundary Jacobian (EBJ) was calculated in Matlab using a modified version of Eidors (Adler and Lionheart 2006). On the 4 million element mesh used for simulations the EBJ was computed with the Parallel EIT Solver (Peits) (Jehl et al 2015) in order to reduce computing time and memory requirements. The implementation, in our case with linear basis functions, was similar in Matlab and Peits: for each electrode a sparse template matrix was constructed to store the contributions of drive and measurement fields. For each injection and measurement electrode pair this template matrix could then be multiplied by the electric potential fields of the drive and measurement currents to obtain the EBJ entry for this electrode for this line of the protocol. For the integral along each electrode boundary edge e with end point indices i 1 and i 2 , a 3 3 [ ] × sub-matrix was added to the template matrix , , , , , 1/3 1/6 1/2 1/6 1/3 1/2 1/2 1/2 1 , c e temp 1 2 elec 1 2 elec elec where the weights were obtained from Gauss-Lobatto quadrature (a detailed description of the procedure is given in appendix A). Once the template matrices were set up for all electrodes, the drive and measurement fields to multiply them with were the same forward solutions that were required to compute the 'traditional' Jacobian matrix and therefore no additional forward solutions had to be computed. The computation of one line of the EBJ took on average 0.03 seconds for a 4 m element mesh on ten processors using Peits, whereas one line of J σ took around 0.1 s to compute. The additional computation of the EBJ therefore had a minimal effect on the computation time. (a) (b) (c)

Simulation parameters
In order to study the characteristics of the voltage changes due to electrode boundary changes, a simulation study was performed on an 800 k element mesh (figure 1(b)). The mesh had refined elements (element size 3 mm instead of 6 mm) in the region where electrodes were located and very fine elements in the vicinity of one electrode (element size 0.4 mm). Current injection was then either simulated to be through a polar electrode pair involving the refined electrode and a coarse electrode on the other side of the head or through adjacent electrodes, again a coarse and the fine one. A representative set of ten different measurement electrode pairs (that were not injecting current) were chosen. The ten measured voltages were simulated for the following permutations: x s coordinate, y s coordinate or electrode diameter change, polar or adjacent injection, fine or coarse electrode. The diameter of 7 mm was altered between ±1.5 mm. Electrode movements were simulated in the range of 0.1-10 mm in both directions and both surface coordinate dimensions. We want to highlight that the mesh was not altered to model the electrode movement. Instead, the assignment of surface facets to the electrodes was used to move the electrodes. This approach was simpler and could account for larger movements of the electrodes. Subsequently, the simulated voltage changes due to the mentioned changes in electrode boundaries were scaled by the changes predicted by the EBJ. This scaling was done in order to be able to show all voltage changes together in one figure, even though they were of different magnitude. The plots of the EBJ accuracy have the simulated electrode change in millimetres on the x-axis and the electrode movement the EBJ would recover given the simulated voltage difference on the y-axis (in millimetres as well). Therefore an optimal prediction of voltage change would result in all 10 curves being perfectly diagonal (black dashed bar in figures 2 and 3). The absolute values of the EBJ with respect to position (both x s and y s ) and the EBJ with respect to radius were taken in order to determine which factor was more important. The accuracy of the three different EBJs for the three studied electrodes for adjacent and polar injection was evaluated by linearly interpolating the voltages from simulated electrode boundary changes and comparing the slope to the value of the EBJ.
The simulations were performed using a current level of 133 μA and contact impedance of the electrodes (with area E ) of z c =1 kΩ E on a realistic head mesh (800 k tetrahedral  . The ten simulated measurements were the first ten lines of the 'EEG31' protocol (Tidswell et al 2001) and for the adjacent injection one of the injecting electrodes was substituted with an adjacent one.

Voltage dependence on electrode characteristics
3.2.1. Comparison of x s , y s and d prediction of EBJ. The EBJ accurately predicted voltage changes caused by electrode movement (figure 2). The changes caused by electrode diameter variations were overestimated by the EBJ, and were significantly smaller than the changes caused by electrode movement (colour bar in figure 2).

3.2.2.
Effect of mesh refinement. The mesh refinement around an electrode did not limit the precision of the EBJ. However, if the electrode positions were corrected on a coarse mesh then the discretisation error was bigger (the big steps in figure 3(b) as opposed to the small steps in figure 3(a)). This is to be expected, since the boundary of the electrode could only change in intervals equal to the element size.

Comparison of polar and adjacent current injection.
The voltage changes caused by electrode movement were only minimally different between adjacent and polar current injection (figures 2(b) and 3(a)). This matches expectations, considering that the formula for the computation of the EBJ only depends on electric potential differences at the electrode boundary.

Precision of electrode boundary jacobian
The slopes of simulated voltage changes due to electrode boundary changes were linearly fitted and compared to the EBJ values, to get a measure of the EBJ precision (first three rows of table 1). The voltage changes due to electrode movement were well approximated by the EBJ, with most values having around 10% mismatch. The predicted changes due to electrode size changes were on average 690% off. The average absolute values of the EBJ were  compared for the different types of electrode perturbation (movement along x s and y s and change radius r) to illustrate that the changes due to electrode size were significantly smaller than the ones caused by movement (last three rows of table 1).

Time-difference image reconstruction algorithm
In both the simulation study and its experimental validation, we were using spherical perturbations in the brain to represent a stroke. We were therefore able to introduce prior information into our reconstruction by using first order Tikhonov regularisation to bias the algorithm towards finding small connected perturbations. All time-difference images were created with a standard least-squares minimisation using generalised singular value decomposition. Many other algorithms have been successfully applied to EIT data (Lionheart et al 2004), and could equally well be used here for simultaneous reconstruction of conductivity changes and electrode movements. In the following, we give a brief outline of our time-difference image reconstruction approach.
In EIT, the image reconstruction problem can be described by the minimisation of the cost functional with x the change in conductivity and electrode positions we wish to reconstruct, v the voltage difference between the two measurements, F x ( ) a non-linear function relating conductivity and electrode position changes to voltage changes, x Γ the expected variance of the reconstructed variables (conductivities Γ σ and positions p Γ ) and the first order Tikhonov regularisation term for the conductivities combined with zero order Tikhonov for the electrode positions Using the substitution D D and the minimum can be found by setting the derivative to zero Using the generalised singular value decomposition (gSVD) of the Jacobian J and scaled regularisation matrix D x this equation can be solved for x where J U X Λ = and D VMX x = . Since both Λ and M are diagonal matrices with a limited number of generalised singular values (Hansen 1994), they can easily be inverted.
To reduce the computational cost of calculating the gSVD and to prevent the 'inverse crime' (Lionheart et al 2004) we used a much smaller hexahedral mesh (2526 elements of 1 1 1 × × cm) for the image reconstructions ( figure 4(a)). The Jacobian matrix J σ which was computed on the fine mesh was summed into geometrically regular cubes J hex and the Laplacian matrix L for the first order Tikhonov regularisation was computed on this hexahedral mesh. In order to simultaneously reconstruct conductivity changes and electrode movements, the 'traditional' Jacobian matrix J hex had to be combined with the EBJ as = J J EBJ Equivalently, if only the conductivity changes or only the electrode movements were reconstructed, only the relevant Jacobian matrix and the corresponding part of D was used.
The expected variance of the conductivity was set to var 0.1 = σ S m −1 and for the electrode positions to var 1 p = mm. For all reconstructed images and electrode positions, the regularisation factor 2.1 10 2 7 λ = ⋅ − was kept constant. The value was chosen based on the shape of the L-curve (Hansen 1994), however the L-curve was not pronounced enough to choose λ in an automated way. The solution norm x 2 ∥ ∥ for our choice of λ was generally around 0.1 and the residual norm Jx v 2 ∥ ∥ − around 1.0 10 4 ⋅ − . The colour bar of all images was scaled according to the largest reconstructed change in the whole mesh. Therefore images of slices sometimes do not contain the maximum value of the colour bar.

Image quantification
In order to evaluate the image quality objectively, three quantification methods were applied. The volume P corresponding to the reconstructed perturbation was identified as the largest connected cluster of voxels with 50% of the maximum of the image (Malone et al 2014). The  region of interest (ROI) was defined as the largest connected cluster of voxels with 50% of the maximum of the simulated conductivity change.
• Localisation error: ratio between the displacement x y z , , P P P ( ) of the centre of mass of the reconstructed perturbation P from the actual perturbation location, and the average dimension of the mesh d d d , , x y z ( ) x y z d d d , , mean , , .
• ROI change: difference of the average value of the reconstructed image (d r σ ) in the ROI and the average value of the actual perturbation (d a σ ), divided by the average value of the perturbation where the Tikhonov smoothing correction factor f tikh corrects for the effect that the reconstructed perturbation is larger in volume and smaller in amplitude by scaling the reconstructed conductivity change within the ROI ( f 10 tikh = was used throughout). • ROI noise: noise-to-signal ratio of the reconstructed image (d r σ ), computed as the ratio between the standard deviation (std) outside the ROI and the average value in the ROI std d mean d .  (Tang et al 2008). The current level was 250 μA, contact impedances of all electrodes were set to z c = 220 E Ω ⋅ and the diameter of the 32 electrodes was 10 mm. The plastic perturbation had a radius of 1.5 cm and a conductivity of 0.0001 S m −1 ( figure 4(b)). The injecting pairs of electrodes were chosen to maximise the distance between electrodes by finding the maximum spanning tree, weighted by the inter-electrode distances. Measurements were made for each injection on all adjacent electrode pairs not involved in delivering the current, giving a total of 869 measured voltages from 31 independent current injections. All voltages, the conductivity Jacobian matrix and the EBJ were computed with Peits (Jehl et al 2015).
The level of noise added to simulated differential voltages was chosen to match the experiments and consisted of 0.006 p ς = % proportional noise and 1 a ς = μV additive noise, such that where rand( ) ς indicates random numbers drawn from a Gaussian distribution with zero mean and standard deviation ς.

Electrode position recovery
To reconstruct conductivities and electrode positions simultaneously, the algorithm outlined in section 4.1 is used. Analogously, if only electrode positions are reconstructed, then the full Jacobian is replaced by only the EBJ. The performance of the EBJ for electrode position recovery was validated with three different recovery modalities: (1) only the EBJ was used to recover electrode positions when the electrodes moved and the conductivities did not change; (2) only EBJ was used when electrode positions and conductivity (plastic ball in back of the head) changed; (3) the full Jacobian matrix was used to reconstruct conductivities and electrode movement at the same time when both electrode positions and conductivity have changed. All these differently recovered electrode positions were plotted together with the actual electrode movement for two cases, one where only one electrode on the back of the head moved by 5 mm ( figure 5(a)) and the other one where all electrodes were moved along both x s and y s by random values with standard deviation 1 mm ( figure 5(b)).
The 2-norm of the difference of recovered movement versus actual movement shows that the simultaneous reconstruction of a conductivity perturbation and electrode movements (last row in table 2) was close to that of the sole recovery of the electrode positions when no perturbation was introduced (first row in table 2). Contrastingly, if the conductivity changed simultaneously with the electrode positions, the recovery of only the electrode positions resulted in a significant over-correction (second row in table 2), especially close to where the perturbation was introduced (numbers 41-64 in figures 5(a) and (b) corresponding to x s and y s coordinates of electrodes 21-32 towards the back of the head).

Images
The best possible images (i.e. in the presence of no electrode movement) that were achievable with time-difference reconstructions without electrode correction (figure 6(a)) and with electrode correction (figure 6(d)) were qualitatively similar. The simultaneous reconstruction of conductivity and electrode positions slightly reduced the contrast and the precision of the imaged perturbation.
If one electrode in the top back of the head was moved by 5 mm, the simple time-difference reconstruction resulted in an unsurprisingly noisy image (figure 6(b)) with large artefacts Figure 5. Recovered electrode movement for three different cases: only EBJ was used and only the electrode positions changed, only EBJ was used and additionally to the electrode movement a perturbation was inserted at the back of the head, and the complete Jacobian J was used when both electrodes and conductivity changed. In (a) an electrode on the back of the head was moved along x s by 5 mm and in (b) all electrodes were moved along x s and y s by a random value. On the x axis of these two plots the entries 1 and 2 correspond to the x s and y s coordinate of electrode 1, then the next two entries correspond to electrode 2 and so on. around the moved electrode. Simultaneous reconstruction of the electrode positions (figure 6(e)) restored the quality of the image almost to the ideal case (figure 6(a)). As with a large movement of one single electrode, smaller movements of all electrodes had a detrimental effect on the image quality for simple time-difference imaging (figure 6(c)). The image reconstructed with the full Jacobian was again largely unaffected by this large degree of electrode movement (figure 6(f)). All these findings were summarised for three different positions of the perturbation by assessing the image quality according to the previously defined error metrics (figures 7(a)-(c)). The larger localisation error (and inherently larger ROI errors) for a perturbation on the side can be explained by the reduced sensitivity of the current protocol used for these simulations.

Experimental setup
Two saline tanks were printed using the 3D printer Makerbot Replicator 2 from Makerbot Ind. The model for the tank was created from the same MRI segmentation used for the mesh creation, whereas the skull was segmented from a corresponding CT scan. The skull model was further edited by introducing small holes, such that the conductivity matched that of a real skull (Avery 2015). One tank was printed with a modified EEG 10-20 electrode placement (Tidswell et al 2001, Avery 2015 and the other tank had perturbed electrode positions with random displacements along x s and y s with standard deviation of 1 mm, matching the electrode movement simulated in section 4. For the recordings, the tanks were filled with 0.4 S m −1 saline. The electrode contact impedances were measured to be 220 E Ω ⋅ . Current was injected at 1.76 kHz with 250 μA amplitude, to approximately match the allowed current level in human measurements (Dybdahl 2009). A slightly modified 32-channel version of the KHU Mark 2.5 system (Wi et al 2014) was used for the recordings and each measurement was repeated 20 times over the course of approximately one minute. Plastic perturbations with 3 cm diameter were placed in approximately the same locations used in the simulations.

Images
With baseline and perturbation measurements in the same saline tank, both using J hex and using J resulted in a reconstruction of the perturbation with only minimally worse quality than for simulated noisy voltages (figures 8(b) and (c) and the first two bars in the quality measures in figure 9). However, when reconstructing a perturbation measured in the tank with the moved electrodes, then large artefacts in the outer layers of the head near the skull were observed (figures 8(d) and (e)). These artefacts strongly influenced the localisation errors (figure 9), where an artefact was interpreted as the reconstructed perturbation. The drop in image quality was most likely caused by small differences in skull placement in the two different saline tanks.
Simultaneous reconstruction of conductivities and electrode movements improved the image visibly (figure 8(e)), but did not remove all artefacts. Since electrode movement cannot completely account for all voltage changes, artefacts caused by other effects (such as skull position in the two tanks) remain in the conductivity image. Parts of the voltage differences caused by these geometrical differences were, however, pulled into the electrode movement recovery. Therefore, the precision of the recovered electrode movement between the two 3D printed tanks was lower than in simulations (error 2-norm 8.0 mm instead of 4.7 mm in simulations).

Absolute reconstruction algorithm
In contrast to time-difference imaging, absolute reconstructions do not require a baseline measurement. In many realistic applications, absolute imaging fails-among other reasons-because of small differences between the model and the imaged object, resulting in large image artefacts. Since the 3D printed saline tank and skull were very close to the corresponding mesh, we wanted to analyse how the correction of electrode positions performed in iterative absolute reconstructions. A generalised minimal residual algorithm with tolerance of 5 10 11 ⋅ − was used to solve the standard Gauss-Newton problem for the search direction d k . The reconstructed change in conductivity and electrode positions was then updated in each iteration where k α is the step size and was optimised using the Brent line search method (Brent 1973) with a gold-section bracketing loop to find the Brent abscissae. Analogous to the time-difference approach, the reconstructions were weighted by the expected variance. In iterative reconstructions it is important to ensure the positivity of the elements' conductivity in the Gauss-Newton step, which was done here by using the substitution Applying this substitution to (21) moved the scaling by the expected standard deviation of the variables from the regularisation term to the Jacobian matrix, making the Jacobian dimensionless. According to the chain rule, the logarithmic part of the substitution results in a multiplication of the Jacobian entries corresponding to conductivities by e σ . After + y k 1 was found using Gauss-Newton, the electrode positions and conductivities were updated according to the corresponding + x k 1 . The regularisation parameter for all absolute reconstructions was set to 10 2 1 1 λ = − and the stopping criterion was set to <10 −7 change between iterations with a maximum of 6 iterations. As a trade-off between computational efficiency and precision of the forward solutions, a tetrahedral mesh with 1 million elements was used for the absolute reconstructions.

Images
All absolute reconstructions from simulated data had a large positive artefact in the front of the head, where the sensitivity is high. When the full Jacobian matrix was used, the algorithm unsuccessfully tried to reduce this artefact by moving the electrodes in the front of the head (left side of figure 10(c)). When evaluating the recovery of electrode movements, we therefore subtracted this 'baseline' movement from the reconstructed movement of the electrodes when they have actually moved. Doing this, we found a good performance of the iterative electrode movement correction (table 3, as compared to the last row in table 2). The large artefact results in a poor image quantification score for all reconstructed images, with or without inclusion of the EBJ. Subsequently, we illustrate the difference the electrode correction makes by showing one reconstruction without and with EBJ, both thresholded to one third of the maximum reconstructed change (figure 10).
We were not able to reconstruct meaningful absolute images from the experimental data. Even on the tank with no electrode movement the reconstructions did not show the perturbation.

Electrode boundary jacobian characteristics
We observed that the measured voltages changed linearly with electrode movement over 1 cm and were very accurately predicted by the EBJ (the linearity was already observed multiple times, as reviewed in Kolehmainen et al (1997)). There were outliers where some voltages did behave highly non-linearly (figure 2(a)), but as seen from the colour coding these outliers were of small amplitude. Therefore they were not expected to adversely influence the reconstruction of electrode movements. Changes in electrode diameter had a less predictable influence on the measured voltages (figure 2(c)). Firstly, the EBJ predicted changes were between one and two orders of magnitude smaller than for electrode movements and, secondly, the observed voltage changes were even smaller than the EBJ predictions. Changes in electrode shape can be viewed as a combination of a change in diameter and the discretisation errors seen in the electrode movement on the coarse electrode ( figure 3(b)). Therefore these changes can be expected to be small, even though we did not analyse this specifically. We concluded that changes in electrode size and shape could be ignored in the presence of electrode movement, and that they were not very well approximated by the EBJ. The focus of the simulation study and the experiments was therefore laid on the correction of electrode movements, which had the strongest impact on the image quality.
The mesh refinement around the electrode had no influence on the precision of the EBJ (figure 3), which indicated that meshes did not have to be altered for the use of electrode movement correction. Electrodes could be moved by simply assigning different surface facets to the electrode, without introducing errors larger than the discretisation errors caused by the element size. This facilitated the electrode movement correction significantly, since the mesh   remained untouched and the electrodes could be moved by more than the size of the finite elements.

Simulation study
Time-difference reconstructions were very stable and gave good images in all cases. Using the EBJ, the reconstruction algorithm could account for very large electrode movements and thereby reduce the image noise notably. In the presence of image artefacts which were not caused by electrode movement, using the EBJ resulted in a less precise reconstruction of the electrode positions. This was because the reconstruction algorithm reduced the cost functional by moving parts of the image artefact into the electrode positions. However, since 64 parameters constituting the positions of the 32 electrodes could not account for these artefacts, they were still visible in the resulting conductivity image.

Experimental validation
Time-difference reconstructions in the printed head tank with the correct electrode positions gave good results. Using the baseline measurement with the correct tank and the perturbation measurement on the tank with changed electrode positions resulted in noisier images than in the corresponding simulations. The additional noise was most likely caused by small geometrical differences in the skull positioning in the two different tanks, slightly different saline levels between the measurements and ambient and system temperature differences. While such errors were not present in conventional time-difference measurements, they manifested themselves when two different tanks were used with a half an hour break in between setting up each tank. Still, the simultaneous reconstruction of conductivity and electrode movement significantly improved the reconstructions and allowed for the detection of the perturbation.

Absolute reconstructions
The iterative absolute recovery of electrode positions was accurate and strongly reduced the electrode movement related artefacts. Still, the absolute conductivity reconstructions with or without electrode correction showed a systematic artefact in the front of the head. This was in a region of strong sensitivity and might have been exacerbated by the comparative imprecision of the skull in the one million element mesh used for the absolute reconstructions. Moving to a finer, more precise mesh would result in very long computation times even for a single image reconstruction. We were not able to get meaningful absolute reconstructions from the experimental data, suggesting that the simulations and experiments were still too different even though the geometry of the tank and inserted skull were printed with a precision of 0.2 mm. However, the placing of the skull in the tank was less precise than this. More advanced absolute reconstruction algorithms might be able to retrieve more information from our experimental data by including geometrical uncertainties around the skull into the imaging method.

Conclusions
We have applied the Fréchet derivative of the EIT forward problem with respect to electrode boundary changes to a realistic 3D model of the human head, using a fast implementation. This allowed us to reconstruct the conductivities and electrode positions simultaneously, and therefore make reconstructions more stable in the presence of electrode movement. For this realistic EIT setting we have found that: (i) The electrode position has a much stronger effect on the boundary voltages than the electrode diameter and shape. The voltages changed linearly for electrode movements up to 1 cm and were accurately predicted by the EBJ. (ii) The simultaneous reconstruction of conductivity and electrode positions worked very well. Image artefacts caused by electrode movements were reliably removed and the electrode positions were accurately corrected. The positive results of the simulation study were confirmed in experiments on a 3D-printed saline tank containing a realistic 3D-printed skull. (iii) While the correction of electrode movements was accurate in iterative absolute imaging as well, the absolute conductivity reconstructions from simulated data had a large artefact in the front of the head. No meaningful absolute images could be reconstructed on the 3D-printed saline tank.
We conclude, that the proposed application of the EBJ works well in time-difference reconstructions and can be used for long term monitoring of physiological changes. For stroke type differentiation, time-difference measurements are not possible. Therefore, existing multifrequency reconstruction algorithms need to be adapted to correct for imprecisely modelled electrode positions using the EBJ. [˜˜˜] = . For the integral along one edge, we kept the dot product of the outward normal and the permutation vector field, v n e ( ) ⋅ , fixed. Therefore, the electrode Jacobian entry for one electrode for one current injection and voltage measurement electrode pair can be written as the sum of the integral approximations along all element boundary edges