Monotonicity-based Electrical Impedance Tomography for Lung Imaging

TL;DR

This paper introduces a novel monotonicity-based EIT method for continuous lung monitoring, leveraging periodic data decomposition and conductivity change constraints to enhance imaging quality, validated through simulations and experiments.

Contribution

The paper proposes a new monotonicity-based approach for lung EIT imaging that improves image quality by exploiting conductivity change properties.

Findings

Enhanced lung imaging quality demonstrated in simulations.

Successful application in phantom and human experiments.

Effective separation of pulmonary, cardiac, and other signals.

Abstract

This paper presents a monotonicity-based spatiotemporal conductivity imaging method for continuous regional lung monitoring using electrical impedance tomography (EIT). The EIT data (i.e., the boundary current-voltage data) can be decomposed into pulmonary, cardiac and other parts using their different periodic natures. The time-differential current-voltage operator corresponding to the lung ventilation can be viewed as either semi-positive or semi-negative definite owing to monotonic conductivity changes within the lung regions. We used this monotonicity constraints to improve the quality of lung EIT imaging. We tested the proposed methods in numerical simulations, phantom experiments and human experiments.

Tables

Table 1. Table 1: Size of two identical cylinders at five different times.

Time	$t_1$	$t_2$	$t_3$	$t_4$	$t_5$
Height	$6cm$	$7cm$	$9cm$	$7cm$	$6cm$
Radius	$1cm$	$2.5cm$	$3.5cm$	$2.5cm$	$1cm$

Table 2. Table 2: Size of radish for experiments at different times.

Time	$t_1$	$t_2$	$t_3$	$t_4$	$t_5$
Height	$9.5cm$	$10cm$	$11cm$	$10cm$	$9.5cm$
Radius	$1.5cm$	$2.2cm$	$3.5cm$	$2.2cm$	$1.5cm$

Equations

$$\nabla \cdot \left({\sigma }^t\nabla u_t^j\right)=0\text{mbox}in\Omega ,$$

$${\sigma }^t\frac{\partial u_t^j}{\partial \mathbf{n}}=0\text{mbox}on\partial \Omega \setminus \bigcup _{k=1}^E{\mathcal{E}}_k,$$

$${\left(u_t^j+z_{j,k}{\sigma }^t\frac{\partial u_t^j}{\partial \mathbf{n}}\right)|}_{{\mathcal{E}}_k}=U^{j,k}(t),k=1,\dots ,E$$

$${\int }_{{\mathcal{E}}_k}{\sigma }^t\frac{\partial u_t^j}{\partial \mathbf{n}}ds=0\text{mbox}ifk\in \{1,\dots ,E\}\backslash \{j,j+1\},$$

$${\int }_{{\mathcal{E}}_j}{\sigma }^t\frac{\partial u_t^j}{\partial \mathbf{n}}ds=I=-{\int }_{{\mathcal{E}}_{j+1}}{\sigma }^t\frac{\partial u_t^j}{\partial \mathbf{n}}ds,$$

$$V^{j,k}(t)=U^{j,k}(t)-U^{j,k+1}(t),$$

$${\mathbb{V}}(t)=\left[\begin{array}[]{ccc}V^{1,1}(t)&\cdots&V^{1,E}(t)\\ V^{2,1}(t)&\cdots&V^{2,E}(t)\\ \vdots&\ddots&\vdots\\ V^{E,1}(t)&\cdots&V^{E,E}(t)\end{array}\right].$$

$$V^{j,k}(t)$$

$$\frac{d}{dt}{V}^{j,k}(t)$$

$$\frac{d}{dt}{V}^{j,k}(t)\approx -{\int }_{\Omega }\frac{\partial {\sigma }^t}{\partial t}\nabla u_0^j\cdot \nabla u_0^{k}dx.$$

$$\frac{\partial }{\partial t}{\sigma }^t(x)=\frac{\partial }{\partial t}{\sigma }_L^t(x)+\frac{\partial }{\partial t}{\sigma }_H^t(x)x\in \Omega ,$$

$$V_L^{j,k}(t)={\int }_{\Omega }{\sigma }_L^t\nabla u_t^j\cdot \nabla u_t^{k}dx,$$

$$\frac{d}{dt}{V}_L^{j,k}(t)=-{\int }_{\Omega }\frac{\partial {\sigma }_L^t}{\partial t}\nabla u_t^j\cdot \nabla u_t^{k}dx,$$

$${\mathbb{V}}_{L}(t)=\left[\begin{array}[]{ccc}V^{1,1}_{L}(t)&\cdots&V^{1,E}_{L}(t)\\ V^{2,1}_{L}(t)&\cdots&V^{2,E}_{L}(t)\\ \vdots&\ddots&\vdots\\ V^{E,1}_{L}(t)&\cdots&V^{E,E}_{L}(t)\end{array}\right].$$

$$\frac{\partial }{\partial t}{\sigma }_L^t(x)\ge 0for\; allx\in \Omega ,$$

$$\frac{\partial }{\partial t}{\sigma }_L^t(x)<0for\; allx\in \Omega .$$

$${\mathbf{a}}^T\frac{d}{dt}{\mathbb{V}}_L(t)\mathbf{a}=-{\int }_{\Omega }\frac{\partial {\sigma }_L^t}{\partial t}\nabla \left(\sum _{j=1}^E{a}_j{u}_t^j\right)\cdot \nabla \left(\sum _{k=1}^E{a}_k{u}_t^k\right)dx$$

$$\underset{x\in \Omega }{\inf }\frac{\partial {\sigma }_L^t}{\partial t}(x)\ge 0\Rightarrow \frac{d}{dt}{\mathbb{V}}_L⪯0,$$

$$\underset{x\in \Omega }{\sup }\frac{\partial {\sigma }_L^t}{\partial t}(x)\le 0\Rightarrow \frac{d}{dt}{\mathbb{V}}_L⪰0.$$

