A MATLAB-Based Boundary Data Simulator for Studying the Resistivity Reconstruction Using Neighbouring Current Pattern

Abstract

Phantoms are essentially required to generate boundary data for studying the inverse solver performance in electrical impedance tomography (EIT). A MATLAB-based boundary data simulator (BDS) is developed to generate accurate boundary data using neighbouring current pattern for assessing the EIT inverse solvers. Domain diameter, inhomogeneity number, inhomogeneity geometry (shape, size, and position), background conductivity, and inhomogeneity conductivity are all set as BDS input variables. Different sets of boundary data are generated by changing the input variables of the BDS, and resistivity images are reconstructed using electrical impedance tomography and diffuse optical tomography reconstruction software (EIDORS). Results show that the BDS generates accurate boundary data for different types of single or multiple objects which are efficient enough to reconstruct the resistivity images for assessing the inverse solver. It is noticed that for the BDS with 2048 elements, the boundary data for all inhomogeneities with a diameter larger than 13.3% of that of the phantom are accurate enough to reconstruct the resistivity images in EIDORS-2D. By comparing the reconstructed image with an original geometry made in BDS, it would be easier to study the inverse solver performance and the origin of the boundary data error can be identified.

1. Introduction

Electrical impedance tomography (EIT) [1, 2] reconstructs the spatial distribution of electrical conductivity or resistivity of a closed conducting domain (Ω) from the surface potentials developed by a constant current injection through the surface electrodes surrounding the domain to be imaged. Before carrying out the practical measurements on patients, it is advised to test an EIT system with a tissue mimicking model of known properties [3] called practical phantoms [4–10]. Hence, phantoms are often required to assess the performance of EIT systems for their validation, calibration, and comparison purposes. Two-dimensional (2D) EIT (2D-EIT) assumes that the electrical current flows in a 2D space which is actually three-dimensional inside real volume conductors. Hence, the development of a perfect 2D practical phantom is a great challenge as the real electrodes always have a definite surface area, and hence the injected current signal cannot be confined in a 2D plane in bathing solution [5]. Researchers have developed a number of practical phantoms which are three-dimensional objects, and those phantoms are designed and developed, generally, for their own EIT systems. Practical phantoms containing electrolyte (or other conducting medium) [4–10] are three-dimensional in shape and hence they will have some data error due to the three dimensional current conduction. Also, the phantoms containing electrolytes (e.g., NaCl solution or saline) [5, 7, 8] are difficult to transport and are prone to errors since the evaporation of the water gives rise to changes in conductivity [9]. In addition, temperature variations have a marked effect on the conductivity because the temperature coefficient is large [11]. Therefore, the practical phantoms will have a poor stability and a gradually increasing data error over time. Network or mesh phantoms [12, 13] are compact, more stable, rugged, portable, easy to move, consistent over time, and less temperature dependent. But these phantoms need a huge number of identical electronic components properly designed in a mesh mimicking the conductivity distribution of a practical biological tissue. Furthermore, for a large tissue structure, a mesh phantom requires a huge number of very precision components. The reproduction of these kinds of phantoms having different properties is often time-consuming [14]. The option for changing the position and property of an inhomogeneity is limited by the phantom structure and the number of elements in mesh phantom but the practical phantoms allow us to put several types of object in different positions in the bathing solution, but they produce several errors contributing to the poor signal to noise ratio (SNR) in boundary data.

Reconstructed image quality in impedance tomography depends on the errors associated with practical phantom, electronic hardware, and inverse solver performance. Image quality is largely affected by the practical phantom design parameters such as phantom geometry, electrode geometry, electrode materials, and the nature and behavior of the inhomogeneity and bathing solution. SNR is also reduced by the error contributed by current injector, data acquisition system, and signal conditioner circuits. In practical phantoms, the voltage data developed by a three-dimensional current conduction are collected form surface electrodes connected to an analog instrumentation. Therefore, it is quite confusing to identify the source of the errors responsible for poor image quality in a 2D-EIT system. In order to overcome the difficulties and limitations of practical and mesh phantoms, a MATLAB-based boundary data simulator (BDS) is developed to generate accurate 2D boundary data for assessing the EIT inverse solvers. BDS is an absolute 2D data simulator which is required to generate the errorless 2D boundary data to study and modify the inverse solver of a 2D EIT system. As the BDS is a computer program, it is free from the instrumentation errors and allows us to generate voltage profile with different types of phantom geometry, inhomogeneity and background conductivity profile, and inhomogeneity geometry (shape, size, and position). Moreover, it is absolutely stable, compact, easy to use, and easy to handle and modify for further development. Boundary data for different phantom geometries are generated in BDS, and resistivity images are reconstructed in standard reconstruction algorithm. BDS is studied to conform its suitability to use for boundary data generation with different phantom configurations which are required to assess the EIT inverse solvers.

2. Methods

2.1. Mathematical Modelling of EIT

EIT image reconstruction is a nonlinear inverse problem [15] in which the electrical conductivity distribution of a closed domain (Ω) in a volume conductor is reconstructed from the surface potential data developed at the boundary (𝜕Ω) by injecting a constant current signal. A low frequency and low magnitude constant sinusoidal current is injected through an array of electrodes attached to the boundary, and the boundary potentials are measured using a data acquisition system. The voltage data collected from surface electrodes are then used by an image reconstruction algorithm [15] which reconstructs the conductivity distribution of the domain under test (DUT). The reconstruction algorithm computes the boundary potential for a known current injection and known conductivity values and tries to compute the conductivity distribution for which the difference between the measured boundary potential (𝑉_𝑚) and the calculated (𝑉_𝑐) is minimum. The reconstruction algorithm is developed with two parts: forward solver (FS) [5, 15–17] and inverse solver (IS) [15–17]. Forward solver calculates the boundary potential data for a known current injection and known conductivity values. Inverse solver computes the conductivity distribution for which the boundary voltage difference (Δ𝑉=𝑉_𝑚−𝑉_𝑐) becomes minimum.

