Anisotropic TV Regularization in Electrical Impedance Tomography: An Experimental Study

Yanying Wang

School of Mathematics and Statistics, Shandong Normal University, Jinan, China.

DOI: 10.4236/eng.2022.143013                

Abstract

Total variation (TV) regularization method is a typical method to preserve the discontinuities structure in EIT. Isotropic TV and anisotropic TV are two well-known variants of TV. The main differences between them are that the latter tends to distort the reconstructed internal inhomogeneities along the coordinate axis. In this article, we adopt the alternating direction method of multipliers (ADMM) to overcome the non-differentiability of the anisotropic TV and verify the characteristics of anisotropic TV regularization by the tank experiments.

Keywords

Electrical Impedance Tomography, Anisotropic Total Variation, Ill-Posedness

1. Introduction

Electrical Impedance Tomography (EIT) is a method of reconstructing interior conductivity distribution of the imaging target. In EIT, a series of low frequency current is injected through the electrodes attached around the boundary and we measure the induced voltages. Compared with the widely used imaging methods, such as computed tomography (CT), magnetic resonance imaging (MRI) and ultrasound imaging, EIT has a promising clinical application due to its advantages of non-invasive, no radicalization, high time resolution, etc. EIT also has very important applications in other fields, for example, to locate high flow drilling [1], to apply in Electrical Cell-substrate Impedance Spectroscopy (ECIS) [2].

However, EIT is a typical ill-posed problem. To be precise, the voltage response due to the conductivity changes decreases rapidly with the distance between the sources of conductivity changes increases [3]. Moreover, EIT measurements often suffer from noises and artifacts especially in clinical environments. Regularization is a widely used technique to deal with the ill-posedness. The main idea of regularization is to approximate the ill-posed problem by a well-posed one.

Depending on the prior information of the target conductivity distribution, the researchers proposed Tikhonov regularization [4], sparsity based regularization [5], total variation (TV) based regularization [6] and so on. Among all these regularizers, except TV based method, the remaining will blur the edges of the internal structure. TV regularizer is widely used in EIT, because it has the ability of preserving the discontinuities structure of the imaging target.

There exists difficulty in using the TV regularization due to its non-differentiability structure [6]. Many methods have been proposed to overcome this difficulty, such as Newton’s method [7] and the Primal Dual-Interior Point Methods (PDIPM) [8]. Both methods, however, are either unstable or too time consuming [9]. Split Bregman method is also used to deal with the non-differentiability, however, the ability of preserving the edges is decreased [10].

There are two common variants of TV, namely the isotropic TV and the anisotropic TV [11]. However, the anisotropic TV may distort the internal inhomogeneities along coordinate axis. Two possible reasons may cause this distortion, one is the characteristics of the model itself, and another is the algorithm to solve the inverse problem. Gonzalez [11] overcame the non-differentiability by introducing an auxiliary smoothing parameter $\beta >0$ and using the Guass-Newton method to solve the TV regularization problem, and verified the distortion characteristics by numerical simulation and tank experiments. However, in EIT, there is still a lack of real tank experiments to solve anisotropic TV regularization problems with ADMM [12] and verify its possible distortion.

In this article, we use an iteration scheme based on the alternating direction method of multipliers (ADMM) to overcome the non-differentiability of anisotropic TV regularization method and do the tank experiment to verify its property of distortion.

2. Forward and Inverse Problems in EIT

To simplify EIT imaging process, we assume that we use 16-channel EIT system and adopt adjacent current injection and measurement pattern. To be specific, low-frequency current with magnitude I is injected through j-th neighboring electrodes pair $\left({\mathcal{E}}^j\mathrm{,}{\mathcal{E}}^{j+1}\right)$ for $j=1,\cdots ,16$ around the boundary of imaging target $\Omega$, where we denote ${\mathcal{E}}^{16+1}={\mathcal{E}}^1$. Then the potential distribution of j-th injection $u_j$ is governed by the following equations [13]:

$\{\begin{array}{l}\nabla \cdot \left(\sigma \nabla u_j\right)=\mathrm{0,}\text{in}\Omega \\ \sigma \nabla u_j\cdot n=\mathrm{0,}\text{on}\partial \Omega \backslash {\displaystyle \underset{i=1}{\overset{16}{\cup }}{\mathcal{E}}^i}\\ {\displaystyle {\int }_{{\mathcal{E}}^j}}\sigma \frac{\partial u_j}{\partial n}\text{d}s=I=-{\displaystyle {\int }_{{\mathcal{E}}^{j+1}}}\sigma \frac{\partial u_j}{\partial n}\text{d}s\\ {u_j|}_{{\mathcal{E}}_i}=\text{constant}\mathrm{,}i=\mathrm{1,2,}\cdots \mathrm{,16}\\ {\displaystyle {\int }_{{\mathcal{E}}^i}}\sigma \frac{\partial u_j}{\partial n}\text{d}s=\mathrm{0,}i\in \left\{\mathrm{1,2,}\cdots \mathrm{,}E\right\}\backslash \left\{j\mathrm{,}j+1\right\}\mathrm{.}\end{array}$

Here, $n$ is the outward unit vector to $\partial \Omega$, ds is the surface element. The voltage between $\left({\mathcal{E}}^i\mathrm{,}{\mathcal{E}}^{i+1}\right)$ subject to the j-th injection can be measured:

$V_j^i\left[\sigma \right]={u_j|}_{{\mathcal{E}}^i}-{u_j|}_{{\mathcal{E}}^{i+1}}\mathrm{.}$

Using EIT scanner we can measure the following datum:

$V={\left[V_1^3\mathrm{,}\cdots \mathrm{,}{V}_1^{15}\mathrm{,}{V}_2^4\mathrm{,}\cdots \mathrm{,}{V}_2^{16}\mathrm{,}\cdots \mathrm{,}{V}_{16}^2\mathrm{,}\cdots \mathrm{,}{V}_{16}^{14}\right]}^{\text{T}}\in R^{16\left(16-3\right)}\mathrm{.}$

where we neglect the measurement near the driving electrode to minimize the measurement error.

The relation between $\sigma$ and $V_j^i$ can be expressed approximately by the reciprocity principle

$V_j^i\left[\sigma \right]=V_i^j\left[\sigma \right]=\frac{1}{I}{\displaystyle {\int }_{\Omega }}\sigma \left(r\right)\nabla u_i\left(r\right)\cdot \nabla u_j\left(r\right)\text{d}r\mathrm{,}$ (1)

where $r=\left(x\mathrm{,}y\mathrm{,}z\right)$ is a position inside $\Omega$. The EIT problem is to reconstruct the conductivity distribution $\sigma$ using the measured voltage datum $V$ and the relation (1). However, the above equation is nonlinear. We consider to linearize the above equation. To be specific, assuming $\sigma \left(r\right)={\sigma }_0+\delta \sigma \left(r\right)$, we can approximate $\delta \sigma$ by replacing $u_j$ by $u_j^0$ :

$\delta V_j^i\left[\sigma \right]\mathrm{<mo>\; :\; </mo>}=V_j^i\left[\sigma \right]-V_j^i\left[{\sigma }_0\right]\approx \frac{1}{I}{\displaystyle {\int }_{\Omega }}\delta \sigma \left(r\right)\nabla u_i^0\left(r\right)\cdot \nabla u_j^0\left(r\right)\text{d}r\mathrm{,}$ (2)

where $u_j^0$ the potential computed under the reference conductivity ${\sigma }_0$. For the purpose of computerized reconstruction, we discretize the imaging target into finite element elements ( ${\Delta }_k,k=1,2,\cdots ,N$ ). Assume that on each element ${\Delta }_k$ the conductivity is a constant. Let ${\sigma }_k$ denote the value of $\sigma$ on the kth element. Then $\sigma$ can be approximated by $\sigma \mathrm{<mo>\; :\; </mo>}={\left({\sigma }_1\mathrm{,}{\sigma }_2\mathrm{,}\cdots \mathrm{,}{\sigma }_N\right)}^{\text{T}}\in R^N$. Thus, (2) can be written as

$S\delta \sigma =\delta V\mathrm{,}$

where $S$ is the sensitivity matrix (or jacobian matrix) given by

$S=\left(\begin{array}{ccc}& ⋮& \\ \cdots & S_k& \cdots \\ & ⋮& \end{array}\right)\in R^{208\times N}$

with $S_k\mathrm{<mo>\; :\; </mo>}=\frac{1}{I}{\displaystyle {\int }_{{\Delta }_k}}s\left(r\right)\text{d}r$, $s\left(r\right)\mathrm{<mo>\; :\; </mo>}={\left[s^{\mathrm{1,3}}\left(r\right)\mathrm{,}{s}^{\mathrm{1,4}}\left(r\right)\mathrm{,}\cdots \mathrm{,}{s}^{\mathrm{16,14}}\left(r\right)\right]}^{\text{T}}$, $s^{i\mathrm{,}j}=\nabla u_i^0\left(r\right)\cdot u_j^0\left(r\right)$.

In practice, the number of measured data $16\times \left(16-13\right)$ is fewer than N (the total number of elements for $\delta \sigma$ ). We can find an estimate of $\delta \sigma$ by minimizing the following least square problem:

$\delta {\sigma }^{\mathrm{<mo>\; *\; </mo>}}=\underset{\delta \sigma \in R^N}{\arg \min }{\Vert S\delta \sigma -\delta V\Vert }_{L^2\left(\Omega \right)}^2\mathrm{.}$

Since the above problem is ill-posed, the regularization technique is widely used to deal with this difficulty. The most common used regularization method is Tikhonov regularization (TR), which is to solve

$\delta {\sigma }_{\lambda }^{\mathrm{<mo>\; *\; </mo>}}=\underset{\delta \sigma \in R^N}{\arg \min }\left\{\frac{1}{2}{\Vert S\delta \sigma -\delta V\Vert }^2+\lambda {\Vert K\left(\delta \sigma -\delta {\sigma }_0\right)\Vert }^2\right\}\mathrm{,}$

where $K$ is the regularization matrix, $\lambda$ is the regularization parameter. However, it has an excessively smooth effect on the solution, which will blur the edge of reconstructed image.

One technique to preserve the discontinuous boundary is the Total Variation (TV) regularization. It was first introduced in [6] for image denoising and applied in EIT inverse problem in [8]. The two basic variants of TV are isotropic TV

$\delta {\sigma }_{\lambda }^{\mathrm{<mo>\; *\; </mo>}}=\underset{\delta \sigma \in R^N}{\arg \min }\left\{\frac{1}{2}{\Vert S\delta \sigma -\delta V\Vert }^2+\lambda {\Vert \sqrt{{\left|D_x\delta \sigma \right|}^2+{\left|D_y\delta \sigma \right|}^2}\Vert }_{l_1}\right\}\mathrm{,}$ (3)

and anisotropic TV

$\delta {\sigma }_{\lambda }^{\mathrm{<mo>\; *\; </mo>}}=\underset{\delta \sigma \in R^N}{\arg \min }\left\{\frac{1}{2}{\Vert S\delta \sigma -\delta V\Vert }^2+\lambda \left({\Vert D_x\delta \sigma \Vert }_{l_1}+{\Vert D_y\delta \sigma \Vert }_{l_1}\right)\right\}\mathrm{,}$ (4)

where $D_x,D_y\in R^{N\times N}$ are the first order discrete partial derivative operators in the horizontal direction and the vertical direction respectively.

To solve the non-differentiability problem (3) and (4), an auxillary smoothing parameter $\beta >0$ was introduced in the time marching method for the corresponding Euler-Lagrange equation [6] and the primal-dual algorithms [8]. However, both of the methods are too time consuming to be used in the online mode. We intend to solve non-differentiability problem (4) using ADMM algorithm. It bypasses the difficulty by using splitting scheme and soft thresholding. To be precise, we first denote $D=\left(D_x\mathrm{<mo>\; ;\; </mo>}{D}_y\right)\in R^{2N\times N}$. By introducing the a new variable $u\in R^{2N}$, (4) can be written into a constrained minimization problem:

$\{\begin{array}{l}\delta {\sigma }_{\lambda }^{\mathrm{<mo>\; *\; </mo>}}=\underset{\delta \sigma }{\arg \min }\left\{\frac{1}{2}{\Vert S\delta \sigma -\delta V\Vert }^2+\lambda {\Vert u\Vert }_{l_1}\right\}\mathrm{,}\\ \text{s}\text{.t}\mathrm{.}u=D\delta \sigma \mathrm{.}\end{array}$ (5)

The augmented Lagrangian functional for (5) is defined as:

$L_{\alpha }\left(\delta \sigma \mathrm{,}u\mathrm{<mo>\; ;\; </mo>}v\right)\mathrm{<mo>\; :\; </mo>}=\frac{1}{2}{\Vert S\delta \sigma -\delta V\Vert }^2+\lambda {\Vert u\Vert }_{l_1}+\frac{\alpha }{2}{\Vert D\delta \sigma -u\Vert }_{l_2}^2+v^{\text{T}}\left(D\delta \sigma -u\right)\mathrm{,}$

where the variable $v\in R^{2N}$ is the Lagrangian multiplier, $\alpha >0$ is a penalty parameter.

By updating $\delta \sigma$ and $u$ in an alternate order, we can get its ADMM iteration scheme:

$\{\begin{array}{l}\delta {\sigma }_{n+1}=\underset{\delta \sigma }{\arg \min }{L}_{\alpha }\left(\delta \sigma \mathrm{,}{u}_n\mathrm{,}{v}_n\right)\mathrm{<mo>\; ;\; </mo>}\\ u_{n+1}=\underset{u}{\arg \min }{L}_{\alpha }\left(\delta {\sigma }_{n+1}\mathrm{,}u\mathrm{,}{v}_n\right)\mathrm{<mo>\; ;\; </mo>}\\ v_{n+1}=v_n+\alpha \left(D\delta {\sigma }_{n+1}-u_{n+1}\right)\mathrm{.}\end{array}$

The detail of solving the above equations can refer to [9].

3. Experimental Studies

In this section, we compare the performance of isotropic TV using PDIPM and anisotropic TV using ADMM by tank experiments. The results in this article are carried out by EIDORS [14].

In order to quantitatively evaluate the regularization method, two indicators are introduced, respectively $RE\left(n\right)$ and $PSNR\left(n\right)$. Since we don’t know the true conductivity in the tank experiments, we use the conductivity reconstructed by TV regularization as a reference value. We also compare the time spent in the reconstruction process.

We did two groups of tank experiments using the Sciospec 16-channel EIT system [15] shown in Figure 1(a). The cylindrical tank with a radius of 10 cm is filled with the tap water. Sixteen electrodes are attached to the perimeter of the tank. We use the adjacent current injection and measurement pattern. The current with the maximum amplitude 1 mA and the frequency 10 kHz is injected using EIT device. The speed of the data acquisition is about 20 frames/s. The sensitivity matrix $S$ is computed under the calculated conductivity of the water using the measured voltage data ( ${\sigma }_0=0.003\text{S}/\text{m}$ ). We set the imaging plane to be the electrode plane. The finite element model of the imaging plane is shown in Figure 1(b).

In the first group of experiments, carrot, cucumber and bean curd are placed inside the tank, separately. In the second group of tests, carrot, cucumber and bean curd are placed inside the tank in pairs. The parameters for the anisotropic TV method are set to $\lambda =1\times {10}^{-8}$, $\alpha =1\times {10}^{-7}$. The parameters of TR and isotropic TV are set to be optimal empirically. The reconstructed images of single inclusion using Tikhonov regularization (TR), isotropic TV and anisotropic TV are shown in Figure 2. Figure 3 shows the reconstructed images of two inclusions cases. Table 1 and Table 2 illustrate the behavior of RE, PSNR of the single and two inclusions tank experiments, respectively. Table 3 compares the computational time of the tank experiments.

There are several observations from the results. All the three regularization methods can well capture the main feature of the inner object, including position and shape. Numerical indicators show that using ADMM to solve the anisotropic

(a) (b)

Figure 1. (a) Tank measurement setup. (b) Finite element model of imaging plane.

Figure 2. Reconstructed images of single inclusion using different regularization methods. The three rows show respectively the reconstructed images of carrot, bean curd and cucumber. The first column is the images of the tank. The last three columns are the reconstructed images using Tikhonov regularization (TR), isotropic TV with PDIPM and anisotropic TV with ADMM.

Figure 3. Reconstructed images of two different inclusions using different regularization methods. From top to bottom: 1) carrot and bean curd, 2) carrot and cucumber, 3) bean curd and cucumber. The first column is the images of the tank. The last three columns are the reconstructed images using Tikhonov regularization (TR), isotropic TV with PDIPM and anisotropic TV with ADMM.