$${\tilde{V}}_L^{j,k}(t):$$

$$\text{mbox}either\frac{\partial {\mathbf{\sigma }}_L^t}{\partial t}\ge 0\text{mbox}or\frac{\partial {\mathbf{\sigma }}_L^t}{\partial t}<0.$$

$$\frac{d}{dt}{\mathbf{V}}_L=\mathbb{S}\frac{\partial {\mathbf{\sigma }}_L^t}{\partial t},$$

$${\mathbf{V}}_L=(V_L^{1,1},\dots ,V_L^{1,E},\dots ,V_L^{E,1},\dots ,V_L^{E,E}{)}^T.$$

$$\text{mbox}min{\Vert \mathbb{S}\frac{\partial {\mathbf{\sigma }}_L^t}{\partial t}-\frac{d}{dt}{\mathbf{V}}_L\Vert }_2^2+\lambda \Vert \frac{\partial {\mathbf{\sigma }}_L^t}{\partial t}{\Vert }_2^2\text{mbox}subjectto\frac{\partial {\mathbf{\sigma }}_L^t}{\partial t}\ge 0,$$

$$\text{mbox}min{\Vert \mathbb{S}\frac{\partial {\mathbf{\sigma }}_L^t}{\partial t}-\frac{d}{dt}{\mathbf{V}}_L\Vert }_2^2+\lambda \Vert \frac{\partial {\mathbf{\sigma }}_L^t}{\partial t}{\Vert }_2^2\text{mbox}subjectto\frac{\partial {\mathbf{\sigma }}_L^t}{\partial t}\le 0,$$

$$t_1<\cdots <t_{n-1}<t_{n+1}<\cdots ,$$

$$\text{mbox}min{\Vert \mathbb{S}({\mathbf{\sigma }}_L^{t_n}-{\mathbf{\sigma }}_L^{t_{n-1}})-({\mathbf{V}}_L^{t_n}-{\mathbf{V}}_L^{t_{n-1}})\Vert }_2^2+\lambda \Vert {\mathbf{\sigma }}_L^{t_n}-{\mathbf{\sigma }}_L^{t_{n-1}}{\Vert }_2^2$$

$$\text{mbox}subjectto{\mathbf{\sigma }}_L^{t_n}-{\mathbf{\sigma }}_L^{t_{n-1}}\ge 0,$$

$$\text{mbox}min{\Vert \mathbb{S}({\mathbf{\sigma }}_L^{t_n}-{\mathbf{\sigma }}_L^{t_{n-1}})-({\mathbf{V}}_L^{t_n}-{\mathbf{V}}_L^{t_{n-1}})\Vert }_2^2+\lambda \Vert {\mathbf{\sigma }}_L^{t_n}-{\mathbf{\sigma }}_L^{t_{n-1}}{\Vert }_2^2$$

$$\text{mbox}subjectto{\mathbf{\sigma }}_L^{t_n}-{\mathbf{\sigma }}_L^{t_{n-1}}\le 0.$$

$$S_{p}^{t}=-\left[\begin{array}[]{ccc}\int_{T_{p}}\nabla u_{t}^{1}\cdot\nabla u_{t}^{1}dx&\cdots&\int_{T_{p}}\nabla u_{t}^{1}\cdot\nabla u_{t}^{E}dx\\ \vdots&\ddots&\vdots\\ \int_{T_{p}}\nabla u_{t}^{E}\cdot\nabla u_{t}^{1}dx&\cdots&\int_{T_{p}}\nabla u_{t}^{E}\cdot\nabla u_{t}^{E}dx\end{array}\right]$$

$$-{\mathbf{a}}^T{S}_p^t\mathbf{a}={\int }_{T_p}{\left|\sum _{j=1}^E{a}_j\nabla u_t^j\right|}^{2}dx.$$

$${\int }_{\bar{\Omega }}({\sigma }_L^{t_{n-1}}-{\sigma }_L^{t_n}){\left|\sum _{j=1}^E{a}_j\nabla u_{t_n}^j\right|}^{2}dx$$

$$-\alpha {\chi }_{T_p}\ge {\sigma }_L^{t_{n-1}}-{\sigma }_L^{t_n},$$

Full text

• January 2018

Keywords: Electrical impedance tomography, continuous lung monitoring, monotonicity, inverse problem, monotonicity-based regularization

1 Introduction

Electrical impedance tomography (EIT) has received much recent attention owing to its unique ability to allow long-term, continuous monitoring of lung ventilation at the bedside [2, 3], which is not possible using other medical imaging techniques such as computerized tomography (CT), magnetic resonance imaging (MRI), ultrasound (US), or single photon emission computed tomography (SPECT). The increasing sophistication of monitoring intensive-care patients might incorporate the use of EIT system to guide strategies for protective lung ventilation via the close monitoring of patient’s lung [39, 43, 37, 52, 8, 14, 15]. Although EIT cannot compete with CT, MRI or ultrasound in terms of spatial resolution or accuracy [32], its ability to provide long-term, continuous monitoring and portability make it clinically useful.

This paper focuses mainly on the lung EIT. In EIT, multiple surface electrodes are attached to an imaging object to inject currents and measure boundary voltages. At frequencies below about 250 kHz, the potential induced by the injection current is dictated approximately by the elliptic partial differential equation $\nabla\cdot(\sigma\nabla u)=0$ inside the object where $\sigma$ is the time-varying conductivity distribution associated with lung ventilation [32, 44]. The measured EIT data can be viewed as a boundary current-voltage map from the injection current to the resulting boundary voltage, which is determined mainly by the effective conductivity distribution, the configuration of the surface electrodes, and the geometry of the imaging object. Lung EIT aims to provide dynamic images of the time-differential conductivity distribution from the time-differential current-voltage map while minimizing the forward modeling errors due to such as electrode position, boundary geometry, and uncertainty of the reference conductivity distribution.

