Machine learning enhanced electrical impedance tomography for 2D materials

In the forward problem, equation (1) is used with suitable boundary conditions to simulate the potential for a given conductivity distribution on the sample. The CEM solves equation (1) using a form of the Robin boundary condition applied to EIT:

$$\begin{equation}{V}_{n}=u+{z}_{n}\sigma \frac{\partial u}{\partial \mathbf{n}}\quad \text{for}\enspace x\in {e}_{n},\end{equation} \tag{ 2 }$$

where z_n is the contact impedance between the sample and the nth electrode and V_n is the voltage reading on the electrode, while e_n represents the position of the nth electrode. x is a spatial variable of the conductive surface, Ω, and n is the normal to the boundary, ∂Ω. The model also enforces the current density to be zero on the boundary except at the electrodes [40]:

$$\begin{equation}\mathbf{j}=-\sigma \nabla u=0\quad \text{for}\enspace x\in \partial {\Omega}\quad \text{and}\enspace x\notin {e}_{n}.\end{equation} \tag{ 3 }$$

The complexity of the second-order differential equation (equation (1)) that is solved in EIT requires numerical methods to approximate the solution. One of the methods for the solution of the forward problem is the finite element method (FEM) that is based on the Galerkin method [41]. Firstly, we consider the surface of interest, Ω, and its boundary, ∂Ω. Ω is discretised by creating a mesh of N triangles which form the elements of the model. The vertices formed by the mesh are known as nodes. Electrode contact points correspond to labelled nodes on the boundary of the mesh, where each contact is modelled using multiple nodes to account for the finite width of the contacts. An example of the square mesh geometry used in this work, with 5 electrodes on each side giving 20 in total, is shown in figure 1. To represent all of the elements of the mesh, an N × N admittance matrix, A, is assembled such that the electric field is the same at the nodes where different triangles are connected. This admittance matrix is then used to form the linear system of equations

$$\begin{equation}\mathbf{A}\mathbf{f}=\mathbf{b},\end{equation} \tag{ 4 }$$

where f is a vector of the potential in each element and b is a vector that gives the total flux of the current at each node. For details of the CEM and a full derivation of the admittance matrix, see [42, 43] and [44]. Once the forward model is established, the Jacobian, J, may be calculated. Its elements are defined as

$$\begin{equation}{J}_{\mathit{ijk}}=\frac{\partial {V}_{ij}}{\partial {\sigma }_{k}}=-{\int }_{{{\Omega}}_{k}}\nabla {u}_{i}\cdot \nabla {u}_{j}\enspace \mathrm{d}s,\end{equation} \tag{ 5 }$$

where V_ij is the voltage measured between electrodes i and j, σ_k is the conductivity of element k in the forward model, and the integral is over the kth element's surface, Ω_k . The Jacobian is a measure which quantifies the sensitivity of the conductivity distribution to the value of the voltages measured along the sample boundary. See the appendix of reference [45] for a full derivation of the Jacobian. This variable is essential for the inverse problem and the electrode selection algorithm (ESA) discussed in the following sections.

Zoom In Zoom Out Reset image size

The inverse problem is non-linear due to the fact that current does not travel in straight lines. More current passes through regions of lower impedance, which in turn makes the current's trajectory non-uniform and non-linear. Accordingly, simple linear models such as back projection (used in CT scans where x-rays are exploited) cannot capture the complexity of the current flow. In contrast to the forward problem, where equation (1) is solved for the potential, u, for an assumed conductivity distribution, the inverse problem serves to fit a conductivity distribution that matches measured boundary voltages. Since physically measuring voltages only samples the value of the potential at discrete points along the boundary ∂Ω, the given input information is incomplete. Subsequently, the inverse problem is by definition ill-posed and ill-conditioned. As such, the EIT inverse problem requires a form of regularisation to obtain stable solutions [46, 47]. One widely used linear numerical method is called GREIT (Graz consensus Reconstruction algorithm for EIT), which relies on Tikhonov regularisation [3]. GREIT and most of the other inverse solver algorithms utilise the Jacobian matrix (equation (5)), derived in [45, 48].

The GREIT algorithm is a difference EIT inverse solver, which reconstructs the vector of conductivity change, x = σ − σ _r , where σ and σ _r are the sample and reference conductivities, respectively, and x is of length m. GREIT is a linear approximation, which makes use of a reconstruction matrix, R, and the vector of voltage measurement differences, Δ V , to estimate x,

$$\begin{equation}\mathbf{x}=\mathbf{R}{\Delta}\boldsymbol{V}=\mathbf{R}(\boldsymbol{V}-{\boldsymbol{V}}_{\mathrm{r}}),\end{equation} \tag{ 6 }$$

