Dbar Algorithm for EIT
Project description
pyDbar for Electrical Impedance Tomography
Thank you for the interest in pyDbar! This package implements the Dbar method for Electrical impedance tomography and follows up on the pyEIT presented in https://github.com/liubenyuan/pyEIT.
1. Overview
The Dbar method can be used for Electrical Impedance Tomography in a 2D framework. This method is based on the theoretical reconstruction procedure of Nachman and a regularizing strategy that was established by Siltanen and provides a natural framework for numerical implementation There is an implementation in Matlab available at https://wiki.helsinki.fi/display/mathstatHenkilokunta/EIT+with+the+Dbar+method%3A+discontinuous+heartandlungs+phantom.
For now we will keep things simple and a proper documentation will follow.
2. Instalation
To install this package it is enough to use pip:
pip install pydbar.
Moreover, dependencies are handled by pip.
3. To start off:
Test cases are present in the tests folder. Soon, there will be a presentation here of the output.
4. Algorithm Review
Here we explain the main steps of the Dbar algorithm with either direct numerical implementation or with theoretical motivation.
DirichlettoNeumann map
The first step in EIT is to apply currents on a finite set of electrodes around the boundary and measure the corresponding voltages obtained thereafter. Let L be the number of electrodes, then we can apply at most L1 linearly independent current patterns (by current pattern we designate a vector of lenght L that describes the values of current applied in the electrodes). We designate by C^{n} the nth current pattern. There are a few choices that can be considered for current patterns, but for now we fix our focus in two. For n=1,..., L1 we have:
 Adjacent current patterns:
 Trignometric current patterns:
The corresponding voltages are denoted by V^{n} and each of them is a vector of dimension L (corresponding to the set of electrodes). One of the required assumptions is that for each current pattern the sum of the measured voltages equals zero. We assume that an EIT device provides the data in this manner.
For our algorithm we create a mapping from this data: DirichlettoNeumman map, which establishes a relation between the Dirichlet boundary conditions (voltages) and the Neumann boundary conditions (currents) in conductivity model.
How to obtain the DirichlettoNeumann map (DtoN map)?

Normalize the measurements:
c^{n}= C^{n}/C^{n}_{2} and v^{n}= V^{n}/C^{n}_{2}

Define the NeumanntoDirichlet map (inverse of DtoN map), which is a L1 x L1 matrix:
 Compute the matrix approximation of DirichlettoNeumann map for electrodes of Area A and a body of radius r:
This mapping is the matrix approximation of the continuum operator for a body with conductivity γ.
Furthermore, the Dbar algorthm requires either the DtoN map of a body with the same geometry but with conductivity 1 or a reference DtoN map that corresponds to the same body with a different physiology (e.g. expiration and inspiration). This last case is designated by tdEIT and is very useful when monitoring health conditions that change with time like air ventilation in the lungs. With this in mind, the read_data class was constructed to compute one DtoN map from one set of current patterns and corresponding measured voltages.
Scattering Transform
The second step of the algorithm is to compute the scattering transform. This transform is of nonphysical nature, i.e. has no physical explanation, but is tightly connect with the the conductivity and can be determined by the DirichlettoNeumann map. It is given in the new spectral parameter k in the following manner:
Now due to the illposedness of the inverse problem, the noise in the measurements highly affects the values of t for large k. The regularization strategy proven by Siltanen shows that we only need to compute the following approximation for the scattering transform:
The beauty of this regularization strategy lies in a finite discretization of k inside a disk of radius R being required. The choice of R dependends in an estimate of the noise present in the measurements.
To completely understand the computation of the scattering transform we need to obtain ψ at the boundary. Nachman's was able to prove that the values of ψ at the boundary can be determined by the solution of a boundary integral equation for k ∈ ℂ \ {0}:
However, this equation takes a lot of time to solve, thus for speedup purposes we choose an approximate solution. We can have two approaches:
 Solve the boundary integral equation with S_{0} substituting S_{k} and with the following definition:
This approach is present in [2].
 Take ψ≃ e^{ ikz}, since this is the leading behavior and do not solve any integral equation. This is the fastest way, but its triviality leads to more errors.