Robust EIT reconstruction algorithms have been sought since the invention of the first EIT devices by Barber and Brown in the early 1980s [5, 6, 1, 42]. However, EIT is not yet suitable for routine clinical use because of its poor sensitivity [48]. Robust reconstruction may rely on incorporating some strong prior information into the algorithm via regularization. Unfortunately, it seems inadequate to achieve robust time difference EIT reconstructions only with regularization, such as Tikhonov regularization [44, 49] and total variation regularization [9], in spatial domain. In other words, temporal regularization should be used with the prior information in time domain, especially for time difference EIT.

Taking into account of the fundamental limitation of EIT, some strong prior information about the individual’s pulmonary function needs to be incorporated into the analysis along with the measured data. We know that the time-varying patterns in lung EIT data depend mainly on the pulmonary and cardiac cycles and diaphragm motion. Given that the pulmonary and cardiac cycles have different frequencies, we can extract patterns corresponding to the ventilation from the measured EIT data by eliminating the signals associated with cardiac motion and other processes [38, 20].

Using the extracted ventilation-induced signals, we apply the characteristic of monotonicity to the regional lung imaging. This is based on the assumption that changes of conductivity distribution induced by lung ventilation are either monotonically non-increasing or monotonically non-decreasing at any fixed time. This assumption can be understood in terms of the effective conductivity, which means the overall conductivity changes induced by ventilation can be regarded as monotonically non-increasing or monotonically non-decreasing. An increase of conductivity decreases the voltage measurements in terms of matrix definiteness. In other words, if the conductivity increases between two time frames, the measured difference data matrix should be negative semi-definite. Conversely, a conductivity decrease would increase the voltage measurements in terms of matrix definiteness. We observe that the time-derivative of the current-voltage data associated with ventilation is a non-negative operator during inhalation and a non-positive operator during exhalation.

We therefore enforce a global non-positivity constraint on the reconstructed conductivity changes during inhalation, and a global non-negativity constraint during exhalation, which is named as global monotonicity-based method (GMM). Moreover, on each pixel, we derive a local lower bound during inhalation (and a local upper bound during exhalation) using a sensitivity-based variant of the linearized monotonicity method for inclusion detection [26, 46, 47], which is named as local monotonicity-based method (LMM). Enforcing the monotonicity constraints in the image reconstruction algorithm can compensate for the inherent ill-conditioned nature of EIT [4]. The effectiveness of the proposed GMM is validated by numerical simulations, phantom experiments and human experiments while LMM is tested only by numerical simulations and phantom experiments.

2 Methods

2.1 Mathematical model

Let an imaging object occupy a two- or three-dimensional region $\Omega\subset{\mathbb{R}}^{n}$ $(n=2,3)$. We denote the time-varying conductivity at position $x$ and time $t$ by $\sigma^{t}(x)$. It is assumed that $\sigma^{t}(x)\in L^{\infty}_{+}(\Omega)$ and $\frac{\partial}{\partial t}\sigma^{t}(x)\in L^{\infty}(\Omega)$, where $L^{\infty}_{+}(\Omega)$ denotes the subspace of $L^{\infty}(\Omega)$, functions with positive essential infima. In $E$-channel EIT system as shown in Fig. 1 (a), we attach $E$ electrodes ${\mathcal{E}}_{1},{\mathcal{E}}_{2},\cdots,{\mathcal{E}}_{E}$ on its boundary $\partial\Omega$ to inject $E$ different currents ($E-1$ linearly independent) using orderly chosen pairs of electrodes. For the ease of explanation, we use the adjacent pair of electrodes ${\mathcal{E}}_{j}$ and ${\mathcal{E}}_{j+1}$ to inject $j$-th current. Here and in the following, we use the convention that ${\mathcal{E}}_{E+1}={\mathcal{E}}_{1}$. Then, the distribution of the voltage subject to the $j$-th injection current, denoted by $u_{t}^{j}$, is governed by the following complete electrode model (CEM) [45, 50]

[TABLE]

where ${\mathbf{n}}$ is the outward unit normal vector on $\partial\Omega$ , $U^{j,k}$ is the voltage on ${\mathcal{E}}_{k}$, and $z_{j,k}$ is the contact impedance of the $k$-th electrode $\mathcal{E}_{k}$. The magnitude of the current driven through the $j$-th and the $(j+1)$-th electrode is assumed to be normalized to $I=1$. The solution $u_{t}^{j}$ is unique up to constant factors. Let’s denote $u_{0}^{j}$ as the reference potential, which is the solution of (1) with the conductivity distribution $\sigma^{t}$ being replaced by the constant 1.

Assuming that the boundary voltages between all adjacent pairs of electrodes are measured, the $k$-th boundary voltage measured between ${\mathcal{E}}_{k}$ and ${\mathcal{E}}_{k+1}$ subject to the $j$-th injection current is the time-varying function

[TABLE]

where $k,j=1,2,\ldots,E$. Please note that $U^{j,E+1}(t)=U^{j,1}(t)$ and $V^{j,k}(t)=V^{k,j}(t)$. Thus, we collect $E^{2}$ number of time-varying boundary data which can be expressed as the following symmetric matrix form

[TABLE]

The inverse problem of lung EIT for monitoring the ventilation is to visualize the time varying distribution of $\sigma^{t}$ in the lung regions from the time-varying data ${\mathbb{V}}(t)$.

2.2 Boundary data separation

It is reasonable to assume $z_{j,k}\sigma^{t}\frac{\partial u_{t}^{j}}{\partial{\mathbf{n}}}\approx 0$ on the voltage sensing electrodes in (1)-(5), so that we have $u_{t}^{j}|_{{\mathcal{E}}_{k}}=U^{j,k}(t)$ for $|k-j|>1$. Therefore, for $|k-j|>1$, the conductivity is related to the measured signals by the following identities

