Mutual information as a measure of image quality for 3D dynamic lung imaging with EIT

An array of 32 electrodes was arranged on the subject as two transverse planes of 16 electrodes equidistantly spaced around the chest at approximately the fourth and sixth intercostal spaces, along with an abdominal reference electrode, as illustrated in figure 1. EIT measurements were acquired using the EIT sub-system of the biomedical fEITER instrument (McCann et al 2011), which offers not only a high signal-to-noise ratio approaching 80 dB on bench-top phantoms, but has also been designed to meet the IEC 60601-1 patient safety requirements. In this study, a nearest-neighbour current injection protocol was used with a total of 20 current patterns which included eight independent horizontal injections per transverse electrode plane. Sinusoidal current injections of 0.5 mA amplitude at 10 kHz frequency were input and the EIT instrument recorded nearest-neighbour voltage pairs at 100 fps. EIT data were recorded whilst the subject carried out basic breathing procedures in 1 min blocks. Breathing procedures involved normal and progressively deeper tidal breathing regimes interspaced with reference conditions of typically 5 s breath-holds at both inspiration and expiration. MRI fiducial markers were substituted for the EIT electrodes at all electrode locations immediately after the last EIT data collection, before carrying out an MRI scan using a T2-weighted half-Fourier acquired single-shot turbo spin-echo (HASTE) protocol on a 1.5 T Philips Achieva scanner (Philips Healthcare, Best, the Netherlands) under the same breathing procedures as used during the EIT tests. The MR scans consisted of 1 cm thick contiguous axial slices, which were cardiac gated and tuned for the conditions of normal breath-holding, at maximum inspiration and expiration. The whole procedure was approved by the local ethics committee and the subjects gave written, informed consent.

Zoom In Zoom Out Reset image size

The external boundary shape of the subject was obtained from a transverse MR image slice at the fifth intercostal spacing at a breath-hold at maximum inspiration. The contour of the external shape was defined by 50 equispaced points and the contour of each lung by approximately 25 equispaced points. A 3D finite element (FE) forward model extruded along the caudal-distal axis was then generated using EIDORS 3.6 (Adler and Lionheart 2006, ⁶) calling the Netgen mesh generator (Schöberl 1997). Ider et al (1990) have previously performed a mathematical analysis of extruded models with a single transverse plane of electrodes. The present methodology is an extension of their approach because multiple planes of electrodes can be incorporated into the model. The method of extrusion has been applied to swine data as described by Grychtol et al (2012). The FE forward model consisted of 32 electrodes with typically 20 000 tetrahedral elements and 5000 vertices, modelled to a depth of 10 cm. The transverse MR image slice and surface mesh of a typical forward model is shown in figure 2. All electrodes were modelled to a diameter of 1 cm, which is consistent with those applied to the subject. A high mesh density was chosen in the vicinity of each electrode where the electric field is largest. The lungs segmented from the MR image were also included in the forward model and used as prior information in the forward problem. A 3D FE model was generated with the same boundary shape as the forward model, but with no electrodes and a lower mesh density, for the reconstruction problem. A reconstruction model typically consisted of 10 000 tetrahedral elements and 2500 vertices. From the forward and reconstruction models a mapping matrix was generated to transform between the coarse and fine discretizations, as described by Adler et al (2008), allowing one to solve the forward problem sufficiently accurately on the fine model as well as representing the reconstructed conductivity on the coarse model.

Zoom In Zoom Out Reset image size

As previously discussed, a major source of artefacts in thoracic EIT reconstructions is due to an inaccurately known boundary shape. Given full knowledge of the Neumann-to-Dirichlet (NtD) on the true domain, with conductivity assumed a priori to be isotropic, it was shown by Lionheart (1998) that the data measured on the true domain will only be consistent with an isotropic conductivity in the model domain if the true domain and model domain are related by a conformal map. Reconstruction algorithms have been proposed that rely on the conformal structure of a deformed conductivity. It has been shown by Kolehmainen et al (2013), that the shape and conductivity in 2D can be determined using techniques from quasi-conformal maps and minimally anisotropic conductivities. Lastly, it has been shown by Kolehmainen et al (2007) that the shape and conductivity of an isotropic conductivity in 3D can be determined up to a rigid motion. Specifically, when the metric corresponding to an isotropic conductivity is pushed forward by a domain distortion it is conformally flat and this is equivalent to the vanishing of the Cotton–York tensor. A reconstruction algorithm including a penalty term on this tensor has been proposed.

