Accurate Measurement of Airway Morphology on Chest CT Images

2.1. Synthetic Modeling of Airways

In order to generate reliable synthetic patches of airways, the main aspects of the structure of interest as well as the characteristics of the CT scanner with regard to resolution, PSF, and imposed noise have to be reproduced. Based on the knowledge that on reformatted axial plane airways have tangent vessels [10], each airway patch consisted of two bright ellipses (inner and outer walls) with a dark central zone (airway lumen) and zero, one or two tangent vessels, represented by bright ellipses rotated around the airway. The parameters to create the synthetic airways were randomly chosen based on physiological values and are reported in Table 1. Although the creation of a synthetic airway presents some limitations, we think that the proposed model represents an appropriate simulation that helps a neural network learn the main features of real airways. Also, using the multi-scale particle extraction method described in [11], 2D patches can be easily extracted along the airway’s main axis, which is given by the first eigenvector of the Hessian matrix. For this reason, we do not consider 3D patches, which due to the different tubular profiles and a wide variation of 3D orientations that should be taken into account would increase the complexity of the modeling.

To reproduce the structure of the parenchyma, a Gaussian smoothing (with a standard deviation of 5) was applied to Gaussian distributed noise, to create some broadly correlated noise, which made a texture of multiple structures that mimicked the parenchyma. Afterward, the correlated noise was altered to have a mean intensity of −900 HU and a standard deviation of 150. All values were empirically chosen.

All patches were created starting at a super-resolution of 0.05 mm/pixel in a sampling grid of 640 × 640 pixels. To obtain the final patch, the obtained images were first down-sampled to a resolution of 0.5 mm/pixel. Then, a PSF was simulated to mimic the blurring caused by the image reconstruction process. To this end, due to the small size of the patch, we assumed that the PSF can be approximated by means of a spatially locally invariant Gaussian function, as demonstrated in [12]. The standard deviation of the Gaussian filter was randomly chosen in an empirically determined range of 0.4 to 0.9 mm to simulate the differences in the PSF across CT scanners and manufacturers. Finally, a spatially correlated Gaussian noise was added to the image based on Gaussian distributed random noise smoothed with a Gaussian filter (with a standard deviation of 2), with the empirically determined mean of zero and standard deviation of 25. As a last step, the image is cropped to a 32 × 32 pixels grid.

2.2. SimGAN Refinement

Although the proposed generative model simulates reasonably well the geometrical aspects of the structure of interest, the generated patches still may present differences to patches extracted from real structures. For this reason, we implemented a SimGAN refinement, similar to the one described in [9], to improve the quality of the synthetic patches. SimGAN makes use of simulated and unsupervised learning by using a generative adversarial network (GAN) that consists of both a generator (refiner) and a discriminator. The purpose of the refining step is to trick the discriminator in deciding whether an image is a synthetic or real image.

For the implementation of this network, we pre-trained the refiner on synthetic images with 1000 steps and a batch size of 256, while the discriminator was pre-trained on real patches (extracted using the multi-resolution particles method described in [11], initialized with the technique of [13]) and refined patches, obtained from the pre-trained refiner, with 100 steps and a batch size of 256. The number of steps was the same as in [9]. Then, the adversarial training of the SimGAN network was trained for 10,000 steps, batch size of 256, and all parameters and loss function set as in [9]. An example of a generated synthetic airway is shown in Fig. 1.

2.3. Measurement of Airway Morphology

To extract both measurements for airways, we implemented a 9-layer 2D network, which consists of seven convolutional layers, five of which had stride 1 and two had stride 2, and two fully-connected layers (see Fig. 2). The network regresses the measure of the central structure in a patch 32×32 pixels, a size chosen to include enough neighborhood information for big structures, without losing specificity for small and thin features. To train the network, we used an Adam update (β₁ = 0.9, β₂ = 0.999, ϵ = 1e⁻⁰⁸, decay = 0.0) with a specifically customized loss function that combines the absolute relative error and the precision of the measure to improve the network performance and stability (see Sect. 2.4).

The network was trained on a NVIDIA Titan X GPU machine, using the deep learning framework Keras [14] on top of TensorFlow [15], for 300 epochs at a learning rate of 0.001 and batch size of 64.

2.4. Customized Loss Function for Airway Morphology Measurement

