Electrical-impedance-tomography imaging based on a new three-dimensional thorax model for assessing the extent of lung injury

Electrical impedance tomography (EIT) is a technique that calculates the distribution of electrical impedance of a living body by measuring the surface voltage of the living body, which is available for continuous monitoring of the lung function to reduce the risk of intensive lung injury. The global inhomogeneity index (GI) is employed to assess the degree of ventilation homogeneity at different levels of lung injury. The GI values calculated in the current research are all based on two-dimensional EIT reconstruction images; however, this method does not correctly detect low levels of lung injury. In this paper, we establish a true 3D thorax model and propose a method for calculating the GI based on 3D EIT reconstruction images to evaluate lung injury. Simulations reveal that this method can accurately reflect the disease state of lung injury compared with the 2D GI calculation method, and even mild damage of lung injury can be adequately detected.Electrical impedance tomography (EIT) is a technique that calculates the distribution of electrical impedance of a living body by measuring the surface voltage of the living body, which is available for continuous monitoring of the lung function to reduce the risk of intensive lung injury. The global inhomogeneity index (GI) is employed to assess the degree of ventilation homogeneity at different levels of lung injury. The GI values calculated in the current research are all based on two-dimensional EIT reconstruction images; however, this method does not correctly detect low levels of lung injury. In this paper, we establish a true 3D thorax model and propose a method for calculating the GI based on 3D EIT reconstruction images to evaluate lung injury. Simulations reveal that this method can accurately reflect the disease state of lung injury compared with the 2D GI calculation method, and even mild damage of lung injury can be adequately detected.


I. INTRODUCTION
Lung injury is a common thoracic surgical disease mainly caused by various internal and external factors, such as severe infection, trauma, shock, inhalation of harmful gases, and poisoning. 1,2 Research suggests that since lung injury can lead to acute respiratory distress syndrome (ARDS) and even acute respiratory failure, its mortality rate is as high as 50-70%. Research 1,2 shows that lung injury develops rapidly, and failure to diagnose in time can lead to serious or even life-threatening complications. Therefore, it is desirable to have a real-time and safe monitoring tool in early diagnosis of patients with lung injury. 3 Common approaches for diagnosis and monitoring of the disease state in these patient groups are based on lung function tests, e.g., spirometry. However, classical lung function tests are only able to provide global parameters representing the entire lung. 4 In recent years, medical imaging technologies, such as computed tomography (CT) and magnetic resonance imaging (MRI), have been employed to obtain high-resolution anatomical images of the thorax so that the regional information of injured lung can be acquired. 5 However, both CT-scan and MRI are costly and depict only a snapshot of the lungs' state. Dynamic characteristics, such as development of respiratory capacity and degree of lung injury, cannot be captured. As a result, these methods are not suitable for monitoring lung injury in real time. 6,7 Electrical impedance tomography (EIT) is an emerging technology for medical imaging. 5,8 Its basic principle is to apply a certain frequency of low-voltage alternating safety current through a set of electrode arrays arranged on the body surface, measuring the voltage data through the scanning electrode to extract the electrical impedance information of the tissue or organ related to the physiological and pathological states of the human body. The advantages of EIT are its low cost, portability, and nonionizing features; thus, it has been investigated by many research groups for clinical applications. 6,7,9 Many algorithms have been developed for EIT imaging, such as the integral equation method 10 and regularization methods. 11 To improve the quality of reconstruction, iterative methods 8,11 have also been used at the expense of reducing imaging speed in modern medicine. Recently, the fusion of machine learning method and EIT is emerging 12 so that artificial-intelligence technology can be used for EIT reconstruction. In addition, EIT can add valuable information through its ability to provide regional information. Based on the reconstructed images, EIT can create quantifiable clinically relevant measures to assess the present disease state of a patient and its trends. [6][7][8] Quantitative EIT measures can be divided into three groups. The first group only calculates the averaged (or summed) pixel values of functional measures across the entire image or its sections. For instance, the sum of tidal impedance changes in all image pixels provides an estimate of regional tidal volume (VT) in the studied thorax slice. The second group aims at characterizing the regional respiratory system (CRS) under quasistatic but mostly dynamic conditions. The third group aims at characterizing the spatial distribution of ventilation. Examples are the global inhomogeneity (GI) or the coefficient of variation (CV). 13 Since ventilation should be evaluated in lung injury diagnoses and further treatment, the parameters of the first two groups study the entire dynamic process of respiration, which is not suitable for evaluating ventilation in lung injury. The GI aims to describe ventilation inhomogeneity within the lungs, which is one of the important factors for the early diagnosis and treatment of lung injury.
GI has been used for the diagnosis of lung injury based on two-dimensional (2D) EIT images, which are not suitable for the diagnosis and monitoring of mild lung injury. [6][7][8]14 In this paper, the GI parameter for 3D EIT reconstruction is proposed to assess the extent of lung injury. The forward and inverse problems of EIT were constructed based on a true 3D thorax model, considering the true shape of the human thorax boundary. Simulation experiments were conducted to verify the accuracy of the proposed method. In this paper, all our results and conclusions are drawn under noiseless condition.
The rest of this study is organized as follows. In Sec. II, we build the model of the forward problem for a true 3D thorax based on sequential CT images. In Sec. III, we demonstrate an efficient subdivision method for the inverse problem according to the true 3D thorax structure. In Sec. IV, we describe the evaluation of the GI for lung injury based on 3D reconstructed images of EIT. Section V provides simulation results of the proposed method. Finally, Sec. VI states the conclusions and describes planned future work.