Alternative reconstruction algorithms to estimate the shape and conductivity from knowledge of the NtD map have recently been proposed in the literature. Firstly, a Bayesian approximation error approach has been proposed by Nissenen et al (2011) to reconstruct the conductivity and a low rank estimate of the shape. This approach treats all sources of modelling error such as the FE discretization and the boundary shape, but not the unknown conductivity, as 'nuisance' parameters. The modelling error covariance is estimated from an ensemble of 150 CT images of thorax shapes and this covariance is then incorporated into a standard regularized Gauss–Newton method to reconstruct the conductivity. Secondly it has been shown by Dardé et al (2013) that the forward operator F(h), now operating on a C¹ perturbation of the boundary given by a function h, with σ fixed, is also Fréchet differentiable at the origin. That is, given $h \in C^{1}(\partial \Omega ,\mathbb {R}^{3})$, and let D_h(x) = x + h(x), x ∈ ∂Ω, represent a perturbation of the boundary, then $F(h) = F(0) + F^{\prime }_{0}h + O(||h||_{C^{1}(\partial \Omega ,\mathbb {R}^{3})}^{2})$ and an explicit form for the derivative operator, $F^{\prime }_{0}$, is calculated. The weak formulation for the derivative falls just short of the conventional H¹ Sobolev space desired for a FE method, but a dual formulation is derived and a sampling method proposed to estimate the derivative and an iterative algorithm outlined to estimate the shape and conductivity distribution.

The approach outlined here is to calculate the Gâteaux derivative of the forward problem with respect to the boundary shape at the electrodes, leaving the rest of the boundary fixed. This method could equally be applied to other boundary locations, but since the largest sensitivity of measured data to shape changes will be at the electrodes, we only calculate this at the electrode positions. In effect we are assuming the shape has been sufficiently well approximated from the MR image and are only accounting for small corrections from electrode movement due to breathing. We denote the domain $\Omega \subset \mathbb {R}^{3}$ with boundary ∂Ω, and the isotropic, real and positive conductivity $\sigma : \Omega \rightarrow \mathbb {R}^{+}$. The boundary has L well separated electrodes, E_l⊂∂Ω, l = 1: L, attached with centre of mass positions $v_{l} \in \mathbb {R}^{3}$. We consider a prescribed electrode current vector, $I=[I_{l}]_{l=1}^{L}$, with ∑_lI_L = 0, with the resulting electrode potential vector, $U=[U_{l}]_{l=1}^{L}$, and internal potential u (defined up to a constant). Let $z=[z_{l}]_{l=1}^{L} \ge 0$, denote the vector of contact impedances present at the electrode-domain interfaces. A finite difference method on the FEM formulation of the forward problem is used, and an approximate derivative with respect to electrode movement, $J^{m} : \mathbb {R}^{3L} \rightarrow \mathbb {R}^{m}$, is calculated, which can be considered as the first partial derivative of the forward problem with respect to each coordinate of each electrode (Soleimani et al 2006, Crabb et al 2012). To calculate this derivative we need some modelling assumptions on the behaviour of electrode movement. Firstly, the shape of the electrode can conceivably deform as it moves along the boundary if flexible 'elastic' electrodes were used. We use an identical perturbation for each node of a given electrode and so we are assuming the electrodes move along the boundary rigidly without changing shape. Secondly, to calculate the normal component of electrode movement in the finite difference approximation, a small piece of conductivity has to be added under each electrode. We assume the conductivity is constant in a neighbourhood of the electrode, and so the conductivity of the small additional piece of tissue added is the same as under the unperturbed electrode. The GREIT reconstruction algorithm (Adler et al 2009), engineered for 2D linear reconstructions, includes a variation on the algorithm discussed in this section to compensate for electrode movement. The array of 32 electrodes, and flexible current injection capability of the fEITER system, allow us to explore the effectiveness of 2.5D reconstructions.