The DUT will have the distinct conductivity values at each points defined by their corresponding coordinates (𝑥,𝑦). Due to a constant current injection, a potential profile is developed within DUT, and its potential profile without any internal energy sources depends on the conductivity profile. Hence, a relationship, called EIT governing equation, between the electrical conductivity (𝜎) of the points within the DUT and their corresponding potential values (Φ) can be established. The governing equation in EIT [1, 2] can be derived from the Maxwell’s equation and can be represented as ∇⋅𝜎∇Φ=0. (1)

To calculate the domain potential developed for a constant current injected to the DUT with a known conductivity distribution, the above equation is essentially to be solved. As the EIT governing equation is a nonlinear partial differential equation, the direct or analytical technique fails to solve it. Therefore, to calculate the domain potential, the equation is solved by developing a mathematical model called “forward model” which is derived from (1) using a numerical technique like finite element method (FEM) [18].

The EIT governing equation has an infinite number of solutions, and hence the FEM formulation of the EIT technique is essentially required to be provided by some boundary conditions [18–20] to restrict its solutions space. The boundary conditions are imposed into the FEM formulation of EIT by specifying the value of certain parameters (voltage or current). The parameters defining the boundary conditions may be either the potentials at the surface or the current density crossing the boundary or mixed conditions.

The boundary conditions, in which the parameters are the potential at the surface, are called the Dirichlet boundary conditions and are represented as [1, 5, 19, 20]Φ=Φ_𝑖, (2a) where 𝑖=1,…,𝑁 are the measured potentials on the electrodes.

The boundary conditions, in which the parameters are current density crossing the boundary, are known as the Neumann boundary conditions [1, 5, 19, 20] which are given by  _𝜕Ω𝜎𝜕Φ=𝜕𝑛+𝐼onthesourceelectrode−𝐼onthesinkelectrode0otherwise, (2b)

where 𝜕Ω is the boundary, and 𝑛 is the outward unit normal vector on an electrode surface.

In EIT, the FEM technique is used to derive the forward model from the governing equation in the form of a matrix equation establishing the relationship between the injected current and the developed potential within a DUT. The relationship can be assumed as the transfer function of the system which is mathematically represented as a matrix called global stiffness matrix (GSM) [18] or transformation matrix constructed with the elemental conductivities (𝜎) and nodal coordinates (𝑥,𝑦). In EIT, FEM discretizes the DUT by a finite element mesh containing finite number of elements of defined geometry and finite number of node. FEM applied on the governing equation to derive the forward model of a DUT in the form of a matrix equation using the 𝜎 and nodal coordinates. In the EIT forward model, the relationship established between the current injection matrix [𝐶] (matrix of the applied signal) and the nodal potential matrix [Φ] (matrix of the developed signal) through the transformation matrix [𝐾(𝜎)] is mathematically represented as [Φ]=[𝐾(𝜎)]⁻¹[𝐶]. (3)

Now, in FEM formulation in EIT, when the current matrices [𝐶] and [𝐾(𝜎)] are known, and the nodal potential matrix [Φ] is unknown, the forward model or the mathematical problem is termed as the “forward problem”. The procedure of calculating the [Φ] by solving the forward problem (3) with known [𝐾(𝜎)] and known [𝐶] is termed as “forward solution”. In EIT, the forward solver first computes the potential distribution with the assumed initial conductivity distribution (𝜎₀) with a known constant current simulation, and then the inverse solver reconstructs the conductivity distribution from the measured boundary potential data for a same constant current injection through surface electrodes. The EIT reconstruction algorithm tries to mathematically find the elemental conductivity values (conductivity distribution) for which the difference between the estimated nodal potentials (𝑉_𝑐) computed in the FS and the potentials measured (𝑉_𝑚) on the surface electrodes (for a same current injection values) becomes minimum.

The inverse solver of the EIT reconstruction algorithm is developed with a mathematical minimization algorithm (MMA) [19–22] such as Gauss-Newton-based mathematical minimization algorithm (GN-MMA). In GN-MMA, the conductivity update vector ([Δ𝜎]) is calculated and the boundary data mismatch vector (Δ𝑉=𝑉_𝑚−𝑉_𝑐) is minimized by an iteration technique like the modified Newton-Raphson iteration technique (NRIT) [19–22]. The [Δ𝜎] matrix is the desired variation in the elemental conductivity values in [𝜎] matrix for which the forward solver calculates the boundary potentials more similar to the measured value in next iteration using NRIT. Therefore, the algorithm starts with an initial elemental conductivity vector ([𝜎₀]), and it is then updated to ([𝜎₁]=[𝜎₀]+[Δ𝜎]) in the next iteration. Using this [𝜎₁], FS calculates a new potential distribution in DUT and a new voltage mismatch vector [Δ𝑉₁] is thus obtained and compared with the previous voltage mismatch vector [Δ𝑉₀]. If the Δ𝑉₁ is not found as the minimum, the iteration process is continued till the kth iteration using the conductivity update vector ([Δ𝜎_𝑘]) developed by GN-MMA. Using, NRIT the [𝜎] matrix is iteratively updated to [𝜎_𝑘+1]=[𝜎_𝑘]+[Δ𝜎_𝑘] and repetitively tries to find out the minimum value of [Δ𝑉].