II. EIT FORWARD PROBLEM FOR TRUE 3D THORAX
A. Model construction of true 3D thorax for EIT forward problem The true 3D thorax model is constructed according to the human thorax structure, which is based on sequential CT images. For accurate extraction of useful information from CT images, each CT image must be processed. 8 First, we use the partial-differentialequation (PDE) method to eliminate edge noise. Then, we binarize the preprocessed CT image, and a threshold segmentation algorithm, namely, maximum interclass variance algorithm, is adopted to extract the contour of the lung and thorax. The detailed segmentation algorithm will be discussed in Sec. IV. Finally, a contour extraction algorithm based on binarized images is used to improve the accuracy of contour extraction. The rule for judging whether the current pixel is a boundary pixel is expressed by According to the contour extraction algorithm for binary images, e(i, j) is the pixel of a CT image after binarization and E(i, j) is the pixel of the extracted boundary for the CT image. When E(i, j) = 1, the current pixel is a boundary pixel and is reserved; otherwise, the current pixel is deleted. As a result, the pixels retained in the image compose the image boundary. 15,16 The specific image processing flow is shown in Fig. 1. Taking the CT image in Fig. 2 as an example, the specific threshold segmentation process is shown in Fig. 3. The extracted outlines of the thorax and lung based on the CT image are shown in Fig. 4(a). The contour sequences of the thorax and lung based on sequential CT images are extracted and shown in Fig. 4(b). We need CT images only when extracting the contour of the thorax and lung. An optical camera can also be used without CT images. The detailed procedure for constructing the true 3D thorax model is shown in Fig. 5. The contours of the thorax and lung are extracted and stacked to obtain the preliminary model of the thorax in Mimics. Then, SolidWorks is used to smooth the 3D model so that the forward problem could be stably solved. Finally, the lung and thorax models are combined by COMSOL Multiphysics. The forward problem of the 3D model is solved by COMSOL. 8 To obtain 3D reconstruction results, multilayer electrodes must be constructed. We design a 3D electrode structure composed of 96 round electrodes, which are equally placed in six sections of the model with 16 electrodes per layer. The interval between each layer is 2.5 cm, and the space between adjacent electrodes per layer is 5 cm. The diameter of each electrode is 1 cm. The true 3D thorax model with electrodes is shown in Fig. 6. The respective imaging positions of different sections are z = 7.5, 10, 12.5, 15, 17.5, and 20. 5,17,18 In this study, injected currents of 10 mA at 3 kHz and the adjacent measure pattern are used for simulation. The finite element model (FEM) is employed for the forward problem. The accuracy of the numerical approximation is largely dependent on the construction of the finite element meshes. If the mesh discretization is too coarse, it may lead to inaccurate solution of the forward model and low spatial resolution of the reconstructed images. If the mesh discretization is ad-hoc and too fine, it requires a significant amount of memory storage and computational power to generate such a mesh and even more to solve the inverse problem due to an increase in the number of unknowns that need to be solved for. The volume of the thorax discussed in this paper is about 460 cm 3

B. Lung injury model
To simulate different severities of lung injury, we established eight lung models based on the true 3D thorax, as shown in Fig. 8. We discuss two static lung conditions: the first condition represents the lung at end-expiration and the second condition represents the lung at end-inspiration. The condition at end-expiration is modeled by setting the conductivity of the lung tissue to σexp = 120 mS/m. Other tissues not belonging to the lung are set to σ bkg = 480 mS/m. End-inspiration is simulated by changing the conductivity of the lung tissue to σ insp = 60 mS/m. The conductivity of lung injury does not change in both static lung conditions and remained at σ injury = 120 mS/m. 6,7 The volume ratio of the lung injury tissue to the overall lung tissue is defined as η. According to different degrees of lung injury, the η values and lung injury levels according to eight different degrees of lung injury models are shown in Table I.

