Regularization of EIT reconstruction based on multi-scales wavelet transforms

1 Introduction

Electrical Impedance Tomography (EIT) is a radiation-free imaging method [1]. A commonly used imaging strategy is called time-difference reconstruction. It attempts to reveal the conductivity distribution changes inside the human body between two time points through the electrical data obtained via electrodes attached on skin. In time-difference lung EIT typically 16 electrodes are placed equidistantly on a horizontal chest plane. For each EIT frame, currents are successively injected into the human body through adjacent electrodes. A conductivity image is then reconstructed by collecting the voltage measurements recorded from the remaining electrodes.

We denote the conductivity changes of the domain between two time steps by Δs and the measured voltage changes on the electrodes by a vector ΔV. Employing the finite element model (FEM) with M elements, the conductivity change Δs is represented by a M × 1 vector. For time-difference imaging, the reconstruction problem can be linearized [2]. Approximately, there exists the following relation:

where J^E denotes the Jacobian matrix of elements calculated at the constant conductivity 1:

This Jacobian matrix of elements could be calculated by EIDORS toolbox [3]. The Jacobian matrix J^E is ill-conditioned because the degree of freedom of this inverse problem is too large. Solving the conductivity changes from the above relation leads to an ill-posed inverse problem. Especially, the solution Δs from equation (eq 1) is unstable. To circumvent this difficulty, additional prior information is introduced to restrict the flexibility of the solution [2]. In this study, we employ a prior illustrating the sparsity hypothesis. Explicitly, we solve Δs from the following regularized optimization problem instead of (eq 1):

where α is a regularization parameter that controls the trade-off between the regularization term and the fidelity term. According to sparse regularization theory, the regularization term ∥R⋅Δ⁢s∥1 promotes a solution Δ⁢s^ with R⋅Δ⁢s^ being sparse.

Triangular finite element meshes are commonly used in EIT reconstructions. Triangular meshes can simulate irregular domains as well as the electrical properties around electrodes without approximation. Wavelet transforms have been widely used in medical imaging such as Computed Tomography or Magnetic Resonance Imaging and shown a great success in these imaging techniques. However, canonical wavelet transforms are not appropriate to apply on general finite element meshes. Instead, in this article we used spectral graph wavelet transforms and view the triangular meshes as graphs.

2 Method

2.2 Spectral graph wavelets

In the proposed method, the finite element mesh was viewed as an undirected connected graph G = {V, E, A} consisting of a set of indexed vertices V with |V| = N vertices, a set of edges E, and the adjacency matrix A. With this graph representation, the spectral graph wavelets for finite element mesh can be constructed. We refer to the original article for the detailed theoretical background as well as the construction steps of the spectral graph wavelet [5]. Indeed, the wavelet transforms are generated from two functions [5], namely, the scaling function φ₀ and the mother wavelet transform φ. For any scale t > 0, a wavelet transform φ_t can be derived from the mother wavelet transform. Intuitively, in the frequency domain, the scaling function φ₀ acts as a low-pass filter, while the wavelet transform φ_t performs like a band-pass filter with the frequency band determined by the scale t. Each of these transforms can be represented as a N × N matrix. The i-th column of a wavelet transform φ_t is the wavelet with scale t centered at the i-th node of the graph. Some selected graph wavelets locate on the same center are plotted in Figure 1. In this paper, we employ the wavelet transform φ₀ and the wavelet transforms φ_t with scales t = [0.5, 0.25, 0.125].

Figure 1

Some selected graph wavelets. These plots represent the wavelets located at the 500-th node with different scales.

2.3 Multi-scales wavelet transform

Multi-scales wavelet transforms are commonly used in image processing. They provide a multi-resolution view of the underlying problem. In canonical 2D image processing, a multi-scales wavelet transform is constructed through down-sampling the original domain and the wavelet transforms scale-by-scale. We adapt this strategy to the graph wavelet transforms developed in Section 2.2. Indeed, such multi-scales wavelet transform can be considered as a filter bank. Given a set of scales t ∈ [2⁻¹, ⋯ 2^−l], this filter bank consists of wavelets selected from the wavelet transforms φ_t and φ₀. We select all those N wavelets with the scale 2^−l and select around 2^−j ⋅ N wavelets with scale 2^{−l + j}. Here, N is the number of nodes in mesh. This process is realized by downsampling the mesh nodes and forming a chain of node sets Vl⊃Vl−1⁢⋯⊃V0. In this article, we fix l = 3. Regarding the mesh in Figure 1, the downsampled node sets are plotted in Figure 2. For j = 0, 1, 2, the wavelets of scale 2^{−l + j} with centers in V_j are selected. For j = 3, the wavelets of φ₀ with centers in V₀ are selected. Indeed, since the wavelets are identified as columns of matrices φ_t, we need only select the corresponding columns in wavelet transform φ_t for each scale. A multi-scales wavelet transform W will be obtained by setting its column as the selected wavelets scale-by-scale.