where V _r is the vector of reference voltages obtained from the forward model using σ _r . GREIT acknowledges that the resolution of EIT is limited and therefore small features in x will be blurred. So, rather than attempting to reconstruct x exactly, GREIT attempts to construct a 'desired image' $\tilde{\mathbf{x}}$,

$$\begin{equation}\tilde{\mathbf{x}}=\mathbf{D}\mathbf{x},\end{equation} \tag{ 7 }$$

where D is a matrix which represents a blurring operation on x. This is a one-to-one mapping, where each $\tilde{\mathbf{x}}$ element is centred on its associated x element. However, each $\tilde{\mathbf{x}}$ is defined to have a larger spatial area than its corresponding x, which gradually decreases to zero at the outer boundary of the $\tilde{\mathbf{x}}$ area. This larger 'desired image' smooths out differences in the resolution between the boundary and the centre, giving a more uniform distribution overall by exploiting the blurry nature of EIT. The original GREIT paper considered circular anomaly targets, where the desired image was a circle with a slightly larger radius [3]. In this work, D was of size m by m and based on a uniform weighting matrix w, which utilised a sigmoid function as the blurring operator; this was done using the GREIT algorithm as it is implemented in pyEIT [1].

2.3.1. Conventional algorithms

2.3.2. Machine learning adaptive electrode selection algorithm (A-ESA)

In this work, a novel A-ESA was developed to use machine learning and representative training data to achieve a more efficient measurement selection procedure. The current implementation of the A-ESA relies solely on the GREIT reconstruction image, and subsequently does not require information on the sample geometry. However, it does make two assumptions about the nature of EIT and the data available:

• (a)  
  Each new voltage measurement provides information gain, but may increase the total error associated with matrix R, which is the parameter to be minimized for reconstruction [3].
• (b)  
  The GREIT algorithm is more likely to neglect anomalies which are close to a previously detected deviation, since the gradient of the conductivity is larger at the deviation.

The A-ESA aims to identify zones of interest which have a high probability of containing a deviation from the current estimate of the conductivity reconstruction image, σ , which is in a n_x × n_y pixel shape. The closer one gets to the boundary of a detected anomaly, the more the conductivity gradient increases. One can establish a deviation map, which is a combination of the gradient map and conductivity map logarithm. The gradient map, G, for a 2D distribution is given by

$$\begin{equation}{G}_{ij}=\sqrt{{\left(\frac{{\sigma }_{i+1,j}-{\sigma }_{i-1,j}}{2}\right)}^{2}+{\left(\frac{{\sigma }_{i,j+1}-{\sigma }_{i,j-1}}{2}\right)}^{2}},\end{equation} \tag{ 8 }$$

where i and j run over the x and y pixel elements of σ_ij , and all contributions of second order and above have been neglected. The logarithmic map, L, is the logarithm of the conductivity

$$\begin{equation}{L}_{ij}=\mathrm{log}({\sigma }_{ij}),\end{equation} \tag{ 9 }$$

which enables the deviation map to be calculated using

$$\begin{equation}\mathbf{D}=\mathbf{G}\odot {\mathbf{L}}^{\diamond {q}_{\mathrm{h}}}.\end{equation} \tag{ 10 }$$

⊙ represents a Hadamard product (element-wise), while the superscript ◊ represents an element-wise power. q_h denotes a hyperparameter intrinsic to the A-ESA, which must be pre-determined by Bayesian optimization and training on simulated data [55].

Next, the influence map S is obtained by applying perturbations to all previously measured voltages, which are stored in a vector V of size m. A constant perturbing matrix P, is used to perturb all voltages at once:

$$\begin{equation}\mathbf{P}=\mathbf{1}+{p}_{\mathrm{h}}\mathbf{I},\end{equation} \tag{ 11 }$$

$$\begin{equation}{\mathbf{V}}_{\text{pert}}=\mathbf{P}\cdot \mathrm{diag}(\boldsymbol{V})\end{equation} \tag{ 12 }$$

where 1 is an m by m matrix of ones, I is an m by m identity matrix, and p_h is the perturbation constant hyperparameter. V_pert (m × m) is a matrix of perturbed voltages, where each row is equal to the original voltage vector, V , except for one perturbed value. That is, in the kth row, ${\mathbf{V}}_{\text{pert}}^{(k)}$, only the kth value is perturbed by a factor of p_h,