A. Mathematical model for EIT
The aim of image reconstruction for EIT is to obtain the conductivity distribution σ by using the boundary voltage vector V and injected current vector I. With the aid of physical modeling and discretization of the FEM, a deterministic observation model of EIT can be written as where U(σ, I) is the forward model mapping σ and I to V and R(σ) is the model mapping to resistance. The model of V = R(σ)I depends nonlinearly on the conductivity σ and linearly on the current I. If changes in conductivity are small, the inverse problem can be solved with sufficient accuracy by considering the linearized equation system,  where δσ ∈ R n×1 (n is the number of pixels of the reconstructed image) is the change in conductivity σ and δU ∈ R m×1 (m is the number of independent voltage measurements) is the perturbation of boundary voltage due to the change in σ. J ∈ R m×n is the Jacobian matrix, which denotes the partial derivatives of voltages with respect to conductivity. 19,20 In the present work, we focus on "time-difference EIT," which is designed to reconstruct relative changes in conductivity distribution rather than absolute values of conductivity. Changes in conductivity are denoted where σ contains the present conductivity values and σ ref is the vector of reference conductivity. Throughout this paper, σ ref is considered to be the conductivity distribution at end-expiration of the patient, while σ represents the conductivity at end-inspiration. 6,7

B. Subdivision method of inverse problem
To obtain a 3D EIT reconstructed image, the conductivity distribution of the layers in which electrodes are located is reconstructed. The boundaries of reconstructed images are related to the location of electrodes and the thorax contour of the human body, which are various based on different parts of the body. In this study, the image boundary of each electrode layer is accurately obtained according to the forward problem of the true 3D thorax model to guarantee the authenticity of the image boundary. We perform inverse problem subdivision within the boundary coordinates for each layer. Taking the fourth layer of the true 3D thorax model as an example, the procedure is as follows:

Obtain boundary of reconstruction
The coordinate sequence X ∈ R k×1 , Y ∈ R k×1 of boundary points is obtained from the true 3D forward problem of the thorax model, where k is the number of boundary points. The maximum and minimum coordinate values for the X-axis of the boundary points are recorded as Xmax and X min , and the maximum and minimum coordinate values for the Y-axis of the boundary points are recorded as Ymax and Y min , respectively. A rectangle, which has the length of Xmax − X min and the width of Ymax − Y min , is constructed with the center point of ( X max −X min 2 , Y max −Y min 2 ). The rectangle is divided into n * n grids, as shown in Fig. 9, according to the resolution of the reconstructed image (in this paper, the resolution is 32 × 32). The side length of each grid is 1 n (Xmax − X min ). The coordinate values of each grid are obtained by Eqs. (5) and (6). Then, the boundary points of the true 3D thorax model are included in the rectangle, We divide the boundary points into upper and lower parts along the red line formed by two points (0, Y max −Y min 2 ) and (X max − X min , Y max −Y min 2 ). The coordinates of the upper and lower boundary points (xu,yu) ∈ R k 1 ×2 and (x d, y d ) ∈ R k 2 ×2 are shown as the blue dot in Fig. 10, where k 1 and k 2 correspond to the number of boundary points, respectively.
Since the square grids divide the thorax model into n units along the X-axis direction, the abscissa set A of boundary mesh points can be written as The coordinate sets of the upper and lower boundary mesh points are shown in Eqs. (8) and (9), according to (7), where yU 1 , yU 2 , . . ., yUn and yD 1 , yD 2 , . . ., yDn are the ordinates of the upper and lower boundary mesh points, respectively, which can be obtained by the interpolation method according to the known A, (xu,yu) ∈ R k 1 ×2 and (x d , y d ) ∈ R k 2 ×2 . Owing to the fact that the boundary points of the thorax model are not necessarily on the boundary mesh grid, the linear interpolation method is adopted to obtain the coordinate of grid points of the meshed boundary. In this AIP Advances where l = 1, 2,. . ., s, s = , xM is the nearest boundary mesh point from the left-hand side of xu (i+1) , xN is the nearest boundary mesh point from the right-hand side of xui, yUm is the ordinate of the mth boundary mesh point, and k = y u(i+1) −y ui x u(i+1) −x ui . Finally, the coordinate sets B and C of the upper and lower boundary mesh points are obtained, and the coordinates after linear interpolation are shown by the red dot in Fig. 11.