[TABLE]

Note that if we further assume the electrodes are perfectly conductive, i.e., $z_{j,k}=0$ for all $j,k=1,2,\cdots,E$, the complete electrode model in (1)-(5) becomes the so-called shunt model [11]. Then (8) and (9) hold for $|k-j|\leq 1$ as well, which are proven in the Appendix A.1. Throughout this paper, the theory development are based on the shunt model.

However, the unknown contact impedances of current-driven electrodes are not ignorable in real application. The measurements on the current-driven electrodes are affected by the unknown contact impedances. Therefore, the identities (8) and (9) fail for $|k-j|\leq 1$.

The potentials $u_{t}^{j}$ and $u_{t}^{k}$ are related to the unknown conductivity $\sigma^{t}$ via (1). Standard linearized reconstruction methods (LM) [44] for lung EIT replace $u_{t}^{j}$ and $u_{t}^{k}$ on the right hand side of (9) by the reference potentials $u_{0}^{j}$ and $u_{0}^{k}$ and solve the resulting linear equation to determine $\frac{\partial\sigma^{t}}{\partial t}$ from $\frac{d}{dt}V^{j,k}(t)$ by relation

[TABLE]

In this work, we will propose two monotonicity-based improvements to this standard linearized reconstruction method in (10). The first improvement is based on the following data separation approach. We decompose the time-derivative of the conductivity $\frac{\partial}{\partial t}\sigma^{t}(x)$ at position $x$ into two parts

[TABLE]

where $\frac{\partial}{\partial t}\sigma_{L}^{t}$ and $\frac{\partial}{\partial t}\sigma_{H}^{t}$ are conductivity change induced by the lung ventilation, cardiac activities and other factors, respectively.

For the application of lung monitoring, we are interested in recovering the ventilation-induced conductivity change $\frac{\partial\sigma_{L}^{t}}{\partial t}$. Hence, it follows from (8)-(11) that we reasonably assume the existence of a ventilation-induced signal

[TABLE]

and aim to extract

[TABLE]

from $\frac{d}{dt}V^{j,k}(t)$. Similar as (7), the ventilation-induced data can be written as an $E\times E$ matrix

[TABLE]

Due to the clear distinction of frequencies for pulmonary (  $12~{}{min}^{-1}$) and cardiac ( $60-80~{}{min}^{-1}$) activities, the ventilation-induced signals $V^{j,k}_{L}(t)$ can be obtained by using Fast Fourier transform and applying band-pass filter to the measured data $V^{j,k}(t)$ with a selected frequency range [20], see Fig. 4 for numerical validation. The cutoff frequencies $8~{}min^{-1}$ and $20~{}min^{-1}$ were used for the band-pass filter to obtain $V^{j,k}_{L}(t)$, which are around the breath rate $12~{}{min}^{-1}$ for healthy adult [16, 20]. Fig. 1 shows a human experiment. We see that the separated ventilation-induced signals in Fig. 1 (c) is periodic and monotonically changing in the period of inhalation and exhalation. Fig. 1 (f) shows the high frequency component $V_{H}^{j,k}$ which is corresponding to the cardiac activities, measurement noise and other factors. The cardiac signal peak is not obviously shown in Fig. 1 (e) because the cardiac signal is significantly smaller in amplitude compare with the ventilation-induced signal.

2.3 Monotonicity constraints

After extracting the ventilation-induced signals $V^{j,k}_{L}(t)$ from the measured signals using a band-pass filter, we can reconstruct the ventilation-induced conductivity change $\frac{\partial\sigma_{L}^{t}}{\partial t}$ by approximating the potentials in (13) by reference potentials and solving the resulting linear equation. We improve this strategy with the following monotonicity-based constraints. By considering the air inhalation and exhalation during breath, it is reasonable to assume that $\sigma^{t}_{L}$ is monotone with respect to $t$ in the sense that the ventilation-induced conductivity is either increasing everywhere during exhalation,

[TABLE]

or decreasing everywhere during inhalation,

[TABLE]

We should mention that the above monotonicity hypothesis may not always be entirely true for all the cases due to the inhomogeneity of lung tissue and the change of outer shape of the domain. However, the monotonicity assumption makes sense when considering pressure changes during inspiration and expiration. During inhalation, air flows into the lungs as a result of decrease in the pressure inside the lungs. As a consequence, the overall effective conductivity in the lungs decreases at inhalation [12, 41]. Similarly, the overall effective conductivity in the lungs increases at exhalation. The use of the monotonicity constraint as a temporal regularization allows to deal with the ill-posed nature of EIT and improve the robustness in reconstruction.

From (13), we have

[TABLE]

for all vectors ${\bf a}=(a_{1},\ldots,a_{E})^{T}\in{\mathbb{R}}^{E}$. For a symmetric matrix $A$, we write $A\succeq 0$ if $A$ is semi-positive definite, that is, ${\bf a}^{T}A{\bf a}\geq 0$ for all vectors ${\bf a}$. We write $A\preceq 0$ if $A$ is semi-negative definite. From (17), we have the following monotonicity relations

[TABLE]

Here, the right sides of (18) and (19) should be understood in the sense of matrix definiteness.

The proposed methods evaluate whether the conductivity of the lungs is in the increasing or decreasing period and reconstruct the lung ventilation-induced conductivity change $\frac{\partial\sigma_{L}^{t}}{\partial t}$ using the additional global monotonicity-based constraint (15), resp., (16). Given the monotonicity assumption in (15) and (16), it has been shown in [26] that testing the matrix definiteness on the right hand sides of (18) and (19) can give the outer shape of a conductivity change for continuous boundary measurements. For electrode measurements one can expect to obtain an upper bound of the conductivity change, cf., [29, 30, 31, 7, 17, 18, 19] for recent works on this monotonicity-based approach in electrical impedance tomography.