Hence, in the EIT inverse solver, it is understood that the desired elemental conductivity matrix is obtained by a minimization algorithm (MMA) which is composed of Gauss-Newton method and Newton-Raphson iteration in which the technique iteratively tries to find out an optimum conductivity distribution [𝜎_𝑘] for which the voltage mismatch vector is minimized [Δ𝑉]. At a particular iteration in this MMA, the elemental conductivity matrix is calculated when the current matrices [𝐶] and [Φ] or [Δ𝑉=𝑉_𝑚−𝑉_𝑐] are known. This process is logically an opposite process to the forward problem. Thus, when the current matrices [𝐶] and [Φ] are known, and the elemental conductivity matrix [𝜎] is unknown, the model or the problem is called the “inverse problem.” The procedure of calculating the [𝜎] or [Δ𝜎] using with known [Δ𝑉] and the known [𝐶] is termed as “inverse solution.”

2.2. Image Reconstruction with GN-MMA and NRIT

Electrical conductivity imaging is a highly nonlinear and ill-posed inverse problem [19–22]. In EIT, a minimization algorithm is used to obtain an optimized elemental conductivity value [𝜎] for which the voltage mismatch vector [Δ𝑉] becomes minimum. In the image reconstruction process, the minimization algorithm [17, 18] first defines an objective function (𝑠) from the computational predicted data [𝑉_𝑐] and the experimental measured data [𝑉_𝑚] and runs iteratively to minimize it. Generally, in the EIT image reconstruction algorithm, the inverse solver searches for a least square solution of the minimized object the function (𝑠) using by a Gauss-Newton method and the NRIT-based iterative approximation techniques.

If 𝑓 is a function mapping a t-dimensional (t is the number of element in the FEM mesh) impedance distribution into a set of M (number of the experimental measurement data ([𝑉_𝑚]) available) approximate measured voltages, then the Gauss-Newton-method-based minimization algorithm [19–26] tries to find a least square solution of the minimized object function (s) [19–26]  which is defined as: 1𝑠=2𝑉_𝑚−𝑓²=12𝑉_𝑚−𝑓^𝑇𝑉_𝑚−𝑓. (4)

Now, differentiating (4) with respect to the conductivity 𝜎, it reduces to 𝑠^′=−𝑓^′𝑇𝑉_𝑚−𝑓=−𝐽^𝑇Δ𝑉, (5) where the matrix 𝐽=𝑓^′ is known as Jacobin matrix [19–22], which may be calculated by a method as described in [19, 22] or by the adjoint method [23] represented by (6) 𝐽=_Ω∇Φ_𝑠⋅∇Φ_𝑑𝑑Ω, (6) where Φ_𝑠 is the forward solution for a particular source location, and Φ_𝑑 is the forward solution for the adjoint source location (source at the detector location and detector at the source location).

Differentiating (5) with respect to 𝜎 again, the equation reduces to 𝑠^′′=𝑓^′𝑇𝑓^′−𝑓^′′𝑇𝑉_𝑚−𝑓. (7)

By Gauss-Newton method, the conductivity update vector [Δ𝜎] is given by 𝑠Δ𝜎=−^′𝑠^′′=𝐽^𝑇Δ𝑉𝑓^′𝑇𝑓^′−𝑓^′′𝑇.Δ𝑉 (8)

Thus, the conductivity update vector is given by Δ𝜎=𝑓^′𝑇𝑓^′−[𝐻]^𝑇Δ𝑉⁻¹𝐽^𝑇Δ𝑉, (9) where the higher-order term 𝐻=[𝑓^′′] is known as the Hessian matrix [24]. In (9) by neglecting 𝐻, the update conductivity vector reduces to Δ𝜎=𝑓^′𝑇𝑓^′−1𝐽^𝑇[∇𝑉]. (10)

In general, using NRIT method, the conductivity update vector expressed as in (10) can be represented for kth iteration (where 𝑘 is a positive integer) as Δ𝜎_𝑘=𝐽_𝑘^𝑇𝐽_𝑘⁻¹𝐽_𝑘^𝑇Δ𝑉_𝑘, (11) where [Δ𝑉_𝑘] and [𝐽_𝑘] are the voltage mismatch matrix and Jacobian matrix, respectively.

The [𝑓^′]^𝑇matrix in (11) is always ill conditioned [19–24], and hence small measurement errors will make the solution of (11) changes greatly. In order to make the system well posed, the regularization method [19–26] is incorporated into the reconstruction algorithm by redefining the object function [19–26] with regularization parameters as 𝑠_𝑟=12𝑉_𝑚−𝑓²+12𝜆‖𝐺𝜎‖², (12) where 𝑠_𝑟 is the constrained least-square error of the regularized reconstructions, 𝐺 is the regularization operator, and 𝜆 (the positive scalar) is called the regularization coefficient [19–26] 𝑠_𝑟=12𝑉_𝑚−𝑓^𝑇𝑉_𝑚+1−𝑓2𝜆(𝐺𝜎)^𝑇(𝐺𝜎). (13)

Differentiating the inject function in (12) with respect to the elemental conductivity: the following relations are obtained 𝑠^′_𝑟=−𝑓^′𝑇𝑉_𝑚−𝑓+𝜆(𝐺)^𝑇(),𝑠𝐺𝜎_𝑟′′=𝑓^′𝑇𝑓^′−𝑓^′′𝑇𝑉_𝑚−𝑓+𝜆𝐺^𝑇𝐺. (14)(15)

Now, using Gauss-Newton- (GN-) method-based minimization process, the conductivity update vector [Δ𝜎] is obtained as 𝑠Δ𝜎=^′𝑟𝑠_𝑟^′′=𝑓^′𝑇𝑉_𝑚−𝑓−𝜆(𝐺)^𝑇(𝐺𝜎)𝑓^′𝑇𝑓^′−𝑓^′′𝑇𝑉_𝑚−𝑓+𝜆𝐺^𝑇𝐺. (16)

