Improved resolution of D-bar images of ventilation using a Schur complement property and an anatomical atlas.

BACKGROUND: Electrical impedance tomography (EIT) is a nonionizing imaging technique for real-time imaging of ventilation of patients with respiratory distress. Cross-sectional dynamic images are formed by reconstructing the conductivity distribution from measured voltage data arising from applied alternating currents on electrodes placed circumferentially around the chest. Since the conductivity of lung tissue depends on air content, blood flow, and the presence of pathology, the dynamic images provide regional information about ventilation, pulsatile perfusion, and abnormalities. However, due to the ill-posedness of the inverse conductivity problem, EIT images have a coarse spatial resolution. One method of improving the resolution is to include prior information in the reconstruction.
PURPOSE: In this work, we propose a technique in which a statistical prior built from an anatomical atlas is used to postprocess EIT reconstructions of human chest data. The effectiveness of the method is demonstrated on data from two patients with cystic fibrosis.
METHODS: A direct reconstruction algorithm known as the D-bar method was used to compute a two-dimensional reconstruction of the conductivity distribution in the plane of the electrodes. Reconstructions using one step in an iterative (regularized) Newton's method were also computed for comparison. An anatomical atlas consisting of 1589 synthetic EIT images computed from X-ray computed tomography (CT) scans of 74 adult male subjects was computed for use as a statistical prior. The resolution of the D-bar images was then improved by maximizing the conditional probability density function of an image that is consistent with the a priori information and the statistical model. A new method to evaluate the accuracy of the EIT images using CT scans of the imaged patient as ground truth is presented. The novel approach is tested on data from two patients with cystic fibrosis. RESULTS AND
CONCLUSIONS: The D-bar images resulted in better structural similarity index measures (SSIM) and multiscale (MS) SSIM measures for both subjects using the mask or amplitude evaluation approach than the one-step (regularized) Newton's method. Further improvement was achieved using the Schur complement (SC) approach, with MS-SSIM values of 0.718 and 0.682 using SC evaluated with the mask and amplitude approach, respectively, for Patient 1, and MS-SSIM values of 0.726 and 0.692 using SC evaluated with the mask and amplitude approach, respectively, for Patient 2. The results from applying an anatomical atlas and statistical prior to EIT data from two patients with cystic fibrosis suggest that the spatial resolution of the EIT image can be improved to reveal pathology that may be difficult to discern in the original EIT image. The novel metric of evaluation is consistent with the appearance of improved spatial resolution and provides a new way to evaluate the accuracy of EIT reconstructions when a CT scan is available.

INTRODUCTION

Electrical impedance tomography (EIT) is a nonionizing imaging technique in which electrodes are placed around the circumference of a patient's chest and a low‐frequency, low‐amplitude alternating current is applied to the electrodes, inducing a voltage distribution on these same electrodes that is dependent on the conductivity distribution in the interior of the body. From this measured voltage data, an inverse problem can be solved to compute the conductivity in the plane of the electrodes and form a two‐dimensional (2D) cross‐sectional image of the chest by plotting the reconstructed conductivity. Since tissues vary in conductivity as the presence of air and fluids changes, dynamic images can be formed. This property has led to the development of real‐time EIT pulmonary imaging systems for bedside use, such as the commercial systems by Dräger, Timpel, Maltron, and SwissTom . Treatment of mechanically ventilated patients with acute respiratory distress syndrome (ARDS) is particularly timely due to the novel coronavirus COVID‐19. The appropriate setting for positive end‐expiratory pressure (PEEP) on the ventilator for ARDS patients is an active area of clinical research. Recent clinical studies have demonstrated the potential of EIT to guide PEEP titration maneuvers in patients with ARDS,. , , , , , , , , , Other applications of EIT for ARDS patients include detection and monitoring of pneumothorax, detection of patient–ventilator asynchrony, and lung perfusion assessment. ,

EIT may also benefit patients with obstructive lung disease, who often require longitudinal monitoring of their disease. Standard tools include pulmonary function tests (PFTs), chest x‐rays, and computed tomography (CT) scans. Cystic fibrosis patients experience pulmonary exacerbations, air trapping, and mucus plugging, and therefore require pulmonary monitoring from a very young age. However, their care comes with the challenges that very young children are unable to perform PFTs and CT scans result in ionizing radiation delivered to a vulnerable group of patients. EIT has the potential to serve as an as‐needed nonionizing monitoring device for these patients due to its ability to assess regional lung function. There are several recent studies , , , , , , exploring the use of EIT for cystic fibrosis patients.