$$\begin{equation}{\mathbf{V}}_{\text{pert}}^{(k)}=\left\{{V}_{1},{V}_{2},\dots ,{V}_{k}+{p}_{\mathrm{h}}{V}_{k},\dots ,{V}_{\mathrm{m}}\right\}.\end{equation} \tag{ 13 }$$

Vectors of the difference between the perturbed and reference voltages are found for each row of V_pert,

$$\begin{equation}{\Delta}{\boldsymbol{V}}_{\text{pert}}^{(k)}={\mathbf{V}}_{\text{pert}}^{(k)}-{\boldsymbol{V}}_{\mathrm{r}}=\left\{{\Delta}{V}_{1},{\Delta}{V}_{2},\dots ,{\Delta}{V}_{k}+{p}_{\mathrm{h}}{V}_{k},\dots ,{\Delta}{V}_{\mathrm{m}}\right\},\end{equation} \tag{ 14 }$$

where ${\Delta}{\boldsymbol{V}}_{\text{pert}}^{(k)}$ is the kth difference vector. Since there are m different ${\Delta}{\boldsymbol{V}}_{\text{pert}}^{(k)}$ vectors, one for each k-value, m conductivity deviation vectors are reconstructed using GREIT,

$$\begin{equation}{\mathbf{\Delta }\mathbf{\sigma }}^{(k)}=\mathbf{R}{\Delta}{\boldsymbol{V}}_{\text{pert}}^{(k)},\end{equation} \tag{ 15 }$$

$$\begin{align}\hfill \,{\Rightarrow}\,{\Delta}{\sigma }_{i}^{(k)}& =\sum\limits _{j=1}^{m}\,{R}_{ij}{\Delta}{V}_{pert,\ j}^{(k)}\hfill \\ \hfill & ={R}_{i1}{\Delta}{V}_{1}+\cdots +{R}_{ik}({\Delta}{V}_{k}+{p}_{\mathrm{h}}{V}_{k})+\cdots +{R}_{im}{\Delta}{V}_{\mathrm{m}},\hfill \end{align} \tag{ 16 }$$

where R is the reconstruction matrix from equation (6). After this stage, the influence map is created by summing all $m\ {\Delta}{\sigma }_{i}^{(k)}$ vectors along the kth axis. Then the resulting vector is reshaped into a matrix of n_x × n_y , which are the horizontal and vertical spatial coordinates of the GREIT pixel grid:

$$\begin{equation}{S}_{i}=\sum\limits _{k=1}^{m}\vert {\Delta}{\sigma }_{i}^{(k)}\vert ,\end{equation} \tag{ 17 }$$

$$\begin{equation}\boldsymbol{S}\in {\mathbb{R}}^{m={n}_{x}{n}_{y}}\stackrel{\text{Reshape}}{\to }\mathbf{S}\in {\mathbb{R}}^{{n}_{x}\times {n}_{y}}.\end{equation} \tag{ 18 }$$

The Hadamard product of the deviation and influence map is taken to produce the n_x × n_y total map,

$$\begin{equation}\mathbf{T}=\mathbf{D}\odot {\mathbf{S}}^{\diamond {r}_{\mathrm{h}}},\end{equation} \tag{ 19 }$$

$$\begin{equation}{T}_{ij}={G}_{ij}{L}_{ij}^{{q}_{\mathrm{h}}}{S}_{ij}^{{r}_{\mathrm{h}}}={G}_{ij}\,\mathrm{log}{({\sigma }_{ij})}^{{q}_{\mathrm{h}}}{\left(\sum\limits _{k=0}^{m}\vert {\Delta}{\sigma }_{ij}^{(k)}\vert \right)}^{{r}_{\mathrm{h}}},\end{equation} \tag{ 20 }$$

where r_h is a third hyperparameter.

Therefore, the total map T is a measure of which pixels are the best for measurements in terms of gaining information about the overall conductivity distribution. T is constructed by combining the deviation map D (the variation of conductivity gradient across the pixels of the previous GREIT reconstruction) and the influence map S (sensitivity of the conductivity of certain pixels to deviations in measured voltage). The total map thus considers both the pixel sensitivity to voltage change and the gradient of pixel conductivity.

Zoom In Zoom Out Reset image size

2.3.3. A-ESA current and voltage pair selection