2.4 EIT data sources consistency

For the monotonicity-based spatiotemporal EIT imaging, it is necessary to determine the conductivity monotone increase and monotone decrease period. According to (18) and (19), the conductivity increasing period can be determined by evaluating the definiteness of $\frac{d}{dt}{\mathbb{V}}_{L}$. However, due to the unknown contact impedances of the current driven electrodes, the tridiagonal elements of ${\mathbb{V}}_{L}$ are not reliable for determining the period. Fortunately, we made an interesting observation about EIT data sources consistency through several human experiments. This allows us to determine the monotone increase and monotone decrease period of conductivity.

This observation is obtained by a simple modification of the data $V_{L}^{j,k}(t)$ based on the equations (12), (13) and the monotonicity assumptions (15) and (16)

[TABLE]

where $sgn$ is the sign function and ${ave}(V_{L}^{j,k})$ indicates the average of $V_{L}^{j,k}$ over one period of time. Fig. 2 shows $V_{L}^{j,k}(t)$ and ${\widetilde{V}_{L}^{j,k}}(t)$, which are obtained from human experiment. As shown in Fig. 2 (b), most of ${\widetilde{V}_{L}^{j,k}}$ are of the same pattern in terms of increasing and decreasing period.

Several human experiments (seven adults in our lab) and numerical experiments showed that the time change of ${\widetilde{V}}_{L}^{j,k}$ (except for a few pair of ($j,k$)) has the same pattern as $\sigma_{L}^{t}$, as shown in Fig. 2 (b), which could help us determine the monotonicity period by observing the pattern of the modified data. It is likely that the few pair of inconsistent data in ${\widetilde{V}}_{L}^{j,k}$ is related to $\frac{\int_{\Omega}\nabla u_{0}^{j}\cdot\nabla u_{0}^{k}dx}{\int_{\Omega}|\nabla u_{0}^{j}\cdot\nabla u_{0}^{k}|dx}\leq c$ for a constant $c<0$, meaning that the sensitivity $\nabla u_{0}^{j}\cdot\nabla u_{0}^{k}$ is negative and overwhelms positive sensitivity by some amount. In Fig. 2, the blue rectangular area shows $\frac{d}{dt}(\sum_{j,k}\widetilde{V}_{L}^{j,k}(t))\geq 0$ which indicates the period of decrease in conductivity ($\frac{\partial\sigma_{L}^{t}}{\partial t}\leq 0$). On the other hand, the pink rectangular area shows $\frac{d}{dt}(\sum_{j,k}\widetilde{V}_{L}^{j,k}(t))\leq 0$ which indicates the period of increase in conductivity ($\frac{\partial\sigma_{L}^{t}}{\partial t}\geq 0$).

The consistency between modified data $\widetilde{V}_{L}^{j,k}$ and conductivity $\sigma_{L}^{t}$ is of great importance to determine the period of increase and decrease in conductivity for imposing temporal regularization constraints. Future studies will be needed to provide a rigorous proof of this observation.

2.5 Imaging with global monotonicity-based constraints

Discretizing the reference imaging domain $\bar{\Omega}$ into $\bar{\Omega}=\cup_{p=1}^{N}T_{p}$ where $T_{p}$ is the $p$-th pixel, the conductivity distribution $\sigma_{L}^{t}$ in domain $\bar{\Omega}$ can be expressed as $N\times 1$ vector ${\boldsymbol{\sigma}}_{L}^{t}$ (bold symbol). For a vector ${\bf b}=(b_{1},\ldots,b_{N})^{T}$, we write ${\bf b}\geq 0$, resp., ${\bf b}\leq 0$ if all components of the vector are non-negative, resp., non-positive. From the monotonicity assumption (15), resp., (16), the time differential $\frac{\partial{\boldsymbol{\sigma}}_{L}^{t}}{\partial t}$ satisfies the following uniform positivity (or negativity) property

[TABLE]

This paper uses this monotonicity constraint (21) to solve the ill-conditioned linear system (called linearized EIT system) that arises from approximating the potentials in (13) by reference potentials

[TABLE]

where $\mathbb{S}$ is the $E^{2}\times N$ sensitivity matrix whose $((j-1)*E+k,p)$ element is $S_{p}^{j,k}=-\int_{T_{p}}\nabla u_{0}^{j}\cdot\nabla u_{0}^{k}dx$, and ${\bf V}_{L}(t)$ is the column concatenated data set of $\mathbb{V}_{L}(t)$ given by

[TABLE]

We obtain the following global monotonicity-based constrained reconstruction method (GMM)

• •

In the conductivity increasing period

[TABLE]

• •

In the conductivity decreasing period

[TABLE]

where $\lambda\geq 0$ is the regularization parameter.

In order to implement the minimization problems of (24) and (25), we discretize time $t$ as

[TABLE]

and reconstruct the time difference ${\boldsymbol{\sigma}}_{L}^{t_{n}}-{\boldsymbol{\sigma}}_{L}^{t_{n-1}}$ by solving the following constrained minimization problems

• •

In the conductivity increasing period

[TABLE]

• •

In the conductivity decreasing period

[TABLE]

2.6 Additional local monotonicity-based constraints

From subsection 2.2 to subsection 2.5, we used a data splitting approach and a monotonicity argument to obtain a global lower (upper) bound for the conductivity change in the conductivity increasing (decreasing) period. In this subsection, we use a refined monotonicity argument to obtain also an upper bound for the conductivity change for the conductivity increasing period(and, analogously, a lower bound for the conductivity decreasing period). The following approach is a new combination of the standard linearized reconstruction method with a sensitivity based variant of the monotonicity method developed in [26]. This is motivated by the recent paper [13] (see also [24]) which uses a sensitivity based variant of the factorization method (see [25, 22, 36]) to regularize the standard linearized reconstruction method.