While EIT has the advantage of being able to provide real‐time functional images at a low cost with EIT systems costing under $100 000, compared to over 1 million dollars for a CT scanner or MRI system, the spatial resolution is on the order of centimeters for EIT compared to millimeters resolution for CT and MRI images. From a physics point of view, the low spatial resolution is due to the long wavelength and high attenuation of the electromagnetic waves used to probe the region. Typically applied currents for EIT are between 10 and 150 kHz with an amplitude between 0.2 and 3.5 mA, which corresponds to wavelengths between 3 km and 30 km. Mathematically, this manifests itself in an ill‐posed inverse problem for determining the conductivity. The data are the current‐to‐voltage mapping, and the conductivity does not have continuous dependence on the data, which is a condition of ill‐posedness. As a result of the discontinuous dependence on the data, large changes in the conductivity distribution in the interior of the body can correspond to very small changes in the measured voltages. For example, in simulations, a 0.35 S/m change in conductivity over one half of a simulated lung resulted in a 49‐mV change in the measured voltages in the first (most sensitive) current pattern and a maximum magnitude change in the current‐to‐voltage matrix of 10−4. This poses challenges for EIT hardware design, and the need for a highly accurate voltmeter in the measurement system. For the reconstruction algorithm, the ill‐posedness manifests itself in ill‐conditioned matrices or blowup of functions computed in the reconstruction algorithm. Regularization is a technique for stabilizing the ill‐posed inverse problem and can be implemented through a penalty term for iterative methods or modification of the scattering transform for D‐bar methods, which are direct (noniterative) solvers. Anatomical priors in the regularization term reduce the ill‐posedness and give the opportunity to provide the reconstruction algorithm with known information about the conductivity distribution. This approach has been followed in iterative reconstruction methods with success. , , , , , , , , , Since direct methods, such as the D‐bar method or Calderón's method, do not have a penalty term in which anatomical information can be provided, spatial priors have been included in the scattering transform for D‐bar , , , , and the Fourier transform for Calderón's method.

Shape reconstruction methods , , , directly incorporate geometry and prior information and preserve sharp edges, reducing the computational burden of the full reconstruction problem. Postprocessing approaches to improve resolution include machine learning methods , and the total variation (TV)‐enhanced D‐bar method.

Previously, a method was proposed to introduce statistical prior information into the D‐bar method based on Schur complement properties. The method presents an improvement of the image obtained by the D‐bar method by maximizing the conditional probability density function of an image that is consistent with a priori information and the model, given a D‐bar image computed from the voltage measurements. Positive results were demonstrated on experimental phantoms using a set of 15 890 simulated samples to generate the joint probability model and covariance matrix.

In the present work, we make use of an archival set of 74 CT scans of adult patients to generate an anatomical atlas of human chest data. The method is then applied to archived data from two cystic fibrosis patients with different levels of disease progression and pathology. To perform the method, an image quality evaluation based on CT scans and the structural similarity index is proposed. The results demonstrate an improvement in resolution, verified through the image quality evaluation method and with CT scans from the two patients at the time of data collection. While the purpose of this paper is to introduce the proposed approach and demonstrate its feasibility, the improved resolution suggests its potential to impact several important clinical applications. For example, detection of regions of air trapping as small as 5 mL and monitoring for their clearance could improve the care of cystic fibrosis patients since clearing mucus plugging, which is evidenced by air trapping, is key to the prevention of lung infection in these patients. , , , Likewise, patients with advanced neuromuscular disorders such as spinal muscular atrophy type 1 and Duchenne muscular dystrophy lose the ability to clear mucus from their lungs and require cough assist systems to aid in their clearance. , , In the intensive care unit, detection of pneumothoraces as small as 5 mL could lead to timely intervention at the bedside. ,

We will refer to the D‐bar method postprocessed with the statistical prior introduced using the Schur complement property as the Schur complement method throughout the text, especially in the figures.

MATERIALS AND METHODS

In this section, we review the mathematical model for the inverse conductivity problem, the D‐bar method, and the Schur complement property.

Governing equations for EIT

Let Ω be a bounded domain in representing a cross‐sectional slice of the chest, the conductivity, and the electric potential. Then the inverse problem of EIT is modeled by the generalized Laplace equation: with the Dirichlet boundary condition where is the measured voltage on the boundary of the region Ω, and the Neumann boundary condition modeling the application of a current density J on the boundary. The map from the current density on the boundary to the voltage that arises on is known as the Neumann‐to‐Dirichlet (ND) map, and it is the discrete version of this mapping that serves as our data. Traditionally in the mathematical literature, the problem is posed with its inverse, the Dirichlet‐to‐Neumann (DN) map, as the data. We briefly describe the construction of the discrete DN map from the measured data. Assume the applied current patterns comprise an orthonormal set consisting of linearly independent vectors with elements , where , , and L is the number of electrodes. Note that the currents must sum to zero over the electrodes in accordance with Kirchhoff's law. Let denote the voltage on the lth electrode corresponding to the kth current pattern . Note that for a unique solution (choice of ground), we require , .

Denote the discrete inner product between vectors u and w by The entries of the (electrode) model approximation ND map are given by where is the area of the ℓth electrode. The discrete DN map .

The D‐bar method and its implementation

