Generative-based Airway and Vessel Morphology Quantification on Chest CT Images
Accurately and precisely characterizing the morphology of small pulmonary structures from Computed Tomography (CT) images, such as airways and vessels, is becoming of great importance for diagnosis of pulmonary diseases. The smaller conducting airways are the major site of increased airflow resistance in chronic obstructive pulmonary disease (COPD), while accurately sizing vessels can help identify arterial and venous changes in lung regions that may determine future disorders. However, traditional methods are often limited due to image resolution and artifacts.
We propose a Convolutional Neural Regressor (CNR) that provides cross-sectional measurement of airway lumen, airway wall thickness, and vessel radius. CNR is trained with data created by a generative model of synthetic structures which is used in combination with Simulated and Unsupervised Generative Adversarial Network (SimGAN) to create simulated and refined airways and vessels with known ground-truth.
For validation, we first use synthetically generated airways and vessels produced by the proposed generative model to compute the relative error and directly evaluate the accuracy of CNR in comparison with traditional methods. Then, in-vivo validation is performed by analyzing the association between the percentage of the predicted forced expiratory volume in one second (FEV1%) and the value of the Pi10 parameter, two well-known measures of lung function and airway disease, for airways. For vessels, we assess the correlation between our estimate of the small-vessel blood volume and the lungs’ diffusing capacity for carbon monoxide (DLCO).
The results demonstrate that Convolutional Neural Networks (CNNs) provide a promising direction for accurately measuring vessels and airways on chest CT images with physiological correlates.
Lung Airway Vessel Deep-learning regression Subvoxel resolution
In the last decades, changes in peripheral lung airways and vessels have been shown to be key elements of the pathophysiological development of lung diseases associated to chronic tobacco exposure, among others [26, 35, 4, 14]. Therefore, the ability to detect small changes in the morphology of those structures is crucial to the early diagnosis of the disease and the development of new therapies that can reserve those early changes. The development of those therapies is hampered by the lack of accurate biomarkers related to small airways and vessels morphological changes that can assess and monitor the therapeutical response [4, 14]. For example, chronic obstructive pulmonary disease (COPD) is mostly affecting peripheral airways  and airway wall thickness predicts airflow obstruction and physical impairment . In asthma, several studies have shown that CT-measured airway wall thickness predicts the severity and duration of attacks [3, 26].
Additionally, several studies have demonstrated that small pulmonary arteries become smaller and shrink at subsegmental levels in patients with COPD [15, 9], and endothelial dysfunction may be caused by both pulmonary and extra-pulmonary vascular alterations in COPD [35, 4]. Therefore, having an automated method for airway and vessel morphology assessment will help precise measurements of the geometrical properties of bronchial and venous trees which, in turn, may lead to improved diagnosis and open the door to new studies on lung disorders.
While in the past several methods have been proposed with the aim of helping physicians accurately locate small pulmonary airways and veins on chest CT images,[31, 5], up to date not much work has been proposed for sub-voxel morphology assessment. Traditional approaches for airway wall thickness detection are based on edge-detection methods, that, although limited by the Nyquist theorem, use the reconstructed CT signal to analyze properties of the structure directly. Among them, the full width at half max (FWHM) , for which the true edge of an ideal step function undergoing low-pass filtering is located at the FWHM location, is one of the most typical algorithms. In , a modified FWHM, as proposed in [22, 23] is used as available in the Virtual Bronchoscopy program (Siemens Medical Solution), to determine the effect of bronchodilation on airway metrics derived from airway wall thickness reflecting airway disease in patients with COPD.
The main drawback of the FWHM method is that if the airway wall is considered as a laminated structure whose scale is close to the scanner scale resolution, the edge principle the FWHM is based on would be violated. Therefore, as reported in , the measurements provided by FWHM are biased towards under- or over-estimation of the inner or outer boundary, respectively. For this reason, in  a model-based fitting method, assuming a Gaussian point spread function (PSF) and an idealized step-like model for the airway, was proposed. However, this method is computationally very expensive and does not take into account deviations from the proposed model, such as areas where the airway wall is in contact with a vessel.
Another popular approach to measure airway walls involves the use of the zero crossings of the second order derivative (ZCSD) and the phase congruency of the local phase [32, 8, 21], which characterizes lumen-to-wall and wall-to-parenchyma transitions. While shown to provide better localization of the airway wall than FWHM and to be less sensitive to different reconstruction kernels and radiation doses, this technique is still computationally expensive and, as admitted by the authors, still sensitive to overestimation.
A method based on 2D dynamic programming approaches was also proposed to detect both the inner and outer boundaries in cross-sectional images, utilizing cost functions combining the first and second derivatives . However, validation was perfomed only on phantom CT images and in-vivo results were not provided.
More recently, a new algorithm for airway wall segmentation that combines coarse airway segmentation and optimal graph construction was proposed . Comparison to manual annotation and correlation of Inner Volume (IV) and Wall Volume Percentage (WV%) with predicted forced expiratory volume in one second (FEV1%) were used to validate the method. However, the proposed technique requires a complex and time consuming parameter tuning and it still computationally expensive. Moreover, no comparison with FWHM and ZCSD is provided. A similar approach was also proposed in [19, 20], where a more accurate validation was presented. Nevertheless, as claimed by the authors, the suggested technique may not always be achievable as it requires the pre-segmentation step to define complete airway trees, which is not always feasible.
While different methods have been proposed in the literature for measuring airway wall thickness and lumen, not many studies have tried to accurately measure the vessels radius with high precision. Some approaches are mainly based on edge detectors and the relation between scale and physical radius based on an idealized Gaussian model for vessels . In a recent work, a technique for vascular morphology quantification based on FWHM was also proposed . However, the method uses the centerline and its distance to the vessel border (previously segmented) to compute the vessel radius. While acceptable for big vessels, this method is prone to high errors for peripheral vessels. Also, no validation of the vessel morphology quantification is provided.
The main drawback of all traditional methods is that they suffer from over- and under-estimation errors, especially for small stuctures at the scanning resolution .
To overcome these issues when measuring small airways and vessels, we propose the development of a convolutional neural regressor (CNR) , that learns the main characteristics of the structures, automatically regress the airway wall thickness, airway lumen, and vessel radius on small 2D patches extracted in the orthogonal plane along the main axis of the structure of interest, can quantify clinically significant morphology, and is robust to variation in CT parameters.
Since manually measuring airways and vessels on clinical CT images is a complicated process and traditional methods are not reliable in providing highly accurate measures of small structures, in this paper we also propose to use a generative approach based on synthetic models for airways and vessels that aim at mimicking their main distinguishing features with known physical dimensions for training the network. A Simulated and Unsupervised Generative Adversarial Network (SimGAN)  is then used to refine the generative model. We use the data created by means of the proposed generator and refined with SimGAN to train the CNR that, through a specifically designed loss function aimed at jointly optimizing both accuracy and precision, to regress the size of the structures of interest.
For validating the proposed approach we used refined synthetic airways and vessels created by the model generator, patches extracted from an specifically-designed airway phantom, and in-vivo cases for an indirect validation on both airways and vessels in comparison to traditional methods.
2 Paper Contributions
In the present paper, we first demonstrate on synthetic generated data that the same technique can be successfully applied both for airway and for vessel assessment with very accurate results, extending our workshop paper , in which our CNR and model generator were introduced, but validation was limited to show only preliminary outcomes.
Then, three experiments are performed on clinical cases that show very encouraging results and demonstrate that the technique is very accurate and precise even when the initial conditions, scanner brand, and scanner protocols are modified.
Finally, we compute an indirect in-vivo validation via a correlative analysis with physiological factors for both airways and vessels. For the bronchial tree, we used 3,038 subjects to assess the association between Pi10 and FEV1% and a comparison to the 3D airway measurement software package Pulmonary Workstation (VIDA Diagnostics, Inc., Iowa City, IA, USA)  is provided. Also, linear models were created to look at the association between our measurements and functional small airway disease (fSAD) using the parametric response mapping (PRM) method, a non-invasive imaging biomarker that identifies small airway loss, narrowing, and obstruction . Then, to carry out a physiological evaluation of our ability to accurately measure vessel lumen radius, we also analyze how three metrics of blood volume distribution (total blood volume, TBV, blood volume of vessels of less than 5 , BV5, and blood volume of vessels of less than 10 , BV10) correlate to lungs’ diffusing capacity for carbon monoxide (DLCO) adjusted by site altitude.
3 Material and Methods
A color-coded scheme of the proposed method for both airway and vessel assessment is presented in Fig. 1, whereas the different parts of the workflow are detailed in the following sections.
3.1 Airway and Vessel Model-based Generator
A model-based generator (MBG, black box in Fig. 1) was developed to synthesize patches that simulate the main characteristics of the structure of interest as well as the CT scanner attributes including resolution, PSF, imposed noise and blurring.
To avoid possible issues in measuring the structures due to their orientation, we simulated airways and vessels on a reformatted axial plane. This is a reasonable simulation, as when considering in-vivo CT images the first eigenvector of the Hessian matrix (or the structure’s centerline) can be used to extract the patches along the airway/vessel’s main axis. An example of a vessel reformatted on its main axis can be seen in Fig. 2.
Based on the structures anatomy and their surrounding morphology , the airway geometrical model 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. Conversely, to simulate vessels, a bright central ellipse was generated, with zero, one or two simulated airways in the surrounding. Since only arteries are tangent to bronchi, the simulated airways might also be separated from the vessel.
Since the structures are simulated as ellipses, for the measurement of airway and vessel lumen we consider the nominal radius of the shape, given by:
where and represent the maximum and minimum diameters of the ellipse.
In Tab. 1, all parameter values (randomly chosen based on the physiological values proposed in ) used to create the structures are reported. To mimic the parenchyma and its texture of multiple structures, we used a Gaussian distributed noise to which a Gaussian smoothing (standard deviation = 5) was applied. The correlated noise was then altered to have a mean intensity of -900 HU and a standard deviation of 150. All values were chosen from nominal values.
|Lumen radius (LR)||[0.5, 6.0] mm||[VR, VR + 0.8] mm|
|Airway wall thickness||[0.1*LR + 0.2, 0.3*LR + 1.5] mm||[0.1*LR + 0.2, 0.3*LR + 1.5] mm|
|Number of vessels||[0, 2]||[1, 3]|
|Number of airways||1||[0, 2]|
|Vessel radius (VR)||[LR, LR + 0.8] mm||[0.5, 4.5] mm|
|Skewness of reconstruction||[-25, 25] degrees||[-25, 25] degrees|
|Airway Lumen Intensity||[-1150,-1050] HU||[-1150, -1050] HU|
|Airway Wall Intensity||[-500, 50] HU||[-500, 50] HU|
|Vessel Intensity||[-50, 50] HU||[-50, 50] HU|
|Noise Level||25 HU||25 HU|
|Smoothness Level||[0.5, 0.875] mm||[0.5, 0.875] mm|
From an accurate analysis of CT images, it is also possible to notice that some peripheral vessels may be located close to the chest wall, the presence of which may affect the measure of the CNR. To deal with this potential issue, we randomly added one or two curved regions at the two opposite corners or borders of the patch of some vessels with radius lower than 1.5 mm. For this added region, we used uniformly distributed values in a range HU.
For training the CNR, we used patches of 3232 pixels, as enough neighborhood information can be included for big structures without losing specificity for thin and small characteristics. A resolution of 0.05 mm/pixel in a sampling grid of 640640 pixels is initially used to generate the images. Then, a down-sampling to 0.5 mm/pixel is applied and a simulated PSF is applied to mimic the CT blurring caused by an image reconstruction process.
Estimating a PSF can be challenging due to its complexity and variance across manufactures . However, experimental measurements of the PSF of a CT scanner demonstrated that it can be approximated by a 3D Gaussian function with the assumption that it is locally space invariant and separable . Due to the small size of the generated patch, this assumption is valid and we approximate the PSF by a spatially invariant Gaussian function with standard deviation randomly chosen in a range that simulates the PSF variation across different CT scanners and manufacturers.
Then, we apply a Gaussian random noise smoothed with a Gaussian filter (standard deviation = 2), with zero mean and a standard deviation , values chosen based on nominal noise characteristics in high dose CT scans .
As a last step, an image cropping is implemented to obtain the final 3232 pixels patch. In Fig 3a-d, an example of airway (top) and vessel (bottom) patch as created by the MBG at each step is shown.
3.2 SimGAN Refinement
The proposed MBG simulates reasonably well the geometrical aspects of the structure of interest. However, although patch values were chosen to be as accurate as possible, differences from real anatomy may still appear. We thus implemented a SimGAN refinement, similar to , that makes use of simulated and unsupervised learning through a generative adversarial network (GAN) and was originally implemented to specifically avoid the introduction of artifacts and preserve annotation information of the refined images. To the best of our knowledge this is the first time a SimGAN approach is applied to a medical image problem.
A refiner, R (purple box in Fig. 1), is trained to create a refined patch that includes the main aspects of real structures, while at the same time a discriminator, D (orange box in Fig. 1), is used discriminate between real and synthetic patches. This starts a minimax game between R and D (with the weights of two networks updated alternately) that continues until D is not any more capable to distinguish between the two different domains.
The network R consists of a pixel-to-pixel fully convolutional neural network with ResNet blocks, while 4 convolutional layers separated by max pooling and a last layer that outputs the probability of the patch of being a refined image are used to build D.
For D, a cross-entropy loss for a two class classification problem is used, while, as in , the cost function of R is given by:
where is the i-th synthetic training patch, Y is the set of real patches,and are the function parameters. In this work, has been empirically set to 0.01.
In function 2, adds realism to the synthetic images forcing D to fail classifying the refined images as synthetic, whereas is a self-regularization loss that is used to minimize per-pixel difference between a feature transform of the synthetic and refined images and thus preserve the annotation information of the MBG and avoid that, for example, the structure change shape or orientation. More details about the cost functions can be found in .
As in , the receptive field of D is limited to local regions, so that multiple local adversarial losses per patch are considered. Moreover, to improve the stability of the network, a mini-batch of refined images (randomly selected from a buffer of refined images generated on previous iterations) are included into the training batch.
For the implementation of this network, R is first pre-trained on synthetic images for 1,000 steps. Then, D is also pre-trained for 200 stepss, as suggested in , using real patches (extracted as described in section 3.2.1) and refined ones, obtained from the pre-trained R. As a final step, an adversarial training of the SimGAN is executed for 10,000 steps. All networks are trained with a constant 0.001 learning rate and 256 batch size.
For the training of the airway/vessel SimGANs, a total of 165,640 in-vivo airway patches and 2,320,000 in-vivo vessel patches were extracted from 30 clinical cases of the COPDGene study (phase 1) . In order to get balanced datasets, we generated 200,000 and 2,500,000 airway and vessel synthetic patches, respectively. The refined output from the SimGANs is depicted in Fig. 3e. It is worth noting the resemblance between a simulated airway/vessel (Fig. 3e) and a real airway/vessel extracted from a CT scan corresponding to the same nominal structure size (Fig. 3f).
In-vivo Airway/Vessel Patches Extraction for SimGAN Training
To generate the dataset of real 2D patches for SimGAN training, 3232 images were extracted around airway and vessel candidates from in-vivo CT images. To this end, the lung region was first segmented using the method in  and a probability map of the different structures of interest was extracted using the algorithm proposed in . A thresholding operation (threshold = 0.7) and a binary skeletonization were then used to define the initial candidate locations.
Since the in-vivo patches have to be extracted with the same axial-oriented appearance as in the simulated images, the main axis of the structure of interest is considered. To this end, the skeletonized mask is used as input to the scale-space particles sampling method presented in , that starts from an airway/vessel mask to identify an “oriented” bronchial/vascular tree centerline by means of the second-order local information of the image (Hessian matrix). The centerline is stored in the form of a collection of points (called particles) that contain information about scale, orientation (Hessian eigenvectors and eigenvalues), and central pixel intensity. This approach capitalizes on the multi-scale self-similarity of the considered tree, making it more robust to noise in the smaller structures than standard methods .
3.3 Airway/Vessel Measurement Regression
The airway/vessel measurements were computed by means of two separate CNRs having the same basic design (light blue box in Fig. 1). A 9-layer 2D network, consisting of 7 convolutional layers, 5 with stride one and 2 with stride two, and 2 fully-connected layers was implemented. While the network for airway measurement has two outputs (lumen radius and wall thickness), for vessels a single output (lumen radius) is used.
The main goal of the network is to regress the size of the structure centered in the 3232 pixels patches. An Adam update (=0.9, =0.999, =) with a loss function defined by a combination of absolute RE and precision of the measure, as described in section 3.3.1, was used to train both networks. To improve invariance of the network, we also applied data augmentation that adds random noise (range HU), randomly inverts intensity values inside the patches, and slightly shifts and flips the images.
Loss Function for Airway and Vessel Measurement
One of the main issues of typical approaches for airway and vessel morphology assessment, is that they tend to under- or over-estimate the measurement of small airways and vessels (especially with size at the image resolution level). Therefore, in this paper we use a new loss function, similar to , that combines the loss of the RE over all images, , and the precision of the measure over a set of replicas of the exactly same structure, :
where is the true measure, the predicted size, and determines the weight of with respect to (in this work, has been empirically chosen). is given by:
where N is the total number of original patches, and M is the number of replicas used (here, we use M=25).
Conversely, the precision loss term, , is calculated over the M replicas of the same geometric model (with fixed physical dimensions) to which varying PSFs as well as a different number of airways and vessels with various locations and rotations are added. This way, the network learns to precisely measure the structure of interest regardless of possible confounding factors. The definition of is given by:
Since airway lumen and wall thickness are measured simultaneously, and for airway assessment are given by the sum of the corresponding loss computed independently for the two measures.
Since we noticed that the standard deviation was higher for the assessment of small structures (radius 1.0 mm), we also assigned a higher weight to of these structures to give them more importance during training. Therefore, equation 3 for vessels becomes:
while for airways, it becomes:
where l indicates the airway lumen, wt stands for wall thickness, and
Training Set Definition
Training data for both structures consisted of 100,000 geometric models, each replicated 25 times by varying PSFs as well as adding a different number of airways and vessels at various locations and rotations. On the other hand, the validation set was generated with 1,000,000 patches (40,000 geometric models, each replicated 25 times) for both structures.
3.4 Experimental Analysis
Both synthetic and in-vivo experiments were performed to evaluate the two algorithms for airway and vessel assessment.
For synthetic validation, we first generated a dataset of 200,000 patches with random values that we used to compute the absolute RE across all cases, compare results to state-of-the-art methods - considering airways (vessels) with wall thickness (vessel radius) of 1.0 mm and 0.5 mm (image resolution) -, and compare wall thickness measurements to FWHM and ZCSD. As a metric, the mean RE was computed for three separate groups, generated based on the known wall thickness value (WT 0.7 mm, 0.7 mm WT 1.5 mm, WT 1.5 mm).
As a second synthetic validation, we computed the accuracy of the algorithm by calculating the RE obtained comparing the CNR measurements and the geometrical model ground truth when varying lumen (wall thickness = 1.2 mm) and wall thickness (lumen size = 2.0 mm) for airways, and the radius for vessels. To this end, we computed the mean RE across 100 patches generated for each size value.
In an attempt to demonstrate the reliability of the method regardless of the presence of noise and smoothness, as a further experiment we generated a dataset by first varying the level of noise ( HU, mm) and generating 100 synthetic patches for each noise value. Then, we created a second dataset by fixing HU and changing the applied smoothness ( mm) to generate 100 synthetic images per smoothing level. We finally computed the mean RE across the 100 patches for each level of noise and smoothness.
For airways, we first fixed the wall thickness at 1.5 mm and used three different airway lumens (small: 0.5 mm; medium: 2.5 mm; large: 4.5 mm), and then we fixed the lumen at 1.5 mm and used three wall thickness values (small: 0.5 mm; medium: 1.2 mm; large: 2.0 mm). For vessels, we fixed three different radius values (small: 0.5 mm, medium: 2.0 mm, and large: 3.5 mm).
In order to compare the proposed algorithm to state-of-the-art methods, we also validated the performance of the algorithm a CT airway phantom consisting of 8 tubes with different wall thickness and lumen diameter. The tubes were constructed using Nylon66 and were inserted into polystyrene, in an attempt to simulate the lung parenchyma.
The CT image of the phantom was taken using a GE Siemens Sensation 64 CT scanner, with a field of view (FOV) of 40 cm and reconstructed with a standard reconstruction kernel to obtain non-overlapping, 0.6 collimation images. An example of a CT slice of the phantom with wall thickness and lumen values of the eight tubes (as measured by a caliper) are presented in Fig. 4.
As additional experiments, since an accurate and reliable in-vivo ground-truth is very complicated to obtain, we performed two indirect validations on in-vivo clinical cases.
First, we evaluated the reliability of the proposed method using 50 inspiratory thoracic cases from the COPDGene Phase 2 study cohort , for which 6 scans per subject acquired with a different dose and reconstructed with varying parameters were available. All scans for a subject were acquired using the same scanner, just varying the dose or using a different reconstruction approach. Of the 50 cases, 25 were taken using a GE scanner and 25 using a Siemens scanner. Each subject was scanned with two high dose reconstructions (standard, vs. sharp kernel), two field of views (bigger FOV vs. smaller FOV, both in standard high dose), and low dose with two reconstructions (standard and iterative). To make sure that patches extracted from one subject are the same for all scans, we computed the scale-space particles from the high dose image reconstructed with standard kernel and big FOV (STD), and we then registered the created particles to the other scans, using the advanced normalization tools (ANTs) technique described in .
For the statistical analysis, we compared the CNR measurements of each scan to those obtained for the reference image (STD) by means of intra-class correlation coefficient (ICC) and concordance correlation coefficient (CCC) as well as analysis of the distribution of the difference with Bland-Altman (BA) analysis, box-plots, and violin-plots.
Finally, as a second in-vivo validation, we computed a physiological evaluation of the measurement. For the bronchial tree, we compared the correlation between Pi10 and the FEV1% predicted using our approach and the 3D airway measurement software package Pulmonary Workstation (VIDA Diagnostics, Inc., Iowa City, IA, USA) on 3,038 clinical cases from the COPDGene Phase 1 study . 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.
Linear models were created to look at the association between our measurements and a metric of functional small airway disease (fSAD) using the PRM method, a non-invasive technique that measures lung density during inhalation and exhalation, processes the resulting images, and classifies each point in three-dimensional space as normal lung parenchyma, functional small airway disease, or emphysema .
For vessel measurement, we analyzed how total blood volume (TBV), blood volume of vessels of less than (BV5), and blood volume of vessels of less than (BV10) correlate to DLCO adjusted by geographical altitude, based on the study in . For this experiment, we compared the correlation obtained using the radius measure provided by our technique and the scale computed by the particle method  on 1,958 clinical cases from the COPDGene Phase 2 study .
4.1 Synthetic Validation
|WT 0.7 mm||0.7 mm WT 1.5 mm||WT 1.5 mm|
The mean absolute RE obtained for wall thickness, airway lumen, and vessel radius across the 200,000 generated patches was 4.9%, 2.02%, and 2.25%, respectively. When considering airways with a wall thickness of 1.0 mm, a mean absolute RE of 6.3% is obtained, while for wall thickness the image resolution level (0.5 mm) the mean absolute RE is 13.09%.
These REs are significantly lower than those previously reported in the literature for structures of similar sizes [29, 32]. Regarding vessels, a mean absolute RE of 6.09% is obtained for structures of 1.0 mm, while considering vessels with a radius of 0.5 mm the mean absolute RE was 11.3%.
In Tab. 2 the mean RE grouped by wall thickness size (0.7 mm, mm, and 1.5 mm) using CNR in comparison to ZCSD and FWHM on 200,000 testing patches is presented. As shown, while traditional methods yield a very high RE, especially for small airways, a small RE and a very high accuracy is obtained with the proposed method.
In Fig. 5(a-c) the tendency of the RE obtained on synthetic patches when varying airway lumen (Fig. 5a), airway wall thickness (Fig. 5b), and the vessel lumen (Fig. 5c) is presented. As expected, a small RE is obtained for big airways, while a tendency to under-estimate the measure (although always below a 10%) appears for small airways (lumen mm).
For wall thickness results (Fig. 5b), at sub-pixels levels (wall thickness mm) a significant under-estimation error is obtained, while for thicker walls (wall thickness mm) the network tends to over-estimate the measurement. Conversely, for vessels (Fig. 5c), the RE obtained is close to zero, with just a small tendency to under-estimation for vessels of 0.5 mm.
Fig. 6 presents the results obtained when varying the noise and smoothing level applied to the generated patch and using three fixed values of airway lumen size (0.5, 2.5, and 4.5 mm), wall thickness (0.5, 1.2, and 2.5 mm), and vessel lumen radius (0.5, 2.0, and 3.5 mm). As shown, for all structures a stable RE across the varying levels of noise and smoothness is obtained. The CNR yields a very high accuracy (RE ) for medium and large structures, while a slightly higher RE with a bigger standard deviation is obtained when airway lumen, airway wall thickness, and vessel lumen size are at the image resolution of 0.5 mm.
For the three structures, CNR yields a stable RE when varying noise and smoothing level with only a small bias introduced determined by a little under-estimation of the small structures, as expected. For airways with a small wall thickness of 0.5 mm, if the level of smoothing is below 0.6 mm a very small RE is obtained, while if higher levels of smoothing are introduced to the patches the RE slightly increases.
One of the key aspects of the proposed algorithm is the usage of an adversarial mechanism that allows a refinement of synthetic patches into more realistic images.
We validate the output of the SimGAN by testing the performance of the trained discriminator to distinguish between real and synthetic patches. To this end, we used the MBG to generate 300,000 new synthetic airway patches and vessel patches with varying parameters and we extracted 31,316 real airway images and 155,000 real vessel patches (randomly selected) from 10 subjects of the COPDGene Phase 2 study that were not used for training the discriminator.
When passed to the discriminator, 98.772% of synthetic airway patches and 98.983% synthetic vessel patches were classified as coming from the real domain, indicating the fidelity of the SimGAN results in terms of discrimination. As for real patches, 99.901% airways and 99.985% vessels were correctly classified by the discriminator.
4.2 Phantom Evaluation
In Tab. 3 the RE obtained when measuring the wall thickness of the eight tubes of the phantom using both the proposed method (CNR) and traditional techniques is presented. Even though the lumen size is not provided by traditional methods, for completness the RE obtained when measuring the airway lumen with our CNR is also reported in Tab. 3.
In general, the proposed CNR yields the lowest RE for all tubes with the exception of tube C that is best measured by FWHM. Also, an important aspect to notice is the small RE obtained for all tubes when measuring the lumen radius with the proposed technique.
|Tube||CNR WT||FWHM WT||ZCSD WT||CNR Lumen|
4.3 In-vivo Indirect Evaluation
For in-vivo evaluation, airway and vessel morphology has to be computed from CT images. In Fig. 7, a scheme of the processing pipeline to obtain bronchial and vascular measurement from the segmented tree is presented. For each particle, the corresponding structure size is computed.
Consistency to CT acquisition parameters
The CCC values obtained measuring airways and vessels when varying a single CT parameter (kernel, FOV, dose, reconstruction) in comparison to the corresponding reference image (STD) are presented in Tab. 4. In general, very high results are obtained, although, as expected, the measurement of wall thickness is the most complicated one.
In accordance to these results, the ICC obtained for all classes in comparison to the reference images was 97.5%. The box-plots comparing the different parameter variations to the STDs are presented in Fig. 8. For completeness, bias and limits of agreement of the Bland-Altman (BA) analysis for all structures in comparison to measurements obtained on the reference image (STD) are shown in Tab. 5.
|Wall Thickness (mm)||Airway Lumen (mm)||Vessel Radius (mm)|
Fig. 9 shows the violin-plot of the measurements made on sharp images in comparison to the reference scan taken with GE or Siemens scanners. For this analysis, only intra-parenchymal airways (with sizes in the range used for training the network) were considered. As shown, the scanner type is not affecting the measure, as similar differences (mean close to zero) between the images taken with standard and sharp kernels are obtained with the two scanners for all structures.
In order to further demonstrate the reliability of the proposed technique, we also considered the difference in the measurement of airway lumen, airway wall thickness, and vessel radius when varying the scanner parameters in comparison to those obtained on STD, sub-divided in four quartiles based on the size of structure of interest. The results are shown in Fig. 10, while a paired Tukey’s t-test analysis of the results showed no significant difference (p0.001) between the groups for all variations.
FEV1% in correlation to Pi10
As shown in Tab. 6a, the Pearson’s correlation coefficients between Pi10 and FEV1%pred for the presented method (CNR) and VIDA Diagnostics in airway patches extracted from clinical CTs were -0.51 (95%CI -0.53, -0.48) and -0.33 (-0.36, -0.29), respectively. The correlation between Pi10 and fSAD was 0.401 (CI: 0.38, 0.439) for our approach and 0.0862 (CI: 0.051, 0.1213) for VIDA. In a multivariate analysis, the association of our method with fSAD was positive (beta=15.93 p<0.001), while the measurements obtained with VIDA and fSAD had a weaker and negative association (beta=-11.28, p0.001).
DLCO in correlation to TBV, BV5, and BV10
Tab. 6b shows the Pearson’s correlation coefficients obtained between DLCO adjusted for site altitude and TBV, BV5 and BV10 when measuring vessels with the proposed method (CNR) and the particle’s scale. Correlations to DLCO of 0.45 (95%CI 0.41, 0.48), 0.45 (95%CI 0.42, 0.49), and 0.47 (95%CI 0.44, 0.50) were obtained with CNR for TBV, BV5, and BV10, respectively, compared to 0.34 (95%CI 0.39, 0.37), 0.43 (95%CI 0.39, 0.47), and 0.40 (95%CI 0.37, 0.44) given by the scale measurements.
In this paper, we presented a novel method for automatic morphology assessment of airways and vessels from CT scans. The use of a neural network in combination with a generative model refined by SimGAN and the customized loss function represent the innovative aspects of this work.
One of the fundamental limitations for the development of subvoxel sizing methods is the ability to generate realistic and high-resolution representations of 3D bronchial and vascular trees. Realistic anthropomorphic phantoms are complex to generate as several parameters need to be adjusted in order to properly define the topology of the structures of interest and their relationship in the 3D space . Moreover, these methods are not ideal for morphology assessment, as an exact size of the structure is difficult to obtain. An important contribution of this paper is given by the introduction of a 2D generative model, which does not require parameter adjustments, and allows for the generation of patches at will with the exact physical dimensions known a priori.
Results obtained with synthetic patches showed a low absolute RE across all measurements. Although a direct comparison is not possible, considering the absolute RE for airways of 1.0 mm the presented CNR for wall thickness assessment outperforms (absolute RE 6%) the method proposed in [22, 23], where the wall thickness measured on plastic tubes of 1.0 mm yield to an absolute RE of around 10%. Also, a comparison to two traditional algorithms for airway measurement shows that our method improves the state-of-the-art, especially for small and complex airways.
One possible explanation to the bad performance of traditional methods is that they are based on sub-voxel detection of an edge, which has been already described in the literature as a tough problem [29, 32]. This is intimately related to the Nyquist limit theorem and the inability to accurately resolve the size of pulse function whose size is closed to the scanner resolution (i.e. the sampling period of the signal). The problem is further complicated by the variation of the spectral response imposed by different reconstruction kernels and the variation of noise due to different radiation doses.
A validation with structures of different sizes, levels of noise, and smoothing showed that the proposed CNR is not affected by those components, and only sizes below image resolution may determine a small increase of the prediction error. This is also confirmed by the results obtained from clinical cases when varying scanner protocols. Although a direct validation is not possible, due to the complexity to obtain a precise ground truth, the CNR provides very similar results over all scanner variations, demonstrating the high precision of CNR for structure measurement.
As a qualitative example, Fig. 11 shows the enlarged part of the CT image of four subjects, each of which taken with different scan parameters (STD, lower FOV, LD, and ITER), overlaying the vessel segmentation scaled based on the radius provided by CNR. To this end, we first computed the vessel particles segmentation using  and for each particle, we generated a cylindrical stencil with the radius scaled according to an equation that relates the CNR vessel size and the CT’s PSF to the actual vessel radius. From this example, the accuracy and precision of the presented algorithm regardless of the scan protocol can be appreciated.
Considering the results on synthetic patches, an increased RE was expected for small wall thicknesses and high smoothing levels (Fig. 6b), due to a blurring effect that may confound the network. Furthermore, a higher RE was also expected for wall thicknesses bigger than 2.0 mm (Fig. 5b), as the network was trained with wall values lower than 1.7 mm, in accordance to the physiological values presented in . Indeed, we would have expected an even higher error, furtherly demonstrating the reliability of the proposed CNR even for structures it was not trained with. An important aspect to take into account is that airways with wall thickness bigger than 2.0 mm are very unlikely in humans and would be an indicator of strong pathological conditions, which were not aimed in this work.
The results obtained from the phantom validation are promising, showing that the proposed technique properly measures even small and thin airways, as in case of tube D. Although a variance is present in the RE obtained for the wall thickness measurement of all tubes, this variance is smaller than the one obtained using traditional methods that for some tubes seem to get really confounded. In addition to showing that the proposed algorithm outperforms traditional methods, these results also confirm the reliability of the proposed CNR in accurately measuring airway wall thickness and lumen size regardless of the starting conditions.
The indirect validation on clinical cases showed very encouraging results, indicating that the CNR can be used to accurately assess the bronchial and vascular system morphology regardless of different starting conditions and scanner brand. The obtained CCC indicates that a standard low dose (LD) and low dose with iterative reconstruction (ITER) are the most confounding factors in comparison to STD. This was an expected result as the lower resolution of the images causes blurring effects and noise artifacts that may affect the measurement, especially for small and sensible structures such as the wall thickness.
Results from indirect physiological validation showed that CNR is robust and outperforms standard methods. While for bronchial assessment the correlation obtained with CNR is significantly different (p 0.001) from VIDA measurements, better explaining lung function decline in smokers, for the venous system the difference to using the particle scale is not statistically significant. However, there is a clear tendency that shows that the CNR approach improves the association between blood volume and DLCO.
As a final demonstration of the reliability of the proposed technique, Fig. 12 shows the rendering of bronchial (top) and vascular (bottom) trees, of clinical cases from the COPDGene study with different levels of Pi10 (low, medium, high, as measured by the VIDA workstation) and BV5 (low and high, measured by using the particle scale technique presented in ). For airways, the physical size of each segment is given by the lumen size provided by CNR, while the color of each point represents the wall thickness. For vessels, the different colors represent the radius size, and a comparison between CNR and the particle scale method [17, 33] is shown. From this example, we can claim that the measurements provided by CNR respect the natural trend of airway lumen and wall to get smaller when going distally into the lung and that CNR also provides a more realistic measurement than traditional approaches.
Finally, an important aspect to take into account is that, although no simulated branching points have been provided during the training phase, the CNR accurately measure airways and vessels in those complicated points. An example is shown in Fig. 13, where two enlarged parts of a CT showing a branching point of a vessel (left) and an airway (right) overlaying the stenciled segmentation, scaled by the CNR measurements, are presented.
In this paper, we present an innovative airway and vessel generative model to train a CNR for bronchial and vascular morphological assessment that represents a novel contribution to an unsolved problem that is limiting the application of image-based techniques for the phenotyping of airway and vessels.
Up-to-date not much work has been proposed in the literature for the assessment of airway and vessel morphology from CT images. The results obtained with the technique proposed in this work on synthetic images, on the airway phantom as well as by indirect in-vivo validation, showed that our method accurately and precisely measures the size of the different structures and may potentially be used for future early diagnosis of lung disorders and to study the level of resistance and obstruction in small airways and vessels.
An important contribution of this paper is given by the introduction of a 2D generative model that, in combination with a generative adversarial network, creates axial-reformatted synthetic structures at will with known physical dimension to train a CNR. The same concept might be extended to other complex problems for which the biological and image-based knowledge is given, such as pulmonary fissures or cerebral vessels.
For future work, we plan to improve the generation of the synthetic model by reducing the level of approximation of the PSF and additive noise. Also, a new adversarial method that may help neural networks to be more domain-invariant when using synthetic data has been recently proposed . We are planning on implement this idea and compare results to those obtained in the present paper. Finally, a method to properly validate image-to-image filtering with GANs will be also investigated.
- (2016) Tensorflow: large-scale machine learning on heterogeneous distributed systems. arXiv preprint arXiv:1603.04467. Cited by: §3.3.
- (2009) Advanced normalization tools (ants). Insight J 2, pp. 1–35. Cited by: §3.4.
- (1998) Airway wall thickness in patients with near fatal asthma and control groups: assessment with high resolution computed tomographic scanning. Thorax 53 (4), pp. 248–253. Cited by: §1.
- (2007) Impaired flow-mediated dilation is associated with low pulmonary function and emphysema in ex-smokers: the emphysema and cancer action project (emcap) study. AJRCCM 176 (12), pp. 1200–1207. Cited by: §1, §1.
- (2018) Small airway segmentation in thoracic computed tomography scans: a machine learning approach. Physics in Medicine & Biology 63 (15), pp. 155024. Cited by: §1.
- (2019) Airway wall thickening on CT: Relation to smoking status and severity of COPD. Respiratory medicine 146, pp. 36–41. Cited by: §1.
- (2015) Keras. Note: \urlhttps://keras.io Cited by: §3.3.
- (2010) Measuring small airways in transverse ct images: correction for partial volume averaging and airway tilt. Academic radiology 17 (12), pp. 1525–1534. Cited by: §1.
- (1968) Newer aspects of the pulmonary vasculature in chronic lung disease: a comparative study. Angiology 19 (7), pp. 399–407. Cited by: §1.
- (2013) Computed tomographic measures of pulmonary vascular morphology in smokers and their clinical implications. AJRCCM 188 (2), pp. 231–239. Cited by: §3.4.
- (2012) Computed tomography–based biomarker provides unique signature for diagnosis of copd phenotypes and disease progression. Nature medicine 18 (11), pp. 1711. Cited by: §2, §3.4.
- (2016) Domain-adversarial training of neural networks. The J of Machine Learning Research 17 (1), pp. 2096–2030. Cited by: §6.
- (2015) Chronic obstructive pulmonary disease: ct quantification of airway dimensions, numbers of airways to measure, and effect of bronchodilation. Radiology 277 (3), pp. 853–862. Cited by: §1.
- (2013) Small airway obstruction in COPD: new insights based on micro-CT imaging and MRI imaging. Chest 143 (5), pp. 1436–1443. Cited by: §1.
- (1967) Vascular changes in pulmonary emphysema: the radiologic evaluation by selective and peripheral pulmonary wedge angiography. American J of Roentgenology 100 (2), pp. 374–396. Cited by: §1.
- (2016) Automatic synthesis of anthropomorphic pulmonary ct phantoms. PloS one 11 (1), pp. e0146060. Cited by: §5.
- (2009) Sampling and visualizing creases with scale-space particles. IEEE TVCG 15 (6). Cited by: §3.2.1, §3.4, §5, §5.
- (2015) Deep learning. nature 521 (7553), pp. 436. Cited by: §1.
- (2005) Optimal surface segmentation in volumetric images-a graph-theoretic approach. IEEE transactions on pattern analysis and machine intelligence 28 (1), pp. 119–134. Cited by: §1.
- (2012) Optimal graph search based segmentation of airway tree double surfaces across bifurcations. IEEE transactions on medical imaging 32 (3), pp. 493–510. Cited by: §1.
- (2013) Accurate measurement of small airways on low-dose thoracic ct scans in smokers. Chest 143 (5), pp. 1321–1329. Cited by: §1.
- (2000) Computed tomographic measurements of airway dimensions and emphysema in smokers: correlation with lung function. AJRCCM 162 (3), pp. 1102–1108. Cited by: §1, §5.
- (2002) Development and validation of human airway analysis algorithm using multidetector row ct. In Medical Imaging 2002: Physiology and Function from Multidimensional Images, Vol. 4683, pp. 460–469. Cited by: §1, §5.
- (2017) CT Image Enhancement for Feature Detection and Localization. In MICCAI, pp. 224–232. Cited by: §3.2.1.
- (2018) Accurate measurement of airway morphology on chest CT images. In Image Analysis for Moving Organ, Breast, and Thoracic Images, pp. 335–347. Cited by: §2, §3.3.1.
- (2000) Airway wall thickness in asthma assessed by computed tomography: relation to clinical indices. AJRCCM 162 (4), pp. 1518–1523. Cited by: §1.
- (2011) Optimal graph based segmentation using flow lines with application to airway wall segmentation. In Bienn. International Conference on Inf. Processing in Medical Imaging, pp. 49–60. Cited by: §1.
- (2011) Genetic epidemiology of copd (copdgene) study design. COPD 7 (1), pp. 32–43. Cited by: §3.2, §3.4, §3.4, §3.4.
- (1997) Accurate measurement of intrathoracic airways. IEEE TMI 16 (6), pp. 820–827. Cited by: §1, §1, §4.1, §5.
- (2009) Lung extraction, lobe segmentation and hierarchical region assessment for quantitative analysis on high resolution computed tomography images. In MICCAI, pp. 690–698. Cited by: §3.2.1.
- (2014) Comparing algorithms for automated vessel segmentation in computed tomography scans of the lung: the vessel12 study. Medical image analysis 18 (7), pp. 1217–1232. Cited by: §1.
- (2006) Accurate airway wall estimation using phase congruency. In MICCAI, pp. 125–134. Cited by: §1, §4.1, §5.
- (2012) Computational vascular morphometry for the assessment of pulmonary vascular disease based on scale-space particles. In ISBI, pp. 1479–1482. Cited by: §1, §3.2.1, Figure 12, §5.
- (2008) Three-dimensional airway measurements and algorithms. Proceedings of the Am Thoracic Soc 5 (9), pp. 905–909. Cited by: §3.1.
- (2003) Enhanced expression of vascular endothelial growth factor in pulmonary arteries of smokers and patients with moderate chronic obstructive pulmonary disease. AJRCCM 167 (9), pp. 1250–1256. Cited by: §1, §1.
- (1993) Dynamic imaging of the upper airway during respiration in normal subjects. J of Applied Physiology 74 (4), pp. 1504–1514. Cited by: §1.
- (2005) The point spread function of spiral ct. Physics in Medicine & Biology 50 (22), pp. 5307. Cited by: §3.1.
- (2017) Learning from simulated and unsupervised images through adversarial training. In IEEE CVPR, pp. 2107–2116. Cited by: §1, §3.2, §3.2, §3.2, §3.2, §3.2.
- (2005) Intrathoracic airway trees: segmentation and airway morphology analysis from low-dose ct scans. IEEE transactions on medical imaging 24 (12), pp. 1529. Cited by: §1, §2.
- (2017-08) Statistical characterization of noise for spatial standardization of CT scans: Enabling comparison with multiple kernels and doses.. MedIA 40, pp. 44–59. Cited by: §3.1.
- (1965) Morphometry of the human lung. Springer. Cited by: §3.1, §3.1, Table 1, §5.
- (2019) Automatic quantitative analysis of pulmonary vascular morphology in ct images. Medical physics. Cited by: §1.