Given the $p$-th pixel $T_{p}\subset\bar{\Omega}$, define the following $E\times E$ matrix

[TABLE]

and note that, for ${\bf a}=(a_{1},\ldots,a_{E})^{T}\in\mathbb{R}^{E}$, we have

[TABLE]

We will use the following quantitative version of the monotonicity relations (18) and (19). For ${\bf a}=(a_{1},\ldots,a_{E})^{T}\in\mathbb{R}^{E}$,

[TABLE]

which is proven in the Appendix A.2.

Consider the conductivity increasing period, $\sigma_{L}^{t_{n-1}}\leq\sigma_{L}^{t_{n}}$. If $\alpha\geq 0$ fulfills

[TABLE]

then, by (30) and (31),

[TABLE]

where $\chi_{T_{p}}$ is the characteristic function of subset $T_{p}\subset\bar{\Omega}$.

Under the assumption that the conductivity change is constant on the pixel $T_{p}$ we hence obtain by contraposition that

[TABLE]

Hence, for $p$-th pixel we aim to find a smallest possible $\alpha\geq 0$ with $\alpha S_{p}^{t_{n}}\not\succeq{\mathbb{V}}_{L}(t_{n})-{\mathbb{V}}_{L}(t_{n-1})$. Collecting the values for all pixels in a vector $\boldsymbol{\alpha}$, we arrive at the following local monotonicity-based constrained reconstruction method (LMM).

• •

In the conductivity increasing period

[TABLE]

Analogously, we obtain for the conductivity decreasing period, $\sigma_{L}^{t_{n}}\leq\sigma_{L}^{t_{n-1}}$, that, for all $\beta\geq 0$,

[TABLE]

so that we aim to find a smallest possible $\beta\geq 0$ with $-\beta S_{p}^{t_{n-1}}\not\preceq{\mathbb{V}}_{L}(t_{n})-{\mathbb{V}}_{L}(t_{n-1})$ for $p$-th pixel and collect the values for all pixels in a vector $\boldsymbol{\beta}$.

• •

In the conductivity decreasing period

[TABLE]

For the implementation of these additional local monotonicity-based constraints, we approximate the potentials $u_{t}^{j}$ in the definition of $S_{p}^{t}$ by reference potentials $u_{0}^{j}$ as this is done in the above-described linearized reconstruction methods. Accordingly, $S_{p}^{t_{n}}$ and $S_{p}^{t_{n-1}}$ are approximated by the matrix $S_{p}$, which is obtained by a rearrangement of the $p$-th column of the sensitivity matrix $\mathbb{S}$ defined in subsection 2.5. The smallest possible values for $\alpha$, resp., $\beta$ on each pixel are determined by a binary search method.

3 Algorithm summary

The constrained minimization problem with the non-negative constraint (26) and (28) can be viewed as non-negative least square problem [51, 10]. Numerous experiments verified that enforcing a non-negativity constraint could lead to more accurate approximate solution[21, 40]. We adopt Algorithm 1 to solve problems (26) and (28).

To solve (33) and (36), we used similar algorithm as solving (26) and (28) with minor modifications of steps ii. 3. (a) and ii. 3. (b) in Algorithm 1. Precisely, the step ii. 3. (a) and ii. 3. (b) in algorithm 1 should be adapted as that in Algorithm 2.

By adding monotonicity constraints of the conductivity change, the inverse problem will get better posed and the solution of (26), (28) and (33), (36) will provide more accurate and stable images compare with the standard LM. Numerical, phantom and human experiments results will be shown in the following sections.

4 Numerical simulations

In this section, we show various numerical experiments to test the proposed GMM and LMM and compare the results with that of the standard LM.

4.1 2D numerical example

Numerical validation for the proposed monotonicity-based reconstruction methods were performed on a chest model occupying two dimensional domain $\bar{\Omega}$. See Fig. 3. Three elliptic perturbations were added, one heart modeled in the region $H=\{(x,y)|(x-3.7)^{2}+(y-1.7)^{2}\leq 0.7^{2}\}$ and two lungs modeled in the regions $L_{1}=\{(x,y)|\frac{(x-1.8)^{2}}{0.7^{2}}+\frac{(y-3)^{2}}{1.25^{2}}\leq 1\}$ and $L_{2}=\{(x,y)|\frac{(x-5.6)^{2}}{0.7^{2}}+\frac{(y-3)^{2}}{1.25^{2}}\leq 1\}$. The conductivity in the regions of $H,L_{1}$ and $L_{2}$ were set to be time dependent; $\sigma^{t}_{H}=2+0.5\cos(3\pi t)$ in $H$ and $\sigma^{t}_{L}=0.5+0.3\cos(\pi t)$ in $L_{1}\cup L_{2}$. The background conductivity is $1$.

Finite element mesh with 766 triangular elements and 452 nodes was generated. Boundary voltages were simulated using 16-channel EIT system. There are a total of 256 measurements each time frame by injecting currents and measuring boundary voltage adjacently. For simplicity, the contact impedances were ignored in the forward model.

With the conductivity shown in Fig. 3 (b), the simulated boundary data ${\bf V}$ is separated into two parts by applying band-pass filter, one part ${\bf V}_{L}$ is related to $\sigma^{t}_{L}$ and the other part ${\bf V}_{H}$ is related to $\sigma^{t}_{H}$, see Fig. 4. The cutoff frequency for ${\bf V}_{L}$ is 0.3 Hz and 0.7 Hz, while the cutoff frequency for ${\bf V}_{H}$ is 0.8 Hz and 1.7 Hz. There are two peak frequency for ${\bf V}_{H}$ in Fig. 4 (b) which may be because of the nonlinearity of conductivity to voltage map. Fig. 4 shows that ventilation-induced data can be filtered out due to the nature of frequencies distinction between pulmonary and cardiac activities.