Now the last step is to solve the Dbar equation, which gives the name to the algorithm.
Dbar equation
Now, having computed the scattering transform or an approximation of it, we are able to solve the following PDE in the k variable:
Afterwards, we can determine the conductivity by:
.
By substituting the scattering transform by its approximations and solving the equation with respect to it we obtain respective approximations to the conductivity.
5. Numerical Implementation
Above, we present an immediate numerical implementation of the DirichlettoNeumann map that trivially translates to code and a theoretical implementation of the regularization strategy for the Dbar and corresponding determination of the conductivity. Accordingly, we present the numerical implementation for both versions of the scattering transform and the solution to the Dbar equation.
texp Scattering transform:
The texp scattering transform corresponds to take the exponential approximation of ψ, which has the following form:
Let us consider that we have equally distributed electrodes around the boundary. By our assumption of Ω being a disc we define, without loss of generality, the positions of the electrodes to be where .
From this, we can take an approximation of the exponentials in these set points as an expansion in the orthonormal basis of current pattern (L1 degrees of freedom):
.
By immediate substitution and approximation of the integral at the electrode points this leads too:
.
tB Scattering Transform:
Here the first step is to solve the boundary integral equation with S_{0}. For such we consider an initial vector approximation ψ^{B} of the solution of the (continuous) boundary integral equation and take the expansion in terms of the orthonormal current patterns:
.
Further, we also take a LxL matrix approximation of G_{0} that defines the S_{0} operator as:
.
Through this matrix approximation, we can translate the boundary integral equation into a linear system of equations. Defining Φ to be the current pattern matrix (L x L1) leads us too:
,
To compute the scattering transform we need to solve the last system of equations for each k in a disk of radius R. Afterwards, we compute it by:
,
where δL is the difference between L_{γ} and L_{1}.A more in depth look at the details can be found in [1] with respect to the solution of the boundary integral equation with S_{k}.
Dbar Equation:
Now we have to solve the Dbar equation in the kvariable. This is the essence of this new spectral parameter, it allows for the reconstruction of the conductivity without iterative solutions of an FEM model of conductivity. Although, we do not solve the PDE directly... we pass it to an integral equation and use Fourier Transform after discretization to solve the linear system of equations.
Integral Equation:
Solve for each z∈Ω:
Recall that, the regularization procedure establishes that for k>R the scattering transform is 0, in both cases. Hence, the integral can be taken in the disk of radius R. To practically solve this equation we need to deal with the singularity properly. A fast approach to solve the equation is to transform it into a periodic setting (initial idea was by Vainikko, [3], and adapted to the Dbar equation in [4]). Let be the fundamental solution of the Dbar operator. First, we write the integral equation in convolution terms:
The periodization is done by tiling the kplane with a square Q=[2R3ε, 2R + 3ε]^{2} for some ε>0, in this manner, Q containes the disk of radius 2R. Now, we take a smooth cutoff function η which fulfils η(k)=1 for k<2R+ε and η(k)=0 for k>2R+2ε. With this we define the periodic approximate green function .
In this framework we have the periodic version of the integral Dbar equation given by:
where ⋆ denotes periodic convolution and T_{R} is the defined as T for k<R and 0 otherwise. The interesting part is that the solution of this equation and of the full equation coincide in k<R, which is what we desire.
By discretizing the kplane into an equally spaced grid with step h, the periodic convolution can be computed efficiently by Fast Fourier transforms and we obtain the following system of equations:
where the boldface notation represents the matrices formed by evaluation of the functions in the kgrid and ⋅ represents matrix multiplication. This equation can be written in a simpler form, which will be helpful to describe our code below:
.
Since μ_{R} is inside the Fourier transform, we can not describe this system through the simple form Ax=b. Hence, the system of equations lends itself well to a solution by a matrixfree iterative solver like the GMRES.
Due to conjugation on μ_{R} the equation is ℝlinear and separation of the real and imaginary parts, which is done by justposition of the real and imaginary parts into a vector. Further details, can be seen in [1].
This is what we need for implementation and it provides a framework for the code we present in this package. A description of each class and function and its role is presented next.
6. Documentation:
Here, we provide an overview of each class and examples of usage. Our package possesses three classes:
sim, read_data, k_grid, scattering, dBar.
The first three classes are independent and concern the input information and framework. The fourth connects the gather data into a single operator, which is transmitted to the dBar class. This class is dependent of both the scattering and k_grid and performs the solver.
sim class
class sim(anomaly, L)
This class is a simulator of voltages obtained through the FEMsolution of the Direct problem. We can create the conductivity of a body in a disk and then solve the direct problem for this conductivity by applying pairwise adjacent current patterns (as described above). Thereafter, we can determine the voltage at the electrode positions and hence create the necessary data to apply an inversion procedure.
This class is based solely on pyEIT.
Parameters:

anomaly: list
By assumption we have that the background conductivity is equal to 1. Then we can create disk anomalies by selecting the center, radius and the conductivity of the region.

L: int
Number of electrodes.
Attributes:

Current: (L1, L)
Matrix with the currents patterns applied to the body.

Voltage: (L1, L)
Matrix with the "measured" voltages obtained to sucessive FEMsolutions of the direct problem.
Methods:

simulate():
This function uses each of the adjacent pair current patterns, solves the direct problem by FEM, and determines the voltages at the electrodes.
read_data class
class read_data(Object, r, AE, L)
Defines the DirichlettoNeumann matrix approximation from the experimental data: electrical measurements, radius of the body, area of electrode and number of electrodes. The current version only holds for circular domains and for objects which have conductivity equal to 1 near the boundary.
Parameters:

Current: (L1, L) 2darray
Matrix of current patterns. Each line is a current pattern.

