Fraction model
A fraction model is a representation of the tissue distribution within an object [21, 27]. It assumes the object consists of a limited number of tissues. After meshing the object, volume fractions or concentration values of different tissues can be determined for each element. Meanwhile, if the conductivity spectra of all tissues are known, the conductivity distribution is also obtained by transforming the fraction model.
The model can be summarized as follows.
- 1.
The domain is composed of T kinds of distinct tissues, \(t_{1} \ldots t_{j} \ldots t_{T}\).
- 2.
The impedance spectra of each tissue are known. For example, the conductivity of the \(j{\text{th}}\) tissue for the \(i{\text{th}}\) frequency is \(\sigma_{ij} = \varepsilon^{{t_{j} }} \left( {\omega_{i} } \right)\).
- 3.
The conductivity of the \(n{\text{th}}\) element is equal to the linear weighted sum of the conductivities of the component tissues [28].
$$\sigma_{n} \left( {\omega_{i} } \right) = \mathop \sum \limits_{j = 1}^{T} f_{nj} *\sigma_{ij}$$
(3)
where \(0\, \le \,f_{nj} \, \le \,1\, {\text{and }}\,\mathop \sum \nolimits_{j = 1}^{T} f_{nj} \, = \,1\). The weighting value \(f_{nj}\) is the volume fraction of the \(j{\text{th}}\) tissue in the \(n{\text{th}}\) element.
Let us take a two-dimensional region \({\varvec{\Omega}}\) that is composed of T kinds of tissues as an example. First, discretization of the domain is performed using the finite-element method (FEM), and \(N\) is the number of elements. As the object contains \(T\) kinds of tissues, there is a fraction value \(f_{nj}\) in each element for every tissue. Subsequently, the fraction matrix \(\varvec{F}_{\text{matrix}} \in \varvec{R}^{T*N}\) is obtained, among which the \(n{\text{th}}\) column is the proportion of each tissue inside the \(n{\text{th}}\) element. Then, the matrix is vectorized by column to obtain the fraction vector, \(\varvec{F} = \varvec{vec}\left( {\varvec{F}_{\text{matrix}} } \right),\varvec{F} \in \varvec{R}^{{\left( {T*N} \right)*1}}\):
Next, the element-based and piecewise constant conductivity distribution is approximated using the fraction model. The conductivity vector can then be represented as \(\varvec{\sigma}\left( {\omega_{i} } \right) \in \varvec{R}^{N*1}\). Finally, currents are injected at the boundary of \(M\) frequencies (\(\omega_{1} \ldots \omega_{i} \ldots \omega_{M}\)), and for each frequency a boundary voltage vector \(\varvec{\upsilon}\left( {\omega_{i} } \right) \in \varvec{R}^{K *1}\) is acquired at a given time, where \(K\) is the number of measurements. By comparing this vector with the background frame, a boundary voltage difference vector \(\Delta\varvec{\upsilon}\left( {\omega_{i} } \right)\) is obtained. The whole difference vector \(\Delta\varvec{\upsilon}\) can be obtained by longitudinally stitching the \(M\) vectors (\(\Delta\varvec{\upsilon}\left( {\omega_{1} } \right) \cdots \Delta\varvec{\upsilon}\left( {\omega_{i} } \right) \cdots \Delta\varvec{\upsilon}\left( {\omega_{M} } \right)\)).
Forward problem
Based on the fraction model described above, the discrete forward problem is deduced as follows:
First, the linear relationship between fraction value and conductivity is obtained according to assumption (3):
$$\varvec{A}\left( {\omega_{i} } \right)\varvec{F} =\varvec{\sigma}\left( {\omega_{i} } \right)$$
(4)
where \(\varvec{A}\left( {\omega_{i} } \right) \in \varvec{R}^{{N*\left( {T*N} \right)}}\) is a coefficient matrix formed by the conductivity spectra of \(T\) kinds of tissues:
$$A(\omega_{i} ) = \left[ {\begin{array}{*{20}c} {\varepsilon^{{t_{1} }} \left( {\omega_{i} } \right)} & {\varepsilon^{{t_{2} }} \left( {\omega_{i} } \right)} & \cdots & {\varepsilon^{{t_{\text{T}} }} \left( {\omega_{i} } \right)} & 0 & 0 & \cdots & 0 & \cdots & 0 & 0 & \cdots & 0 \\ 0 & 0 & \cdots & 0 & {\varepsilon^{{t_{1} }} \left( {\omega_{i} } \right)} & {\varepsilon^{{t_{2} }} \left( {\omega_{i} } \right)} & \cdots & {\varepsilon^{{t_{\text{T}} }} \left( {\omega_{i} } \right)} & \cdots & 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 & 0 & 0 & \cdots & 0 & \cdots & {\varepsilon^{{t_{1} }} \left( {\omega_{i} } \right)} & {\varepsilon^{{t_{2} }} \left( {\omega_{i} } \right)} & \cdots & {\varepsilon^{{t_{\text{T}} }} \left( {\omega_{i} } \right)} \\ \end{array} } \right]$$
Second, in the traditional tdEIT algorithm, the discrete forward problem without noise is as follows [29]:
$$\varvec{J}\left( {\omega_{i} } \right)\Delta\varvec{\sigma}\left( {\omega_{i} } \right) = \Delta\varvec{\upsilon}\left( {\omega_{i} } \right)$$
(5)
where the Jacobian matrix \(\varvec{J}\left( {\omega_{i} } \right) \in \varvec{R}^{K*N}\) is the first derivative of \(\varvec{\upsilon}\left( {\omega_{i} } \right)\) with respect to \(\varvec{\sigma}\left( {\omega_{i} } \right),\) and can be calculated using the standard derivation method or the compensation theorem method [1, 4].
By substituting Eq. (4) into Eq. (5), the ideal discrete forward problem of fraction change reconstruction is obtained:
$$\varvec{J}\left( {\omega_{i} } \right)\varvec{A}\left( {\omega_{i} } \right)\Delta \varvec{F} = \Delta\varvec{\upsilon}\left( {\omega_{i} } \right)$$
(6)
where both \(\varvec{J}\) and \(\varvec{A}\) vary with frequency. Contrarily, \(\Delta \varvec{F}\) is independent of frequency, which makes it possible to simultaneously employ multi-frequency measurements to reconstruct one-frame image by viewing \(\Delta \varvec{F}\) as the unknown parameter.
Inverse problem
To simultaneously utilize the multi-frequency information, the following objective function is constructed to minimize the norm of the data mismatch under all frequencies:
$${\varvec{\Phi}} = \frac{1}{2}||\varvec{S}\Delta \varvec{F} - \Delta\varvec{\upsilon}||^{2}$$
(7)
\(\varvec{S}\) is an assembly matrix that includes \(\varvec{M}\) element matrices, each of which is obtained by multiplying the \(\varvec{J}\) and \(\varvec{A}\) under the same frequency:
$$\varvec{S} = \left[ {\begin{array}{*{20}c} {J\left( {\omega_{1} } \right)A\left( {\omega_{1} } \right)} \\ {J\left( {\omega_{2} } \right)A\left( {\omega_{2} } \right)} \\ { \ldots \ldots } \\ {J\left( {\omega_{M} } \right)A\left( {\omega_{M} } \right)} \\ \end{array} } \right]$$
(8)
The objective function can also be written in an equivalent form [30] (Eq. 9):
$${\varvec{\Phi}} = \frac{1}{2}\mathop \sum \limits_{i = 1}^{M}|| \varvec{J}\left( {\omega_{i} } \right)\varvec{A}\left( {\omega_{i} } \right)\Delta \varvec{F} - \Delta\varvec{\upsilon}\left( {\omega_{i} }\right)||^{2}$$
(9)
Using the standard Tikhonov regularization method [3], the objective function becomes as follows:
$${\varvec{\Phi}} = \frac{1}{2}\left[||{\varvec{S}\Delta \varvec{F} - \Delta\varvec{\upsilon}||^{2} + \lambda ||\varvec{R}\Delta \varvec{F}||^{2} } \right]$$
(10)
The regularization parameter \(\lambda\) controls the trade-off between fidelity and robustness [31, 32]. The L-curve method is chosen for both SC and DLS to determine the regularization parameter \(\lambda\) because it can be implemented efficiently and easily [33,34,35]. \(\varvec{R}\) is the regularization matrix and represents the inverse of the covariance of the expected image. The Standard Tikhonov method is chosen for both SC and DLS, which sets \(\varvec{R} = {\text{diag}}\left( {\varvec{S}^{T} \varvec{S}} \right)\) in SC and \(\varvec{R} = {\text{diag}}\left( {\varvec{J}^{T} \varvec{J}} \right)\) in DLS.
To satisfy the constraint \(\mathop \sum \nolimits_{j = 1}^{T} f_{nj} = 1 \forall n\), a substitution is made in the objective function (\(f_{n1} = 1 - \mathop \sum \nolimits_{j = 2}^{T} f_{nj}\)). Then, the unknown parameter becomes \(\Delta \varvec{F}_{{\varvec{T} - 1}}\) (\(t_{2} \ldots t_{j} \ldots t_{T}\)). Correspondingly, matrix \(\varvec{S}\) will change to \(\varvec{S}^{\prime}\). The final objective is as follows:
$${\mathbf{\varPhi^{\prime}}} = \frac{1}{2}\left[ {\varvec{S}^{\prime}\Delta \varvec{F}_{{{\text{T}} - 1}} - \Delta\varvec{\upsilon}^{2} + \lambda \varvec{R}\Delta \varvec{F}_{{{\text{T}} - 1}}^{2} } \right]$$
(11)
The solution is obtained by setting the first derivative of \({\varvec{\Phi}} '\) to zero [36]:
$$\Delta \varvec{F}_{{T - 1}} = \left( {\varvec{{S}^{\prime}}^{T} {\text{*}}\varvec{S}^{\prime} + \lambda \varvec{R}} \right)^{{ - 1}} {\text{*}}\varvec{{S}^{\prime}}^{T} {\text{*}}\Delta \varvec{\upsilon }$$
(12)
The temporary fraction vector \(\tilde{\varvec{F}}\) of \(\varvec{T}\) tissues is obtained by combining initial fraction distribution \(\varvec{F}^{0}\) and \(\Delta \varvec{F}_{T - 1}\). The final \(\varvec{F}\) is obtained by adding the constraints of Eq. (13) to \(\tilde{\varvec{F}}\):
$$f_{\text{ni}} \, = \,\left\{ \begin{aligned} & 0\quad {\text{if}}\;\widetilde{{f_{{{\text{ni}}}} }} \le 0 \\ & 1\quad {\text{if}}\;\widetilde{{f_{{{\text{ni}}}} }} \ge 1 \\ & \widetilde{{f_{{{\text{ni}}}} }} \quad {\text{otherwise}} \\ \end{aligned} \right.$$
(13)
The fraction change \(\Delta \varvec{F}\) is derived from Eq. (14):
$$\Delta \varvec{F} = \varvec{ F} - \varvec{F}^{0}$$
(14)
Under the condition of satisfying a certain imaging speed, the performance of SC can be improved by linear iteration of limited steps (\({\text{NUM}}\)). Then, the reconstruction process is summarized as follows:
- Step 0::
Initialize fraction vector \(\varvec{F}^{0}\) and the number of iteration steps to \(k = 1\)
- Step 1::
If \(k \le {\text{NUM}},\) repeat Steps 2–5; otherwise, output \(\Delta \varvec{F}\)
- Step 2::
Set the first derivative of the objective function to zero and get \(\Delta \varvec{F}_{T - 1}^{k}\)
- Step 3::
Combine \(\varvec{F}^{k - 1}\) and \(\Delta \varvec{F}_{T - 1}^{k}\) to get \(\tilde{\varvec{F}}^{k}\)
- Step 4::
Impose the constraints in Eq. (13) on \(\tilde{\varvec{F}}^{k}\) to get \(\varvec{F}^{k}\) and update the objective function and \(\Delta\varvec{\upsilon}\)
- Step 5::
Set \(\Delta \varvec{F} = \varvec{F}^{k} - \varvec{F}^{0}\) and \(k = k + 1\); return to Step 1
Finally, to compare with DLS, the final fraction change needs to be transformed into conductivity change using the linear relationship in Eq. (4).
Matrix \(\varvec{S}\) evaluation
In this section, the reasons why SC can improve image quality are theorized and two indicators of \(\varvec{S}\) are proposed to verify them.
First, SC needs to impose the constraints in Eq. (13) on the calculated solution \(\tilde{\varvec{F}}\). According to the linear relationship in Eq. (4), the corresponding conductivity will also be limited. This inherent characteristic of SC limits image artifacts to a certain extent and improves image quality.
Second, SC simultaneously uses multi-frequency data to reconstruct one shot image. This is in direct contrast to the DLS, which uses the data of only one frequency. The number of equations is \(M\) times that of the DLS algorithm; that is, \(M*K\). The unknowns are a multiple of \(\left( {T - 1} \right);\) that is, \(\left( {T - 1} \right)*N\). When \(M*K \ge \left( {T - 1} \right)*N\) is satisfied, SC is expected to increase the number of independent equations, reduce the degrees of freedom and morbidity of the inverse problem, and improve image quality. To verify this conjecture, the following two indicators are proposed for \(\varvec{S}\), which is inversed in the reconstruction process [37]:
Rank
Rank corresponds to the maximum number of linearly independent columns or rows of a matrix. It is thus a measure of the “nondegenerateness” of the system encoded by this matrix. There are numerous equivalent solutions for rank. In this study, the Singular Value Decomposition (SVD) method is adopted.
Condition number
In the field of numerical analysis, the condition number of a matrix signifies how much the output vector of the system encoded by this matrix can change for a small change in the input vector. A problem with a high condition number is said to be ill-conditioned [36], which means that a small amount of noise at its input will have a significant impact on its output. It is defined as the ratio between the maximum and the minimum singular value of this matrix.
Image evaluation
We considered the case of two tissues and one target occupied by \(t_{2}\), while the background is occupied by \(t_{1}\). In accordance with the objective quantification method proposed by Adler [5], four indicators—position error, deformation error, image noise, and total image error—are used to evaluate the performance of the two algorithm. First, the region of perturbation (RP) is defined, which satisfies the following two conditions [21]:
- a.
Values larger than 50% of the maximum displacement from the mean value of the image [21].
- b.
The largest connected cluster of voxels in (a).
Position error (PE)
PE is defined as the ratio of the displacement distance of the centroid of the reconstructed perturbation from the real position to the diameter of the mesh \(d^{\text{MESH}}\):
$${\text{PE}} = \frac{{\left| {d^{\text{RP}} - d^{\text{REAL}} } \right|}}{{d^{\text{MESH}} }}$$
(15)
where \(d^{\text{RP}}\) and \(d^{\text{REAL}}\), respectively, represent the distance from the centroid of the reconstructed perturbation and that of the real perturbation to the center of the model.
Deformation error (DE)
DE is defined as the averaged difference in X–Y dimensions between the real and reconstructed perturbations.
$${\text{SD}} = \frac{1}{2}\left( {\frac{{\left| {l_{x}^{\text{RP}} - l_{x}^{\text{REAL}} } \right| + \left| {l_{y}^{\text{RP}} - l_{y}^{\text{REAL}} } \right|}}{{d^{\text{MESH}} }}} \right)$$
(16)
Among them, \(\left( {l_{x}^{\text{RP}} ,l_{y}^{\text{RP}} } \right)\) and \(\left( {l_{x}^{\text{REAL}} ,l_{y}^{\text{REAL}} } \right)\), respectively, represent the width on the X–Y dimensions of the reconstructed and real perturbations.
Image noise (IN)
IN is defined as the inverse of the contrast-to-noise ratio (CNR) between the real perturbation and the background:
$${\text{IN}} = \frac{{\sqrt {\frac{1}{{N^{B} - 1}}\mathop \sum \nolimits_{n \in B} \left( {\Delta f_{n2} - \Delta \tilde{f}_{2}^{B} } \right)^{2} } }}{{\left| {\Delta \tilde{f}_{2}^{\text{REAL}} - \Delta \tilde{f}_{2}^{B} } \right|}}$$
(17)
where \(\Delta \tilde{f}_{2}^{\text{REAL}}\) and \(\Delta \tilde{f}_{2}^{B}\) are the mean fraction change of the real perturbation and the background, and \(N^{B}\) is the number of elements of the background.
Total image error (TE)
Overall, the error of an image is described by adding the three errors described above:
$${\text{TE}}\,{ = }\,{\text{PE}}\,{ + }\,{\text{SD}}\,{ + }\,{\text{IN}}$$
(18)
Experimental setup
Tissue impedance spectra
In numerical validation, the conductivity spectra of human brain tissues taken from the literature were used [38]. The background and perturbation conductivities were set to the values for normal brain and ischemic brain.
In the phantom experiment, fruit and vegetable objects with frequency-dependent conductivities were used to mimic the properties of live tissues. The background medium was a mixture of 0.1% NaCl solution and pomelo granules, and the perturbation was a cucumber segment of diameter 2.5 cm and height 7 cm. Conductivity measurements were acquired with a Solartron 1294 impedance analyzer for 25 frequencies in the range 1–200 kHz using Ag–AgCl electrodes. This frequency range is consistent with our imaging system, FMMU-EIT5 [39]. The cucumber was cut into 3 mm × 3 mm cubes for placement in the measurement box, which had a diameter of 1 cm and a height of 1.2 cm and was connected to the analyzer via a conductor. Three different samples were selected for each tissue, and measurements were made five times for each sample at room temperature. The conductivity was calculated using the specific calculation formula for this measurement box, and the final conductivity spectra were obtained by averaging.
Three frequency points were intercepted from the broad spectra obtained via the two different methods mentioned above to generate the final spectra, because when \(M = 3\) the number of equations is greater than the number of unknowns for our reconstruction model.
Design of the numerical validation
In the numerical validation, SC and DLS were applied to synthetic data, and the reconstructed images were compared. The generation of synthetic data is based on a circular mesh model with a radius of 300 pixels divided into 800 elements and 441 nodes. Further, 16 electrodes are evenly placed on the edge (Fig. 11a). A current of peak amplitude 1 mA was injected into two electrodes placed polarly, and the difference between the voltages on all adjacent pairs of electrodes not involved in delivering the current was measured, for a total of 192 measurements per frequency. The ground point was fixed at the center of the mesh. After generating the boundary voltage data, two levels of Gaussian noise were added, with SNRs of 60 dB and 80 dB, respectively.
An elliptical perturbation with an average length of six pixels was set at five different locations (Fig. 1a). The disturbance centers were 0, 60, 120, 180, and 240 pixels from the center of the circle. The conductivity of the background and perturbation at different frequencies were set to the values shown in Table 2.
Each target generated two frames of measurement data, namely, background frame without target and foreground frame with target. Then, SC and DLS were applied to reconstruct the images. The reconstruction was based on a circular mesh model divided into 512 elements and 289 nodes, which differed from the generation model to avoid the inverse crime (Fig. 11b). The regularization parameter was selected via the L-curve method and the two indexes of matrix \(\varvec{S}\) were calculated simultaneously in the reconstruction process.
Design of the phantom study
A phantom study was designed to reproduce the experimental setup introduced previously in the simulation. The boundary voltage data were generated across an acrylic tank filled with biological tissues. The tank was 17 cm in diameter and 7 cm in height, and 16 electrodes are evenly placed on its edge.
The cucumber was placed in three different positions, as shown in Fig. 12a, and immersed in the saline–pomelo mixture. The three perturbations, Targets 0–2, were 0 cm, 4.25 cm, and 6.5 cm away from the center of the tank, respectively.
Data acquisition was conducted using the FMMU-EIT5 system [39] shown in Fig. 12b. This EIT data acquisition system can produce a programmable current with SNR greater than 89 dB and can also measure the voltage difference precisely with CMRR higher than 75 dB [39]. The three excitation frequencies were set at 20 kHz, 50 kHz, and 100 kHz. The peak amplitude of the injected current was set at 1 mA. The used background frame and the foreground frame were averaged over 10 frames, respectively. Images were reconstructed using the same mesh employed in numerical validation. In the following, unless otherwise specified, the regularization parameter was selected using the L-curve method.
In addition to the experiments above, which were designed to compare the performance of SC and DLS, another two preliminary simulation experiments were conducted to explore the possibility of multi-tissue imaging and the influence of spectral errors. The specific implementation plans were presented in the results section.