After the total map has been calculated, the A-ESA uses it to identify all possible current pairs and corresponding voltage pairs that offer the most information gain on the conductivity distribution. The A-ESA first traverses all combinations of current pairs, then selects an optimal subset using the total map value sum for each pair. Then for each optimal current pair, it traverses all voltage pairs to select an optimal subset of voltage pairs for that current pair. This is done by considering a straight line between every possible pair of current electrodes and evaluating the sum of the total map values for the pixels crossed by the line, as shown by the illustration of a pixel grid in figure 2(a). The current pairs are ranked according to their 'crossed pixel total map sum', current pairs with the highest sum of total map values are most optimal. Subsequently, the single best current pair, or a subset of the top pairs from the ranking, are the first to be considered. After an optimal current pair has been identified, all possible voltage pairs associated with the selected current pair will be evaluated in turn to obtain a subset of optimal voltage pairs. This stage utilises the Jacobian matrix, which was calculated when solving the forward problem. The Jacobian is a measure of how sensitive each voltage pair is to the conductivity of each pixel, whereas the total map is the metric of how much information can be gained by measuring the conductivity over each pixel. Each voltage pair has a unique set of Jacobian elements, where each element corresponds to a particular conductivity pixel and associated total map value. As shown in figure 2(b), regions of interest are identified by considering pixels whose total map value is above a given percentage of the maximum total map value;

$$\begin{equation}{T}_{ij}\enspace > \enspace {t}_{\mathrm{h}}\enspace \times \enspace \mathrm{max}(\mathbf{T}),\end{equation} \tag{ 21 }$$

where t_h is a hyperparameter bound between 0 and 1. For each voltage measurement, V_a , its corresponding Jacobian elements associated with the regions of interest are summed

$$\begin{equation}\sum\limits _{b}\,{J}_{ab}=\sum\limits _{b}\frac{\partial {V}_{a}}{\partial {\sigma }_{b}},\end{equation} \tag{ 22 }$$

where σ_b is the conductivity value corresponding to pixel b of the pixels which make up the region of interest. Voltage pairs that are the most sensitive to the conductivity in the region of interest offer the most information and potential improvement in the next reconstruction. The subset of voltages with the highest sum of Jacobian elements are considered to be the most sensitive, hence the most effective measuring pairs. The described current and voltage pair selection offers an optimised set of four-terminal measurements to be performed on the sample, which will subsequently give a better GREIT reconstruction. The process repeats iteratively until some user-defined stopping condition, such as after 10 iterations of 50 measurements each, or a desired accuracy is achieved.

One far-reaching application of machine learning is the use of CNNs for image post-processing. CNNs have proven to be effective at filtering images for feature extraction and noise suppression, resulting in higher clarity images [56], and have recently seen application in EIT [16, 18, 20, 24]. One such CNN commonly used for image segmentation is known as U-Net, which was designed for use in neuron imaging and cell tracking [39]. The deep D-Bar CNN architecture is shown in figure 3, which was adapted from [9] for this work to use GREIT instead of D-Bar reconstructions. CNNs utilise a series of convolutional layers that transform the input into a latent space, which achieves a more compact representation of features in the data. This is performed using element-wise multiplication between a 3 × 3 pixel convolutional filter and a region of the image of the same size. By scanning the filter sequentially across the image, the filters act to average regions of pixels together, which are then pooled and combined. A typical feature of this model is the skipping of multiple layers that gives the architecture's floor-like structure [9]. U-Net utilises the common rectified linear activation function to perform discrete convolutions on an input image. For a full description of how the U-Net CNN works, see [9, 39].

Zoom In Zoom Out Reset image size

In this work, we used the quadratic loss function, $\mathcal{L}$, as the measure of error between the true conductivity map and the output produced by the deep learning model via the back propagation process:

$$\begin{equation}\mathcal{L}=\sum\limits _{i,j}{\left({\sigma }_{ij}^{\text{True}}-{\sigma }_{ij}^{\text{Pred}}\right)}^{2},\end{equation} \tag{ 23 }$$

where σ ^True is the true conductivity map and σ ^Pred is the prediction output of the CNN. We modified the standard U-Net architecture by introducing zero padding, achieved by adding a frame of zeros to the result of the convolution at each layer. This ensured the preservation of the size of the input image as it propagated through the network to the output. The neural network structure was written using the Keras library and Tensorflow back-end [57, 58]. Note that the same loss measure was also used to evaluate the forward solvers and ESAs, in which case σ ^Pred was the raw reconstruction of the GREIT algorithm.

The newly implemented forward solver is the first forward solver that utilises the CEM and is written in Python using the external Numpy and CuPy (for GPU acceleration) libraries [59, 60]. The code is completely vectorised to reduce the execution time and uses the 2D FEM to assemble an admittance matrix from three sub-matrices [44]. The sub-matrices were calculated prior to the assembly by using integrating matrices that deliver the corresponding integral result to each matrix element. This CEM-enhanced forward solver models the width of the electrodes on the sample and the contact impedance, which was the main difference between this solver and the simplified model used in the original pyEIT package [1]. To evaluate their relative performance, both solvers were used to generate GREIT reconstructions for the same test set of conductivity distributions.

