A Point-Matching Method of Moment with Sparse Bayesian Learning Applied and Evaluated in Dynamic Lung Electrical Impedance Tomography

Dynamic lung imaging is a major application of Electrical Impedance Tomography (EIT) due to EIT’s exceptional temporal resolution, low cost and absence of radiation. EIT however lacks in spatial resolution and the image reconstruction is very sensitive to mismatches between the actual object’s and the reconstruction domain’s geometries, as well as to the signal noise. The non-linear nature of the reconstruction problem may also be a concern, since the lungs’ significant conductivity changes due to inhalation and exhalation. In this paper, a recently introduced method of moment is combined with a sparse Bayesian learning approach to address the non-linearity issue, provide robustness to the reconstruction problem and reduce image artefacts. To evaluate the proposed methodology, we construct three CT-based time-variant 3D thoracic structures including the basic thoracic tissues and considering 5 different breath states from end-expiration to end-inspiration. The Graz consensus reconstruction algorithm for EIT (GREIT), the correlation coefficient (CC), the root mean square error (RMSE) and the full-reference (FR) metrics are applied for the image quality assessment. Qualitative and quantitative comparison with traditional and more advanced reconstruction techniques reveals that the proposed method shows improved performance in the majority of cases and metrics. Finally, the approach is applied to single-breath online in-vivo data to qualitatively verify its applicability.


Introduction
Electrical impedance tomography (EIT) is a medical imaging technique which reveals the conductivity or admittance distribution of a subject under test (SUT) [1]. In EIT, an alternating, low-amplitude current usually up to 1MHz is induced into a cluster of electrodes, while the measured electrode potentials are used as raw data for the image reconstruction. Unlike other medical imaging modalities, EIT is characterized by the absence of ionizing radiation, its low cost and its notable temporal resolution. This makes EIT a useful tool for real-time lung function monitoring. Many studies have shown the importance of EIT for revealing vital signs related to ventilation properties, such as tidal volume (TV) [2], or pathological situations, such as acute respiratory distress syndrome (ARDS) [3].
Despite its potential advantages, EIT is lacking in spatial resolution, something that still keeps its application in medical equipment limited. In addition, EIT images often present artefacts that, in some cases, may degrade their clinical value and diagnostic efficacy. Such artefacts are related to the highly ill-posed and ill-conditioned nature of the EIT inverse reconstruction problem. This means that the image quality presents with high sensitivity to voltage signal noise. In real-time EIT imaging, the signal-to-noise ratio (SNR) is often limited, because higher frequency currents need to be injected in order to achieve a high temporal resolution [4,5]. Indeed, in practice, lower frame rate EIT hardware systems [6,7] usually present better voltage SNR levels than higher frame rate systems [8,9].
Furthermore, EIT is susceptible to modeling errors. Particularly in dynamic thoracic imaging, the lack of a homogeneous background, the patient's unknown chest boundary shape, as well as its changes over time due to the patient's breathing cycles and the high chance of the electrodes' displacement during signal acquisition also introduce significant modeling errors [10][11][12]. Despite the fact that a percentage of modeling errors can be compensated when time-difference EIT is applied, their impact on image quality may be still noticeable. Moreover, EIT is a highly non-linear problem, which means that highly inhomogeneous admittance inclusions cannot be accurately estimated with simple linear methods [12]. The variations in the lungs' admittance due to the continuous change in air volume is a typical case [13].
Over the years, many approaches have been proposed for EIT image reconstruction. The most simple approaches assume small conductivity or admittance inclusions and linearize the problem around a predefined homogeneous value. Then, the problem is treated directly using either truncated singular value decomposition (TSVD) or standard Tikhonov regularization (STR), where an identity matrix prior is considered. Generalized Tikhonov regularization (GTR) schemes [14] such as the Laplace prior, the NOSER prior [15] and the high-pass filter Gaussian prior [16] can be also applied, improving the image reconstruction performance compared to TSVD and STR. To compensate boundary movement effects, [17,18] introduced electrode movement priors, combining them with the previous schemes. Another single-step approach which makes use of a figure-of-merit (FoM) framework to optimize both parameters and performance is the Graz consensus reconstruction algorithm for EIT (GREIT) [19]. GREIT uses dual-mesh schemes, i.e., a coarse 2D domain for the admittance reconstruction and a fine 3D domain, extruded from the coarse 2D one, for accurate forward calculations. This greatly improves accuracy, but at the cost of additional time and complexity.
Despite the benefits of single-step approaches in dynamic EIT imaging speed, EIT is a non-linear problem. Therefore, most linear approaches cannot accurately capture significant conductivity changes. In such cases, iterative approaches have been developed to deal with the non-linearity. A common iterative framework is based on the Gauss-Newton (GN) algorithm, which makes use of the GTR schemes mentioned above (L 2 -norm regularization) [20]. Another popular iterative approach is the total variation (TV) approach, where L 1 -norm regularization priors are used [21][22][23]. A hybrid non-linear difference EIT imaging approach was proposed in [12,24] in an effort to effectively deal with both the models' mismatches and the problem's non-linearity. Another L 1 -norm approach uses the Bregman distance scheme, showing improved performance in lung imaging compared to the traditional TV scheme [25].
Sparse Bayesian learning (SBL) was firstly proposed as a mathematical formulation in [26][27][28] but was only recently applied in EIT [29][30][31][32]. Instead of the traditional regularization schemes, it treats the inverse problem as a log-likelihood optimization procedure, assuming a sparse conductivity distribution. SBL approaches show robustness to signal noise, while the non-trivial hyperparameter selection needed in regularization techniques is avoided. Although some SBL methods, such as structure-aware SBL (SA-SBL) [30] and time-sequence learning SBL [33] have been proposed for EIT, their evaluation is limited to simple circular and cylindrical structures. Hence, SBL has not been applied in dynamic thoracic imaging, where the structures have more complex geometries, usually unknown, and are highly non-homogeneous.
The point-matching method of moment (PM-MoM) for EIT was also recently proposed in [34]. It uses a global integral equation approach with Green's functions. The logarithm of conductivity is expressed as a linear combination of modified radial basis functions (RBFs). Contrary to the traditional finite element (F.E.) approach, which treats the problem as a weak form that does not hold for significant conductivity changes, the PM-MoM is formulated globally, decreasing the problem's non-linearity. The PM-MoM has therefore shown to converge faster than the traditional F.E. approaches both in L 1 and L 2 -norm inverse problem schemes. Despite the fact that PM-MoM has been tested on circular and cylindrical structures, it has not been quantitatively evaluated in dynamic thoracic imaging.
Motivated by the benefits of the PM-MoM and the SBL optimization scheme, in this work we introduce an approach that combines these two methods. In particular, the proposed PM-MoM SBL approach undertakes the image reconstruction problem's non-linearity, offering robustness to noise and reduced susceptibility in modeling errors. To evaluate the PM-MoM SBL approach, we apply it to 3D F.E. thoracic structures (cases) based on 3 male subjects' CT images, available online. For each case, five sub-structures are built, considering five corresponding breath-cycle states from expiration to the end-inspiration. Each structure includes the lungs, heart, vertebrae, muscle and skin tissues to avoid the assumption of a uniform background [11]. It is noted that most previous studies for EIT approaches in dynamic thoracic imaging consider only the inspiration and expiration ends and only the lung tissues in their models for quantitative evaluation. The proposed approach is compared with traditional (GN, TV, difference of absolute images) and more advanced (prior movement, hybrid non-linear imaging) F.E.-based regularization approahces, as well as the regularized MoM approach, showing increased noise robustness and improved performance both qualitatively and quantitatively. Finally, the proposed method is tested in in vivo human breath data which is available online, verifying its proper applicability.
The rest of this paper is organized as follows. In Section 2, the EIT problem's principle, as well as state-of-the-art regularization-based methods, are outlined. In Section 3, the proposed PM-MoM SBL approach is presented, while in Section 4, the 3D thoracic structures, the evaluation FoM, the method adopted to extract the reference images and the in vivo data are described. In Section 5, the image reconstruction results, as well as the quantitative results, are demonstrated and discussed. Finally, Section 6 concludes this work.