Fig. 5 shows reconstructed images at ten time fames: $t_{1}=0.1s,0.3s,\ldots,1.9s$. The first column of Fig. 5 shows the true conductivity distributions. The red color in the images indicates the increasing conductivity or positive values while blue color indicates the decreasing conductivity or negative values. The second column shows EIT images using the standard LM with the unfiltered data ${\bf V}$. The third column shows EIT reconstructions using the standard LM with the filtered data ${\bf V}_{L}$. The fourth row shows EIT images using GMM. The fifth row displays images using LMM. Fig. 6 shows images of $\boldsymbol{\alpha}$ and $-\boldsymbol{\beta}$ which are used for LMM in Fig. 5. It shows clearly that $\boldsymbol{\alpha}$ and $-\boldsymbol{\beta}$ provides the upper and lower bound of the conductivity change, respectively. Moreover, numerical experiments with noisy data were conducted. Whiten noise was added to the voltage data to make the signal-to-noise ratio as 50dB. Fig. 7 is the results with the noisy data.

As shown in Fig. 5 and Fig. 7, results from the standard LM using unfiltered and filtered data have many artifacts. On the other hand, the proposed GMM and LMM provide reconstructions with less artifacts. We chose best regularization $\lambda$ for the standard LM from many tests by varying $\lambda$. For the proposed GMM and LMM, very small $\lambda$ was chosen for these simulations and similarly for the subsequent simulations and experiments. Moreover, the performance of GMM and LMM are insensitive to the choice of $\lambda$ provided $\lambda$ is sufficiently small. However, the performance of the standard LM is highly depending on the choice of $\lambda$. We think the monotonicity constraints play an important role as a temporal regularization in the proposed methods.

4.2 3D numerical example

For 3D numerical simulation, thorax model was used as displayed in Fig. 8. To compute the forward solution, 4843 tetrahedra elements and 1195 nodes were used. The thorax model $\bar{\Omega}$ was contained in the cuboid $26cm\times 17cm\times 12cm$, and 32 electrodes were placed as shown in Fig. 8. For each pair of current injection through adjacent electrodes, there are 32 boundary voltage measurements between all pair of adjacent electrodes which lead to a total of 1024 measurements each time frame. Inside the domain $\bar{\Omega}$, we put two identical cylinders whose size change with time. The size of cylinders at five different times $t_{1},\cdots,t_{5}$ are given in TABLE 1. And the vertical center of all the cylinders are at $z=6cm$. The conductivity value of these two cylinders was set to be $0.5$ and background conductivity was $1$.

Fig. 9 and Fig. 10 show EIT reconstructions of $xy$-slice images at $z=9cm$ and $z=3cm$, respectively. The first rows of Fig. 9 and Fig. 10 show the true distribution of conductivity changes. The second rows of Fig. 9 and Fig. 10 show the reconstructed images using the standard LM. The third rows of Fig. 9 and Fig. 10 show reconstructions using GMM. The fourth rows of Fig. 9 and Fig. 10 display the results using LMM.

Fig. 9 and Fig. 10 show that the reconstructions from the standard LM have many artifacts on the background while GMM and LMM have less artifacts. Moreover, the size of reconstructions from LMM are more accurate than that from LM and GMM.

For a quantitative analysis, we evaluated the performance figures of merit including the position error (PE, the smaller and less variable the better), resolution (RES, the smaller and more uniform the better ) and shape deformation (SD, the lower and more uniform the better) of the reconstructions introduced in GREIT [3]. Fig. 11 shows PE, RES and SD of the results in Fig. 10. The dashed blue lines, red lines and green lines show the results using the standard LM, GMM and LMM, respectively. These quantitative evaluations show that GMM and LMM have better performance than LM since the results from GMM and LMM have less position error, resolution and shape deformation. LMM performs slightly better than GMM as we see from Fig. 9, Fig. 10 and Fig. 11.

5 Phantom and human experiments

The proposed GMM and LMM were tested with phantom experiments. And for human experiments, results of GMM are presented. It is difficult to apply the propose LMM to phantom and human experiments, because the tridiagonal elements of the data matrix in (7) as well as in (14) are not reliable due to the unknown contact impedances of the current-driven electrodes. However, we include the reconstructions using LMM of phantom experiments to see how the results are affected by the unknown contact impedances of the current-driven electrodes.

5.1 Phantom experiments

Experiments were conducted with thorax shape phantom using 32-channel Swisstom pioneer EIT system [54]. To simulate the monotonically conductivity change of the lungs, five well trimmed radishes were used as shown in Fig.12. The detailed information of the size of radishes are shown in TABLE 2. The size of thorax phantom is bounded by cuboid $26cm\times 17cm\times 12cm$ with two rings of electrodes. The phantom was filled with saline solution of conductivity $0.32s/m$. 5mA current was injected adjacently at frequency 50kHz and the adjacent boundary voltage data was measured.

The Swisstom EIT system provides a full set of measurements, $V^{j,k}(t)$, $j,k=1,\ldots,32$ for each time frame. However, the voltage data on current-driven electrodes are known to be prone to errors and affected by unknown contact impedances, so that only $32\times 29$ measurements ($V^{j,k}(t)$ with $|j-k|>1$) are used for LM and GMM. We applied the full data of $32\times 32$ measurements to LMM and find that the erroneous data from the current-driven electrodes can still provides acceptable results.

Fig.13 and Fig.14 display the time difference conductivity images for phantom experiments at slices $z=9cm$ and $z=3cm$, respectively. The first rows of Fig.13 and Fig.14 show the reconstructed images using the standard LM. The second rows of Fig.13 and Fig.14 show the reconstructed images using GMM. The third rows of Fig.13 and Fig.14 show the reconstructed images using LMM with full data.