When trying to accurately measure small airways with sizes at image resolution level, typical approaches usually have problems of under- or over-estimation. For this reason, in this paper we suggest the usage of a new loss function that combines the loss of the relative error over all images, ${\mathcal{L}}_{\mu }$, and the precision of the measure over 25 replicas of the same structure,${\mathcal{L}}_{\sigma }$:

where y is the true measure of a synthetic patch, $\widehat{y}$ is the measure predicted by the CNR, and λ defines the weight of ${\mathcal{L}}_{\sigma }$ with respect to ${\mathcal{L}}_{\mu }$. The definition of ${\mathcal{L}}_{\mu }$ is given by:

where N indicates the total number of patches. On the other hand, the loss term for the precision, ${\mathcal{L}}_{\sigma }$, is computed over a number of replicas of the same geometric model (with fixed physical dimensions) to which varying PSFs are applied and a different number of airways and vessels are added with varying locations and rotations. This way, the network learns to accurately measure the structures of interest regardless of possible confounding factors inside the patch. The definition of ${\mathcal{L}}_{\sigma }$ is given by:

where N represents the total number of images, and M indicates the number of replicas considered. In this work, we used M = 25.

Since for airways lumen radius and wall thickness are measured simultaneously, for this structure the two terms of the loss, ${\mathcal{L}}_{\mu }$ and ${\mathcal{L}}_{\sigma }$, are given by the sum of the corresponding loss computed independently for the two measures.

Since we noticed that the measurement of small airways (lumen less than 1.0 mm) was the most affected by a high standard deviation, we also empirically assigned a higher weight to the precision term of these structures so that they acquire more importance when computing the loss. Therefore, Eq. 1 becomes:

where λ = 2.0 has been empirically selected, l indicates the airway lumen, wt stands for wall thickness, and

and

2.5. Training Set Deftnition

The training dataset consisted of 100,000 × 25 replicas of the same geometric model, to which varying PSFs were applied, and different additional vessels were added at varying locations and rotations. Therefore, a total of 2,500,000 training patches were used. Conversely, for the validation set we generated 1,000,000 patches (40,000 × 25 replicas).

The values of the parameters used for the creation of the images were randomly chosen in ranges that were empirically defined based on physiological measures of the structures of interest, as shown in Table 1. We trained the network using all images refined by SimGAN.

Finally, in order to help the network focus more on geometry than intensity values, during training, we applied a data augmentation that in addition to adding random noise it also randomly inverts intensity values inside the patches. Furthermore, we introduce a small random shift and random axes flipping to the patch to improve the learning of the network.

2.6. Experimental Setup

We evaluated the proposed approach for airway measurements on both synthetic and in-vivo cases. For the synthetic validation, we first generated a dataset of 200,000 patches (with random values chosen in the range of Table 1) that were used in three different experiments. First, we evaluated the accuracy of the algorithm by calculating the relative error (RE) between the CNR measurement and the ground truth defined by our geometrical model when varying lumen and the wall thickness size. To compare our results to the state-of-the-art methods, we also computed the absolute error obtained for airways with a wall thickness of 1.0 mm at the image resolution level (0.5 mm).

In order to demonstrate the ability of the method to accurately measure the structures of interest regardless of presence of noise and smoothness, as a second experiment we generated 100 images for each level of noise (σ_n ∈ [0, 40] HU) and for each level of Gaussian smoothing (σ_s ∈ [0.4, 0.9] mm) and computed the mean RE (in percentage) across the 100 patches. We repeated the same experiment first fixing the wall thickness at 1.5 mm and considering three values of airway lumen (small: 0.5 mm; medium: 2.5 mm; large: 4.5 mm), and then fixing the airway lumen at 1.5 mm and using three wall thickness values (small: 0.5 mm; medium: 1.2 mm; large: 2.0 mm).

As a final test on synthetic images, we compared the proposed method for airway measurement to FWHM and ZCSD computing the mean RE (in percentage) on patches of different sizes.

As a further validation, we tested the performance of the algorithm on a CT airway phantom of known lumen size and wall thickness. The phantom was constructed using Nylon66 tubing inserted into polystyrene to simulate lung parenchyma surrounding the airways. Non-overlapping, 0.6 mm collimation images, 40 cm FOV, were acquired using a GE Siemens Sensation 64 CT scanner and reconstructed with a standard reconstruction kernel. Eight tubes with varying wall thickness and lumen diameter (reported in Table 2), as measured by a caliper, were studied. An image taken from the CT scan of the phantom showing the eight tubes is presented in Fig. 3.