The GREIT algorithm was evaluated on its ability to produce reconstructions which identify the main features of the conductivity distribution. Here, sample conductivity distributions were generated artificially, including both continuous and discontinuous anomalies, as well as their combinations. The continuous anomalies were Gaussian distributions, which imitates the presence of nearby point charges, while the discontinuous anomalies were triangles, ellipses and lines of uniform conductivity. 1 to 5 anomalies of random size were used per sample and were randomly positioned on the sample. The triangles and ellipses had a normally distributed conductivity value with a mean at the reference conductivity, whereas the lines had a conductivity at 0.001% of the reference, to simulate an electrical discontinuity or cut through the sample. Synthetic Gaussian noise with a standard deviation of 3% was added to each simulated voltage measurement. Since the samples were simulated, their true conductivity distributions were known and were subsequently used to calculate the $\mathcal{L}$ losses of the corresponding GREIT reconstructions.

The A-ESA was written in Python, in a fully vectorised form, using the NumPy library [59] and the pyEIT implementation of the GREIT algorithm [1, 3]. Before use, the A-ESA's four hyperparameters were optimised on simulated training data, where 1000 samples with different conductivity distributions were used. Each conductivity distribution was reconstructed using 15 iterations of 10 voltage measurements and their losses calculated. 150 measurements per sample for 1000 samples resulted in significant computation time for a single evaluation of the objective function. Bayesian optimisation is a technique commonly used to sample time-expensive objective functions [55]. As such, this method was applied to the A-ESA for 100 calls, where each evaluation used different values of the hyperparameters. The A-ESA can be represented in an objective function form by

$$\begin{equation}{f}_{i}^{(k)}=\frac{1}{{n}_{i}^{(k)}}\frac{{\mathcal{L}}_{i+1}^{(k)}-{\mathcal{L}}_{i}^{(k)}}{{n}_{i+1}^{(k)}-{n}_{i}^{(k)}},\end{equation} \tag{ 24 }$$

where ${f}_{i}^{(k)}$ represents the ith component of the objective function vector for the kth training sample, which measures the speed of convergence towards the true conductivity map. ${n}_{i}^{(k)}$ is the number of performed voltage measurements before the ith loss value, ${\mathcal{L}}_{i}^{(k)}$, is calculated. The index i runs from 1 to N_m/10, where N_m is the total number of measurements (150 in this case). The sample number k runs from 1 to 1000 and to get the objective function value for a sample k, the elements of f^(k) are summed. Bayesian optimisation evaluates the error loss for different values of the hyperparameters, quickly finding an optimal set with a minimised loss. The final normalised loss value for the optimised A-ESA is given by

$$\begin{equation}\mathbf{F}=\frac{1}{1000}\sum\limits _{k=1}^{1000}\,\sum\limits _{i=1}^{{N}_{\mathrm{m}}/10}\,{f}_{i}^{(k)}.\end{equation} \tag{ 25 }$$

The Bayesian search was performed on continuous conductivity samples. The bounds on the hyperparameters were chosen according to the mathematical limits of the parameters and some initial intuitions gained from small, trial optimisation runs. After several Bayesian searches over 1000 different samples, a final set of optimal hyperparameters was identified, as shown in table 1.

After the training process, these optimised hyperparameters were used in the A-ESA for its evaluation against two standard measuring techniques: the opposite–adjacent ESA (opposite electrode current driver with adjacent electrodes as the voltage measurement pairs) and the adjacent–adjacent ESA.

One of the motivations of this research was to determine whether EIT is suitable to detect charges close to the surface of an examined sample. As such, point charges were modelled using Gaussian-distributed spots of conductivity with added uncorrelated noise. The U-Net CNN was used for post-processing of the GREIT conductivity reconstructions, to see whether the neural network was able to improve upon the reconstruction accuracy. To enable this the CNN was trained and tested on a total of 130 000 samples, where each of their losses were minimised during the training process using the Adam optimiser [61]. To train, sample reconstructions were input into the network, processed, and the network output evaluated against the true conductivity map of the sample using equation (23).

