Are patient specific meshes required for EIT head imaging?

Head imaging with electrical impedance tomography (EIT) is usually done with time-differential measurements, to reduce time-invariant modelling errors. Previous research suggested that more accurate head models improved image quality, but no thorough analysis has been done on the required accuracy. We propose a novel pipeline for creation of precise head meshes from magnetic resonance imaging and computed tomography scans, which was applied to four different heads. Voltages were simulated on all four heads for perturbations of different magnitude, haemorrhage and ischaemia, in five different positions and for three levels of instrumentation noise. Statistical analysis showed that reconstructions on the correct mesh were on average 25% better than on the other meshes. However, the stroke detection rates were not improved. We conclude that a generic head mesh is sufficient for monitoring patients for secondary strokes following head trauma.


Background
Electrical impedance tomography (EIT) is a non-invasive three-dimensional medical imaging technology relying on a small, low cost device and is therefore of interest in many medical diagnostic fields. Because of the ill-posed image reconstruction problem, small measurement noise and modelling errors introduce large artefacts into the images and can make detection of physiological changes impossible. In static imaging, where no baseline measurement is available, it has been shown that inaccurate head models have a detrimental effect on image quality (Kolehmainen et al 1997, Jehl et al 2015b. Therefore, most EIT applications use time-difference (TD) data, based on the observation that most geometric and system related errors cancel out when a reference measurement is subtracted from the data measurement (Brown 2003). TD applications of head EIT include the analysis of epileptic seizures (Bagshaw et al 2003) and imaging of cortical activity (Tidswell et al 2001), but the focus of this paper is the monitoring of patients with traumatic brain injury (TBI) for secondary bleeding (Xu et al 2010).
It has previously been found that even in TD, a head-shaped finite element model produces better images than spherical models (Bagshaw et al 2003). Based on this finding, ways have been investigated to create patient-specific head models based on either a computed tomography (CT) or a magnetic resonance imaging (MRI) scan of the head (Tizzard et al 2005, Vonach et al 2012. However, it has not yet been ascertained if differences between heads introduce enough errors into reconstructed images to justify the use of patient specific head meshes. If not, then images for studies of EIT of the adult head could all be undertaken with one anatomically accurate model of a generic head, which would represent a considerable saving in time and effort.

Purpose
The purpose of this work was to assess how accurately a human head has to be modelled to enable EIT monitoring and imaging. The questions to be answered were: (1) how large are voltage and image errors introduced by using a generic head mesh, as opposed to a subject specific one? (2) Is a subject specific head mesh required for the detection of localised conductivity changes?
These were addressed in computer simulation on accurate head meshes from MRI and CT scans of four patients' heads, created using a novel method. Two conductivity perturbations with different magnitude, haemorrhage and ischaemia, were chosen to test the reconstruction quality. Boundary voltages and Jacobian matrices were computed on all four meshes, for ischaemia or haemorrhage at five different locations. Time-difference images were then reconstructed from the simulated noisy boundary voltages, using the four different meshes and a coarse, homogeneous head model. Analysis of the resulting TD image quality-evaluated as the inverse of objective measures of image errors-showed that even though the image quality was generally better with the correct mesh, the number of correct perturbation detections was not higher than with the other head meshes and even the homogeneous mesh.

Mesh creation
The mesh generation pipeline used required a CT scan and an MRI scan of the same head

Segmentation
Previous head segmentations for EIT were based on either a CT scan and then morphing a brain into the shape of the skull (Vonach et al 2012), or on MRI scans with poor sensitivity to the skull shape (Tizzard et al 2005, Sadleir et al 2010, Vonach et al 2012. The method proposed by Vonach et al (2012) did not allow for meshing of white matter or other tissues apart from brain, cerebrospinal fluid (CSF), skull and scalp. Also, it required that one tissue completely enclosed the inner one, which is not anatomically correct. The greatest difficulty in head segmentation for EIT is, however, the accurate modelling of the skull with respect to the soft tissues. The combination of the conductive scalp and highly resistive skull underneath means that a significant proportion of the injected current is shunted around the brain (Abascal et al 2008). For this reason, CT scans were used in this study for accurate skull segmentations and MRI scans of the same heads for soft tissue segmentation. The segmentation was mostly manual and took around two working days per head. In the following, the approximate procedure is described, with segmented layers printed in italics and Seg3D tools in capitals (figure 1).
• CT and MRI head scans of the same patient were loaded into Seg3D. The rotational mismatch was estimated and the MRI scan was then rotated to align with the CT scan using MIPAV. Once the two scans were aligned, all segmentation work was done manually in Seg3D. The resolution of the CT and MRI scans was increased to 512 512 512 × × , using Gaussian interpolation (resample tool). Increasing the number of pixels while smoothing the images was required to generate smooth surfaces in the resulting finite element meshes.
• The skull was segmented first by thresholding the CT scan at an appropriate intensity level (around one third of the maximum). The diploë was then found by filling holes in the skull and then removing the skull from this newly created layer. The modelling of the diploë included anisotropy into the overall skull conductivity. • The soft tissues were then extracted from the MRI scan. The first soft tissue to be segmented was the white matter, because it had a clearly visible contrast to the surrounding grey matter and could therefore be accurately found with the threshold tool (between around 40% and 50% intensity). Areas with similar intensity to the white matter were mostly outside the skull cavity, and could therefore be removed from the white matter by running remove with the skull segmentation as the mask, and then finding connected components on the white matter. • Next, the grey matter was found in a similar manner to the white matter. Because the grey matter MRI intensity (around 30-40%) was similar to tissues in the lower part of the head, it was generally not possible to disconnect the actual grey matter from the tissues with similar intensity by simply removing the skull. Therefore, this was done manually using the paint brush tool to remove the connections between the grey matter and other tissues, based on physiological atlases. This was the most time-consuming step in the head segmentation. • The superior sagittal sinus had a high intensity on the MRI scans (50-70%) and was easily segmented by thresholding and finding connected components. It was added to the diploë layer. • By thresholding the MRI scan at very low intensity (below 10%), the background was segmented. At the nose and ears, the air filled cavities were manually removed (paint brush) from the background and assigned to a new segmentation layer air. This was done, in order not to have highly irregular surfaces in the resulting finite element meshes. • The cerebrospinal fluid (CSF) and eyes were combined to one layer, and found by thresholding the MRI scan at an intensity level of 10-20% and subsequent subtraction of grey matter, white matter, skull and diploë and superior sagittal sinus. Any unassigned gaps between grey matter and skull were filled with CSF by iteratively dilating it and removing the other layers to prevent overlap.
• All unassigned voxels were finally assigned to a scalp layer, including most of the lower head, nose and throat. This simplification was deemed acceptable, because these areas were far from the electrodes and the region of interest (i.e. the brain). The scalp was found by thresholding the MRI to include all relevant tissues (∼15-70%) and then subtracting all other segmented layers. dilation and removal of the background ensured that the scalp covered the surface of the segmented head. • When all tissue layers were segmented, they were smoothed by running smooth binary dilate -> erode and overlaps were prevented with boolean remove. Finally, the segmentation was saved in .nrrd format and, using MedInria, translated into the .inr format required by the mesher.

Meshing
The high quality 3D mesh generator of the computational geometry algorithms library Cgal (Cgal 2015) was used to create tetrahedral meshes directly from the joint MRI&CT segmentation. The default 3D mesh generator of Cgal (written in C++) was adapted, by defining a problem specific element sizing field. This sizing field defined the desired finite element size throughout the head and was used to refine elements near the electrodes and towards the surface of the head. Specifically, all elements within a radius of 10 mm of the centre of an electrode were required to be 0.5 mm small, while the size of the other elements was assigned linearly from the surface to the centre of an ellipsoid fitted into the head (2 mm on the surface to 4 mm in the centre of the head). The effect of the mesh refinement around electrodes can be seen on the left of the slice taken of one head mesh, which cuts through an electrode located in posterior position (figure 2). The resulting meshes contained 3.3-3.5 million tetrahedral elements and their quality was controlled with the Joe-Liu quality measure (Liu and Joe 1994) where q is the quality, V the element volume and l the length of the edge between the ith and jth element vertex. All meshes had an average quality above 0.89 and not more than 7 elements with quality below 0.1 (with the worst element having a quality measure of 0.073).
After the manual location of nasion, inion and left and right ear tragus, the 32 electrode positions were computed automatically for each mesh using Matlab (MATLAB, The MathWorks, USA). The electrode positions were defined according to the EEG 10-10 system (Nuwer et al 1998) and matched the electrode assignment of the EasyCap EC40 (Brain Products GmbH, Germany). The head dimensions were analysed based on the distances between EEG 10-10 electrode positions Fpz and Oz (length), T7 and T8 (width) and the height was found as the difference between maximum and minimum z dimensions of nodes in the meshes (table 1). These head dimensions agree with anthropometric findings of an average head length of 188 8.1 ± mm and an average head width of 145 6.2 ± mm (Motmans and Ceriez 2005).

Voltage simulation
The conductivities σ of the used tissues were set to values compiled from literature for frequencies around 10 kHz (table 2; Horesh 2006). The current level was 250 μA and contact impedances of all 32 electrodes were z E 1000 c = Ω⋅ , with E being the contact area of the electrodes with diameter 10 mm. The current injection pairs of electrodes were chosen to maximise the distance between electrodes by finding their maximum spanning tree. 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 and the Jacobian matrix were computed with Peits (Jehl et al 2015a), resulting in computation times of less than 11.5 min for voltages and Jacobian matrix on each healthy head model, and less than 3.5 min for the computation of only the voltages for each simulated stroke (on 12 processors with 20 MB cache each).
The five positions of the simulated perturbations were based on EEG 10-10 nomenclature electrode positions. Anterior, central and posterior were 25%, 50% and 75% along a line between Fpz and Oz, respectively. The lateral and superior positions were halfway between the central position and T8 and Cz, respectively. The simulated strokes had a radius of 1.5 cm and conductivity 0.7 S m −1 for haemorrhages and 90% of the healthy conductivity for ischaemias (Horesh 2006, Dowrick et al 2015. Three levels of noise were added to simulated voltages and were chosen to match, in ascending order, tank experiments (Jehl et al 2015b), measurements on humans with relatively low noise and with higher noise (Goren et al 2015). This noise is referred to as system noise, and accounts for noise from the instrumentation. The noise levels were 0.006% p ς = proportional noise and 1 a ς = μV additive noise, 0.01% where rand( ) ς indicates random numbers drawn from a Gaussian distribution with zero mean and standard deviation ς.

Image reconstruction
First order Tikhonov regularisation was used to bias the algorithm towards finding small connected perturbations (Lionheart et al 2004). All images were created with a standard leastsquares minimisation using generalised singular value decomposition (gSVD). The advantage of using gSVD was that it only had to be computed once for each reconstruction mesh and  could then be used for all image reconstruction by simple matrix multiplication (Jehl et al 2015b).
To reduce the computational cost of calculating the gSVD and to prevent the 'inverse crime' (Lionheart et al 2004), much smaller hexahedral meshes (3-4 thousand elements of 1 1 1 × × cm) were used for the image reconstructions. The Jacobian matrices were computed on the fine meshes and then projected onto the geometrically regular cubes of the hexahedral mesh. The Laplacian matrices for the first order Tikhonov regularisation were computed on the hexahedral meshes.
For all reconstructed images, the regularisation factor was kept constant for each noise level, and was λ = × − 5 10 4 for low noise, λ = × − 1 10 3 for medium noise and λ = × − 2 10 3 for high noise. The value for the lowest noise was chosen manually based on the shape of the L-curve (Hansen 1994). As the noise was roughly doubled, the regularisation factor was also doubled to account for the higher noise level. The colour bars of the images were scaled according to the largest reconstructed change in the whole mesh. Therefore images of slices sometimes do not contain the maximum value indicated in the colour bar.

Image and voltage error measures
Three metrics were used to objectively evaluate the image errors. These were error in location, and the ratio of background to perturbation, measured first as the average and then as the absolute sum of changes. All these were 0% in ideal images. The overall error was expressed as the sum of these three metrics in percent, and image quality as its reciprocal. Acceptable quality images were empirically found to have <100% image error or >0.01 image quality, respectively.
The volume P corresponding to the reconstructed perturbation was identified as the largest connected cluster of voxels with at least 75% of the maximum absolute change in the image (Jehl et al 2015b). The region of interest (ROI) was defined as the largest connected cluster of voxels with 50% of the maximum of the simulated conductivity change. The used error metrics were then • Location error: ratio between the distance x y z , , The overall image quality was then defined as the inverse of the sum of the three error measures in percent (e.g. the image quality for 10% location error, 5% ROI contrast error and 25% ROI noise would be 1/(10 + 5 + 25) = 0.025). Examples of image reconstructions with differently assessed quality illustrate that up to 100% image errors it was possible to detect the stroke visually (figure 3). Between 100-120% some strokes could still be seen, while others could not. If the image errors were larger than 120%, it was no longer possible to detect strokes. 100% image errors were therefore defined as the limit of stroke detectability. The occurrence of a reconstruction with more than 120% was highlighted in tables 3 and 4 by colouring the entry red.
Voltage errors between two measurements v v 1 2 − were assessed with two metrics.

Analysis of the image quality
4.1.1. Was the image quality better with the correct mesh? For reconstruction on the correct mesh, the image quality metric was on average 25% better, than when different meshes were used for simulation and reconstruction. Sometimes, image reconstructions on the wrong mesh were better than reconstructions on the correct mesh. The image quality of reconstructions on the correct mesh was plotted against the average image quality of reconstructions on the other three head meshes (figure 4). For low noise, the improved quality using the correct mesh is clearly visible (blue crosses). The difference in image quality between correct mesh and other mesh reconstructions reduces with increasing noise. When filtering out reconstructions that were unsuccessful in both modalities (i.e. quality below 1/100% 0.01 = for both), then the following ratios of better reconstruction with the correct mesh versus better reconstruction with the other meshes were obtained: 32 : 8 at low noise, 22 : 16 at medium noise and 20 : 10 at high noise. An example of a reconstructed image on the correct mesh and the same reconstruction on a different mesh is shown in figure 5. 4.1.2. Were stroke detection rates higher with the correct mesh? While using the correct mesh for reconstruction increased image quality on average, the total amount of successful stroke detections was not improved (figure 4). The number of cases where reconstructions with the correct mesh were successful, while reconstructions with the wrong meshes were not, was 3 (all at the highest noise level). On the other hand, there were 5 cases where only the reconstructions on the wrong meshes gave an acceptable image quality (2 at the highest noise level, 1 at medium noise and 2 at low noise).

Which factors did stroke detection rates depend on?
In this section, we compare the influence of mesh differences, stroke position, stroke type and system noise level on stroke detection rates. The system noise was found to have the strongest impact on image quality, and mostly affected reconstructions of ischaemia.
No clear differences were found between the meshes used for voltage simulation and image reconstruction (table 3). The only surprise was, that even at the highest noise level, almost all reconstructions from voltages simulated on mesh 1 were successful. No explanation was found for this observation. The least successful reconstructions of haemorrhages were observed for the posterior position, while the lateral position was the least successful for ischaemia recovery (table 4). The added system noise had a very strong influence on the successful reconstructions of ischaemic strokes. At the medium noise level, 66 out of 80 ischaemias were detected, while with the highest noise this number dropped to 33 (table 4).

How large were voltage changes caused by the perturbations?
The average 2-norm of the voltage change caused by haemorrhages was 301 μV (1.3% proportional) and for ischaemias 27 μV (0.1% proportional). These values were obtained by averaging the changes of the respective stroke types for all five positions. The position of the stroke had only a small influence on the signal. The largest changes were caused by posterior strokes, and were approximately 1.5 times larger than the smallest changes, which were caused by central strokes. Note: the first value corresponds to haemorrhage, the second to ischaemia. The maximum score for each stroke type is 16, which means that all combinations of voltage simulation and image reconstruction meshes resulted in a good image. 0-100% image error were considered successful and everything above unsuccessful. If all images had an error of less than 120% the number is black, otherwise red.

How large were voltage errors introduced by wrong meshes?
Baseline voltage differences were compared between all four head meshes, and the average for these six comparisons was computed. The average voltage differences were large with 2-norms between 16-35 mV (59-100% proportional). Time-difference voltage changes were compared between all four head meshes for perturbations simulated in posterior position. The average TD voltage differences were 146 μV (0.63% proportional) for haemorrhage and 13 μV (0.07% proportional) for ischaemia.

How did this compare to the introduced system noise?
The average 2-norm of the added proportional and additive noise was 28 μV (5.6% proportional) for low noise, 59 μV (11.3% proportional) for medium noise and 148 μV (28.3% proportional) for the highest used noise level. These values were computed based on the noise added to the baseline. Consequently, time-difference data contained a combination of the noise on the baseline and equivalent noise on the perturbation voltages.

Comparison of the voltage differences.
In time-difference imaging, the geometrical errors were in the same range as the signal. On the other hand, the proportional and additive noise added to the simulated voltages was 2-10 times smaller than the haemorrhagic signal. This means that reconstructions of haemorrhagic strokes had a ten times better signal to system noise ratio than reconstructions of ischaemic strokes, while the signal to mesh difference ratios were comparable (figure 6).
Successful reconstructions of ischaemic strokes reduced strongly between medium and high system noise. The 2-norm of the signal was two times smaller than the 2-norm of the medium noise, suggesting this ratio as a good rule of thumb for estimating feasibility of head EIT applications.

Mesh creation
The created head meshes were anatomically more accurate than previous EIT head models, due to the use of both CT and MRI scans, instead of just one of the two (Tizzard et al 2005, Vonach et al 2012). The quality of the meshes was further improved by using the established open source Cgal mesh generator directly on the segmentation, thereby skipping previously required separate steps of enclosed 3D surface generation and meshing (Vonach et al 2012). Problem specific element sizing fields were defined to reduce element size near the electrodes, where electric potential gradients were large. Once the manual CT and MRI segmentations were performed, meshes with variable element sizes could be efficiently created by simply varying the Cgal settings in a text file.

Impact of mesh differences
Voltage differences between meshes, voltage changes caused by the simulated perturbations and the degree of added system noise were quantified by the 2-norm. The image quality of the reconstructions was assessed objectively by predefined metrics and a limit of detectability was set at 100% image errors. From these objective measures, the following answers were found to the questions posed at the beginning: (1) How large are voltage and image errors introduced by using a generic head mesh, as opposed to a subject specific one? Time-difference voltage errors between head meshes were of the same order than voltage changes caused by the perturbations. The image quality was on average 25% better, when the reconstructions were made on the correct mesh. When the added system noise was more than twice as high as the signal, then the stroke detection rate decreased significantly.
(2) Is a subject specific head mesh required for the detection of localised conductivity changes? While using the correct mesh for reconstructions improved the overall image The plot on the left shows the effects of using a generic head mesh for reconstructions. The signal to mesh difference ratio was computed as the 2-norm of the voltage changes caused by the simulated stroke, divided by the 2-norm of the difference in time-difference voltage changes on the mesh used for voltage simulation and the mesh used for reconstruction. The plot on the right shows the image quality with respect to signal to system noise ratio for both types of stroke. For each stroke type there are 240 data points, corresponding to four meshes used for simulation, four meshes used for reconstruction, five stroke positions and three noise levels. However, the signal to mesh difference ratio was infinite if the same mesh was used for simulation and reconstruction, which is why only 180 data points for each stroke type are visible on the left plot. quality, stroke detection rates were not improved. In three cases, only the correct mesh lead to a correct stroke identification, while there were five cases, where the reconstruction on the correct mesh failed and the average image errors on the other meshes were below 100%.
Finally, it was observed that even a coarse, homogeneous reconstruction mesh resulted in an acceptable image quality in most cases.

Conclusions
For monitoring, it was found not to be essential to use patient specific head meshes for image reconstructions. The results of this study, however, suggest that the noise on the measured voltages has to be kept very low to enable reliable detection of small conductivity changes in the head, such as ischaemia, epileptic seizures or even fast neural activity. Ways in achieving this may include long averaging, correction of small electrode movements (Jehl et al 2015b) and reduction of system errors (Fitzgerald et al 2002).
An EIT head imaging application in which it is very unlikely to have access to patient specific head meshes, is the differentiation of acute ischaemic and haemorrhagic strokes (Holder and Tidswell 2004). Since a baseline measurement is not available in this case, time-difference (TD) reconstructions are not available and multi-frequency (MF) methods have to be used (Malone et al 2014a). Previous results and preliminary work suggest that MF reconstructions are more sensitive to geometrical errors than TD (Malone et al 2014b), but more analysis needs to be done to give a definitive answer, whether patient-specific meshes are required in MF imaging.
This work suggests that monitoring the development of cerebral ischaemia or haemorrhage over time may be possible. This could be of clinical benefit following stroke thrombolysis, sub-arachnoid haemorrhage or head injury. For such studies, our recommendation based on this work is: construction of the Jacobian matrix with an accurate three million element mesh of a generic head and reconstruction on either the same mesh or a coarser hexahedral mesh.