Neglecting the Hessian matrix [24] in (15) 𝑠Δ𝜎=^′𝑟𝑠_𝑟^′′=𝑓^′𝑇𝑉_𝑚−𝑓−𝜆(𝐺)^𝑇(𝐺𝜎)𝑓^′𝑇𝑓^′+𝜆𝐺^𝑇𝐺. (17)

Replacing 𝑓^′ by 𝐽 and 𝐺^𝑇𝐺 by𝐼 (identity matrix) (21) reduces to 𝐽Δ𝜎=^𝑇𝑉_𝑚−𝑓−𝜆𝐼𝜎𝐽^𝑇,𝐽+𝜆𝐼 (18)

where the matrix 𝐽=𝑓^′ is the Jacobin as stated earlier.

Thus, the conductivity update vector ([Δ𝜎]) is found as Δ𝜎=𝐽^𝑇𝐽+𝜆𝐼⁻¹𝐽^𝑇𝑉_𝑚−𝑓−𝜆𝐼𝜎. (19)

Sometimes, the last term (𝜆𝐼𝜎) is neglected [22], and the conductivity update vector [Δ𝜎] is calculated as Δ𝜎=𝐽^𝑇𝐽+𝜆𝐼⁻¹𝐽^𝑇𝑉_𝑚−𝑓. (20)

In general, the EIT image reconstruction algorithm provides a solution of the conductivity distribution of the DUT for the kth iteration as 𝜎_𝑘+1=𝜎_𝑘+𝐽^𝑇𝐽+𝜆𝐼⁻¹𝐽^𝑇𝑉_𝑚−𝑓−𝜆𝐼𝜎_𝑘. (21)

The EIT algorithm starts with the solution of FP obtained from the EIT governing equation, and the [𝑉_𝑐] is calculated for a known current injection matrix [𝐶] and an initial guess (known or assumed) conductivity matrix [𝜎₀]. The voltage mismatch matrix [Δ𝑉] is estimated, and then it is used to calculate the conductivity update matrix [Δ𝜎] using GN-MMA and is added to the initial conductivity matrix ([𝜎_𝑜]) to update it to a new conductivity matrix [𝜎₁=𝜎_𝑜+Δ𝜎] using NRIT. New update matrix [𝜎₁] is used in forward solver to obtain a new calculated boundary data matrix [𝑉_𝑐1] which provides a new voltage mismatch matrix [Δ𝑉₁]. Therefore, the NRIT algorithms iteratively calculate the [Δ𝜎] using GN-MMA to find out an optimized [𝜎] matrix for which the [Δ𝑉] reaches its minimum value. Thus, the EIT reconstruction algorithm is found to work in the following sequences:(1)forward solver calculates the boundary potential matrix [𝑉_𝑐] for a known current injection matrix [𝐶] and an initial guess (known) conductivity matrix [𝜎₀],(2)measured voltage data matrix  [𝑉_𝑚] is compared with [𝑉_𝑐] to estimate the [Δ𝑉] as [Δ𝑉=𝑉_𝑐−𝑉_𝑚],(3)Jacobian (𝐽) is computed,(4)conductivity update vector [Δ𝜎] is calculated by Gauss-Newton-based minimization algorithm,(5)[𝜎_𝑜] matrix is updated to a new conductivity matrix [𝜎₁=𝜎_𝑜+Δ𝜎] by adding [Δ𝜎] to [𝜎] using Newton-Raphson iteration technique (NRIT),(6)new update matrix [𝜎₁] is used in forward solver to calculate the new voltage mismatch matrix [Δ𝑉₁],(7)check whether the [Δ𝑉₁] is minimum or not or compare the [Δ𝑉] with a specified error limit (𝜀) if provided,(8)stop the algorithm if Δ𝑉≤𝜀 condition is achieved, otherwise repeat the steps 1 to 7 until the specified stopping criteria (Δ𝑉≤𝜀) is achieved.

2.3. Boundary Data Simulator (BDS)

A two-dimensional boundary data simulator (BDS) is developed in MATLAB R2010a [27] using finite element method (FEM) [15] to generate accurate boundary data for studying the EIT reconstruction algorithms. The MATLAB-based BDS is developed as an absolute 2D data simulator for EIT image reconstruction studies, and it is used suitably to generate the errorless 2D boundary data to study and modify the inverse solver of a 2D EIT system. As BDS is developed in a computer software, it is found free from errors produced by the EIT instrumentation and phantom. BDS also allows us to generate boundary potential data for different type of phantom geometry, inhomogeneity geometry (shape, size, and position), inhomogeneity conductivity profiles, and background conductivity profiles. Moreover, it is developed as a compact, absolutely stable, and easy to use and handle for EIT studies. It is developed in such a way that it can be modified for further modifications.

BDS is developed with MATLAB-based computer program consisting of four-part imaging domain simulator (IDS), EIT model developer (EMD), current injection simulator (CIS), and boundary data calculator (BDC). Imaging domain simulator (IDS) in BDS simulates a domain with inhomogeneity with their corresponding conductivity distributions. EIT model developer (EMD) derives a mathematical model of the forward solver by applying FEM on the governing equation of the DUT in the form of a matrix equation. Current injection simulator (CIS) simulates a constant current injection through the definite points at the domain boundary with neighbouring current injection protocol [1, 2, 28–30]. The boundary data calculator (BDC) solves the governing equation by solving the forward model and calculates the potentials at all electrodes at the domain boundary.