Background
In this section, a brief review of the EIT mathematical formulation and the state-ofthe-art inverse reconstruction approaches used for the comparisons is performed.

EIT Principle
Assume a N-electrode EIT setup and a d-dimensional domain (d ∈ {2, 3}) Ω, where the current is injected. The problem can be described according to the following Laplace equation: that implies the following boundary conditions, according to the complete electrode model (CEM) [35]: where r ∈ R d is the position vector, σ(r) is the conductivity, u(r) is the potential, n is the normal outward-pointing vector, e l is the electrodes' positions set, z l is the electrodes' contact impedances, U l is the lth electrode voltage and I l is the current injected on the lth electrode.

Time-Difference EIT
In time-difference (dynamic) EIT imaging, we assume two consecutive states. The corresponding computed boundary voltages' vectors can be written as follows: and where m is the total voltage measurements acquired for each current injection electrode pair, according to the selected measurement pattern [36]. Accordingly, we assume that two voltage data measurement frames V (1) ∈ R Nm×1 and V (2) ∈ R Nm×1 are acquired from the EIT system. We then set the differential voltage frames and Considering Gaussian noise e n ∈ R Nm×1 between the simulated boundary voltages and the measurements, we have δV = δU + e n .
Furthermore, we assume σ (1) (r) and σ (2) (r) as the corresponding conductivity distributions, setting the difference If the finite element method (F.E.M.) discretization scheme is applied in Ω, we assume a number of L elements and that each element i presents a constant conductivity σ i .
Hence, we can write the following conductivity vectors: The general approach is to formulate the inverse problem as a weighted-least squares (WLS) minimization problem between δV and δU, adding a regularization term P(δσ) to stabilize the problem's ill-conditioned nature, such that where W ∈ R Nm×Nm is a diagonal, noise covariance matrix and λ is the regularization hyperparameter. The problem is to find the optimal δσ that minimizes F(δσ) in (14). From hereon, we assume that all the measurement channels have the same noise; hence, W = I Nm .