Denote the forward operator $F_{\sigma ,v} : \mathbb {R}^{L} \rightarrow \mathbb {R}^{L}/\mathbb {R}$, mapping a given current pattern applied to the electrodes to the measured voltages on the electrodes. The vector of measured voltages over all current patterns can then be written as V = Z_IF_{σ, v}I, where $Z \in \mathbb {R}^{m \times L}$ is a linear measurement operator mapping the boundary voltages over all current patterns, $I \in \mathbb {R}^{L}$, to the specific measurement protocol. The forward problem is modelled through the complete electrode model, considered the most physically accurate forward model for EIT (Somersalo et al 1992). We discretize the domain into N elements and V nodes and choose a piecewise constant and piecewise linear basis to represent the conductivity and interior potential respectively. The conventional EIT reconstruction problem is to determine σ stably from knowledge of V. It is well known that the forward problem in EIT is Fréchet differentiable with respect to L^∞ conductivity perturbations (Kaipio et al 2000). If there are a total of m measurements over all input current stimulations, and using a piecewise constant representation of the conductivity, the Fréchet derivative simplifies to a Jacobian matrix $J^{c} : \mathbb {R}^{N} \rightarrow \mathbb {R}^{m}$, readily calculated via a piecewise linear FEM (see Polydorides (2002) for details.) The forward problem can now be linearized through a Taylor series expansion as

$$\begin{equation} {\bf V}_{\sigma + \delta \sigma ,{\bf v+\delta v}} = {\bf V}_{\sigma ,v} + J_{{\bf v},\sigma }^{c} (\delta \sigma ) + J_{{\bf v},\sigma }^{m} (\delta {\bf v}) + O(||(\delta \sigma ,\delta v)||^2), \end{equation} \tag{ 1 }$$

where ||(σ, v)|| := ||σ||_∞ + ||v||_∞. The estimation of conductivity and electrode positions x := (σ, v) from the data is approached from a probabilistic viewpoint. The voltages are assumed to be related by V = Z_IF(σ, v)I + n, where $n \in \mathbb {R}^{m}$ represents the measurement channel noise. We assume mean zero Gaussian noise with covariance matrix $\Gamma _{e} \in \mathbb {R}^{m \times m}$ and a Gaussian smoothness prior on the conductivity and electrode positions, with covariance matrix $\Gamma _{\sigma ,v} \in \mathbb {R}^{3L+ N \times 3L + N}$ and mean x_r := (σ_r, v_r). The measurement noise is also assumed to be independent and identically distributed (i.i.d) and the conductivity and electrode position changes are assumed to be independent from one another and separately i.i.d. These assumptions imply that $\Gamma _{e}^{-1} = I \in \mathbb {R}^{m \times m}$ and $\Gamma _{\sigma ,v}^{-1}$ is a diagonal matrix with entries $\Gamma _{\sigma ,v}^{-1}(i,i) = \alpha ^{2}$ if i ⩽ N and $\Gamma _{\sigma ,v}^{-1}(i,i) = \beta ^{2}$ if N < i ⩽ N + 3L. Physically α and β can be interpreted as the ratios of expected changes of conductivity and electrode position to the measurement noise standard deviation.

The maximum a posteriori (MAP) estimate of xⁱ, at the ith temporal data frame, V_i, is equivalent to

$$\begin{equation} {\bf x}_{i} = \mbox{arg min}_{{\bf y}} \big\lbrace (Z_{I}F({\bf y})I - {\bf V}_{i})^{T} \Gamma _{e}^{-1}(Z_{I}F({\bf y})I - {\bf V}_{i}) + ({\bf y} - {\bf x}_{r})^{T}\Gamma _{\sigma ,v}^{-1}({\bf y} - {\bf x}_{r}) \big \rbrace . \end{equation} \tag{ 2 }$$

The minimizer of the functional is approached via a linearized Gauss–Newton method to compute an update ${\bf x}_{i}^{(1)}$ at the ith measured temporal data frame. The initial guess is chosen here as the prior distribution, x_r, and we denote V_r = Z_IF(x_r)I and δV_ir = V_i − V_r. The gradient, ${\bf g}_{ri} \in \mathbb {R}^{N + 3L}$, of the functional (2), at the point x_r and ith temporal data frame is