Imaging domain simulator (IDS) first defines a DUT with a desired area (𝐴_𝐷) defined by a required diameter and defined with a particular coordinate system. Imaging domain simulator applies the FEM to discretize the domain with a 2D finite element mesh containing finite element of triangular elements (𝑡) and finite number of nodes (𝑛). In IDS, a circular domain (Ω) to be imaged is defined with a required radius (𝑅_𝑝) using the Cartesian coordinate system (Figure 1(a)), and the domain is discretized with a finite element (FE) mesh (Figure 1(b)). The mesh is symmetrically composed of the first-order triangular elements with linear shape functions [18, 31]. The FE mesh is generated with the pdetool of MATLAB R2010a in such a way that it can be refined further to increase the number of elements as per the requirement. All the coordinates and parameters assigned to the finite elements and the nodes are stored in corresponding matrices. Boundary nodes are identified, and the sixteen nodes among the boundary nodes are assigned as the electrodes called the electrode nodes. Inside the domain one (or more) smaller region (regions) is (are) defined as the inhomogeneity (inhomogeneities) positioned at a particular place. The center point (𝑃) of the inhomogeneity with the required shape and size is positioned inside the phantom domain by defining its center with a polar coordinate (𝑟,𝜃) as shown in Figure 1(a). Single or multiple inhomogeneities are defined with their desired areas (𝐴_𝐼) inside the DUT, and elements within the inhomogeneity and the background are identified. The background area is defined as the area of the domain surrounding the inhomogeneity (𝐴_𝐵=𝐴_𝐷−𝐴_𝐼), and the elements within the background area (𝐴_𝐵) are identified. The elements within the inhomogeneity are assigned with a particular conductivity called inhomogeneity conductivity, (𝜎_𝑖) while the rest of the elements are assigned with a different conductivity called background conductivity (𝜎_𝑏) as shown in Figure 1(b). The assigned conductivity values of all the elements are assumed to be featured at their corresponding centroids.

(a)

(b)

Figure 1 

(a) A circular phantom domain (Ω) with an inhomogeneity defined by polar coordinate (𝑟 and 𝜃); (b) a phantom domain (discretized by an FE mesh with 2048 elements and 1089 nodes) with a circular inhomogeneity (𝑅𝑝=75mm,𝑟𝑖=25mm,𝑟=37.5mm,𝜃=45∘,𝜎𝑖=0.005S/m, and 𝜎𝑏=0.21S/m).

EIT model developer (EMD) develops the mathematical model of the forward solver by applying FEM on the governing equation and derive the forward model of a DUT in the form of a matrix equation (3) using the elemental conductivities and nodal coordinates. The EMD establishes a relationship between the current injection matrix, [𝐶] (matrix of the applied signal), and the nodal potential matrix, [Φ] (matrix of the developed signal), through the transformation matrix [𝐾(𝜎)] which is mathematically represented by (3). The global stiffness matrix [𝐾(𝜎)] in EIT is actually an admittance matrix [23] that is formed [16] using the nodal coordinates of all the elements with their corresponding conductivities. Thus, the [𝐾(𝜎)] inforward model represents the transfer function of the EIT system obtained from the governing equation by FEM formulation [19].

The current injection simulator (CIS) is used to simulate a constant current injection through the sixteen nodes called simulated electrodes (SE) on the domain boundary with neighbouring current injection protocol. The CIS works in a “for” loop to execute all the projections [1, 28, 30, 32] of current injection process. In BDS, a constant current injection is simulated into the DUT surrounded by the sixteen simulated current electrodes (SE_𝐼) with all the possible combination of SE_I pairs, and the potential data are calculated on all the electrodes called voltage electrodes (SE_𝑉) in BDC. The current injection through a particular current electrode pair (say SE_𝐼1 and SE_𝐼2) and corresponding voltage data collection from all the possible voltage electrodes (SE_𝑉1, SE_𝑉2, SE_𝑉3, SE_𝑉4, SE_𝑉15,SE_𝑉16, SE_𝑉7, SE_𝑉8, SE_𝑉9, SE_𝑉10, SE_𝑉11, SE_𝑉12, SE_𝑉13, SE_𝑉14, SE_𝑉15 and SE_𝑉16) is known as a simulated current projection (SCP). Hence, in an N-electrode EIT system, there will be N-different current projections each of which will inject current through a particular current electrode pair and collect m voltage (differential/grounded) data where m may be either equal to N or less than N depending on the EIT data collection strategy called the current pattern [1, 28, 30, 32]. Therefore, a complete scan (containing all the current projections) conducted on the DUT yields 𝑁×𝑚 voltage data. As the BDS is studied for sixteen electrode system, the CIS runs for sixteen times and provides sixteen current projections (SCP_𝑉1, SCP_𝑉2, SCP_𝑉3, SCP_𝑉4, SCP_𝑉15, SCP_𝑉16, SCP_𝑉7, SCP_𝑉8, SCP_𝑉9, SCP_𝑉10, SCP_𝑉11, SCP_𝑉12, SCP_𝑉13, SCP_𝑉14, SCP_𝑉15, and SCP_𝑉16). Therefore, a complete data collection procedure (called a complete scan) in the BDS collects m voltage data from the voltage electrodes or voltage electrode pairs in all the sixteen current projections and computes 16×𝑚 voltage data.

Boundary data calculator (BDC) calculates the potentials (developed for a constant current injection by CIS) at all electrode points (electrode nodes) at the domain boundary in each current projection for a particular current pattern. The current injection matrix [32] is formed in CIS using the Neumann type boundary conditions, and the potential matrix is calculated from (3) using the matrix inversion technique working on L-U factorization [33] process. The BDS is developed to run in an another “for” loop for m times to calculate the m electrode potentials from voltage electrodes or voltage electrode pairs at each of the steps of the loop. This second “for” loop runs within the first “for” loop for m times and collects m voltage data for each step of first “for” loop and hence collects 16×𝑚 voltage data as first “for” loop runs for sixteen times. Moreover, as the EIT reconstruction process needs a complete scan, the BDS runs in each current projection and computes sixteen electrode potentials at each projection. The domain potential is calculated from the forward model (3), and the potential values of all the nodes are stored in a nodal potential matrix [33, 34] denoted by [𝑀_NP]. Boundary potential data are separated from [𝑀_NP] and stored in a different matrix called boundary potential matrix [𝑀_BP]. The electrode potential data are extracted from the nodal potential matrix [𝑀_NP] and are stored in a separate matrix called electrode potential matrix [𝑀_EP]. In sixteen electrode EIT system, the [𝑀_EP] is formed as a column matrix and contains the 16×𝑚 electrode potentials (differential or grounded) obtained for all the projections.