Subdivision inside the boundary
The grid contained inside the boundary points of the thorax after linear interpolation is defined as the valid grid. Yi is the ordinate of the ith grid. The effective value of a grid is shown in If H(Xi, Yi) = 1, the current grid is marked as a valid grid; otherwise, the current grid is an invalid grid. Figure 12 shows the valid grids within the interpolated boundary points, where the green grids are valid grids, and the white grids are invalid grids. The coordinate data, number, and quantity of the valid grids are saved for later image reconstruction. To display the subdivision process more intuitively, the flowchart of the subdivision process is shown in Fig. 13. In this paper, the subdivision of each layer is stacked to obtain multilayer 3D subdivision. This provides visualization of several tomographic cross sections ("slices"), which are on the plane of the electrode structure. 21 Multilayer subdivision is shown in Fig. 14.

IV. CALCULATION METHOD OF GLOBAL INHOMOGENEITY INDEX
According to the GI calculated based on 2D images, [4][5][6] in this paper, we propose a method for calculating the GI based on 3D EIT reconstructed images. The GI is calculated based on the following formula: where δσij is the conductivity change of the ith pixel in the jth layer, δσ lungij is the conductivity change of the ith pixel in the jth layer belonging to the lung, Nj is the number of pixels belonging to the lung in the jth layer, and M is the number of image layers, which is six in the present study. The threshold for distinguishing between the lung region and the nonlung region is determined by the maximum interclass variance algorithm. The gray value of the reconstructed image is denoted f (x, y, z), with L gray levels (L is chosen as 20 in this paper). The gray level is represented by the variable bi, i = 0, 1, . . ., L − 1, and b 1 < b 2 <⋯< bL− 1 . p(i) is the probability of gray level bi, which is calculated as follows: where fi is the number of all pixels with gray level bi and N is the number of pixels for reconstructed images, which is defined as follows: Nj.
The segmentation threshold bt of the thorax and lung tissue can be expressed by where ∂ 2 (t) is the interclass variance, which is defined as The average gray value λ 0 , λ 1 of the thorax and lungs and the total average gray value λ of the image are as shown in Eqs. (18)- (20), where are the probabilities for the appearance of the thorax and lung in the reconstructed image, respectively, and ω 0 + ω 1 = 1. The set of lung pixels is defined as The pixel values greater than bt are assigned to the lung region and the rest to the thorax region. 22

A. Optimal excitation method selection
Excitation mode refers to how the exciting current is injected into the test field. For the same system, different excitation and measurement methods directly impact the sensitivity, and the amount of measurement information and antinoise performance of the test field ultimately affect the imagery quality and accuracy of measurement parameters. 23 Four excitation modes are compared in this paper. To facilitate the comparison of imaging results, for all excitation modes, the measurement method is normalized to the same layer adjacent electrode pair voltage measurement (except for the excitation electrodes). 24 The details of the four excitation modes are as follows: (1) AD: The current is applied through adjacent electrodes, and the voltage is measured sequentially from all other adjacent electrode pairs without the pairs containing one or both of the current electrodes. (2) QSD: There are six electrodes between the two excitation electrodes. A more uniform current distribution is obtained in this method. To evaluate the impact of different excitation modes on the result of 3D reconstructed images, the imaging relative error (RE) is adopted, where σij is the estimated conductivity distribution of the ith pixel in the jth layer and σ * ij is the corresponding true conductivity. Nj is the number of pixels in the jth layer. M is the number of reconstructed image layers, according to the true 3D forward problem thorax model, which is shown in Fig. 6. M is taken as 6. Figure 15 shows the variation of the imaging RE with the number of iterations in different excitation modes based on the 3D reconstructed images by the conjugate gradient (CG) method.
As illustrated in Fig. 15, the imaging RE of the AD method is less than that of the other three excitation modes. To simplify the calculation in real reconstruction, the iterative number is empirically selected as 85 for all of the methods. As a result, the 3D reconstructed images based on the four excitation methods are shown in Fig. 16, where t is the time of image reconstruction. For the convenience of comparison, the gray levels are normalized to the range from 0 to 1. 25 It can be seen from Fig. 16 that the imaging error of the AD mode is smaller than that of the other excitation modes, which means that using the AD mode we can obtain relatively high precision image reconstruction results. Therefore, the AD mode is selected as the optimal excitation method. The 3D reconstructed images with eight different severities of the lung injury are established based on the optimal excitation method.