Single-Step Linear Reconstruction
As mentioned in the introduction, EIT image reconstruction is a non-linear problem, i.e., the relation between U and σ(r) is non-linear. However, assuming relatively small conductivity changes, and using Taylor approximation around a linearization point σ o , we can write where J ∈ R Nm×L is the Jacobian matrix around σ o . In the simplest case (dimensionless electrodes), J is computed according to the following formula [37]: where Ω i denotes the i-th element's domain, d denotes the current injection differential channel and m denotes the voltage measurement differential channel. The minimization function is then written according to the following form: In this case, smooth L 2 priors, such as the standard Tikhonov, the Laplace, the NOSER or the Gaussian priors, are commonly utilized [14][15][16]. Hence, we can write where Q ∈ R L×L is the prior matrix. The linearized problem has the following closed-form solution: This is the most commonly used linear EIT scheme.

Regularized Reconstruction Approaches
In this section, we perform a brief review of the regularization-based linear and non-linear approaches used for the comparisons.
(I) Non-Linear Gauss-Newton (GN): The traditional non-linear GN approach makes an iterative estimation of conductivity change distribution to optimize An initial solution is taken from (20). Then, J, as well as δσ, are re-estimated in each iteration until convergence.
(II) Total Variation (TV): This non-linear, iterative approach assumes that intense conductivity changes occur between neighbouring elements of Ω by applying L 1 -norm priors. The functional to be minimized is written as follows: where N ed is the total number of edges between the domain's elements, and L ∈ R N ed ×L is a sparse matrix which shows the relation between the elements and their edges. L i refers to the ith row of the matrix L, and β > 0 is a parameter that prevents the non-differentiability of the regularization term [22]. In terms of this paper, the primal-dual interior point (PD-IPM) TV method is adopted [21][22][23].
(III) Movement Prior: This approach, which was firstly proposed by [17] and furtherly developed in [18], performs linear difference-EIT reconstruction while considering the electrodes' movement effect. The electrode movement δx ∈ R dN×1 is also estimated along with the conductivity change. To this end, Tikhonov and Laplace priors have been proposed for the simultaneous estimation of δσ and δx, properly modifying the Jacobian matrix J and the prior matrix Q [17,18]. Apart from λ, a µ > 0 regularization hyperparameter for the electrode movement prior is needed. In terms of this paper, both δσ and δx estimations are performed using the Laplace prior, in the way described in [18]. This approach has proven to reduce the artefacts caused by the electrodes' movement and boundary changes and has been exclusively developed for applications in dynamic lung imaging.
(IV) Difference of Absolute Images: In this approach, absolute, instead of difference, EIT reconstruction is applied particularly for each measurement frame [38]. The problem's objective functions are the following: and Defining σ (1) * ∈ R L and σ (2) * ∈ R L as the obtained solutions from (23) and (24), respectively, the final estimated conductivity change is simply obtained by In this paper's reconstructions, the minimization of F 1 and F 2 is performed by using the absolute GN non-linear approach with a Laplace prior.
(V) Multiple Priors (Non-Linear Difference Imaging-N.L.D.): This approach, proposed in [12,24], concatenates the voltage data measurement frames as follows: and the computed boundary voltages as follows: while it is assumed that conductivity changes occur in a particular region of interest (Ω ROI ⊆ Ω), which is discretized in L ROI ≤ L elements, such that where M is an operator that maps δσ ROI with the domain's elements. The following optimization problem is defined: where λ 1 and λ 2 are regularization hyperparameters, N ed,ROI is the total number of edges between the ROI elements and L i,ROI ∈ R N ed,ROI ×L ROI shows the relation between the ROI's elements and their edges. A smooth L 2 -norm prior (prior matrix Q 1 ) is used to estimate σ (1) , while a L 1 -norm prior is used to estimate the local conductivity change δσ ROI . The solution σ * that minimizes (29) is achieved via an iterational process. A linesearch to perform the updates is essential due to the problem's high non-linearity. It is worthwhile to mention that apart from these methods, many other approaches have been developed for difference-EIT imaging. For example, the TSVD, a traditional linear approach, gives similar results to (20) by performing thresholding SVD instead of regularization. Finally, a very well-known single-step approach is the D-Bar, which uses a non-linear scattering transform and low-pass filtering [39][40][41]. A comparison between D-Bar and regularized methods can be found in [42].