D‐bar methods are a class of direct reconstruction methods based on inverse scattering theory and complex geometrical optics (CGO) solutions to related partial differential equations (PDEs). The reader is referred to Refs. [26] and [63] for details and a survey of D‐bar methods. Here we employ the D‐bar method for 2D conductivity reconstruction. , , This method has a real‐time implementation and has been used to image the lungs of patients with cystic fibrosis. , Below we summarize the key equations in the method and their numerical solution and refer the reader to the references above for further details.

The CGO solutions in this work are special solutions to the Schrödinger equation that is obtained from the conductivity equation (1) through the change of variables , : To define the CGO solutions, let z be the spatial point identified with a point in the complex plane and be a nonphysical complex frequency that is introduced as an auxiliary parameter.

Since q is zero outside Ω, Equation (6) can be extended to the plane . It is shown in Ref. [64] that there exists a unique solution to the equation where Equation (8) imposes the asymptotic condition for large or large , and is a Sobolev space. The related CGO solution is then defined by It can be readily seen from (7) and (9) that can be directly computed from by which makes the computation of from the measured data the key procedure to the method.

The function (and hence ) is related to the measured data through the nonlinear Fourier transform of q, known as the scattering transform . The scattering transform can be computed from the integral where Λ1 denotes the DN map corresponding to in Ω. Note that Equation (11) requires knowledge of on the boundary of Ω. While Refs. [64, 65, 66] provide equations and numerical methods for the computation of ψ on from the data, in this work we take the exp‐approximation in which ψ is approximated by its asymptotic behavior on and Equation (11) becomes In this work, we compute images that represent a change from a reference conductivity distribution. Denoting the DN map for this reference conductivity distribution by , the difference scattering transform is defined in Refs. [67, 69] and [67] by

A numerical approximation to can be computed as follows. Let denote the center of the ℓth electrode and write the expansion Then since annihilates the constant function 1, discretizing equation (13) results in where Φ is the matrix of current patterns with entries and is the discretized .

The second step in the exp‐form of the D‐bar method is to compute from at each point of the region of interest (ROI) in Ω, so that σ can be computed from . For this, we must solve the D‐bar equation for numerically. The function satisfies the following D‐bar equation : This equation can also be written in the integral form using the generalized Cauchy integral formula, as follows, and simultaneously we replace by the approximation and denote the corresponding function by :

The integral in (18) is also tantamount to a convolution of the right‐hand side of (17) with the Green's function for the operator, . Defining , Equation (17) can be written as where ⋆ denotes the convolution.

The inverse problem of determining σ is ill‐posed due to the discontinuous dependence of the conductivity on the measured data. This can be seen, for example, in the classic example by Alessandrini. In this example, it is shown that a homogeneous conductivity distribution in a disk compared to one containing a concentric disk of radius r 0 with constant conductivity of contrast A yields DN data that differs in the operator norm of the DN map with an upper bound that can be controlled by r 0, while the amplitude of the conductivity difference can be made as large as desired by increasing A. This ill‐posedness manifests itself in the computation of , which will blow up for large values of , and the blowup radius in the k‐plane decreases as the noise level increases. The signal‐to‐noise ratio (SNR) for the active complex electrode 1 (ACE1) system was measured in the lab and found to range from 32 to 90 dB, per electrode, with the highest SNR on the injection electrode and the lowest SNR on electrodes distant from the injection electrodes. For the clinical situation under investigation, namely monitoring lung function and structure of patients with cystic fibrosis, imaging takes place in a private clinic room, which did not add a measurable amount of additional environmental noise to the signals compared to the laboratory setting.

To handle the noise, the domain of computation of is truncated in the k‐plane. For simplicity, we truncate to a disk of radius R and will refer to R as the truncation radius, setting for , and including the subscript R as a reminder of the truncation, . Theoretical bounds for the truncation radius R based on the noise level are presented in Ref. [66]. Here, a truncation radius of was used, chosen by inspection. Equation (19) can then be solved numerically at each point in the ROI, using a grid in the k‐plane with equal spacing h. The convolution is performed using 2D discrete Fourier transforms: where and are matrices formed by evaluating the functions on the k‐grid, ⊙ is the Hadamard product, and represents the 2D discrete Fourier transform. The resulting matrix equation is then solved using GMRES. (See Ref. [73] for details.) The overall reconstruction procedure is illustrated with a flowchart in Figure 1.

FIGURE 1

Block diagram of the human data reconstruction using the D‐bar method with the approximation

Schur complement property

In this section, we review a relevant property of the Schur complement for computing the conditional probability density between two Gaussian random vectors, particularly when the elements of the two vectors have significant cross‐covariance.

Let be a symmetric positive definite matrix in a block form where , , and . The Schur complements of Γ are defined by If the matrices Γ11 and Γ22 are invertible, the inverse of Γ is given by