Voltage: (L1, L) 2darray
Matrix of voltage measurements. Each line corresponds to the current pattern on the same line of the current matrix.

r: float
Radius of the body;

AE: float
Area of one electrode;

L: int
Number of electrodes;
Attributes:

Current: (L, L1) 2darray
Matrix of the normalized current patterns, each column represents one normalized current pattern;

Voltage: (L, L1) 2darray
Matrix of the normalized voltages, where the ith column represents the voltages measured for the ith current pattern;

DNmap: (L1, L1) 2darray
Matrix approximation of DirichlettoNeumann map for the measured data.
Methods:

load():
This function starts by performing the normalization of the current and voltage measurements and thereafter computes the DirichlettoNeumann as described above. It is called directly from the constructor and the user does need to use it directly.
Examples:
k_grid class
class k_grid(R, m)
Definition of the grid discretization of the kplane with respect to the parameters R and m. As required by the periodic discrete convolution, the grid is defined in [2.3R, 2.3R]^2 with step size h=2(2.3R)/2^{m}.
Parameters:

R: float
Defines the size of the k grid

m: int
Defines the step size.
Attributes:

R: float
As the parameter.

m: int
As the parameter.

s: float
It is defined as 2.3R, to simplify definitions further ahead.

h: float
Step size of the k grid.

N: int
It is defined as 2^{m. Therefore, it is the number of elements on each line of the grid.}

k: (N, N) complex 2darray
Complex values of the k_grid. Essential for complexvalued computations. This grid contains the value 0.

pos_x, pos_y: array_like
Arrays that store the indexes of the elements of k which are inside the disk of radius R.

index: int
Represents the position of zero.

FG: (N, N) complex 1darray
The fundamental solution of the Dbar operator, with the same structure of k.
Methods:

generate():
Defines the complex kgrid k and determines the positions inside the disk of radius R.

Green_FS():
Defines the FFT of fundamental solution FG. Recall that we perform a periodization of the convolution equation and therefore we need to establish a smooth decay to 0 close to the square limits of the grid.
scattering class
class scattering(scat_type, k_grid, Now, Ref):
Parameters:

scat_type: string
Allows to choose between the approximations of the scattering transform. Choices: "partial", "exp"!.

k_grid: Object
An object of k_grid type.

Now, Ref: Objects
Objects of read_data type (we require the DN map for our body of choice and an homogeneous/reference body).
Attributes:

tK: (k_grid.N, k_grid.N) complex 2darray
Matrix containing the values of an intermediare operator containing information about the DN map of our body of choice. Essential for the Dbar method.
Methods:

load_scattering(scat_type, k_grid, Now, Ref):
Given the choice of approximation for the scattering transform this function selects between two functions: one compute the exponential aproximation, another computes the partial approximation;

exp_scattering(Now, Ref, k_grid):
Computes the exponential approximation of the scattering transform.

partial_scattering(Now, Ref, k_grid):
Computes the partial approximation of the scattering transform.
dBar class
class dBar(R_z, m_z)
Parameters:

R_z: float
Defines the size of the zgrid.

m_z: int
Defines the step size of the zgrid.
Attributes:

Z: (2^{m_z}, 2^{m_z}) complex 2darray
Planar region of the body.

sigma: (2^{m_z}, 2^{m_z}) complex 2darray
Conductivity values at the points of the Zplane grid.
Methods:

load_mesh(R, m):
Defines the complex zgrid Z.

dBar(mu, k_grid, tK, zz):
Defines the operator [I + P] for a zz element of the zgrid, which is the lefthand side of the Dbar system. It takes as input the vector μ which is a 1darray that contains the real part concatenated with the imaginary part. We assume μ is 0 outside the disk of radius R, for speedup purposes, however this is not a constraint since the scattering transform would cutoff this terms when defining the operator.

solve(k_grid, tK):
Solve for each z on the zgrid the Dbar system with GMRES. We use the dBar function to define the [I+P] operator for each z and use the initial approximation of mu being close to 1 to start off. After solving the system, we store the value of the conductivity through γ(z)=(μ(z,0))^2.

plot():
Plot the obtained conductivity on the respective Z grid.
6. Bibliography:
[1] Mueller, J. L., & Siltanen, S. (2020). The dbar method for electrical impedance tomographyâ€”demystified. Inverse problems, 36(9), 093001.
[2] El Arwadi, T., & Sayah, T. (2021). A new regularization of the Dbar method with complex conductivity. Complex Variables and Elliptic Equations, 66(5), 826842.
[3] Vainikko, G. (2000). Fast solvers of the LippmannSchwinger equation. In Direct and inverse problems of mathematical physics (pp. 423440). Springer, Boston, MA.
[4] Knudsen, K., Mueller, J., & Siltanen, S. (2004). Numerical solution method for the dbarequation in the plane. Journal of Computational Physics, 198(2), 500517.
Project details
Release history Release notifications  RSS feed
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.