Method-of-Moment with Sparse Bayesian Learning (Pm-Mom SBL)
This section presents the proposed method as a combination of the PM-MoM system matrix formulation and an SBL approach for the occurring inverse problem.
The PM-MoM formulates a global integral equation (holding in the whole Ω) instead of using the weak form of (1). It expresses the voltages and fields as Green's functions and their gradients, respectively. The conductivity is non-linearized in the integral equation, and, unlike the conventional F.E.M. formulation, is not assumed to be a piecewise constant at each element. Instead, its logarithm is expressed as a summary of radial basis functions (RBFs) [34].
The governing integral equation takes the following form: where G(r, r ) denotes the Green's function set between an observation point r and a source point r , u o denotes the domain's voltage distribution when the conductivity is constant (homogeneous with a value σ o ) and r + , r − denote the electrode coordinates from where the current is sourced and sinked, respectively. It is also found that We then express the logarithm of conductivity as follows: where r j is the jth pixel's center point. The RBF θ is selected according to the following generalized form: where p 1 > 0 and p 2 > 0 are integers (p 1 even), and D is a parameter which adjusts the RBF's width. By discretizing (30), replacing the conductivity logarithm with (32) and taking the electrodes voltages' differences, according to the measurement pattern adopted and the method described in [34], we are led to a linear system of equations where M o ∈ R Nm×L is the system matrix, c = [c j ] L j=1 ∈ R L×1 is the weighting coefficients' vector and δU ∈ R Nm×1 denotes the numerically expected electrode potentials' differences. Adopting the measurement model described by (9), we treat the inverse problem as a minimization of the following WLS objective function: The minimization of (35) can be performed with the traditional L 2 or L 1 -norm regularized approaches [34]. Unlike J, the matrix M o occurs directly from the discretization of (30) without the assumption of ln(σ(r )) σ(r ) − 1 near σ o . Hence, the expression of the boundary voltages as a function of conductivity is more accurate, leading to a faster convergence of the inverse solution. We note that the Green's function G and its gradient can be separately precomputed either analytically for canonical geometries or by using the F.E.M. or the finite difference method (F.D.M.) (at the same discretization mesh as MoM) to solve the Laplace equation for the potential and the field for non-canonical geometries.
In this particular work, we make use of an SBL formulation [30] to minimize (35). To this point, we interpret the objective function in a Bayes log-likelihood context where Θ is a set of hyperparameters. Considering that c is a superposition of some clusters overlapping each other with an equal size h, and g = L − h + 1 is the total number of clusters, the following factorization is performed: where (9) is approximated as We also define Φ = M o Ψ ∈ R Nm×gh . Furthermore, we assume that the weight vector x ∈ R gh×1 obeys the following Gaussian distribution: with zero mean value and Σ 0 ∈ R gh×gh covariance matrix.
} and adopting the expectationminimization (EM) method, as in [30], we get an a posteriori estimation of the weight vector's x mean values vector µ x ∈ R gh×1 and covariance matrix Σ x ∈ R gh×gh . Hence, we get a maximum a posteriori (MAP) estimation of x. For clarity, we summarize the SBL process in Algorithm 1.
The updating rule for γ i is based on the majoration-minimization method [26,43]. In addition, instead of the linearized Jacobian matrix, we use the PM-MoM M o system matrix as an input in order to limit the problem's non-linearity effect and avoid the necessity of recalculating J.
The SBL method shows increased robustness to noise and modeling errors, while the reconstruction artefacts are minimized. Furthermore, unlike the traditional regularized schemes, the choice of the hyperparameter h (number of clusters) slightly affects the reconstruction quality [30]. Hence, the cumbersome and non-trivial process of hyperparameter selection is avoided. Moreover, the regular "moment" grids used in PM-MoM are suitable for performing sparse-based reconstructions. However, despite recent developments, the SBL methods are overall characterized by high complexity. For instance, the SBL approach applied presents a complexity of O(N 2 m 2 gh) per iteration. Nevertheless, avoiding the search process for optimal hyperparameters partially reduces the increased complexity effects. A simplified flow chart of the whole PM-MoM SBL, as well as an example of a domain's clustering, are depicted in Figure 1. The SBL approach adopted in this particular paper is in the form used in [30]. However, some modified SBL approaches using approximate message passing (AMP) to accelerate the E-step of the algorithm have been researched for 3D imaging [31], frequency-difference EIT [32] and multiple measurement vector (MMV) time-sequence measurements [33]. Such techniques can also be appropriately combined with the PM-MoM in a similar manner.

Evaluation Methods
In this section, the thoracic structures' extraction and the corresponding tissues' electrical properties are demonstrated. In addition, the evaluation metrics, including GREIT FoMs with minor modifications, the Pearson CC, the RMSE and the recently proposed FR are briefly explained. Finally, an in vivo online available dataset demonstrating a subject's full-breath cycle is briefly described.