The training samples were generated in the same way as outlined in section 3.1. GREIT reconstructions were made for each of the samples using an ordered electrode selection. First, the number of source-sink pairs was set to 15. 15 sources were taken from a list, ordered by the remainder of their division by 5 (i.e. 0, 5, 10, 15, 1, 6, 11, 16, 2, 7, 12, 17, 3, 8, 13, 18, 4, 9, 14, 19). The 0th electrode was at the leftmost edge of one of the square sample sides and the rest were ordered clockwise. Next, a random fixed electrode distance between the source-sink pairs was used to select the corresponding sink electrode. This process was used to spread the source-sink pairs around the sample, enabling a more consistent reconstruction. A limited number of voltage pairs were selected randomly, using a random distance between the voltage pair electrodes for each sample. Any measurement lines where sources or sinks crossed over were discarded. Subsequently, for each sample, an array of current and voltage pairs without repeated lines was created. This selection algorithm allowed for randomisation in the selection of current–voltage pairs, while retaining consistency of reconstructions by spreading measurements around the training sample. In order to improve the similarity to real EIT data, white Gaussian noise with a mean of zero was added to the simulated voltage readings. 60 000 samples of mixed geometric and Gaussian type anomalies were produced for training, 10 000 more were used for testing of the CNN. A further 40 000 reconstructions of samples with only a single Gaussian anomaly were also used to train the network.

An initial test was performed only on the geometric type anomalies. The CNN processing was applied to qualitatively evaluate how well the CNN could reproduce the sharp, discontinuous edges of these anomalies. As EIT is inherently blurry, there is some utility in post-processing procedures that can reduce this. Next, to quantify the performance of the CNN, sensitivity tests were devised, which featured uncorrelated Gaussian anomalies of 20 different amplitudes and 10 different standard deviations, giving 200 different types of anomalies in total. The amplitudes were uniformly distributed from −1 to 1, at intervals of 0.1, while the standard deviations varied between 5% to 50% of the side of the square sample, at intervals of 5%, where all values were fractional deviations from the reference conductivity. Each type of Gaussian anomaly was simulated with their means at differing spatial positions within a 10 × 10 grid format, where the x–y coordinates varied between ±0.9, as shown in figure 4(c).

Zoom In Zoom Out Reset image size

This resulted in a total of 20 000 conductivity distributions for the CNN's testing set. For each anomaly of particular amplitude and standard deviation, their reconstruction losses were averaged over the 100 different mean locations to obtain an average $\mathcal{L}$ loss. The deviation in position between the Gaussian anomaly mean for the true, raw and processed maps was also calculated, which gave a direct measure of how similar the raw and processed maps were to the true data. Once the raw and processed reconstructions were made for each anomaly type, their averaged losses and errors in mean position were plotted against the Gaussian's amplitude, with different plots for each of the 10 different standard deviations. At the end of simulated data generation there was a total of 130 000 true maps and GREIT reconstructions that were used for training and testing of the CNN, which included geometrically shaped and Gaussian type anomalies. This training data allowed evaluation of the CNN on its ability to process the GREIT reconstructions and obtain an image closer to the true conductivity map.

The comparison between the new, CEM-enhanced FEM solver and the original simplified version found in the pyEIT package involved evaluation of the losses $(\mathcal{L})$ of the conductivity reconstructions for 10 000 different samples. A visual comparison of the reconstructions between the new and the old solvers for three different samples is shown in figure 5. From left to right, the first three images show the true conductivity map (a), the old forward solver reconstruction (b), and the new CEM-enhanced forward solver reconstruction (c). Panel (a) shows that the sample was composed of a triangular, elliptical, and linear anomaly. Visually comparing (b) and (c) shows that both the old and new solvers performed to the same degree. However, the CEM-solver was able to better capture the conductivity with more consistency than the old solver. In images (d)–(f) the new solver performed better than the old solver, resolving more of the true triangle anomaly shape. In images (g)–(i), both solvers distinguished the two anomalies, capturing their approximate shapes and the surrounding conductivity variation. However, in this case the old solver performed slightly better at capturing the shape of the elliptical and triangular anomalies than the new solver. Overall, the old and new solver reconstructions only show minor differences under visual comparison. However, the calculated loss of the new solver when averaged over 10 000 samples was 3% lower than that of the old simplified forward solver.

Zoom In Zoom Out Reset image size

To evaluate the A-ESA, its reconstruction error loss was compared to the losses achieved by two commonly-used ESAs: the opposite–adjacent and adjacent–adjacent measurement protocols. Initially a direct comparison between all three methods was performed on two different samples, where their losses were recorded over 120 measurements in steps of 10, as shown in figure 6. Panels (a) and (b) show the reconstruction loss comparison for the two samples, where the losses are normalised by the lowest error achieved by the A-ESA. Both plots show that the A-ESA achieved a consistently lower loss than the opposite–adjacent method before both methods converged beyond 100 measurements. The adjacent–adjacent method clearly performed very poorly, with a final loss of approximately three times that of the A-ESA in either case.