Figure 2

The down-sampled node sets from the mesh used in Figure 1. With the decreasing of wavelet scales t, the denser node sets are used to select the wavelets from the corresponding φ_t.

2.4 Regularization with multi-scales wavelet transform

Combining the above tools, we introduce a sparse solver using the multi-scales wavelet transform as regularization matrix. Given the multi-scales wavelet transform W constructed in Section 2.3, the regularized EIT inverse problem (eq 3) has been rewritten as:

(4)Δ⁢s^=argminΔ⁢s⁢12⁢∥Δ⁢V−JN⋅Δ⁢s∥22+α⁢∥WT⋅Δ⁢s∥1

Here we abused Δs to denote the conductivity changes on nodes. This optimization problem can be solved by Split Bregman iterative solver [6].

We summarize the proposed sparse solver with wavelet-based regularization (SSWR) by flowcharts in Figure 3.

Figure 3

Workflow of the proposed framework.

3 Simulation results

The performance of this wavelet based regularization framework has been evaluated through simulation. The simulations were executed on a circular finite element domain with radius 1 (see. Figure 4A). On the boundary, 16 electrodes were attached equidistantly. The background conductivity is fixed to s₀ = 1.0 S ⋅ m⁻¹. Alternative currents of 10 mA were injected successively and adjacently through the electrodes, the voltages between adjacent electrodes were recorded and denoted by V_h. Two contrasts with conductivities 1.6 S ⋅ m⁻¹ (blue contrast) and 0.6 S ⋅ m⁻¹ (red contrast) are embedded into the background (see Figure 4B). After embedding, another voltage measurement (denoted by V_ih) was simulated. The voltage difference was calculated by ΔV = V_ih − V_h . In addition, 10% white noise was added to ΔV.

Figure 4

Simulation results. (A) shows the background of the simulated phantom. (B) shows the ground truth. Two contrasts are embedded into the background. (C) and (D) are the reconstructed images from the proposed SSWR solver and the one step GN solver with Tikhonov regularization respectively. The reconstructions are performed on a coarser mesh (C and D).

The reconstructions were applied on the coarser mesh shown in Figure 4C and D to avoid the so-called “inverse crime” problem. This coarser mesh is the same as the meshes used for Figures 1 and 2. Indeed, the multi-scale wavelet transform matrix was obtained as described in Section 2.3 using the down-sampled node sets depicted in Figure 2. In Split Bregman iterative solver, the number of iterations was fixed to be 5. The regularization parameter is set to be α = 10⁻⁵. In simulation, computation time for SSWR solver takes around 0.4 s to accomplish the reconstruction. Please insert your manuscript here.

The proposed SSWR solver has been compared with the standard one step Gauss-Newton method (GN) [3]. GN was combined with standard Tikhonov regularizer which employs the identity matrix as regularization matrix. The regularization parameter used for Tikhonov regularization is set to be 5 ⋅ 10⁻³. The reconstructed images are plotted in Figure 4D. All the parameters are heuristically chosen for different solvers to get the best performance.

4 Discussion

The reconstructed images (see Figure 4C and D) demonstrated that the proposed wavelet-based regularization generates reconstruction with fewer artefacts. In particular, the image reconstructed by SSWR is more smooth and homogenous in contrast value. Moreover, the shape of the contrasts is better detected by the SSWR solver.

In this study, all parameters, such as the regularization parameters and the number of iteration for SSWR solver, are determined heuristically. A strategy to find the optimal parameter setting is preferred. Another shortcoming is the computational time. The proposed SSWR solver is an iterative solver. It is very efficient, but not enough for real-time imaging.

5 Conclusion

In this preliminary study, a multi-scales wavelet transform is constructed in the context of irregular finite element mesh. This is realized by applying the spectral graph wavelet transforms and down-sampling the nodes. Sparse regularization using this multi-scales wavelet transform was applied in simulation. Simulation results indicated that the proposed wavelet-based sparse solver obtains better images. In particular, the image reconstructed by the wavelet-based sparse solver avoids small turbulences and the shape of contrasts can be better recovered.