Impact of ISTA and FISTA iterative optimization algorithms on electrical impedance tomography image reconstruction

In EIT, the reconstruction of conductivity maps within a medium is achieved by injecting alternating current (I) through strategically placed electrodes and measuring the resulting voltage (V) at the surface [3]. This process involves employing the Finite Element Method (FEM) to model the medium, account for electrode placements, and solve the associated inverse problem using various regularization techniques [39]. A commonly adopted practical approach is the adjacent measurement strategy, in which voltage measurements are not performed at the current injection electrodes. This strategy reduces redundant data by systematically rotating the current injection and measurement points, thereby improving the data acquisition efficiency.

EIT image reconstruction involves solving both forward and inverse problems simultaneously. The complexity of the forward problem is influenced by the FEM mesh and the number of electrodes [14]. However, the inverse problem requires the construction of a numerical model that discretizes the domain into fine elements to approximate the conductivity distribution. A widely used approach to solving the inverse problem is time-difference imaging, which aims to find a stable resistivity distribution (ρ) that minimizes the difference between the measured voltages (V) and the model’s predicted voltages, V (ρ) as shown in (1). (1) minV−V(ρ)2 \min \,{\left\| {V - V(\rho )} \right\|^2} where V(ρ) represents the baseline voltage measured under homogeneous conditions or at a distinct temporal point.

Because the inverse problem in EIT is inherently nonlinear and ill-posed, estimating the impedance distribution requires careful regularization to obtain a stable and meaningful solution. Regularization incorporates prior information to constrain the solution, as shown in (2). (2) minV−V(ρ)2+αLρ2 \min \,{\left\| {V - V(\rho )} \right\|^2} + \alpha {\left\| {L\rho } \right\|^2} where α is the regularization parameter, and L is the regularization operator, which can enforce smoothness or other constraints on the solution. Regularization enhances the agreement between the estimated and measured voltages, and its effectiveness is strongly dependent on the chosen method. Iterative algorithms, such as ISTA and FISTA, were employed to solve the regularized cost function. By optimally adjusting α and iteratively refining the solution, these methods play a critical role in achieving accurate and reliable IT-image reconstruction.

ISTA and FISTA are iterative optimization algorithms designed to solve regularized inverse problems. Although these algorithms are not regularization methods, they efficiently compute solutions to problems formulated with regularization terms. ISTA and FISTA are particularly effective for L1-norm regularization, such as sparse signal reconstruction. These algorithms iteratively refine the solution, progressively reducing the error between the measured and reconstructed data, guided by the chosen regularization method [40, 41, 35].

Proximal operators and the proximal gradient method are fundamental tools in optimization, for solving problems involving non-smooth functions [42, 43]. The proximal operator for a given function f, denoted as prox_λf(v), is obtained by determining the argument x that minimizes the expression in (3). (3) f(x)+12λx−v2 f(x) + {1 \over {2\lambda }}{\left\| {x - v} \right\|^2} where λ > 0 is the parameter. A notable property of this operator is that a point x^* minimizes the function f if and only if it is a fixed point of prox_λf [44], which can be expressed as in (4). (4) x∗=proxλfx∗ {x^ * } = {\rm{pro}}{{\rm{x}}_{{\rm{\lambda }}f}}\left( {{x^ * }} \right)

This result is particularly useful when deriving the proximal operator for the L1-norm, f (x) = ∥x∥₁. The minimum function can be determined as shown in (5). (5) ϕ(x,v)=x1+12λx−v2 \phi (x,v) = {\left\| x \right\|_1} + {1 \over {2\lambda }}{\left\| {x - v} \right\|^2}

Differentiating this function element-wise for x and setting the result to zero yields the result shown in (6). (6) ∂ϕ∂xi=sign(xi)+λ(xi−vi)=0 {{\partial \phi } \over {\partial {x_i}}} = {\rm{sign}}({x_i}) + \lambda ({x_i} - {v_i}) = 0