$$\begin{equation} {\bf g}_{ri} = 2 [ J^{c}({\bf x}_{r}) | J^{m}({\bf x}_{r}) ]^{T} \delta {\bf V}_{ir} =: 2J_{r}^{T} \delta {\bf V}_{ir}. \end{equation} \tag{ 3 }$$

The Gauss–Newton method assumes all second order partial derivatives of the forward problem F_I with respect to σ and v can be neglected. This, and the assumption that the prior covariance matrix is diagonal, allows us to write the Hessian matrix, $H_{r} \in \mathbb {R}^{N + 3L \times N + 3L}$, of the functional (2) at the point x_r, as a 2 × 2 block matrix of the form

$$\begin{eqnarray} H_{r} = 2 \left[\begin{array}{@{}cc@{}}J^{c}({\bf x}_{r})^{T}J^{c}({\bf x}_{r}) + \alpha ^{2}I & J^{c}({\bf x}_{r})^{T}J^{m}({\bf x}_{r}) \\ J^{m}({\bf x}_{r})^{T}J^{c}({\bf x}_{r}) & J^{m}({\bf x}_{r})^{T}J^{m}({\bf x}_{r}) + \beta ^{2}I \end{array} \right] \end{eqnarray} \tag{ 4 }$$

The first update of the Gauss–Newton scheme is then

$$\begin{equation} {\bf x}_{i}^{(1)} = {\bf x}_{r} - \tau ^{(1)} H_{r}^{-1} {\bf g}_{ri}, \end{equation} \tag{ 5 }$$

where τ⁽¹⁾ is a linesearch parameter. We are only interested in the linearized difference between two data frames V_i and V_j, and so just set τ⁽¹⁾ = 1. The linearized reconstruction, δx_ij, between two data frames V_i and V_j is then

$$\begin{equation} \delta {\bf x}_{ij} := {\bf x}_{i}^{(1)} - {\bf x}_{j}^{(1)} = -H_{r}^{-1} J_{r}^{T} ({\bf V}_{i} - {\bf V}_{j}). \end{equation} \tag{ 6 }$$

This solution was computed in EIDORS via the following procedure. Firstly the extruded mesh was generated as described previously and a constant conductivity assigned to the model. Simulated boundary voltages were acquired from this model using the same measurement protocol as the experiment, and in order to get a consistent scaling between the simulated and measured voltages, V_s and V_m respectively, a best fit homogeneous reference conductivity, σ_r, was computed using a formula derived by Kaipio et al (2000). The prior information of the lungs was then included by assigning the conductivity in the lung regions as 0.3σ_r. Using this reference conductivity the movement Jacobian, J^m, was calculated in EIDORS 3.6⁷ along with the standard conductivity Jacobian, J^c, using the default Jacobian function⁸. Image reconstructions were performed using two linear difference imaging methods. Firstly, using the described shape derivative method and secondly, a standard EIDORS method using a one-step linearized Gauss–Newton technique with standard Tikhonov regularization. In the following section we discuss how optimal regularization parameters α and β have been chosen using a mutual information technique.

In this section we describe the EIT and MR image co-registration process, and the method used to calculate the mutual information between the EIT and MR images. The image co-registration process and mutual information calculation provides both a qualitative and quantitative measure of EIT reconstruction performance. Although mutual information techniques have been used in previous multimodal imaging studies (see Pluim et al (2003) for a review) such techniques have yet to be applied in the area of lung imaging using EIT.

The Manchester Confeitir software (McCormick et al 2007) was used to convert the reconstructed EIT data from the irregular tetrahedral mesh into a matrix with 2 mm-cubic isotropic voxels each representing a conductivity change and the coordinates of each voxel were then transformed so that they were in the same space as the MRI data. The software efficiently identifies the tetrahedron whose centroid is closest to the centre of any given voxel. Additionally, the application of a small threshold to the EIT data aided improved clarity of the subsequent visualizations. This typically involved discarding 5% of extreme points from the visualization. The EIT images were transformed to the Nifti data format for subsequent importing into 3D Slicer (Pieper et al 2004). The MRI data were imported into 3D Slicer separately and co-registered using methods validated with bench-top phantoms by Davidson et al (2012). This involved marking up a set of control points around the outside of the torso on both the EIT and MRI data using the fiducial registration module of 3D Slicer, followed by least squared minimization of the distance between the two sets of data in space.