As a final experiment, since an accurate and reliable in-vivo ground-truth is very complicated to obtain, we performed an indirect validation by means of a physiological evaluation. To this end, we computed the Pi10 parameter with our approach and with ZCSD, and analyzed its correlation to FEV1% on 590 clinical cases, with airway particles extracted using [11]. Pi10 is a metric of airway thickness that is computed measuring the square root of the wall area across the whole airway tree and regressing the value at a hypothetical airway with an internal perimeter of 10 mm. The wall area is found by subtracting the area of the lumen from the airway area, while Pi is computed from the lumen radius.

3.1. Synthetic Evaluation

Figure 4 shows the tendency of the RE for predictions obtained on the synthetic data when varying the lumen radius (Fig. 4a) and the wall thickness (Fig. 4b) of the airway. As expected, the error is small for airways with a large lumen (Fig. 4a), while it increases (with a tendency to under-estimate the measure) for lumens smaller than 1.0 mm, although it is always below a 10% RE. Regarding the wall thickness (Fig. 4b), a significant under-estimation error is obtained at sub-voxel levels (below the image resolution of 0.5 mm), while a tendency to over-estimation is obtained when the wall thickness is bigger than 2.0 mm.

On average, an absolute RE of 6.3% is obtained for airways with a wall thickness of 1.0 mm, while when the airway wall thickness is at the image resolution (0.5 mm) the absolute RE is at 13.09%. These REs are significantly lower than those previously reported in the literature for structures of similar sizes [5,7].

Results obtained when fixing three values of airway lumen (0.5, 2.5, and 4.5 mm) and three values of wall thickness (0.5, 1.5, and 2.5 mm) and varying the level of noise and smoothing are presented in Fig. 5. As shown, for both measurements the RE is stable across the different levels of noise and smoothness. While for medium and large structures a very high accuracy is obtained (RE close to 0), the smallest structures (generated with airway lumen or wall thickness at the image resolution of 0.5 mm) are the one confusing the network the most determining also a bigger standard deviation. In all cases, the RE is stable when varying noise and smoothness, and the bias introduced by the CNR is small, with a little under-estimation for small structures, as expected. For small wall thicknesses (0.5 mm), when the smoothing level is low (<0.6 mm) a very small RE is obtained, while this error increases when applying higher levels of smoothing (>0.6 mm).

Finally, Table 3 shows the mean relative error (in percentage) obtained for different sizes of wall thickness using the proposed method in comparison to ZCSD and FWHM on the 200,000 testing patches. Three wall thickness intervals were chosen: lower than 0.7 mm, between 0.7 mm and 1.5 mm, and bigger than 1.5 mm. As shown, while traditional methods tend to have a very high relative error, especially for small airways, the proposed method yields a very high accuracy and outperforms them. Similar results were obtained for the airway lumen.

3.2. Phantom Evaluation

The relative error obtained measuring the wall thickness of the eight tubes of the phantom using the proposed method (CNR) in comparison with traditional techniques are presented in Table 4. For completeness, the relative error obtained when measuring the lumen with CNR (not measured by traditional methods) is also reported. The proposed CNR has the lowest RE for all considered tubes, with the exception of tube C where FWHM gives the best result, and in general is able to well measure the wall thickness even for small and thin tubes, as in case of tube D. Although there is variance in the RE for the measurement of the wall thickness of all tubes, this variance is smaller than the one obtained using traditional methods, that for some tubes seem to really confounded. An important aspect to notice is the small RE obtained for all tubes when measuring the lumen radius with the proposed CNR.

3.3. In-Vivo Indirect Evaluation (FEV1% in Correlation to Pi10)

Table 5(a) shows the Pearson’s correlation coefficient between FEV1% and the Pi10 metric computed with our approach and ZCSD in airway patches extracted from a real CT. The correlation coefficient between FEV1% and Pi10 calculated by ZCSD and CNR was −0.38 and −0.54, respectively, indicating a significantly higher correlation of the Pi10 computed by the CNR with FEV1%. This result suggests that the proposed method could potentially be used to accurately measure FEV1% in patients with COPD.