Let and be the expected values of the Gaussian multivariate random variables X and Y, and its covariance matrices, and and its cross‐covariance matrices. If the random variable is known (it has been measured), the best estimate for the random variable X (in the mean square sense) is the conditional expectation which we denote by . The conditional covariance is denoted by and their computation reduces to linear algebra. We state the related theorem below for reference.

(Theorem 3.5 [ Assume the joint probability density function of the Gaussian multivariate random variables X and Y is of the form Then the conditional probability density function is of the form where

In such case, we use the notation where In the context of this paper, X corresponds to the actual conductivity and Y to the D‐bar reconstruction.

Postprocessing of D‐bar reconstructions

In this section, we consider the use of the anatomical atlas and Schur complement property to postprocess the D‐bar reconstructions. The anatomical atlas is used to model the marginal distribution . The construction of the model for joint distribution is then carried out via first drawing from and computing the D‐bar reconstructions and then computing the second‐order sample statistics as described in more detail below.

Approximate joint normal model based on conductivity samples and its D‐bar reconstructions

Let be a prior model of the conductivity distribution based on a three‐dimensional thorax, where the mean and covariance are computed from a set of samples of the conductivity distribution. Let refer to the ith sample of this set of samples. Let be the conductivity in the 2D section of in the plane of the electrodes, and let be the vectorization of the discretized . Let the vector represent the vectorized EIT reconstruction of .

To compute the vector , first a finite element mesh that represents the sample is built. Then, the finite element method is applied to solve the forward problem to compute simulated voltages . Additional noise representing the measurement noise of the ACE1 EIT system was then added to the voltages, and an EIT reconstruction algorithm is used to obtain an image from the voltages. In this work, the D‐bar method is used to compute the conductivity distribution , and represents the vectorization of this D‐bar image.

Define the vector to be the stacked vector of and , . The diagram in Figure 2 shows how to compute the vector for each  sample.

FIGURE 2

Block diagram of the computation of from samples

Following the diagram in Figure 3, the second‐order sample statistics of the random variable Z are then computed. The joint sample mean and covariance are computed by taking into account all computed samples of vector .

FIGURE 3

Block diagram of the Schur complement matrices computation

As seen from the form of the conditional mean, Equation (28), the estimate (in the context of this paper) is an affine map. By using the same EIT reconstruction algorithm applied to compute the images , the conductivity distribution may be estimated through the measured data . The postprocessed estimate of the actual σ given the may be computed by The approximate conditional covariance can be computed by

Due to the small number of thorax samples, the marginal sample covariances (diagonal blocks of ) are typically not strictly positive definite and thus not invertible. Therefore, we use a diagonal ridge regression type stabilization () as follows: Then, the matrix and the vector are computed by

Human data

EIT data collected in previous experiments as part of a larger study conducted in accordance with the amended Declaration of Helsinki–Ethical Principles for Medical Research Involving Human Subjects were used. These data were collected at Children's Hospital Colorado, Aurora, Colorado, under the approval of the Colorado Multiple Institutional Review Board (COMIRB) (approval number COMIRB 14‐0652) with informed written parental consent and children's informed assent. EIT data were collected using the ACE1 EIT System during pulmonary function tests (spirometry maneuvers) at the patient's regular clinic visit.

Patient 1 was a 17‐year‐old male in stable condition at the time of data collection, with no air trapping, consolidation, or bronchiectasis noted in the radiologist's report. Inspiratory CT scans in the axial and coronal planes are shown in Figure 4a and 4b, respectively. Patient 2 was an 18‐year‐old male. The radiology report indicated that air trapping was suggested through the majority of the right lung as well as the presence of peripheral mucous plugging most prominent within the right lobe, and cylindrical bronchiectasis throughout the chest bilaterally greater on the right. Inspiratory CT scans in the axial and coronal planes for Patient 2 are shown in Figure 4c, and 4d, respectively.

FIGURE 4

CT scans of Patients 1 and 2. Axial plane scans are one slice in the plane of the EIT electrodes. The images are shown in digital imaging and communications in medicine (DICOM) orientation, with the subject's right at the viewer's left

Anatomical atlas

This prior model considers anatomical information, through the use of segmented images of thoraces, and physiological information, by taking into account the conductivity distribution of the segmented tissues. The computation of sample‐based densities obtained via other imaging modalities and surgical experiments are sometimes called anatomical atlases. [ , Sec. 3.3.5] To build the sample‐based prior, CT scans of 74 adult male subjects were segmented in bones, lungs, heart, and muscles. The CT scans were collected for a previous study approved by the Research Ethics Committee of the University of São Paulo Medical School (CAAE 52619216.2.0000.0068). The choice of patient characteristics can limit the quality of the prior model since the atlas is intended to be representative of the patients under consideration, such as of the same age group and gender, yet rich enough to include examples of pathologies that may be encountered. Thus, increasing the number of samples without making the set too diverse (for instance, including very different age groups) would be expected to improve the performance of the proposed approach.

Each CT slice had 512 × 512 pixels, and the CT slice thickness was between and , depending on the set. The CT images were downsampled to 169 × 169 pixels in the axial plane (X–Y plane) using a nearest‐neighbor interpolation, and 10 slices were used to describe a total of in the cranial‐caudal direction (Z‐axis) centered at the T7 thoracic vertebra. To take into account the variability of the electrode belt position, all possible sets of 10 equidistant slices in a range of were considered, resulting in 1589 datasets. For example, a CT scan set tall with a slice thickness of (44 slices total) provided 4 datasets with 10 slices distant from each other.

The images were rotated by using the center of the T7 thoracic vertebra and the center of the sternum as reference. Also, the images were scaled so that all images have the same rib cage dimension in the T7 thoracic vertebra height, using an affine transformation.

After this transformation, the average contour was computed as the pixels from the body that appear in at least 50% of all three‐dimensional images. Then, the muscle tissue from all images was expanded or cropped to fit this average contour.

Table 1 shows the conductivity distribution of the tissues used to compute the anatomical atlas. The conductivity distribution of the lungs, heart, and muscle was obtained from in vivo experiments on an animal model. Since the conductivity distribution of the bones was not available from in vivo experiments, statistics of the bone conductivity distribution were obtained from in vitro experiments.

TABLE 1

Normal distribution parameters of the segmented tissues

		Mean (S/m)	Std. (S/m)
In vivo	Lungs	0.272	0.240
	Heart	0.460	0.192
	Muscle	0.447	0.140
In vitro	Bone	0.052	0.011

Computation of the approximate joint normal model

To obtain the samples, 10 samples with different tissue conductivities were built from each of the 1589 three‐dimensional images described in Section 2.6. For each sample, the conductivity values of the segmented tissues (bones, lungs, heart, and muscles) were randomly sampled according to their conductivity distribution from Table 1. When the sample value of the bones was smaller than , and when the samples values of lungs, heart, and muscles were smaller than , the values were disregarded and new values were sampled. To take into account variation in conductivity between the lungs such as when one or both lungs have a pathology, the conductivity of the left and right lungs were sampled independently in 50% of the samples.

A finite element mesh was built using the average contour of the anatomical atlas, and two sets of voltages were computed for each of the 15 890 conductivity datasets. One set of voltages was computed with the conductivity distribution of all segmented tissues (bones, lungs, heart, and muscles), and other set of voltages (reference voltage) was computed with heart and lungs segmented tissues replaced by the conductivity of the muscle. The voltages were computed with the same current patterns and number of electrodes used in the human data collection. To simulate the noise measurement, Gaussian noise with the same mean and standard deviation of the ACE1 EIT system was added to the voltages.

For each set of electrical potentials , the D‐bar method presented in Section 2.1 was applied, and the images were obtained. The reference voltage of the dataset was used to compute the DN map , and the image represents the D‐bar method reconstruction of the variation of heart and lungs tissues in the domain, using as reference the set described above in which the heart and lungs are replaced with the conductivity of muscle. The D‐bar method settings were the same as those used to reconstruct the human data.

To compute , a Gaussian interpolation was applied to each of the 15 890 conductivity datasets and 2D images were computed. For difference images, is defined to be the difference between the conductivity dataset with all segmented tissues and the reference image where heart and lungs segmented tissues are replaced by the conductivity of the muscle. The reference dataset simulates a homogeneous conductivity distribution with the conductivity of muscle, and so each reconstructed image represents the difference of the sample from a background reference of muscle. The reference is used for the computation of the DN map and scattering transform . ,

Figure 5 shows 2D images that represent the statistical distribution of the conductivity variation for the anatomical atlas. Figure 6 shows the diagonal of the conditional covariance matrix for the anatomical atlas for two different values of κ.

FIGURE 5

Mean and standard deviation of the conductivity variation distribution for the anatomical atlas. The images are shown in DICOM orientation, with the subject's right at the viewer's left

FIGURE 6

Standard deviation (square root of the diagonal) of the conditional covariance matrix for the anatomical atlas of two different values of κ. The images are shown in DICOM orientation, with the subject's right at the viewer's left

One‐step Newton's method

To compare the results obtained by the D‐bar method and the proposed method, reconstructions using the one step of a (regularized) Newton's Method (OSNM) were computed. This method is based on the NOSER code, and a simplified description is presented below. The rationale for the choice of a one‐step Newton's method algorithm for comparison is that it is the most commonly used method with a real‐time implementation for EIT reconstruction, and the work uses this method for reconstruction in its comparison of functional EIT lung images with CT scans.

Let be the vector of voltages that represents the solution of the forward problem of the conductivity distribution σ, and σ0 is the linearization point where Newton's method is applied. By truncating the Taylor series expansion of , we have where denotes the computed voltages of the conductivity distribution σ0 and represents the Jacobian matrix that maps the conductivity distribution to the computed voltages at the conductivity distribution σ0 (linearization point). The computation of the Jacobian matrix is described in the literature [ , Sec. 4.1].

From Equation (34), the difference image can be computed by The term is not invertible, and the generalized Tikhonov regularization is applied, where Q is the Tikhonov matrix and α is the regularization parameter that adjusts the weight of the term .

For the human thorax data, the is the data averaged over the 100 first frames, is computed using the finite element method and a high‐pass Gaussian filter is applied to regularize the one‐step Newton's method. The difference image represents the conductivity distribution variation due to the difference between a measured voltage frame and the average voltage represented by .

For a proper comparison, the finite element mesh used to estimate the difference image may have the same contour of the domain and position of the electrodes used to compute the D‐bar reconstructions. The contact impedance of the electrodes is modeled using the complete electrode model with a prismatic element with a triangular base to represent the electrode face.

Image quality evaluation

The CT scans were prepared for use as ground truth for comparison with structural similarity index measures (SSIM) using the following steps. First, from the segmented CT scans, the intersection of the two outer contours and the union of the lung regions were selected to form a mask of the body and lungs. The extracted lung region was then assigned a Hounsfield unit of −1. This is the lung mask that is used as ground truth comparison for the “mask” approach for image evaluation. To take relative conductivity values into account, the aligned inspiratory and expiratory CT scans were subtracted inspiratory minus expiratory, and the lung mask was applied to extract the lung region. The Hounsfield units of the subtracted image were then normalized to lie between 0 and −1. These images are used as ground truth comparison for the “amplitude” approach for image evaluation.

The SSIM and multiscale structural similarity (MS‐SSIM) are methods used to quantify the quality of digital images when a reference image is known. The SSIM is defined as follows. Given an image x and a reference image y, consider three comparison measures: luminance ℓ, contrast c, and structure s defined by where and are the mean values of x and y, respectively, and are the variances of x and y, respectively, is the covariance of x and y, , , and . Here L is the dynamic range of the image x. Then the SSIM is defined by More, generally, weights can be applied to each of the factors in (40).

The MS‐SSIM provides a convenient way to incorporate image details at different resolutions. The main idea of the MS‐SSIM is the following. A low‐pass filter is iteratively applied to the input image x, and the filtered image is then downsampled by a factor of 2. The original image is indexed as Scale 1, and iterations are applied to obtain the final Mth image (Scale M). At the jth scale, the contrast comparison (38) and the structure comparison 39 are calculated and denoted by and , respectively. The luminance comparison (37) is computed only at Scale M and is denoted as . Then the MS‐SSIM is defined by Here the exponents in (41) are all set to 1.

The SSIM method has been used as a metric to evaluate EIT algorithms improvements. , , By using simulated and experimental phantoms, EIT images are estimated and the results are compared to images that represent these phantoms. However, the way that the SSIM metric is used may not be appropriate for human data since the reference image that describes the conductivity (or resistivity) distribution of an in vivo human thorax are not feasible to be obtained.

To quantify the improvement of the EIT images of ventilation human data, an image quality evaluation process based on CT scans and the SSIM is proposed. This procedure aims to use the SSIM metric by considering aspects of the SSIM analysis as luminance, contrast, and structure even when the electrical property of the reference is unknown.

To compute the reference image for the ventilation human data, CT scans of the patient at the end of the inhalation and expiration are required. For each CT scan, the slice representing the plane of the electrodes is selected, and four reference points are determined. Then, the slices are registered based on these set of reference points, and the contour of the body and the lungs regions are identified in each slice.

To determine the reference points the center of the vertebral column and the center of the sternum are identified. The first point (P1) is placed in the center of the vertebral column and the second point (P2) in the center of the sternum, the other two points are determined by drawing a perpendicular line that crosses the middle distance between P1 and P2 and ends at the edge of the lung and the rib cage in both directions. The last two points (P3 and P4) are set at the end of the drawn line, P3 is placed in the edge of the left lung, and P4 is placed in the edge of the right lung.

Later two new images are generated by thresholding the CT scans to identify the contours and lungs regions. The first image (Ref. Image 1) is composed of the intersection of the contours and the union of the inspiratory and expiratory lungs regions. This reference image represents the area where the ventilation may occur during the respiratory cycle. The second reference image (Ref. Image 2) is composed of the intersection of the contours and the normalized difference between the inhalation and expiration slices in the region where the lungs are identified. This reference image represents the area where the ventilation may occur and at the same time underestimates the area where the heart movements occur during the ventilation.

As an example, Figure 7 shows the steps of the process to register the images for CT scans of Patient 1. Figure 7a,b shows the reference points of the inspiratory and expiratory slices. Figure 7c shows the registered images overlaid, and the white marks in the image show the four reference points. The reference points in the inspiratory slice are rotated to make points P1 and P2 vertically aligned, and they are used as reference to register the two slices by applying a 2D projective geometric transformation to the images.

FIGURE 7

Reference points and the geometric transformation for Patient 1

Figure 8a shows the contour and lungs regions overlaid. Then the regions are shifted to vertically centralize the contours in the image, and Figure 8b is created by computing the intersection of the contours and union of the lungs regions. The contour and background are set to 0.0, and the lungs region is set to −1.0.

FIGURE 8

Reference images, Patient 1

To compute Figure 8c, an image representing the Hounsfield unit value only in the lungs regions of the inspiratory slice is subtracted by an image representing the Hounsfield unit value only in the lungs regions of the expiratory slice, the values smaller than zero are set to 0.0, the result is normalized and all values are shifted by −1.0 resulting in the image that represents Ref. Image 2.

To perform the SSIM and MS‐SSIM, an EIT image formed by the subtraction of the reconstruction of the full inspiration from the reconstruction of full expiration is computed. Then all positive values are set to zero, and the image is registered by applying a 2D projective geometric transformation to adjust the contour of the EIT image and the contour of the reference images. Comparisons to Ref. Image 1 are referred to as “mask” evaluations, and comparisons to Ref. Image 2 are referred to as “amplitude” evaluation.

RESULTS

In this section, we show the results of the Schur complement method applied to the two patient datasets. Also shown are reconstructions from one step of a (regularized) Newton's method and the D‐bar reconstructions prior to postprocessing with the Schur complement. Difference images and SSIM measures using the “mask” and “amplitude” image evaluation approaches are provided.

Difference images

Difference images for Patients 1 and 2 were computed using the D‐bar method as described in Section 2.2. The reference voltages for the reconstructions were computed from data averaged over the 100 first frames. The normalized mean of the raw D‐bar images for Patients 1 and 2 are shown in Figure 9. To compute this parameter, the mean of all pixels in each raw D‐bar image is computed and normalized taking into account the average of the first 100 raw D‐bar images. After the computation of the mean of the images, the curve is normalized by subtracting the mean of the first 100 values and dividing the result by the maximum range of the first 100 values. In the spirometry maneuver performed by each subject, several normal breaths are taken, followed by a deep inhale prior to a forceful, rapid exhale with the goal of expelling all air from the lungs. Following the long, forceful exhale, a deep breath is taken. Movies of the reconstructions are shown in Ref. [84].

FIGURE 9

Normalized mean of the raw D‐bar images for Patients 1 and 2

Figures 10 and 11 show reconstructions at maximum inhalation (first column), maximum expiration (second column), and the difference between these, which was computed by subtracting the image in column 2 from the image in column 1. The frame number (time) corresponding to maximum inhalation is indicated with the first marker (marker symbol +) in Figure 9a,b. The second marker (marker symbol x) corresponds roughly to the point at which as much air as possible had been expelled from the lungs, which is the frame number used for the expiratory images in Figures 10 and 11.

FIGURE 10

Reconstructions for Patient 1. First row: Reconstructions using one‐step Newton's method. Second row: Reconstructions using the D‐bar method. Third row: Postprocessed D‐bar images with the Schur complement method using . Fourth row: Postprocessed D‐bar images with the Schur complement method using . The first column contains images of full inspiration. The second column contains images of full expiration. The third column is the difference between full inspiration and full expiration, formed by subtracting the image in Column 2 from the image in Column 1. The images are shown in DICOM orientation, with the subject's right at the viewer's left

FIGURE 11

Reconstructions for Patient 2. First row: Reconstructions using one‐step Newton's method. Second row: Reconstructions using the D‐bar method. Third row: Postprocessed D‐bar images with the Schur complement method using . Fourth row: Postprocessed D‐bar images with the Schur complement method using . The first column contains images at full inspiration. The second column contains images at full expiration. The third column is the difference between full inspiration and full expiration, formed by subtracting the image in column 2 from the image in column 1. The images are shown in DICOM orientation, with the subject's right at the viewer's left

SSIM measures

The reference points and geometric transformations of the CT scan preparatory to creating a reference image for the SSIM measures are depicted in Figure 12 for Patient 1 and Figure 13 for Patient 2. The segmented scans for each patient with the contour and lungs regions overlaid are shown in Figure 14. Since the relationship between the Hounsfield units in the CT scan and conductivity values in the EIT difference images are unknown, the values of each were normalized to a dynamic range of −1 to 0 before computing the SSIM and MS‐SSIM measures of comparison. The adjusted contours and reference images for the “mask” and “amplitude” measures for Patients 1 and 2 are given in Figures 15 and 16, respectively. The segmented lung regions from the reconstructions representing the difference between inspiration and expiration by the one‐step Newton's method, D‐bar method, and Schur complement method with and were also normalized to have a dynamic range between 0 and −1 and are plotted in Figure 17 for Patient 1 and Figure 18 for Patient 2. SSIM and MS‐SSIM measures using the “mask” and “amplitude” approach of image evaluation are given in Tables 2 and 3 for Patients 1 and 2, respectively.

FIGURE 12

Reference points and geometric transformation for Patient 1

FIGURE 13

Reference points and geometric transformation for Patient 2

FIGURE 14

Contour and lungs regions overlaid

FIGURE 15

Contour and Reference images, Patient 1

FIGURE 16

Contour and Reference images, Patient 2

FIGURE 17

Segmented lung region from EIT reconstructions, inspiration minus expiration, Patient 1

FIGURE 18

Segmented lung region from EIT reconstructions, inspiration minus expiration, Patient 2

TABLE 2

SSIM measures, Patient 1

	Evaluation method	OSNM	D‐bar	Schur 10−5	Schur 10−12
SSIM	Mask	0.60825	0.72216	0.73790	0.72772
	Amplitude	0.49720	0.63623	0.64382	0.63105
MS‐SSIM	Mask	0.61776	0.68823	0.71813	0.71173
	Amplitude	0.57669	0.66269	0.68189	0.67509

Abbreviations: MS‐SSIM, multiscale structural similarity index measures; OSNM, one‐step Newton's method; SSIM, structural similarity index measures.

TABLE 3

SSIM measures, Patient 2

	Evaluation method	OSNM	D‐bar	Schur 10−5	Schur 10−12
SSIM	Mask	0.60475	0.68119	0.72349	0.74175
	Amplitude	0.54916	0.64633	0.66403	0.68327
MS‐SSIM	Mask	0.67426	0.67503	0.71159	0.72594
	Amplitude	0.64429	0.64889	0.67682	0.69155

Abbreviations: MS‐SSIM, multiscale structural similarity index measures; OSNM, one‐step Newton's method; SSIM, structural similarity index measures.

Discussion

Application of the proposed algorithm to a specific clinical patient consists of the following steps. First, an anatomical atlas is required. In practice, this would be on file and previously computed from a collection of CT scans of a representative group of patients. The atlas consists of previously computed D‐bar reconstructions of each conductivity distribution in the atlas as well as the mean, variance, and standard deviation of the set of images. From the EIT data, a preliminary D‐bar reconstruction is computed for each frame. The D‐bar algorithm can be implemented in real time , and so this step can be running during the EIT data collection. To postprocess the images, the Schur complement matrices are computed as in Figure 3 and applied to each image frame.

The normalized mean of the raw D‐bar images in Figure 9 clearly shows a relative volume of air in the lung as the patient performs the spirometric maneuver. While the raw D‐bar images in Figures 10 and 11 clearly show conductivity differences in the lungs between inspiration and expiration, the Patient 2 right lung ventilation deficiency is not easily interpreted. The one‐step Newton's method reconstructions have some fundamental differences from the D‐bar images. While neither method achieves good lung separation, Newton's method reconstructions have more connected lung regions near the heart than the D‐bar reconstructions, which can be seen in the segmented images in Figures 17 and 18. The D‐bar reconstructions have some artifacts near the underarms which compress the lung region toward the center of the images. When compared with the SSIM measures, the D‐bar reconstructions have a closer agreement with the CT scan‐derived reference image using the mask or amplitude evaluation and SSIM or MS‐SSIM for both patients.

The postprocessed images yield an improvement in the spatial resolution, which can be seen qualitatively by comparing the postprocessed EIT images to the axial CT scans in Figure 4. It can be seen from the axial CT scan for Patient 2 that the right lung can accommodate less air than the left, and the presence of air trapping, as noted by the radiologist, reduces the difference between inspiration and expiration in the right lung. This is clearly evident in the postprocessed reconstructions for Patient 2. For Patient 1, the SSIM measures indicate improved agreement with the CT scan‐derived reference image using the Schur complement method with using both the mask and amplitude evaluation measures. However, for , an improvement in the SSIM measures over D‐bar only holds for the mask evaluation method. For Patient 2, the use of the Schur complement method results in an improvement of SSIM measures for both evaluation methods and both values of κ.

Related work includes Ref. [78] in which functional EIT images were compared to thoracic CT scans in pigs. In this work, normalized time difference EIT images of conductivity change were reconstructed using a one‐step Newton's method with the mean over the dataset serving as reference. Heart and lung regions were identified using an unsupervised learning method and used to extract the corresponding signals. It was found that the EIT pixels with the strongest heart and lung signals were located inside the anatomical boundaries of their respective organs in the CT scans. In Ref. [50], D‐bar images were postprocessed using a convolutional‐neural network. Results on data collected on saline‐filled tanks with two different EIT systems demonstrated improved resolution measured by SSIM.

CONCLUSION

The use of a statistical prior built from an anatomical atlas to postprocess 2D cross‐sectional EIT images of the human chest computed by the D‐bar method may be a practical approach to improve the resolution of real‐time images. This may help in the real‐time detection of pathology in the bedside imaging of patients with ARDS and chronic lung disease since EIT is a nonionizing imaging modality that can be performed continuously at the bedside or as needed in a clinical setting.

CONFLICT OF INTEREST

The authors have no conflicts to disclose.