We denote two images A and B, each with N cubic voxels, and each voxel having a positive greyscale value. The probability distribution of A, p_A(a), is defined as the number of voxels in image A that have pixel value a, normalized by N. The joint probability distribution of A and B, denoted p_{A, B}(a, b), is then calculated as the number of times out of N that pixel in A contains a and the same pixel in B contains b normalized by N. The Shannon entropy, in the imaging context, is a measure of the information content of an image, measured in bits. The information content of a single event, that is a particular greyscale value of an image, is proportional to the log  of the inverse of the probability of an event. The total information content of an image, or entropy, is the information content of a single event, weighted by the probability that the event occurs, summed over all events. The total entropy of A, and the joint entropy of A and B, are thus expressed as

$$\begin{eqnarray} &&H_{A} = -\sum _{a} p_{A}(a) \log ({p_{A}(a)}), \nonumber\\ &&H_{A,B} = -\sum _{a} \sum _{b} p_{A,B}(a,b) \log (p_{A,B}(a,b)), \end{eqnarray} \tag{ 7 }$$

respectively. The mutual information of A and B, I_{A, B}, is defined as the total entropy of A and B minus the joint entropy, I_{A, B} = H_A + H_B − H_{A, B}, which equates to

$$\begin{equation} I_{A,B} =\sum _{a}\sum _{b} p_{A,B}(a,b) \log {\frac{p_{A,B}(a,b)}{p_{A}(a)p_{B}(b)}}. \end{equation} \tag{ 8 }$$

It can be shown that 0 ⩽ H_A ⩽ log (N), where H_A = 0 when the image conveys no information i.e. it is featureless or homogeneous, and H_A = log (N) for white noise. It is also true that 0 ⩽ I_{A, B} ⩽ H_A, where I_{A, B} = 0 when A and B have no features in common and I_{A, B} = H_A = H_B = H_{A, B} when A and B are the same. Thus, the larger the mutual information between two images, the more similarities the two images share.

To calculate the mutual information, a threshold is first applied to the MR image to generate a binary image representing just lungs and chest, similar to the segmentation show in figure 2(a). The intensity values of the cubic voxel based EIT image are then scaled linearly so they lie within the same range as the binary MR image. Both images have a finite number of voxels, and the probabilities in (7) are interpreted as histograms with 256 equispaced bins. If a bin is empty, 0 × log (0) is interpreted as 0. The MR and EIT images are both sampled at their voxel centres to estimate H_A, H_B and H_{A, B}, and hence I_{A, B} via (8). This calculation is therefore solely measuring how well a linearized EIT reconstruction is able to find the overall shape of the lungs. We perform the mutual information calculation over a range of regularization parameters, (α,β), for both a standard (β = 0) and shape corrected reconstruction algorithm. This effectively results in a numerical parameter study to maximize the mutual information between MRI and EIT as a function of the parameters α and β.

Figure 3(a) shows typical voltage measurement data for two measurement sites, where i–j–k–l means current injection between electrodes i and j, and voltage measurement between electrodes k and l; site 13–14–15–16 refers to voltages at an anterior right-hand side (RHS) location at the fourth intercostal space, whilst site 23–24–25–26 refers to a posterior left-hand side (LHS) location at the sixth intercostal space. The effect of a brief breath-hold at approximately t = 20 s can be clearly seen along with periodic changes consistent with the breathing rate of the subject. Additionally, the significant differences in the phase of the voltage changes for the two different sites suggests spatial differences of the underlying conductivity distribution within the lungs. Figure 3(b) illustrates raw measurement data from a fourth intercostal LHS posterior site (7-8-9-10) over a 10 s period with data points marked every 10 ms. The observations illustrate the excellent low-noise characteristics (<6 μV rms in recorded human lung data) of the raw human EIT data from the fEITER system at 100 fps. (Note that, despite its visual prominence, the approximate 1.5 μV rms quantization noise level evident in figure 3(b) contributes less than 10% of the total noise power.)