Thoracic Structures
To examine and compare the previously described algorithms' performance in dynamic imaging, we have created 3 3D fine F.E. thoracic structures based on 3 CT-images of 3 corresponding different healthy adult male subjects. The CT images were taken between the third and the fourth intercostal levels and are included in a large medical database which is available online in [44]. The 3D models have been created in MATLAB using the EIDORS and the NETGEN software [45,46] and include the following tissues: left lung, right lung, heart, vertebra and skin, while muscle is assumed to be the background.
For each structure, 5 breath-cycle states have been considered from end-expiration (deflated) to end-inspiration (inflated). Hence, a total number of 15 F.E. models have been created, demonstrating 3 subjects' cases in 5 breath-cycle states. Each state presents chest boundary changes (a total change of 5-8% of the chest's width) and lung shape changes (total expansion 10-15% of the lungs' width at the inflated state) [12]. The lungs' admittance changes between the states have also been considered.
The tissues' admittance values are loaded from an open-source database, demonstrated in [47][48][49]. All admittance values depend on the selected current signal frequency f in which the EIT measurements are performed, while the lungs' admittances also depend on the breathing state. For this particular work, we assume f = 100 kHz, which is in the range of the current frequencies used for dynamic lung imaging. This is actually a common frequency choice for such applications [7,50,51]. Nevertheless, higher frequencies, such as those applied in high-performance modern EIT systems [7,9], can be also considered. Furthermore, to take into account each tissue's inhomogeneity, the following standard deviations (std) of the admittance values assigned to each tissues' elements have been taken into account: 1% for the skin, 2% for the heart and the muscle background and 3% for both lungs. These values fall within the range depicted in the mentioned database.
The conductivity and permittivity values assigned to each tissue at 100 kHz, as well as their std are shown in Table 1. For the lungs, the deflated and inflated states' values were taken from [47][48][49]. To find the intermediate states' values, we firstly assumed that the lungs' volumes increase linearly over time during the inhalation process [52]. A relative (arbitrary unit-A.U.) volume has been defined as follows: where P is the total number of states from end-expiration to the end-inspiration. In our case, P = 5. In actuality, the lungs' volume change is more complex and heavily depends on each particular breath. The main changes in the lungs' admittance occur due to the air-flow. The lungs' conductivity as a function of their volume can be expressed by [52,53] where s b , w and s i are lung morphological parameters described in [53], with their values selected to be s b = 0.5, w = 1.5 and s i = 2. In addition, K 1 and K 2 are coefficients used to scale the lungs' conductivity between the known values and end-inspiration and end-expiration. The permittivity of the lungs is correspondingly defined as [53] where e rb and e rm are also lung morphological parameters described in [53], with their values selected at e rb = 10 4 and 10, respectively. Furthermore, L 1 and L 2 are scaling coefficients. We note that blood-cycle related changes have not been taken into consideration, since the HR frequency is 3-6 times higher than the breath frequency. In each case (1-3), the boundary extracted from the corresponding CT image is used as a cross-section to create the 3D F.E. structure. For each structure, a height of h = 1 A.U. has been considered, while the x-axis limits have been normalized between −1 and 1 A.U. A number of N = 16 circular electrodes of radius R el = 0.05 A.U. have been placed at the z = 1/2 A.U. level. In addition, an electrode position error has been added: 5% height std and 3% angle std, since this is a more realistic case.
The thoracic structures are demonstrated in Figures 2-4. Their boundary and lungs' shape changes are demonstrated at the cross-section level in Figure 5. Finally, the numbers of each model's tetrahedral elements and nodes are shown in Table 2.
To simulate the measurements, the adjacent (skip-0) current and voltage measurement pattern was considered [36,54], while a Gaussian noise of −50dB SNR was added to the extracted raw signals. The EIDORS library tool in MATLAB was used to perform the simulations [45].

Reconstruction Domain
When the EIT image reconstruction is performed using simulated models, we need to avoid inverse crime. This occurs when the simulated model's and the reconstruction domain's mesh or boundary is equal [55]. Instead, the reconstruction needs to be performed on a significantly different mesh, usually coarser than the simulated model's one.
In this work, all the image reconstructions are performed on a 2D coarse thoracicshaped domain, called Ω, which presents a different boundary than any of the original model's boundaries. For the FEM-based reconstruction approaches, the domain contains L = 1024 triangular elements and n e = 545 nodes. The shunt electrode model has been assumed to simulate the electrodes' effects [35,56]. Furthermore, for the PM-MoM reconstructions, a L = 1060 uniform pixel grid has been used, considering the electrodes as points [57]. The reconstruction domain for the two discretizations is shown in Figure 6.

Reference Image Extraction
In order to perform a quantitative evaluation of EIT imaging, a corresponding "ground truth" reference image has to be defined on the reconstruction domain Ω. However, the simulated models have completely different shapes and discretization meshes compared to Ω in both the standard FEM and MoM reconstruction cases. In order to "match" the simulated models with the reconstruction domain, the approach presented in [58] is adopted.
This approach considers that the simulated models' shape in all three cases is not constant, as well as that difference-EIT imaging is performed. Hence, we firstly get five absolute reference images, each one representing a particular state. Then we take the differences between the 2nd-5th images and the 1st "reference frame", resulting in 4 reference images.
Each "true boundary" domainΩ k,l , k = {1, 2, 3}, l = {1, 2, ..., 5} is extracted from the 3D models' electrodes' cross-section plane. Then, we scale Ω and eachΩ i,j in the x-axis by normalizing its limits between −1 and 1. We secondly define A i as the Ω ithelement's/pixel's area, with i = {1, 2, ..., L}. Then, the percentage of A i which is within the curves defined by the following six tissues, left lung, right lung, vertebra, heart, skin and muscle, is expressed as a weight vector, for the kth case and lth state. At this point, we define the vector which represents the mean admittances for each one of the 6 mentioned tissues at the lth state. Then, the ith element's or pixel's reference admittance is estimated as follows: The kth case, lth state (absolute) reference admittance vector is then defined as For each image frame, we get the following difference reference admittance vector: where δγ r,k,l+1 ∈ C L×1 for l = {1, 2, 3, 4}, assuming the kth case.
A simple example of this process and the F.E.M. reference images for k = 1 are demonstrated in Figure 7.

Figures of Merit
To quantitatively evaluate the EIT reconstructions, we use the following five GREIT FoM: target amplitude-TA, position error-PE, shape deformation-SD, resolution-RES and ringing-RNG [19]. These have been properly adapted to the examined thoracic cases. Furthermore, the CC, the RMSE and the FR metrics are applied.

Target Amplitude-TA
Assume δσ * ∈ R L×1 isthe conductivity difference estimated from the EIT reconstructions. The TA can be defined as the normalized summary of elements'/pixels' amplitudes in the image, TA is a FoM similar to the amplitude response-AR, which is considered to be the most important GREIT FoM [19]. Its absolute value should be relatively low and stable during the breath process. When admittance values are reconstructed, TA can be estimated by taking only the real values.

Position Error-PE
The position error-PE-shows the precision of the reconstructed inclusions' center of gravity. In our case, we define the right and the left lung as the corresponding inclusions. Then, we get two PE values: for the left lung, where r tLL is the true center of the left lung and r iLL is the reconstucted left lung's center, and PE RL = |r tRL − r iRL | (50) for the right lung, where r tRL is the true center of the right lung, and r iRL is the reconstucted right lung's center. The total PE is given by To detect an inclusion as "lung", we first filter the reconstructed image, setting all the elements'/pixels' absolute values that are below a selected threshold (−1/4 of the maximum absolute value) to zero and all non-zero values to 1. We denote x f ∈ R L×1 as the filtered image conductivity distribution, where Secondly, the left and right lung inclusions are separated in the reconstructed image with a y-axis line, as shown in Figure 8.

Shape Deformation-SD
Shape deformation-SD-denotes the percentage of the reconstructed and filtered inclusion which is not within the "true lung's" boundary. Assume LL and RL are the "true" left and right lungs' domains, respectively. If x f is the filtered reconstructed image conductivity distribution, we get for the left lung and for each element/pixel left from the y-axis with an area A i (see Figure 8). For the right lung we have for each element/pixel right from the y-axis with an area A i . The total SD is given by SD should also be low and stable.

Resolution-RES
If A LL represents the reconstructed left lung inclusion's area, A RL is the reconstructed right lung inclusion's area and A o is the Ω area, the resolution-RES-is given by We can estimate A LL and A RL from the following expressions and RES should be low and uniform [19].

Ringing-RNG
Ringing-RNG-demonstrates whether the reconstructed inclusion causes areas of opposite sign near the target inclusion. It is given by RNG is an important GREIT FoM, since in dynamic lung image reconstructions, conductive areas often appear between the lungs that are sometimes wrongly recognized as "heart" [19]. It should also be low and as stable as possible. Apart from the above GREIT parameters, we also apply the following FoM.

Pearson Correlation Coefficient-CC
The Pearson correlation coefficient-CC-is one of the most common metrics that quantify an image's quality. It indicates the similarity between a "ground truth" reference image and the reconstructed image. It is given by where δσ r = Re{δγ r } for each particular case k and state l.

Root Mean Square Error-RMSE
An additional FoM used for the image evaluation is the well-known RMSE, which is estimated according to the following formula:

Full Reference-FR
This metric was recently proposed as a universal FoM for EIT systems' evaluation on the reconstructed images [59]. It has been extensively presented and applied in phantom experimental setups. However, this is the first time that FR is applied in dynamic thoracic models.
To estimate FR, the normalization of the reference images δσ r and the reconstructed images δσ * (element/pixel) data between −1 and 1 needs to be performed. We define as EDre f ∈ R L×1 the normalized reference image data and EDtest ∈ R L×1 the normalized reconstructed image data.The global FR (GFR) is defined as follows: We also define the local FRs for the left and the right lungs, respectively, as and Both local and global FR indicate a high-quality reconstruction when they take low values.

In Vivo Data
A qualitative comparison is attempted using online available in vivo EIT data [45]. This data consists of 34 data frames of a single breath cycle captured by the 16-electrode serial-data EIT Scanner [60] using the adjacent current and voltage measurement pattern. This system performs demodulation of the input signal with an AM signal of a higher order of magnitude frequency than that of the electrode voltage signal. The injected current frequency (carrier) signal was set at 65 kHz. Since difference-EIT imaging is performed, the first frame is used as reference, resulting in 33 image reconstructions.

Results and Discussion
Image reconstructions were performed for the 3 structures presented in Section 4.1, considering each one of the 5 mentioned breathing states and resulting 4 images per structure. Reconstructions were also performed for the in vivo data described in Section 4.5, resulting in 33 EIT images. The regularization scheme-based approaches described in Section 2.4, the MoM-regularized approach, as well as the proposed MoM SBL approach described in Section 3 were used. Particularly for the multiple priors difference non-linear approach (N.L.D.), we consider that the ROI where δσ ROI occurs is equal to Ω, since the "lungs" area covers a significant part of Ω [12].
For all cases, the reconstruction hyperparameter λ value, as well as the µ, β and h parameters' values (for the movement-prior, TV and PM-MoM SBL reconstructions, respectively) were heuristically selected, as shown in Table 3. The selection of λ was performed in such a way that CC is maximized. Although some methods for the λ selection, such as the L-curve, the noise figure (NF) and the BestRes calibration methods [19,61], have been proposed, this process is beyond this work's scope. Finally, for the PM-MoM SBL, we set the maximum number of iterations κ max to 5 and the minimum tolerance min to 10 −5 . Table 3. Selection of reconstruction parameters per algorithm.

Simulation Results
The resulting image reconstructions and the FoM values demonstrating the simulated cases are depicted in Figures 9-14. Specifically, the image reconstructions that represent each one of the structures 1-3 are shown in Figures 9, 11 and 13, respectively. We define the reference image extracted according to the process described in Section 4.3 as the "true" image. Furthermore, the corresponding FoM values obtained are demonstrated in Figures 10, 12 and 14, respectively.
A visual inspection of the images which resulted from the 1st structure ( Figure 9) shows that the air-filled lungs are successfully detected from all the approaches. However, the lungs' shape and area is deformed, while "pseudo-heart" and boundary artefacts often appear. Such effects are less intense in the MoM Laplace and the MoM SBL cases. As we proceed to the full-inhalation state, the conductivity contrast increases, resulting overall in increased absolute TA, PE and sometimes RNG (Figure 10). At the same time, both local and global FR decrease, while non-significant changes occur at the RES, SD and CC. The RMSE value remains almost constant for all the approaches, except for the TV, where it decreases. Comparing the metric values for each algorithm, better results are obtained by the PM-MoM SBL approach, which shows the lowest and most uniform absolute TA, the lowest PE, RES, SD, RMSE and GFR and the highest CC, followed by the PM-MoM regularization approach.
The images obtained from the simulations of the second structure ( Figure 11), which is characterized by closer distance between the lungs, show almost all the artefacts near the boundary instead of between the lungs. The PM-MoM SBL method also shows less intense artefacts, while, along with the regularized PM-MoM and the multiple priors N.L.D. approach, achieving the best CC and lower RMSE, local and GFR values ( Figure 12). The proposed method also achieves the lowest absolute and most constant TA, RES and SD. However, the best local FR levels are extracted from the N.L.D, while the PM-MoM shows an increased RNG metric.
The third case results in Figure 13 are characterized by overestimation of the air-related conductivity change near the chest. This occurs due to the presence of lung tissue very close to the chest boundary, as the EIT measurements are sensitive to conductivity changes near the boundary [62]. A visual comparison of the images in Figure 13 indicates that the PM-MoM SBL method has the best performance, an absence of "positive conductivity change" artefacts and less lung deformation. Considering the quantitative results, the PM-MoM SBL approach achieves the lowest PE, RES, SD and RMSE ( Figure 14). Although the best TA, RNG, CC and GFR are demonstrated by the GN, N.L.D, N.L.D. and the regularized PM-MoM methods, respectively, the PM-MoM SBL shows the most constant TA, acceptable levels of RNG, a CC which is close to the best one, and the second-lowest GFR.

In Vivo Results
The in vivo EIT reconstructed images that demonstrate a subject's single breath, as described in Section 4.5, are shown in Figures 15-18. In particular, Figure 15 shows the reconstructed images using the GN and TV approaches, reviewed in Section 2.4. Figure 16 shows the reconstructed images using the movement prior linear approach, Figure 17 depicts the reconstructed images using the difference of absolute images and the multiple priors non-linear difference imaging (N.L.D.) approaches and Figure 18       A qualitative observation of the images leads to the outcome that all the approaches, except for GN, are able to detect the full-inspiration state. However, the presence of ringing is significant near the centre (between the lungs) in the movement prior, GN and N.L.D. approaches. This might lead to misleading conclusions about the presence of "heart" tissue between the lungs, as mentioned above. In fact, the heart tissue is not directly detectable in difference EIT imaging, since its conductivity does not significantly change during the breath cycle. In addition, any blood-cycle-related conductivity changes are synchronized with the heart rate HR, while the "pseudo-heart" inclusion is synchronized with the breath cycle. Meanwhile, the TV, difference of absolute images and both PM-MoM approaches demonstrate boundary "positive conductivity change" artefacts which are related to the mismatch between the patients' thoracic shape and ∂Ω. This effect is less intensive in the regularized PM-MoM, PM-MoM SBL, the difference of absolute images and the TV methods, which overall perform better than GN, N.L.D. and movement prior. However, the TV method appears to underestimate the lungs' area in relation to the total thoracic area. We also observe that the inequality between the lungs' volumes is successfully detected by most of the approaches (except for TV), but is more clear when enacting PM-MoM (regularized or SBL), difference of absolute images and N.L.D.

Discussion
In this work, a proposed EIT reconstruction method, which combines PM-MoM for the system matrix formulation and SBL for the inverse problem solution, is applied to dynamic thoracic imaging. A number of evaluation criteria is adopted, and an extensive comparison is performed with numerous state-of-the art approaches. Qualitative and quantitative studies have been performed on 3D time-variant non-homogeneous thoracic models. In vivo imaging of a patient's full-breath cycle has also been applied.
The qualitative results both in the simulation (Figures 9, 11 and 13) and the in vivo studies (Figures 15-18) reveal that, overall, the proposed PM-MoM SBL approach shows a better spatial resolution than both the traditional and some more advanced linear and non-linear EIT reconstruction methods, as well as the regularized PM-MoM method. This is confirmed quantitatively in Figures 10, 12 and 14 for the simulated cases. In all three subject-cases, the PM-MoM SBL outperforms the other approaches in most of the FoMs.
It is worthwhile to mention that, of the other approaches, the regularized PM-MoM appears to be the most efficient. In addition, the more recently proposed movement prior, difference of absolute images and multiple priors N.L.D. apaproaches appear to outperform the traditional GN and TV methods.
Considering the time needed for the reconstructions, the best performance (about 50 ms per frame) is achieved by the regularized PM-MoM when using the Laplace L 2 -norm prior with a single step. The movement prior linear approach needs about the same amount of time per reconstruction, while the non-linear GN and TV methods need significantly more time (about 4 to 6 seconds, depending on the number of iterations needed). Due to the relatively high complexity of SBL (O(N 2 m 2 gh) per iteration), the PM-MoM SBL approach needs, on average, 5.6 s per image frame reconstruction, when h = 4, N = 16, m = 13, g = 1055 and the number of iterations is 5. However, despite the time needed for PM-MoM SBL, the hyperparameter selection process, which is usually time-consuming, is avoided, contrary to the regularization approaches. Finally, the difference of absolute images as well as the multiple priors N.L.D. approaches require significantly longer times to reconstruct the images. The times mentioned above have been achieved using an AMD Ryzen 5 3600 system.
In conclusion, the PM-MoM SBL approach outperforms the regularized MoM one regarding the images' quality and spatial resolution. Additionally, there is no need for hyperparameter selection, which partially reduces the SBL process complexity effect on execution time. The PM-MoM SBL can be directly applied either for offline imaging (after collecting the measurements) or online on particular breath states where the lung conductivity change is significant. Another choice for faster online imaging is to reduce either the maximum number of iterations κ max or tolerance min (see Algorithm 1). It is worthwhile to mention that further research on optimizing the SBL approaches' complexity, as well as evolution in hardware, might improve the total time needed per image reconstruction.

Conclusions
An EIT reconstruction approach based on a method of moment which expresses the conductivity logarithm with radial basis functions and an SBL method was applied to dynamic EIT lung imaging. In this study, 3D CT-based F.E. thoracic cavities considering 5 breath cycle states and the basic thoracic tissues were developed to simulate the measuring process. Quantitative evaluation was performed using a variety of metrics, and an extensive comparison with other reconstruction approaches took place. In vivo imaging using online available data was also carried out. The results show that the proposed (PM-MoM SBL) approach appears to improve spatial resolution in the reconstructed EIT images. Furthermore, despite the method's high time complexity, the lack of necessity for hyperparameter selection reduces the effects of this disadvantage. Future research must be performed in the following directions: A) Optimization of the SBL algorithm in terms of imaging quality and time complexity; and B) Application in 3D multi-layer EIT imaging.