By solving for x_i, we obtain the result as shown in (7). (7) xi=vi−λif vi<−λ0if −λ≤vi≤λvi+λif vi>λ {x_i} = \left\{ {\matrix{{{v_i} - \lambda } \hfill & {{\rm{if}}\,{v_i} < - \lambda } \hfill \cr 0 \hfill & {{\rm{if}} - \lambda \le {v_i} \le \lambda } \hfill \cr {{v_i} + \lambda } \hfill & {{\rm{if}}\,{v_i} > \lambda } \hfill \cr } } \right.

The soft-thresholding operator, denoted by S_λ(u), can be more succinctly expressed as shown in (8). (8) Sλ(u)=sign(u)maxu−λ,0 {S_\lambda }(u) = {\rm sign}(u)\max (\left| u \right|- \lambda ,0)

A mapping is considered Lipschitz continuous if a constant L exists, satisfying the condition [45] as shown in (9). (9) f(x)−f(y)2≤Lx−y \left\| {f(x) - f(y)} \right\|_2 \le L\left\| {x - y} \right\|

If the Lipschitz constant L is less than 1 (L < 1), the mapping is termed as contractive. Conceptually, a contractive map reduces the distance between points upon application, potentially leading to convergence to a fixed point through an iterative application. The Lipschitz constant L is particularly significant in discussions related to iterative soft-thresholding algorithms (ISTA), as it allows the verification of a function’s contractivity by adjusting the scale of the function with 1/L.

Fixed-point iterations (FPIs) are a fundamental technique where successive iterations are computed according to (10). (10) xk+1 :=g(xk) {x_{k + 1}}\,:=\,g({x_k}) where g : ℝⁿ → ℝⁿ is a given function. The contractive mapping theorem ensures convergence to the unique fixed point of g provided g is a contractive mapping (discussed in detail in the following section).

To optimize differentiable functions f, the structure of the FPIs can be leveraged by expressing g as shown in (11). (11) g(x)=x−Λ∇f(x) g(x) = x - \Lambda \nabla f(x) where ∇f = 0 at the extrema of f, rendering them as the fixed points of g. If Λ∇f satisfies the contractive property, FPIs can lead to the minima or maxima of f.

The Newton method exemplifies this by substituting Λ with the Jacobian of f at x_k. Similarly, the Gradient Descent algorithm is employed, as shown in (12). (12) Λ=η \Lambda = \eta where η is a predefined learning rate set by the user before the commencement, which requires careful selection to avoid divergence.

In signal and image processing, the ISTA has proven to be an effective tool for addressing sparse optimization problems. ISTA iteratively applies soft-thresholding steps to optimize the loss function and is commonly used for image reconstruction and signal compression. A specific loss function must be defined to address an optimization problem with ISTA, and appropriate tuning parameters must be set. The general formula for ISTA can be expressed as shown in (13). (13) minimize Ax−y22+λx1 {\rm{minimize}}\,\left\| {Ax - y} \right\|_2^2 + \lambda {\left\| x \right\|_1} where A ∈ ℝ^n×m, m ≫ n (i.e., a fat or over-complete dictionary), y ∈ ℝⁿ is the signal, and x ∈ ℝ^m is the sparse code to be determined, respectively. The proximal gradient method can be used to obtain the update step to label the quadratic term as f and the relaxed sparsity term as g shown in (14). (14) xk+1:=Sλη(xk−ηAT(Axk−y)) {x^{k + 1}}: = {S_{\lambda \eta}}({x^k} - \eta {A^T}(A{x^k} - y)) where S_λ is the proximal operator for the L1-norm and soft-thresholding with parameter λ, as previously derived. The update step for ISTA is expressed in (15). (15) xk+1:=Sλ/Lxk−1LAT(Axk−y) {x^{k + 1}}: = {S_{\lambda /L}}\left( {{x^k} - {1 \over L}{A^T}(A{x^k} - y)} \right) where L = σ_max(A)², representing the square of the maximum singular value of A, which is essentially the largest eigenvalue of A^T A. This acts as a Lipschitz constant for ∇f and ∇g. This version of the ISTA uses a fixed maximum step size to ensure the contraction of ∇f and ∇g.

FISTA is an efficient optimization method widely used in non-smooth optimization problems, particularly in sparse coding and ℓ₁-regularized regression. This method combines the accuracy of gradient optimization algorithms with the efficiency of shrinkage thresholding techniques, while leveraging acceleration techniques for faster convergence. The updated formulas for FISTA are shown in (16) and (17). (16) tk+1=1+1+4(tk)22 {t^{k + 1}} = {{1 + \sqrt {1 + 4{{({t^k})}^2}} } \over 2} (17) yk+1=xk+1+tk−1tk+1(xk+1−xk) {y^{k + 1}} = {x^{k + 1}} + {{{t^k} - 1} \over {t^{k + 1}}}({x^{k + 1}} - {x^k}) where t^k is the acceleration parameter used to accelerate convergence. FISTA provides a flexible and robust approach for addressing the image reconstruction problem in EIT, particularly in scenarios involving large datasets and complex data models. By integrating gradient optimization algorithms with thresholding techniques, FISTA optimizes the image reconstruction process and mitigates the effects of noise and measurement errors.

The Newton-Raphson method is an iterative approach used for unconstrained minimization when solving nonlinear functions [35, 36]. The core concept of this method involves approximating the objective function with a quadratic function and finding the minimum of this approximation, which serves as an estimate of the minimum of the original objective function. The iterative formula for this method is expressed as shown in (18). (18) x(k+1)=x(k)−(ATA+rl)−1AT(Ax(k)−y) {x^{(k + 1)}} = {x^{(k)}} - {({A^T}A + rl)^{ - 1}}{A^T}(A{x^{(k)}} - y) where r represents the regularization parameter, and I is the identity matrix. In this study, the value of r was set to 0.1.

The NOSER method is a fast, static EIT algorithm derived from the Newton method. However, it operates without the need for iterations [35, 34]. Instead, it only uses the first step of the Newton method. This approach can be summarized as shown in (19). (19) E(x(0))=y−Ax(0)2,F(x(0))=∂E(x(0))∂xi, i=1,2,…,N,x=x(0)−[A(x(0))]−1×F(x(0)), \left\{ \matrix{E({x^{(0)}}) = {\left\| {y - A{x^{(0)}}} \right\|^2}, \hfill \cr F({x^{(0)}}) = {{\partial E({x^{(0)}})} \over {\partial {x_i}}},\;\;\;\,i = 1,2, \ldots ,N, \hfill \cr x = {x^{(0)}} - {\left[ {A\left( {{x^{\left( 0 \right)}}} \right)} \right]^{ - 1}} \times F({x^{(0)}}), \hfill \cr} \right.

To evaluate the efficacy of the iterative methods in enhancing the quality and accuracy of signal or image reconstruction, the following metrics were employed: Mean Absolute Error (MAE), Mean Squared Error (MSE), Signal-to-Noise Ratio (SNR), and Peak Signal-to-Noise Ratio (PSNR) [46].

The MAE metric quantifies the average magnitude of errors between the reconstructed (or predicted) values and the true values, providing a straightforward measure of reconstruction accuracy, as shown in (20). (20) MAE=1N∑i=1NYi−Y^i {\rm{MAE}} = {1 \over N}\sum\limits_{i = 1}^N {\left| {{Y_i} - {{\hat Y}_i}} \right|} where Y_i represents the actual values, Ŷ_i denotes the values predicted by the model, and N denotes the total number of observations.

MSE measures the average squared difference between the estimated and actual values, highlighting the variance of the reconstruction errors, as shown in (21). (21) MSE=1N∑i=1N(Yi−Y^i)2 {\rm{MSE}} = {1 \over N}\sum\limits_{i = 1}^N {{{({Y_i} - {{\hat Y}_i})}^2}}

The SNR evaluates the quality of the reconstructed signal by comparing its power with the power of the noise present, providing an indicator of the clarity of the signal as shown in (22). (22) SNR=10⋅log10σsignal2σnoise2 {\rm{SNR}} = 10 \cdot {\log _{10}}\left( {{{\sigma _{signal}^2} \over {\sigma _{noise}^2}}} \right) where σsignal2 \sigma _{signal}^2 and σnoise2 \sigma _{noise}^2 are the variances of the signal and noise, respectively.

The PSNR is relevant in image processing and assesses the quality of a reconstructed image by measuring the ratio of the maximum possible pixel value to the mean squared error between the reconstructed and original images, as shown in (23). (23) PSNR=10⋅log10MAXI2MSE {\rm{PSNR}} = 10 \cdot {\log _{10}}\left( {{{MAX_I^2} \over {{\rm{MSE}}}}} \right) where MAX_I represents the maximum possible intensity of the image pixel.

Assessing the quality and fidelity of reconstructed images is paramount for accurate diagnosis and imaging. Various metrics have been used to gauge different aspects of the reconstruction process, shedding light on factors such as image accuracy, resolution, and shape fidelity [47].

The Amplitude Response (AR) is a metric that quantifies the proportion between the amplitudes of pixels within an image’s target area and those in its reconstructed counterpart. It is relevant when evaluating the reconstruction of a spherical target occupying a volume V_target located in the electrode plane and having a conductivity of σ_t within a medium with a uniform reference conductivity of σ_baseline. The formula for computing the AR is shown in (24). (24) AR=Sk[x˜]kVtargetΔσσbaseline AR = {{{S_k}[\tilde x]k} \over {V_{\rm{target}}{{\Delta \sigma } \over {\sigma_{\rm{baseline}}}}}} where ∆σ represents the difference between the target conductivity σ_target and reference conductivity σ_baseline, indicating the degree of deviation in the conductivity within the target from that of its surrounding medium.

The position error (PE) gauges the accuracy of aligning the target location within a reconstructed image. This metric contrasts the actual target coordinates, r_target, with the centroid of the reconstructed target area, r_recon, formulated as shown in (25). (25) PE=rtarget−rrecon PE = {r_{target}} - {r_{recon}}

A positive SD value suggests a tendency for the reconstructed targets to shift towards the center of the imaging medium.

The resolution (RES) evaluates the dimensions of the delineated targets relative to the overall imaging field and serves as an analog to the breadth of the point spread function (PSF), as shown in (26). (26) RES=rreconAreaAtotal RES = {{{r_{reconArea}}} \over {{A_{total}}}} where r_reconArea represents the aggregate pixel count within the identified target zone, x̂_q, and A_total is the pixel count spanning the entire area subject to the reconstruction.

Shape distortion (SD) is a reconstruction algorithm commonly yield circular representations of centrally positioned targets. However, deviations from the expected circular shape often occur, particularly for targets near the imaging medium’s boundary. The SD metric quantifies the proportion of the reconstructed image that deviates from a circular shape with equal area as shown in (27). (27) SD=∑k∉C[x^q]k∑k[x^q]k SD = {{\sum\nolimits_{k \notin C} {{{[{{\hat x}_q}]}_k}} } \over {\sum\nolimits_k {{{[{{\hat x}_q}]}_k}} }} where C denotes a circle centered at the centroid of x̂_q with an area equivalent to that of the target.

The Ringing (RNG) metric evaluates whether the reconstructed images exhibit areas of contrasting polarity surrounding the primary reconstructed target area. RNG quantifies the ratio of image amplitudes of opposite polarity outside the circle C to those within C as shown in (28). (28) RNG=∑i∉CX&[x^]i<0[x^]i∑i∈C[x^]i RNG = {{\sum\nolimits_{i\not \in C} {{X_\& }{{[\hat x]}_i} < 0{{[\hat x]}_i}} } \over {\sum\nolimits_{i \in C} {{{[\hat x]}_i}} }}

In this study, simulation data were generated to evaluate the initial performance of the iterative algorithms. Data were generated from three sine waves with frequencies of 10, 50, and 100 Hz, amplitude of 1 V, number of samples N_sample of 200, and added noise of 12 dB with subsampling locations of 0.2, as shown in Fig. 1.

Figure 1:

Simulate the wave signals from three frequencies of 10, 50, and 100 kHz with (A) the sum of the signal and (B) the discrete point in the signal with a threshold of 0.2.

The results were compared between the two iterative methods, ISTA and FISTA, in terms of their performance in signal reconstruction for EIT applications. The results revealed that the FISTA method (R² = 0.9978) improved ISTA (R² = 0.9925), particularly by reducing noise and measurement errors, as shown in Fig. 2.

Figure 2:

Reconstructed signal from ISTA and FISTA iterative method.

The further analysis considered additional evaluation metrics such as the MSE, SNR, and PSNR. The results revealed that FISTA exhibited superior convergence and stability compared to ISTA, particularly between iterations of 80 and 100. Furthermore, the stable values of MSE, SNR, and PSNR for FISTA were lower than those for ISTA, with values of 0.0015, 30, and 40, respectively, compared to 0.002, 25, and 30, as shown in Fig. 3. This highlights the enhanced signal reconstruction capabilities of FISTA over ISTA. The adaptability and efficiency of FISTA stem from the integration of thresholding and gradient optimization techniques. Through this amalgamation, FISTA improves signal reconstruction and alleviates the impact of noise and measurement errors in the data.

Figure 3:

Evaluation of ISTA and FISTA performance based on MAE, MSE, SNR, and PSNR metrics over iterations.

In the application of EIT, simulation data were performed using the EIDORS library and MATLAB software, with a phantom with 16 electrodes, and a circular object assumed to be the conductor was placed into the phantom with the FEM model, as shown in Figure 4 [48, 49].

Figure 4:

16-electrode model simulates an abnormal object with a radius of 0.3 in the phantom mesh (left) and intensity conductivity (right).

The influence of regularization on the reconstruction of conductivity images is commonly addressed as an inverse problem. The results of the experiments with simulated EIT data are shown in Figure 5, using two iterative methods from previous studies—Newton-Raphson and NOSER—alongside the two methods investigated in this study, ISTA, and FISTA.

Figure 5:

Reconstructed EIT image (top) and conductivity intensity (bottom) using Newton-Raphson, NOSER, ISTA, and FISTA regularization methods, respectively

The results demonstrated that all methods successfully identified the presence of objects within the phantom. However, the FISTA regularization approach shows the highest convergence and highest contrast in conductivity between the object and background, followed by ISTA, Newton-Raphson, and NOSER. This comparative analysis highlights the significance of selecting an appropriate iterative method for enhancing the image quality in EIT. The performance of FISTA can be attributed to its advanced algorithmic structure that efficiently balances data fidelity and regularization, leading to more accurate and visually discernible reconstructions. Furthermore, the ability of the FISTA method to rapidly converge while maintaining a high contrast underscores its potential in applications requiring precise imaging, such as medical diagnostics and industrial process monitoring.

Figure 6 shows the normalized conductivity profiles according to the vertical and horizontal axes that intersect the center of the object in Figure 5 for the four iterative methods. ISTA and FISTA show high contrasts and sensitivities at the interfaces between the object and the background. Furthermore, the FISTA reconstructs conductivity with an accuracy of over 0.8. Similarly, ISTA achieves an accuracy of more than 0.6, and Newton-Raphson and NOSER achieve an accuracy of over 0.4.

Figure 6:

Normalized conductivity profile through the center of the object vertically (top) and horizontally (bottom).

In the experiment with the 3D EIT model, the 32 electrodes were divided into two rings. A spherical object with a radius of 0.3 is established with conductivity twice that of the background. The hyperparameter (α) was set to 1e–3 for regularization. The adjacent current injection strategy was used with a current of 2.0 mA, and 512 voltage measurements. Figure 7 shows the simulated 3D model with sand slices along three axes, as shown in Fig. 7(a), and the 3D mesh observations as shown in Fig. 7(b).

Figure 7:

3D simulation model 32 electrodes (2 rings) of an irregular spherical object with radius 0.3 in phantom with (a) in three axes and (b) in three-dimensional space.

Figure 8 shows the reconstruction of the 3D image using different iterative methods from the original structure shown in Fig. 8(a). The internal structure, which was difficult to observe in Fig. 8(b) when using Newton-Raphson, and Fig. 8(c), when using NOSER, becomes more obvious with ISTA. Fig. 8(d) and FISTA Fig. 8(e), respectively.

Figure 8:

3D EIT reconstruction using iterative methods with (a) original object, (b) Newton-Raphson, (c) NOSER, (d) ISTA, and (e) FISTA with corresponding conductivity scale.

The effectiveness of the proposed technique was validated through cross-sections of the 3D structure, with the original structure shown in Fig. 9(a). The blur effect is better suppressed by the iterative methods FISTA and ISTA than by the Newton-Raphson and NOSER methods, as shown in Figs. 9(b)–(e), respectively.

Figure 9:

Cross-section of the 3D EIT image shown in Fig. 8 corresponds to the iterative methods from left-to-right original object, Newton-Raphson, NOSER, ISTA, and FISTA, respectively.

Validation experiments with a heterogeneous tissue environment. The experiment used a background environment of ground meat filled with phantom material. The ground pork meat used was commercially prepared and purchased from a supermarket, ensuring compliance with ethical standards. Two objects made of acrylic resin were placed, with shapes that mimic the shape of the left and right lungs. A total of 32 electrodes were used with the adjacent current injection strategy (I = 5.0 mA) at a frequency of 50 kHz, and 1024 voltage values were acquired. The 2D FEM model is set up with a high vertex density and is two-dimensional. The tuning hyperparameter was set to 1e–3 and the time-difference rendering method was used. The impedance-signal acquisition device used in this study is shown in Fig. 10. Fig. 11 shows an experimental phantom without an object, as shown in Fig. 11(a) and with the object, as shown in Fig. 11(b):

Figure 10:

Process of signal acquisition and image reconstruction of the EIT device.

Figure 11:

Phantom model with ground pork meat and acrylic resin objects simulating the lung shape with (a) without objects and (b) with objects.

Figure 12 shows a reconstructed EIT image of the experimental phantom. The blurred and unrecognizable edges of the lung simulation object when using Newton-Raphson in Fig. 12(a), and NOSER in Fig. 12(b) are improved by using ISTA in Fig. 12(c), and FISTA in Fig. 12(d). In addition, FISTA exhibits high convergence. Therefore, the shape of the reconstructed object closely approximated its actual shape.

Figure 12:

EIT image reconstructed from the lung simulation phantom using iterative methods with (a) Newton-Raphson, (b) NOSER, (c) ISTA, and (d) FISTA.

The conducted research is not related to either human or animal use.