2.4. Neighbouring or Adjacent Current Injection Method

In neighbouring or adjacent current injection method, first reported by Brown and Segar [35], the current is applied through two neighbouring or adjacent electrodes, and the differential voltages is measured successively from all other adjacent electrode pairs excluding the pairs containing one or both of the current electrodes. For a sixteen electrode EIT system with domain under test surrounded by equally spaced sixteen electrodes (E₁, E₂, E₃, E₄, E₅, E₆, E₇, E₈, E₉, E₁₀, E₁₁, E₁₂, E₁₃, E₁₄, E₁₅, and E₁₆), the neighbouring method injects current through the current electrode pairs for sixteen current projections (Figure 2), and the differential voltages are measured across the voltage electrode pairs using four electrode method in each projection.

(a)

(b)

Figure 2 

Current injection and boundary data collection in neighbouring current injection method; (a) data collection method as suggested by Brown and Segar and (b) data collection strategy as suggested by Cheng et al.

As shown in Figure 2(a) in the first current projection (P₁) of adjacent method, the current is injected through electrode 1 (E₁) and electrode 2 (E₂), and the thirteen differential voltage data (𝑉₁,𝑉₂,𝑉₃,…,𝑉₁₃) are measured successively between the thirteen electrode pairs E₃-E₄, E₄-E₅,…, and E₁₅-E₁₆, respectively (Figure 2(a)). As reported by Brown and Segar, in neighbouring current injection method, the current density within the DUT is found highest between the current electrodes (E₁ and E₂ for P₁); the current density then decreases rapidly as a function of distance [35]. Similarly, in current projection 2 (P₂), the current signal is injected through electrodes 2 (E₂) and 3 (E₃), and an another set of thirteen differential voltage data (𝑉₁,𝑉₂,𝑉₃,…,𝑉₁₃) are collected between the thirteen electrode pairs E₄-E₅, E₅-E₆,…, E₁₆-E₁, and so on. Lastly, in the current projection 16 (P₁₆), the last set of thirteen differential voltage data (𝑉₁,𝑉₂,𝑉₃,…,𝑉₁₃) are collected between the thirteen-electrode pairs E₂-E₃, E₃-E₄,…, and E₁₄-E₁₅ by injecting the current through the electrodes E₁₆ and E₁. Thus, the neighbouring current injection method in a sixteen electrode EIT system data collection procedure consists of sixteen current projections (P₁,P₂,P₃,…,P₁₅, and P₁₆), and each of the current projection yields thirteen differential voltage data (𝑉₁,𝑉₂,𝑉₃,…,𝑉₁₃). Therefore, a complete data collection scan with the neighbouring current injection method in a sixteen electrode EIT system yields 16×13=208 voltage measurements.

Though in neighbouring method, EIT boundary data are not collected across the electrode pairs containing one or two current electrode for contact impedance problem [35], but sometimes it is advantageous to collect the boundary data from all the electrodes including the current electrodes to obtain the greatest sensitivity to the resistivity changes in the domain as reported by Cheng et al. [36]. In the present study, the boundary potentials are calculated at all the electrodes (Figure 2(b)) with respect to a virtual ground point selected within the DUT. Hence, in a complete data collection scan, the potentials on all the electrodes are collected in all the sixteen current projection and are stored in [𝑀_EP]. Therefore, the [𝑀_EP] is found as a column matrix containing 16×16 voltage data all collected with respect to the virtual ground point of the DUT. Hence, in the present study, with neighbouring current injection method, the [𝑀_EP] is found as a 256×1 matrix containing 256 electrode potentials. In the present study, 1 mA current injection is simulated through the electrodes of the simulated domain containing sixteen nodal electrodes using adjacent or neighboring current injection protocol (Figure 2(b)). The potentials on all the sixteen electrodes are calculated using boundary data calculator (BDC) for all the current projections, and the electrode potential matrix [𝑀_EP] is used as the calculate boundary potential matrix [𝑉_𝑐] to reconstruct the conductivity distribution of DUT.

The BDS is designed in such a way that a huge number of voltage data sets can be generated using different types of phantoms with their different design parameters. Boundary potential data [𝑉_𝑐] are generated for different type of phantom configurations, and the boundary data have been tested with electrical impedance tomography and diffuse optical tomography reconstruction software (EIDORS) [37, 38] for 2D-EIT. A large number of data sets are generated by changing the values of one or more phantom parameters like: phantom diameter (𝐷=2𝑅_𝑝), inhomogeneity radius (𝑟_𝑖), inhomogeneity geometry (shape, size, and position), inhomogeneity number (𝑁_𝑖), bathing solution conductivity (𝜎_𝑏), and inhomogeneity conductivity (𝜎_𝑖). 1 mA current injection is simulated to the domain boundary, and corresponding boundary data sets are used for image reconstruction in EIDORS. Data generation in BDS and image reconstruction in EIDORS are studied for different inhomogeneity geometries in DUT. Reconstruction is also studied for different iterations and for multiple inhomogeneity reconstruction to evaluate the BDS.

3. Results and Discussion