Fig.13 and Fig.14 show that even the results from GMM are not perfect, they are relatively good compare with the results from the standard LM considering the shape and position of the objects.

Similar as the quantitative analysis for 3D numerical simulations in Fig. 11, the quantitative evaluation for the phantom experiments are shown in Fig. 15. We see that results from GMM have less position error and shape deformation. The results from LMM and LM have comparative position error and shape deformation.

Please note that the data measured on the current-driven electrodes are removed for the reconstruction of the standard LM and proposed GMM. Full data should be used for the proposed LMM since the LMM requires measurements on the current-driven electrodes. This is a preliminary result to indicate that LMM will also be applicable to real data, but it is not comparable to GMM or LM as the voltage data on driven electrodes are being used for LMM.

We would like to emphasize that it is still an open problem to get reliable full measurements without being affected by the unknown contact impedances of current-driven electrodes. Solving this contact impedance problem is the topic of ongoing research which may rely on two approaches. One approach is the data interpolation using the measured data on the non-current-driven electrodes which is shown in a very recent work [28]. The other approach is to use the compound electrode [53, 33] with which different parts of electrode are used for the current injection and voltage measuring.

5.2 Human experiments

32-channel Swisstom BB2 EIT system [55] was used to conduct human experiments. A ring of 32 electrodes were attached to the thorax as shown in Fig. 1 (a). We used the Sense 3D scanner (Cubify 3D System) to scan the 3D geometry of the body. This 3D image with our Matlab GUI software allows us to extract the boundary geometry and electrode positions. See Fig. 16 (a) and (b). Current of 1mA at 150kHz was injected with skipping 4 electrodes, and boundary voltages were measured at 10 frames per second. We applied GMM using a band-pass filtered data surrounding respiratory rate [20]. The source consistency observation in section 2.4 was used to determine the monotonic increase period and decrease period of conductivity.

Fig. 17 shows the results using the standard LM and the proposed GMM. The first row of Fig. 17 shows the transformed voltage. The second row shows the reconstruction using the standard LM with unfiltered data ${\bf V}$ and the third row shows the results from the standard LM with filtered data ${\bf V}_{L}$. The fourth row of Fig. 17 shows the results from GMM. The dark blue color indicates decease of conductivity while light blue color implies increase of conductivity in the second, third and fourth rows in Fig. 17. We compared the profiles of results from the standard LM and the proposed GMM in the last row of Fig. 17.

Fig. 17 shows that results from GMM have less artifacts than that from the standard LM. The major advantage of the proposed method is its robustness; the reconstructed images are continuously depending on the measured data.

6 Conclusions

The proposed reconstruction methods allow lung EIT to visualize the monotonically varying conductivity distribution during inhalation and exhalation. It is based on the assumption that lung conductivity monotonically decreases during inhalation (due to the air flowing into the lungs) and monotonically increases during exhalation (due to the air leaving the lungs). The periodicity of lung ventilation is used to extract its associated current-voltage data; all voltage differences between electrodes increase during inhalation and decrease during exhalation in terms of matrix definiteness regardless of the injection currents. This correlation between the time-differential of the current-voltage map and the changes of conductivity can be enforced in the reconstruction algorithm as a constraint.

We know that the inverse problem of EIT is ill-posed, and therefore any least square method by data-fitting alone may not be able to provide useful images. This means that the boundary current-voltage data alone are insufficient to achieve robust reconstructions for making clinically useful images. Due to the inherent methodological limitation, the EIT reconstruction algorithm requires a strategy that balances data fitting and a suitable regularization by imposing a certain constraint on the expected image. The proposed methods use monotonicity as such a regularization.

Although EIT has limited resolution, its unique advantage lies in its capability for continuous monitoring at the bedside. Given the drawbacks of EIT (e.g., technical difficulties of the ill-posedness related to EIT data being insufficient to probe local conductivity changes), we need to focus on a robust reconstruction method to allow this technique to provide indispensable information in clinical medicine.

Acknowledgment

J. K. Seo and L. Zhou were supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MEST) (No. 2011-0028868, 2012R1A2A1A03670512).

Appendix A.

A.1. Proof of the identities (8) and (9).

Assume $\sigma^{t}(x)\in L^{\infty}_{+}(\Omega)$ and $\frac{\partial}{\partial t}\sigma^{t}(x)\in L^{\infty}(\Omega)$, where $L^{\infty}_{+}(\Omega)$ denotes the subspace of $L^{\infty}(\Omega)$, functions with positive essential infima. By considering the shunt model which ignores the contact impedance between electrodes and the imaging domain, we prove the identities (8) and (9) for all the $j,k=1,2,\cdots,E$.

First note that

[TABLE]

which shows (8) and $V^{j,k}(t)=V^{k,j}(t)$. With the same argument we obtain that

[TABLE]

and

[TABLE]

Hence,

[TABLE]

and it follows from $\nabla u_{t+\delta t}^{k}\approx\frac{\sigma^{t}}{\sigma^{t+\delta t}}\nabla u_{t}^{k}$ [24] and $\lim_{\delta t\rightarrow 0}\frac{\sigma^{t}}{\sigma^{t+\delta t}}=1$ that

[TABLE]

which is the asserted identity (9). $\Box$

A.2. Proof of estimate (31).

This type of monotonicity estimate goes back to [35, 34], see also [23, 26, 27] for recent applications.

Let $a_{1},\ldots,a_{E}\in{\mathbb{R}}$. For brevity, we write $u_{t}^{a}:=\sum_{j=1}^{E}a_{j}u_{t}^{j}$. As in the proof of (8) and (9) we have that

[TABLE]

From

[TABLE]

it follows that

[TABLE]

which is the first inequality in (31). The second inequality in (31) follows from interchanging $t_{n}$ and $t_{n-1}$.