Zoom In Zoom Out Reset image size

Next, the different ESAs were evaluated over 2000 different samples. This was repeated over three runs using the same samples but re-initialising the forward solver mesh each time. For each run the 2000 sample losses were summed and normalised, after which the three data sets were averaged and plotted, as shown in figure 7. For all 120 measurements the averaged $\mathcal{L}$ losses of the A-ESA reconstructions were lower than the losses of the opposite–adjacent protocol and significantly lower than the adjacent–adjacent method. This difference in performance was even more apparent in the region of 40 to 60 voltage measurements, where the loss of the A-ESA was approximately 75% the loss of the opposite–adjacent method. The decrease in the loss was steeper at fewer measurements for both the A-ESA and Opp–Adj method. However, the A-ESA loss gradient was slightly larger up until the 80th measurement, after which both methods became comparable and the final loss of the Opp–Adj was only 12% higher than the A-ESA. The Adj–Adj method performed poorly, never achieving anything lower than 3.5 times the final loss of the A-ESA. Despite the eventual loss convergence the of Opp–Adj and A-ESA methods, by the 60th measurement the A-ESA loss was within 20% of the final loss, whereas the Opp–Adj was only within 60%. The goal of the A-ESA development was not to get a lower loss overall, but to achieve some desired loss level faster than conventional ESAs. Figure 7 clearly shows this goal was achieved, with the extra benefit of lower overall losses as well.

Zoom In Zoom Out Reset image size

The CNN was applied to a variety of GREIT reconstructions, examples of which are shown in figure 8. In general, the CNN post-processing was able to offer significant improvement over raw GREIT reconstructions. The extremely high clarity and accuracy may partially be the result of some over-fitting of the training data. Nevertheless, the analysis demonstrates the utility of machine learning neural networks for post-processing of the reconstruction images: whether for improving image clarity, as in this work; or for classification of specific anomalies that may occur. A focused and thorough analysis would enable a more rigorous training and testing procedure to be developed, along with a more developed CNN for use in EIT applications.

Zoom In Zoom Out Reset image size

The sensitivity of the GREIT algorithm was tested on samples with continuous conductivity distributions by simulating Gaussian anomalies at 100 different positions on the sample, and averaging the losses of their reconstructions. The losses of the reconstruction and the errors in the estimation of the mean position were plotted against the amplitudes of the anomalies, as a function of deviation from the reference conductivity. The graphs from before and after CNN post-processing for a Gaussian anomaly with a standard deviation of 0.5 (25% of the square side length) are illustrated in figure 9. An example of the true image, reconstruction and processed reconstruction are shown in (a)–(c) for the positive amplitude of 0.25, with a mean at the slightly off-centre position of (0.1, 0.1). These show that the raw reconstruction reduced the true Gaussian anomaly peak to a much smaller region of high intensity, failing to retain the broad profile of the true conductivity deviation. The CNN improved upon the raw reconstruction significantly, managing to reinstate the broad nature of the true Gaussian, but skewing it to the right slightly. Three images depicting an anomaly of the same amplitude but centered on (0.7, −0.7) are shown in (d)–(f). Again, the CNN processing retrieved the smooth conductivity gradient, although not without incorrectly drawing more intensity into the bottom right corner at (1.0, −1.0).

Zoom In Zoom Out Reset image size

The $\mathcal{L}$ losses of the raw and post-processed reconstructions shown in (g) are comparable at small amplitudes. However, for large amplitudes the raw reconstructions had a significant amount of error loss. The CNN processing improved significantly upon this, reducing the loss by almost 500% in the case of the extreme +0.5 amplitude point. As shown in (h), the deviation in Gaussian mean positions of the raw and processed reconstructions are almost identical, with very small deviation. However, the CNN data points are slightly lower by up to 2%. Despite the differences between the reconstructions and the true image, knowledge of the mean position was retained to the same degree between the raw GREIT maps and CNN processed maps.

The new, CEM-enhanced forward solver performed better than the original PyEIT forward solver, achieving a 3% lower averaged squared loss. The loss decrease with the new forward solver could be due to the addition of electrode width to the model, as the electrodes in the new solver encompassed a greater proportion of the sample boundary, thus capturing more information about the conductivity of the interior elements. The GPU acceleration from the CuPy package was implemented as part of the new forward solver enhancement [60], which resulted in an increased computation speed by up to five times that of the old solver, despite the increased model complexity.