B. Reconstructed images with eight different severities of the lung injury
To compare the accuracy of calculating GI values based on 2D and 3D reconstructed images, differential images are reconstructed between end-expiration and end-inspiration lungs with different severities of lung injury. Figure 17(a) shows the 2D reconstructed images at the third layer of electrodes placed in the thorax model. Figure 17(b) shows the 3D reconstructed images with eight different severities of lung injury. All the images are reconstructed by the CG method with 85 iterations. The average reconstruction times of 2D and 3D images are 0.603 and 0.916 s, respectively. For the convenience of comparison, the gray levels are normalized to the range from 0 to 1. 25 It can be seen from the results that 3D reconstructed images reflect the conductivity distribution in a certain space instead of a plane. As the degree of lung injury increases, the amount of ventilation gradually decreases.

C. Calculation of GI based on reconstructed images
The GI values based on 3D and 2D reconstructed images for the eight severities of lung injury are calculated according to (12) and the methods in Refs. 4, 6, and 7, respectively. To facilitate the comparison between the different GI calculation methods, the GI is normalized to range from 0 (mild lung injury) to 1 (severe lung injury). The normalized GI values for different severities of lung injury are shown in Fig. 18. Figure 18 shows that the normalized GI values increase almost linearly from 0 to 1. The increasing GI based on the 2D reconstructed images is only 10% for the first two severity levels, which corresponds to 15% of the lung tissue affected by the injury. However, there is more than a 25% increase in GI corresponding to the same severity level based on 3D reconstructed images. are more sensitive than those calculated based on 2D reconstructed images, especially for a low or moderate injury level. The 3D imaging contains multiple layers of body physiological information; however, 2D imaging contains only one layer of body physiological information. Lung EIT is a truly 3D imaging problem. In addition, since the human thorax has the characteristic of three dimensions, the induced currents for EIT are not constrained to the plane within the excitation electrodes. Therefore, 3D imaging is in line with the mechanism of electrical measurement. 3D EIT imaging provides more information of the thorax and lung so that more comprehensive physiological information can be observed from different image layers, which is dependent on the measurements from different electrode planes. The GI calculated based on 3D EIT imaging is more accurate than that based on 2D imaging, considering the 3D characteristics. As a result, the GI based on 3D imaging is more objective. Compared with the traditional method, the GI results based on 3D reconstructed images are closer to the ground truth and thus can more accurately reflect the increasing disease state of a lung than 2D EIT, even for low levels of lung injury.

D. Statistical analysis
In this paper, we studied 30 samples of human CT sequence images with mild lung injury. To demonstrate the effectiveness of the new method, the paired sample test was separately conducted on the two GI calculation methods and the true GI values of lung injury. Taking the GI calculation method based on the 3D EIT reconstruction image as an example, the specific steps of the significance test are as follows:

ARTICLE
scitation.org/journal/adv H 0 represents that there is no significant difference between the true GI and the GI calculated based on the 3D EIT reconstructed images. H 1 represents that there is a significant difference between the true GI and the GI calculated based on the 3D EIT reconstructed images. The data were compared with the results of a two-tailed test and the level of significance α is 0.05. (2) Calculated the statistics from where d = are the mean and standard deviation of the paired sample difference, respectively. n is the number of paired samples and μ 0 is the difference of the paired-sample mean. In this paper, μ 0 = 0.
(3) Calculate the P value using SPSS statistical software. 26,27 If P > α, the hypothesis H 0 is established; otherwise, the hypothesis H 1 is established. The statistical analysis of calculating GI based on different reconstructed images methods is shown in Table II. Table II shows that the GI calculated based on the 3D EIT reconstructed images does not differ from the true GI (P > 0.05). In contrast, the GI calculated based on the 2D EIT reconstructed images and the true GI differed significantly (P < 0.05). We can thus see that the GI values from the calculation method based on the 3D EIT reconstruction images are closer to the true GI values. This method can more accurately reflect the state of mild human lung injury, and it is therefore of great significance for the timely diagnosis and treatment of lung injury.

VI. CONCLUSIONS
In this paper, we propose a new 3D modeling of the thoracic EIT method for assessing the extent of lung injury, which can effectively increase the accurate diagnosis of lung injury in small areas. We built the model of the forward problem for the true 3D thorax based on prior information, which can obtain the real thorax and lung boundary. 3D images of different degrees of lung injury models were reconstructed by increasing the number of electrodes and selecting a better excitation method. Simulation results revealed that the GI calculation method based on 3D EIT reconstruction images can improve imaging accuracy and better reflect the disease state of lung injury. Even low and moderate-level lung injuries could be adequately detected. In this work, the subdivision of each layer is stacked to obtain multilayer 3D subdivision, and the 3D characteristics of the electric field were not fully utilized. Further studies should focus on direct 3D imaging methods based on multilayer electrode measurements to obtain more accurate conductivity distribution information.