Zoom In Zoom Out Reset image size

Figures 4 and 5 are typical of the 3D EIT image reconstructions acquired at 100 fps using the fEITER system. All images displayed in this paper are difference images relative to a maximum exhalation data frame, and blue and red indicate negative and positive conductivity changes respectively. Figures 4(a) and (b) illustrate difference images for a subject sitting upright at maximum (max) inhalation with the standard and shape correction reconstruction algorithms respectively. It can be clearly seen that the shape correction algorithm yields fewer boundary artefacts and sharper contrast between the lungs and background compared with the standard algorithm. These effects are also clearly visible in the 2D slices from the 3D image reconstructions for a second data set shown in figure 5. Figures 5(a) and (b) correspond to the standard method, α = 10⁻², and the bottom figures 5(c) and (d) correspond to the perturbed method, (α, β) = (10⁻², 4 × 10⁻²). Compensating for movement results in largely lateral changes in the electrode positions as well as a large reduction in boundary artefacts. This is promising as we do not expect any physiological mechanisms to produce large conductivity changes near the boundary.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 shows an image co-registration example of MRI and EIT using 3D Slicer, for the subject in the sitting position at maximum inhalation. The left and right column are superior and inferior transverse planes respectively. The top row is the original MR image, the middle row the EIT reconstruction and the bottom row is the resultant image co-registration. From the images, it can be seen that the EIT reconstruction typically underestimates the physiological size of the actual cross-sectional areas of the lung. Inspection of the complete 3D data along the caudal-distal axis showed this to be true for the entire scanned lung volume.

Zoom In Zoom Out Reset image size

Another observation to note is that the conductivity distribution exhibits a significant twisting effect in the superior transverse slice consistent with the 2D linear image reconstructions obtained by Bikker et al (2011) using the Dräger system. This is likely to be an artefact resulting from the inadequacies of the simple 2.5D extruded model used in the reconstruction. These inadequacies include the lack of inclusion of the heart in the model and possibly the model cut-off height along the caudal-distal axis being too small.

Figure 7 shows EIT and MRI co-registration examples for the subject in the seated and supine positions for the two reconstruction algorithms. The images clearly show that the shape correction algorithm decreases the boundary artefacts, and the EIT reconstructions again are typically found to underestimate the total lung volume. The shape correction algorithm also has an increased tolerance to positional changes of the subject in the supine position. Additionally, moving towards the superior position, the reconstructed lung regions tend to decrease in contrast to that of inferior slices, and we believe this could be due to the natural tapering of the lungs as we move in a superior direction. Three-dimensional 'image fusion' video clips created in 3D Slicer show that the shape correction algorithm provides a more physiologically meaningful reconstruction compared with the standard algorithm for all the data along the caudal-distal axis (see supplementary data, available at stacks.iop.org/PM/35/863/mmedia).

Zoom In Zoom Out Reset image size

Figure 8 shows the mutual information with a lung-segmented MR image as a function of the parameter α, for the shape corrected and standard reconstruction in both the sitting and supine positions at maximum inspiration. For comparison, the maximum possible mutual information, that is the mutual information of the lung-segmented MR image with itself, was 1.51 bits. From figure 8 it can be seen for the subject in the sitting position, there is a maximum of mutual information at approximately 1.11 bits for the standard reconstruction algorithm, with α = 3 × 10⁻², and 1.21 bits for the shape corrected reconstruction, with (α, β) = (10⁻², 4 × 10⁻²). In the supine position, there is a maximum of mutual information of approximately 1.09 bits for the standard reconstruction algorithm, with α = 3 × 10⁻², and 1.20 bits for the shape corrected reconstruction algorithm, with (α, β) = (2 × 10⁻², 4 × 10⁻²). This corresponds to an approximate 10% increase in the maximal mutual information with the shape correction algorithm. It can also be seen at this particular value of β that the shape corrected reconstruction generally has increased mutual information with the standard reconstruction over a wide range of α. The increase in mutual information from a standard to shape correction reconstruction is also visually evident in the images of figures 4 and 5 (which are reconstructions at (α, β) corresponding to the maximum mutual information), when compared to figure 2(a).

Zoom In Zoom Out Reset image size