Table 1. The behaviors of RE and PSNR for single inclusion tank model.

Table 2. The behaviors of RE and PSNR for two inclusions tank model.

Table 3. Comparison of computational time for tank experiments.

TV regularization problem can obtain an accurate image and take less time than using PDIPM. However, the images reconstructed by isotropic TV have obvious ladder effect. As a result, there exists pseudo edge. As expected, the edges of the reconstructed images using the anisotropic TV distort along the coordinate axes, whether it has one or two inclusions.

4. Conclusion and Future Work

In this article, we use experimental data to verify the difference between isotropic TV and anisotropic TV. Experimental results demonstrate that the reconstructed image with anisotropic TV regularization will cause geometric distortions along the coordinate axis. This proves that it is the model itself that causes the distortion. In the future work, we will focus on a method that can avoid the distortions along the coordinate axis and do not depend on the selection of regularization parameters.

Acknowledgements

The author would like to thank editor and referees for their valuable advice for the improvement of this article.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	kokponhoué, N.Y., Yalo, N., Akokponhoué, B.H., Houngue, R. and Agbahoungba, G. (2019) Contribution of Electrical Resistivity Tomography and Boring Technique in the Realization of Ten (10) Large Boreholes in a Crystalline Basement Rocks in the Centre-West of Benin. Journal of Geoscience and Environment Protection, 7, 114-130. https://doi.org/10.4236/gep.2019.79009
[2]	Castro, J.A., Olmo, A., Pérez, P. and Yúfera, A. (2016) Microcontroller-Based Sinusoidal Voltage Generation for Electrical Bio-Impedance Spectroscopy Applications. Journal of Computer and Communications, 4, 51-58. https://doi.org/10.4236/jcc.2016.417003
[3]	Barber, D.C. and Brown, B.H. (1988) Errors in Reconstruction of Resistivity Images Using a Linear Reconstruction Technique. Clinical Physics and Physiological Measurement, 9, 101-104. https://doi.org/10.1088/0143-0815/9/4A/017
[4]	Vauhkonen, M., Vadász, D., Karjalainen, P.A., Somersalo, E. and Kaipio, J.P. (1998) Tikhonov Regularization and Prior Information in Electrical Impedance Tomography. IEEE Transactions on Medical Imaging, 17, 285-293. https://doi.org/10.1109/42.700740
[5]	Wang, J. (2021) Non-Convex Regularization for Sparse Reconstruction of Electrical Impedance Tomography. Inverse Problems in Science and Engineering, 29, 1032-1053. https://doi.org/10.1080/17415977.2020.1820001
[6]	Rudin, L.I., Osher, S. and Fatemi, E. (1992) Nonlinear Total Variation Based Noise Removal Algorithms. Physica D: Nonlinear Phenomena, 60, 259-268. https://doi.org/10.1016/0167-2789(92)90242-F
[7]	Vogel, C.R. and Oman, M.E. (1996) Iterative Methods for Total Variation Denoising. SIAM Journal on Scientific Computing, 17, 227-238. https://doi.org/10.1137/0917016
[8]	Borsic, A., Graham, B.M., Adler, A. and Lionheart, W.R. (2010) In Vivo Impedance Imaging with Total Variation Regularization. IEEE transactions on medical imaging, 29, 44-54. https://doi.org/10.1109/TMI.2009.2022540
[9]	Wu, C. and Tai, X.C. (2010) Augmented Lagrangian Method, Dual Methods, and Split Bregman Iteration for ROF, Vectorial TV, and High Order Models. SIAM Journal on Imaging Sciences, 3, 300-339. https://doi.org/10.1137/090767558
[10]	Zhou, Z., dos Santos, G.S., Dowrick, T., Avery, J., Sun, Z., Xu, H. and Holder, D.S. (2015) Comparison of Total Variation Algorithms for Electrical Impedance Tomography. Physiological Measurement, 36, 1193-1209. https://doi.org/10.1088/0967-3334/36/6/1193
[11]	González, G., Kolehmainen, V. and Seppänen, A. (2017) Isotropic and Anisotropic Total Variation Regularization in Electrical Impedance Tomography. Computers and Mathematics with Applications, 74, 564-576. https://doi.org/10.1016/j.camwa.2017.05.004
[12]	Boyd, S., Parikh, N., Chu, E., Peleato, B. and Eckstein, J. (2011) Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers. Foundations and Trends in Machine Learning, 3, 1-122. https://doi.org/10.1561/2200000016
[13]	Lee, K., Woo, E.J. and Seo, J.K. (2017) A Fidelity-Embedded Regularization Method for Robust Electrical Impedance Tomography. IEEE Transactions on Medical Imaging, 37, 1970-1977. https://doi.org/10.1109/TMI.2017.2762741
[14]	Adler, A. and Lionheart, W.R. (2006) Uses and Abuses of EIDORS: An Extensible Software Base for EIT. Physiological Measurement, 27, S25-S42. https://doi.org/10.1088/0967-3334/27/5/S03
[15]	Sciospec Company. https://www.sciospec.com/