Image reconstruction quality in EIT depends on the boundary data accuracy which is dependent on the geometric accuracy of the inhomogeneity developed in BDS. Dimensional accuracy of the inhomogeneity depends on the number of finite elements in the FE mesh or mesh refinement number (𝑁_mr) as shown in Figure 3. As the 𝑁_mr increases, the number of elements in the FE mesh is increased, and hence the geometric accuracy of the inhomogeneity increases which gives more accurate boundary data and better image reconstruction (Figure 3). But the BDS with a highly refined mesh needs a high PC memory and large computation time. In this paper, the mesh refinement is found suitable as 𝑁_mr=4 as per the configuration of the PC (2.4 GHz/1.5 GBRAM/ P-IV) used. It is observed that the FE mesh with 𝑁_mr=4 (containing 2048 elements and 1089 nodes) gives almost an accurate geometry (Figure 3) to the desired inhomogeneity and generates a reconstructible data set in less than 10 seconds. EIDORS reconstructs the resistivity images from the BDS data sets using regularized image reconstruction technique.

(a)

(b)

(c)

(d)

Figure 3 

Circular inhomogeneity (𝑅𝑝=75mm,𝑟=0,𝑟𝑖=37.5mm,𝜎𝑖=0.005S/m,and𝜎𝑏=0.21S/m) with FEM mesh with different number of finite elements: (a) 512 elements and 289 nodes, (b) 2048 elements and 1089 nodes, (c) 8192 elements and 4225 nodes, and (d) 32768 elements and 16641 nodes.

Results show that the resistivity or conductivity can be successfully reconstructed from the boundary data generated by our BDS using a circular domain (𝑅_𝑝=75mm) with a circular inhomogeneity (𝑟=37.5 mm, 𝑟_𝑖=25mm, 𝜃=45^∘, 𝜎_𝑖=0.005S/m, and 𝜎_𝑏=0.21S/m) in the 9th iteration (Figure 4). It is also observed that the reconstructed shape of the inhomogeneity is similar to that of the original one (Figure 4(a)), and the reconstructed conductivity profile in Figure 4(b) is almost similar to that of the original object in Figure 4(a).

(a)

(b)

Figure 4 

(a) Simulated domain with a circular object (𝑅𝑝=75mm, 𝑟=37.5, 𝜃=45∘, 𝑟𝑖=25mm, 𝜎𝑖=0.005S/m, and 𝜎𝑏=0.21S/m); (b) reconstructed image of (a).

Iteration studies shows that in different reconstruction steps called iterations (Figure 5), the reconstructed images become more localized from iteration to iteration and the reconstruction errors (appeared by the red color at phantom periphery) are gradually reduced (Figure 5).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

Figure 5 

Reconstructed images of a simulated domain with a circular inhomogeneity (𝑅𝑝=75mm, 𝑟=37.5, 𝜃=45∘, 𝑟𝑖=25mm, 𝜎𝑖=0.005S/m, and 𝜎𝑏=0.21S/m) for different number of iterations in inverse solver in EIDORS-2D: (a) 1st iteration, (b) 2nd iteration, (c) 3rd iteration, (d) 4th iteration, (e) 5th iteration, (f) 6th iteration, (g) 7th iteration, (h) 8th iteration, (i) 9th iteration, (j) 10th iteration, (k) 11th iteration, (l) 12th iteration, (m) 13th iteration, and (n) 14th iteration.

It is observed that the resistivity is successfully reconstructed from the boundary data in the 9th iteration (Figures 5(i) and 5(j)), though the shape of all the reconstructed images in 9th–12th iterations is almost similar to that of the original one (shown by dotted circles in Figure 5). As the reconstructed resistivity profile similar to that of the original is obtained only in the 9th iteration, the 9th iteration is taken as the optimum reconstruction. In 13th and 14th iterations, the resistivity is overestimated, and the images are lost. The optimum iteration number depends on the data accuracy and reconstruction algorithm, and hence the BDS can be used to generate the boundary data sets required for assessing the inverse solver in EIT.

Voltage data are also generated for a domain (𝑅_𝑝=75mm) with the circular inhomogeneities (𝑟_𝑖=25mm, 𝜎_𝑖=0.005,S/m, and 𝜎_𝑏=0.21S/m) positioned at different places using the BDS (Figure 6). It is observed that the reconstructed image is more circular for an inhomogeneity positioned at the phantom centre where 𝑟=0 and 𝜃=0^∘ (Figure 6(a)). On the other hand, for 𝑟≠0, that is, for the inhomogeneities near domain boundary (Figure 6(b)), reconstructed images are not perfectly circular because of the comparatively less accurate shape of the original object obtained for 𝑟≠0. For a less number of mesh refinements, the geometry of the original side objects is not exactly circular itself (Figure 4), and hence the corresponding boundary data have lower accuracy. An FE mesh with large 𝑁_mr can easily produce an accurate geometry for the boundary objects (objects near domain boundary) with proper shape, which gives a boundary data without geometric error and automatically improves the image shape.

(a)

(b)

(c)

(d)

Figure 6 

Image reconstruction of a simulated domain with a circular inhomogeneity (𝑅𝑝=75mm, 𝑟𝑖=30mm, 𝜎𝑖=0.005S/m, and 𝜎𝑏=0.21S/m) at different positions (𝑟,𝜃): (a) simulated domain with inhomogeneity at 𝑟=0mm and 𝜃=0∘, (b) reconstructed image of the domain shown in Figures (a) and (c) simulated domain with inhomogeneity at 𝑟=37.5mm and 𝜃=0∘, and (d) reconstructed image of the domain shown in Figure (c).