The Bayesian optimisation successfully found an optimal set of hyperparameters that enabled the A-ESA to outperform the opposite–adjacent and adjacent–adjacent electrode selection methods. The simulation studies showed that the A-ESA achieved approximately 20% lower overall loss than the Opp–Adj method, and faster convergence to lower losses.

The A-ESA requires time to perform analysis of the total map between iterations of measurements. For example, for 10 sets of 100 measurements, the A-ESA will reconstruct the GREIT image 10 times. If the time taken to perform a set of electrical measurements is significantly faster than the operation of the A-ESA, then this will cause a computational bottleneck. The time taken for a single iteration of the A-ESA scales on the order of $\mathcal{O}({n}^{2})$ for the number of pixels used in the GREIT image. Generally for an increasing number of contacts, a greater number of potential measurements are available. Standard ESAs will increase in the number of required measurements, but are not flexible in which ones can be taken, whereas the A-ESA will also increase in total number of measurements, but is not restricted to any particular pattern. In the A-ESA the user can define the number of initial measurements to take before activating the A-ESA, the number of iterations (N_iter), and the number of measurements to perform each iteration (N_m). An appropriate choice of N_iter and N_m is situation- and application-specific; a rule-of-thumb applied in this work was to perform no more than half of the measurements that would otherwise be taken if the opposite–adjacent technique was used. The time increase to perform a single iteration of the A-ESA against contact number and number of pixels per row is shown in figure 13. Note, this is only the time taken for the total map analysis and pair selection, the time taken for electrical measurements was omitted. The overall complexity of the A-ESA can be taken as the pixel scaling and the product of N_iter and N_m, $\mathcal{O}({n}^{2})$ + N_iter N_m, where any extra scaling with contact number is accounted for by deciding the total number of measurements to perform.

Zoom In Zoom Out Reset image size

As this is the first attempt to apply ML to the selection of electrodes for EIT measurements, we foresee potential improvements that can be implemented in the future: (1) for samples with little variation, or if the conductivity deviations are comparable to the noise floor, the A-ESA is unlikely to select the same set of measurements each run. This may result in different reconstructions, limiting reconstruction reproducibility. This could be addressed with more varied samples or different types of meshes in the training data. (2) Training and optimisation of the A-ESA hyperparameters is computationally expensive. There are many parameters within EIT, such as the lambda regularisation, contact number and contact width, which could influence the reconstruction. On top of these, the A-ESA has four hyperparameters that determine certain thresholds of variables and raises matrices to powers element-wise. Introducing new variables, more rigorous training and better tuning of the hyperparameters would be beneficial, enabling a higher performance of the A-ESA. (3) The underlying deviation and total map theory was developed specifically for application to GREIT reconstructions. The relative simplicity of the operations underlying the A-ESA makes its improved performance over the opposite–adjacent method even more apparent. However, additional signal processing steps or the use of static inverse solvers, such as the Gauss–Newton method, can be applied to obtain more consistent solutions.

For most clinical or industrial applications, the number of electrodes used is small, hence the current speed of EIT is sufficient. In this work, we are extending the application of EIT to the domain of microfabricated devices based on 2D materials such as graphene. Using conventional microfabrication techniques, devices with hundreds of electrodes can be routinely made, which could enable novel types of microscopic charge or current imaging sensors. For these types of applications in nanotechnology and biomedicine, the time required for electrical measurements themselves are the limiting factor. This is particularly important for the time consuming electrical measurements of small signals in nano- or microscale devices, where the measurement time dominates over typical EIT computation times. Our work offers a solution to this problem: we keep the electrode count high (to allow for flexibility in measurements) but reduce the reconstruction time by only measuring a small subset (selected by the A-ESA) of all available electrode configurations. Therefore, by adopting the A-ESA method, reconstruction performance comparable to or better than that of the conventional methods can be achieved, but at a fraction of the measurement time.

The CNN post-processing of both the Gaussian and geometric anomalies showed an appreciable improvement over the raw reconstructions. The comparison of geometric anomalies in figure 8 shows that the CNN achieved an impressive degree of improvement over the raw GREIT reconstruction, demonstrating the benefits of machine learning and post-processing. However, this drastic improvement may partially be due to the similarity between the training data and testing data, despite introduction of synthetic noise and anomaly placement randomisation. In the Gaussian anomaly sensitivity tests, the CNN did not achieve such a radical improvement over the raw reconstructions. However, it retained more consistency as the Gaussian amplitudes increased, whereas the raw reconstruction deviated significantly. This is somewhat apparent in figures 9(a)–(f), where the CNN processing smoothed out some of the anomalous features introduced by EIT.