Boundary data sets are also generated with a circular domain (𝑅_𝑝=75mm and 𝜎_𝑏=0.21S/m) with a circular inhomogeneity (𝜎_𝑖=0.005S/m) with different diameters (2𝑟_𝑖) and all positioned at the phantom center (𝑟=0). The boundary data are calculated and used for reconstructing the resistivity images. Results show that for the domain discretized with 𝑁_mr=4, the data sets, generated with a diameter larger than 13.3% of the phantom diameter, are accurate enough (Figures 7(a)–7(f)) to reconstruct the resistivity images in EIDORS-2D. It is clearly observed that for 𝑁_mr=4, the triangular elements within the inhomogeneity with smaller 𝑟_𝑖 are unable to shape themselves into a proper circle (Figure 7(g)). Hence, the data obtained for the inhomogeneity with a diameter of 20 mm has low accuracy (Figure 7(g)), and hence the resistivity image (Figure 7(h)) is found with low resolution showed and some reconstruction error (appeared in the red color at phantom periphery). Increasing the FE elements in BDS, the boundary data error can be minimized, and the improved resistivity image can be achieved even for smaller inhomogeneities with a diameter less than 13.3% of 𝑅_𝑝.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 7 

Image reconstruction of circular inhomogeneities (𝑅𝑝=75mm, 𝑟=0, 𝜎𝑖=0.005S/m, and 𝜎𝑏=0.21S/m) with different diameters: (a) original object with 𝑟𝑖=40mm, (b) reconstructed image of the object shown in (a), (c) original object with 𝑟𝑖=30mm, (d) reconstructed image of the object shown in (c), (e) original object with 𝑟𝑖=20mm, (f) reconstructed image of the object shown in (e), (g) original object with 𝑟𝑖=10mm, and (h) reconstructed image of the object shown in (g).

Boundary potential data are also generated for domains (𝑅_𝑝=75mm) containing multiple circular inhomogeneities (𝑟_𝑖=25mm, 𝑟=37.5mm, 𝜎_𝑖=0.005S/m, and 𝜎_𝑏=0.21S/m) placed at different positions inside the domain (Figure 8). Figure 8(a) shows a domain with two circular inhomogeneities (180^∘ apart from each other) which are placed at a central distance (𝑟) of 37.5 mm. Similarly, another domain with three circular inhomogeneities (120° apart from each other) placed inside the phantom domain is shown in Figure 8(c). All the inhomogeneities in both the domains are positioned at a central distance (𝑟) of 37.5 mm. 1 mA current is simulated with the neighbouring current pattern, and the boundary data are collected for resistivity reconstruction. It is noticed that the resistivity images (Figures 8(b) and 8(d)) of inhomogeneities in both the domains are reconstructed successfully.

(a)

(b)

(c)

(d)

Figure 8 

Image reconstruction of multiple inhomogeneities (𝜎𝑏=0.21S/m and 𝜎𝑖=0.005S/m): (a) simulated domain with two circular objects (𝑟𝑖=25mm, 180° apart), (b) reconstructed image of the domain shown in (a), (c) simulated domain with three circular objects (𝑟𝑖=25mm, 120° apart), and (d) reconstructed image of the domain shown in (c).

Results show that the boundary data simulator can be efficiently used to generate boundary potential data for a huge number of phantom configurations in less than 10 seconds. BDS is software-based virtual EIT phantom, and hence it has a number of advantages over the practical and mesh phantoms. The literatures [39–41] presenting the phantom simulations are limited, and they only discuss the software phantoms developed for their own systems. BDS is a software-based versatile boundary data simulator which generates boundary data suitable for studying the reconstruction algorithm required for several EIT systems, and hence it is better suited for assessing the performance of the inverse solver of 2D electrical impedance tomography.

4. Conclusions

A MATLAB boundary data simulator (BDS) is developed for studying the resistivity reconstruction in inverse solvers of 2D-EIT. BDS is developed with four parts: imaging domain simulator (IDS), EIT model developer (EMD), current injection simulator (CIS), and boundary data calculator (BDC). Imaging domain simulator (IDS) simulates a domain with single or multiple inhomogeneities of different geometries defined with their corresponding conductivity distributions, whereas the EIT model developer (EMD) derives a forward model using FEM to solve the governing equation of the DUT. Current injection simulator (CIS) simulates a constant current injection through the simulated electrodes positioned at the domain boundary with the neighbouring current injection protocol. The boundary data calculator (BDC) solves the forward model to solve the governing equation and calculates the potentials at all the simulated electrodes. Boundary data are generated with different type of domains simulated in BDS by changing its input parameters. Resistivity images are reconstructed from the boundary data using standard EIT reconstruction software called EIDORS, and the BDS is evaluated. It is observed that the BDS with FE mesh with 2048 elements can simulate an inhomogeneity of desired geometry with suitable accuracy. The BDS with 2048 elements suitably generates the boundary data for simulated domains containing the objects with different geometries which are found efficient for image reconstruction in EIDORS. Results also show that the conductivity or resistivity profiles of the domains simulated in BDS are successfully reconstructed from their corresponding boundary data generated for different type of single and multiple inhomogeneities. By changing the inhomogeneity position, diameter, and number in BDS, boundary data are successfully generated as well as the resistivity images are reconstructed successfully. Multiple inhomogeneity imaging shows that the BDS suitably generates boundary data with the desired accuracy, and the boundary data are found efficient for resistivity reconstruction in EIDORS. Results also show that for the simulated domains discretized with 𝑁_mr=4, the boundary data sets generated for circular inhomogeneity with a diameter larger than 13.3% of the phantom diameter are accurate enough to reconstruct the resistivity images in EIDORS. Increasing the FE elements in BDS, the boundary data error can further be minimized, and the improved resistivity image reconstruction can be obtained even for smaller inhomogeneities. Hence, it is concluded that the BDS generated a number of boundary data sets which can suitably be used for inverse solver assessment in EIT.