Identification of COPD From Multi-View Snapshots of 3D Lung Airway Tree via Deep CNN

Chronic obstructive pulmonary disease (COPD) is associated with morphologic abnormalities of airways with various patterns and severities. However, the way of effectively representing these abnormalities is lacking and whether these abnormalities enable to distinguish COPD from healthy controls is unknown. We propose to use deep convolutional neural network (CNN) to assess 3D lung airway tree from the perspective of computer vision, thereby constructing models of identifying COPD. After extracting airway trees from CT images, snapshots of their 3D visualizations are obtained from ventral, dorsal and isometric views. Using snapshots of each view, one deep CNN model is constructed and further optimized by Bayesian optimization algorithm to indentify COPD. The majority voting of three views presents the final prediction. Finally, the class-discriminative localization maps have been drawn to visually explain the CNNs’ decisions. The models trained with single view (ventral, dorsal and isometric) of colorful snapshots present the similar accuracy (ACC) (86.8%, 87.5% and 86.7%) and the model after voting achieves the ACC of 88.2%. The ACC of the final voting model using gray and binary snapshots achieves 88.6% and 86.4%, respectively. Our specially designed CNNs outperform the typical off-the-shelf CNNs and the pre-trained CNNs with fine tuning. The class-discriminative regions of COPD are mainly located at central airways; however, regions in HC are scattering and located at peripheral airways. It is feasible to identify COPD using snapshots of 3D lung airway tree extracted from CT images via deep CNN. The CNNs can represent the abnormalities of airway tree in COPD and make accurate CT-based diagnosis of COPD.


I. INTRODUCTION
Chronic obstructive pulmonary disease (COPD) is a common respiratory disease with a trend of growing prevalence globally [1]. The estimated global prevalence of COPD in 2015 is about 174 million [2]. Currently, COPD has been the third The associate editor coordinating the review of this manuscript and approving it for publication was Chao Zuo . leading cause of death (3 million every year) and a major cause of disability [3]. Extensive researches are urgently demanded to improve COPD outcomes through early identification and exploration of biomarkers for phenotype diagnosis and personalized therapy [4], [5].
COPD is characterized by persistent and incompletely reversible airflow limitation and gas trapping caused by multiple pathological alternations including emphysematous VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ lung tissue destruction, gross airway disease and functional small airway disease [3]. FEV 1 (forced expiratory volume in 1 second), measured using spirometry, has been the standard and dominated indictor in the diagnosis and therapies of COPD due to the straightforward and inexpensive characteristics [6]. However, FEV 1 does not allow for differentiation between subtypes and is weakly correlated with the exacerbations and mortality of COPD. In 2017, the GOLD report has clearly stated that to FEV 1 should not be used for assessing the degree of reversibility of airflow limitation or making therapeutic decisions [7]. Low dose and high resolution computer tomography (CT) and subsequently quantitative analysis in COPD have emerged as a new approach to visualize and quantitatively measure regional airflow limitation, gas trapping, emphysema, and airway remodeling [8]- [12]. Nambu et al. [13] found CT measures have significant correlations with COPD physiological indices, where CT measures include of % low attenuation areas (%LAA), square root of the wall area at an internal airway perimeter of 10 mm (Pi10), and airway wall thickness, luminal diameter and airway wall area percent (WA%) in the segmental, subsegmental and subsubsegmental bronchi. Emphysema in COPD can be divided into centrilobular, panlobular, and paraseptal types according to its appearance [14]. Using expiration and inspiration CT scans and a voxel-wise image analysis technique, the Parametric Response Map (PRM) of functional small airways disease (fSAD) and emphysema have been found be able to differentiate between specific COPD phenotypes accurately [5].
Airway trees can be automatically segmented from CT images via advanced imaging processing methods [15]. However, no consensus on the most appropriate methods for quantifying the airways has been achieved though many measures have been proposed. The abnormalities in COPD airway can be thicker wall thickness, smaller airway intraluminal area (Ai) of third-to sixth-generation bronchi [16], more branch variation [17], smaller child-to-parent diameter ratios, larger length-to-diameter ratios, and smaller fractal dimensions [18]. They can also be Pi10 [13], tracheal section [19], CT total airway count [20], and tracheal shape [21].
Visual assessment of CT image-based airways of COPD can provide complementary information to quantitative measurements [6]. Kim et al. [22] have demonstrated that visual assessment by radiologists can lead to the reproducible, physiologically substantial knowledge on airways of COPD. The Fleischner society has also recommended using visual assessment of CT airway to define the subtypes of COPD [14]. However, visual assessment has been reported to be of the relatively low interreader agreement, highly subjective and poorly reproducible, though the use of standard images as references ameliorates the situation [22].
Deep learning, deep convolutional neural network (CNN) in particular, as a kind of computer vision method, has been developed to understand the objects or images. Deep CNN can learn representative features from the data with multiple levels of abstraction and thus the design, extraction and selection of handcrafted features are unnecessary [23]. Via easy practicability and excellent performance, deep CNN has made great breakthroughs in processing both natural and medical images [24]- [29].
Studies on CT imaging analysis of COPD by deep learning or machine learning are not many. Using four canonical CT slices at predefined anatomic landmarks as inputs, González et al. [30] trained deep CNN models with an accuracy (ACC) of 0.773 for the detection of COPD in the COPDGene testing cohort (n = 1000). Meanwhile, the trained CNN model can also accurately stage COPD within one stage (74.95% for COPDGene and 74.6% for ECLIPSE). The prediction ACC of acute respiratory disease events reaches 0.604 and 0.606 in the COPDGene testing cohort (n = 1000) and the ECLIPSE cohort (n = 1062), respectively. Cheplygina et al. [31] utilized instance-transfer leaning to classify COPD from different centers, scanners, or subject distributions. The only study using airway tree to classify COPD is done by Petersen et al. [32], where geometric tree kernels have been employed to extract the features of airways, support vector machine (SVM) is used further and the highest ACC of 64.9% is achieved for 1966 individuals including 893 COPD. To the best of our knowledge, our work is the first to use deep CNN and CT based airway trees for identification of COPD.
In this paper, we propose to use deep CNN to automatically assess 3D lung airway tree from the perspective of computer vision and construct models of identifying COPD. It is noted that, beside the abnormalities of lung airway tree, COPD will also lead to other abnormalities of lung parenchyma, blood vessel, even extra-pulmonary organs (e.g., osteoporosis). We investigate the feasibility of identifying COPD by deep CNN models which only use the abnormalities of lung airway tree with the aim to emphasize both the importance of abnormalities and effectiveness of specifically designed deep CNN models. If one aims to build up one model for accurately identify COPD from CT images, comprehensive features related to lung airway tree, lung parenchyma, blood vessel, even extra-pulmonary organs should be considered.
The main contributions of this study are threefold. First, inspired by the face recognition method [33], the 2D snapshots of 3D lung airway tree, rather than 3D airways themselves or original CT slices, are employed. To utilize 3D airways directly is unfeasible due to the limited number of training samples in the current study and also has the higher computational complexity. The role of the abnormalities of airway tree for COPD identification cannot be clarified while using the original CT slices because these slices include all chest components (e.g., the chest wall, bone, lung parenchyma, heart, pulmonary vessels, airway tree, pulmonary interstitium, etc.). Second, Bayesian optimization algorithm (BOA) is utilized to obtain the effective convolution network structure and hyper-parameters [34]. The specially designed and BOA optimized CNN models achieve the highest identification accuracy of 88.6% and outperform the typical off-the-shelf CNN models, the CNN model using the original CT slices [30], and other state-ofthe-art models using transfer learning [31] and geometric tree kernels [32]. Third, Gradient-weighted Class Activation Mapping (Grad-CAM) is used to highlight the important regions in the snapshots for predicting healthy control (HC) and COPD airways [35]. The class-discriminative regions are found to be located at the central airways in COPD, but to be scattering and with many edges of peripheral airways in HC. This study has demonstrated that the CNNs can represent the abnormalities of airway tree in COPD and make accurate CT-based diagnosis of COPD.

A. PARTICIPANTS AND DATASET
Totally 280 participants (190 COPD patients, 90 healthy controls (HC)) are enlisted from subjects who underwent CT scan in Central Hospital Affiliated to Shenyang Medical College during the period of 2016-2018. COPD is diagnosed by clinicians in respiratory department according to the spirometry measurements of pulmonary function tests. If the subjects fulfill the Global Initiative for Chronic Obstructive Lung Disease (GOLD) criteria, i.e., FEV 1 /FVC is less than 0.7 (FEV 1 is the postbronchodilator forced expiratory volume in 1 second and FVC is the forced vital capacity).
The characteristics of the participants and CT image acquisition parameters are given in Table 1. By the Chi-square test, there is no significant difference for the gender between COPD (male/female: 110/80) and HC (male/female: 54/36) (p = 0.738). According to the two-sample t-test, there is no significant difference for the age (years) between SD and HC groups (mean ± S.D.: 76.6 ± 9.2 vs. 68.7 ± 9.3, p = 0.999).
As listed in Table 1, the tube voltage is set as 120 kV and the slice thickness is 1.5 mm for all the acquisitions of CT images. Two CT scanners (Philips iCT 256 / Neusoft NeuViz 128) are utilized. The pixel size, the X-ray tube current, and the exposure vary across the subjects. The study has been approved by the Medical Ethics Committee of Central Hospital Affiliated to Shenyang Medical College and is in accordance with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. All participants signed a written informed consent in accordance with the Declaration of Helsinki (2000).

B. OVERVIEW OF STUDY PROCEDURES
As shown in Figure 1, there are two main steps in our study. In the first step, the CT images are input into the medical imaging process software of Mimics (Materialise Corp, Belgium) and the airway tree is extracted using the tool ''Deep airway segmentation'' (Pulmonary Module) embed- ded in Mimics [36]. No any further user interaction is used. Colorful snapshots of 3D visualization of surface rendering of the extracted airway trees are obtained from ventral, dorsal and isometric views. The size of all snapshots is set as 224 × 224 by a fix-sized box.
In the second step, the multi-view CNN models are built to classify the participants into COPD and HC. Three schemes are proposed and implemented for the classification. In Scheme I. Colorful snapshots, three deep CNN models are constructed using colorful snapshots from ventral, dorsal and isometric views. According different views, three models are named as CNN-CV, CNN-CD and CNN-CI, where ''C'' indicates the colorful snapshots. Majority voting is used to fuse the predictions of three models and yield the final classification.
Scheme II. Gray snapshots is designed to examine whether the rendering pseudo-color of airway tree influences the performance of identification of COPD. Through the gray-level transformation, we change the input snapshots into gray ones. Similarly three models of CNN-GV, CNN-GD and CNN-GI are constructed, where ''G'' means the gray snapshots.
Scheme III. Binary snapshots is to study if the gray contrast has effect on the classification performance. Actually the gray contrast distributions of airway tree characterize the smoothness of airway tree. The gray levels will be small at the smooth places of surface of airway and the drastic changes of gray distribution will appear if the surface of airway is irregular. Binarization of gray snapshots yields binary ones which are employed to build three models of CNN-BV, CNN-BD and CNN-BI. Herein, ''B'' represents the binary snapshots.

C. CNN MODEL DESIGN AND PARAMETERS FOR OPTIMIZATION
We have designed CNN models by the widely used strategy of ''stage-wise building-block'' [37]. As shown in Figure 2, the CNN model consists of several ''RANGE''s, each RANGE has some convolutional layers, and each convolutional layer is followed by a batch normalization layer and a leaky rectified linear activation layer (Leaky ReLU). One MaxPooling layer follows each RANGE. After the final Max-Pooling layer, one AvgPooling, two fully connected layers, a dropout layer and Softmax layer are linked successively.
Bayesian optimization algorithm (BOA) is employed to find the effective convolution network structure and hyperparameters [34]. As give in Table 2, the encoder depth, network depth and filter size are optimized within the ranges of [1,3], [1,4] and [2,5], respectively. The encoder depth is defined as the number of RANGEs and network depth is the number of convolutional layers in a RANGE. The filter size indicates the height and width of the kernel and it can be 2 × 2, 3 × 3, 4 × 4, or 5 × 5. Meanwhile, the leaky factor, the dropout factor, the initial learning rate, the momentum and regularization factor are optimized by BOA within the range of values listed in Table 2. The leaky factor is a constant of the leaky-ReLU activation function and the dropout factor is the probability of drop neuron nodes. Correspondingly, the momentum is defined as the contribution of the gradient step from the previous iteration to the current iteration of training and regularization factor is the factor for the L2 regularizer. The number of filters in RANGE1 and RANGE2 is determined as 64 and 128 according to experiences.
BOA is a sequential model-based approach to find a global maximizer of an unknown objective function  Q (f |x i , y i for i = 1, . . . , t); (2) Find the new point x that maximizes the acquisition function a (x). Here, the expected improvement (EI) is used as a (x) and defined as here x best is the location of the lowest posterior mean and µ Q (x best ) is the lowest value of the posterior mean. For detailed information on definitions of parameter and implementation of BOA, one can refer to the publications [34], [38]. In current study, BOA is realized using the function embedded in Deep learning Toolbox in Matlab 2018a.

D. CNN TRAINING AND CROSS-VALIDATION EXPERIMENTS
The performance of CNN models is evaluated by using a stratified ten-fold cross-validation. Specifically, snapshots of each single view are divided into the training, validation and testing datasets with the ratio of 8:1:1. In each fold of crossvalidation, the testing and validation datasets are completely different. The method of focal loss is used to tackle the class imbalance problem, where the weighting factor α is set as 0.32 (90/280) and the focusing parameter γ is set as 2.0 [39].
The snapshots in the training and validation datasets have been augmented through the random rotation, reflection, scale, shear and translation with a defined range, where the reflection, scale, shear and translation have been done in both horizontal and vertical directions. In summary, the number of training and validation data is increased to 10 times.
To alleviate the overfitting, several ways such as the data augmentation, dropout, L2 regularization and early stopping have been adopted. The early stopping is triggered if the accuracy of the validation dataset does not increase within 5 iterations.
The measures of accuracy (ACC), confusion matrix, receiver operating characteristic (ROC) curve, area under ROC curve (AUC) are calculated to evaluate the classification performance. Finally, a parameter named the ratio of fulfill area (RFA) of airway tree is defined in order to investigate whether the volume or surface area of airway tree has a great influence on the identification of COPD. It is calculated as the proportion of the airway tree area in one snapshot. The two-sample t-test is employed to examine whether there is a significant difference for RFA between COPD and HC.

E. COMPARATIVE CNN MODELS AND EXPERIMENTS
To compare the classification performance of our specifically designed and BOA optimized CNNs with that of other available CNNs, three typical off-the-shelf CNN models including AlexNet [40], VGG-16 [41] and Inception-V1 [42] are employed using the gray snapshots of dorsal view. Inception-V1 module consists of three convolutional layers (1 × 1 convolutions, 3 × 3 convolutions and 5 × 5 convolutions) and one MaxPooling layer, where 1 × 1 convolutions are used to compute reductions before the expensive 3 × 3 and 5 × 5 convolutions.
We have also tried the way of using the pre-trained CNN with fine tuning. For the pre-trained AlexNet and VGG-16 VOLUME 8, 2020 with fine tuning, we replaced the final three layers of the pre-trained CNNs by three new layers of a fully connected layer, a softmax layer and a classification output layer. Using the snapshots of 3D lung airway tree and the corresponding labels, the new layers are trained. To accelerate the learning of the new fully connected layer, its learning rate factors of weight and bias are set as 50. Similarly, for the pre-trained GoogLeNet with fine tuning, the final three layers are replaced and trained by our new dataset. Since the depth of GoogLeNet is larger than that of AlexNet and VGG-16, the weights in all the remaining 110 layers of GoogLeNet are frozen (i.e., the learning rate is set as zero) while training the new fully connected layer. BOA is also be utilized for the fine tuning of pre-trained CNNs. Three optimized parameters are Dropout factor, Initial learning rate and Momentum.
Moreover, we have implemented a 3D CNN model for comparison. In this specifically designed CNN model, the input is the matrix of 64 × 64 × 64 and there are 6 convolutional layers (C1-C6), 3 dilation convolutional layers (DC1-DC3), 1 global average pooling layer (GAP), and 2 fully connected layers (FC1-FC2). A batch normalization layer follows each convolutional layer and the activation function is ReLU. Dropout is used between FC1 and FC2 with a dropout factor of 0.7.
The experiments were implemented in Matlab 2018a (Deep Learning Toolbox) under a Windows 7 operating system on a workstation with CPU Intel Xenon E5-2650 v4 @ 2.20 GHz, GPU NVIDIA Quadro P4000 and 64G RAM. Recently, Deep Learning Toolbox has been one alternative framework for designing, optimizing and implementing deep CNNs. It has the great features of the training flexibility, data preprocessing, training acceleration, and visualization. One also can easily exchange models between this toolbox and other frameworks including TensorFlow, PyTorch, TensorFlow-Keras and Caffe.

F. CLASS-DISCRIMINATIVE LOCALIZATION MAP
With the aim of visually explaining the decisions made by deep CNN model, Gradient-weighted Class Activation Mapping (Grad-CAM) is used to highlight the important regions in the snapshots for predicting HC and COPD airways [35]. The class-discriminative localization (CDL) map is calculated by Grad-CAM here A k is the kth feature map in the last convolutional layer, K is the number of filters or feature maps (K = 128 in current study) and α c k represents the importance of A k for a target class c (i.e., COPD or HC). α c k is defined as here y c is the output (before the softmax) of deep CNN for one snapshot from class c, i and j indicate the dimensions of A k and Z = I × J . In our deep CNN models, I = 100 and J = 100. Finally, L c Grad−CAM with the dimension of I × J is interpolated into the one with the same dimension as the input snapshot, i.e., 224 × 224.

A. NETWORK ARCHITECTURE AND PARAMETER OPTIMIZATION
After BOA, the final optimized architecture of the proposed CNN model is shown in Figure 2. The number of kernel (or feature maps) is 64 and 128 in RANGE1 and RANGE2, respectively. The kernel size is 3 × 3 and the batch size is 16. The encoder depth is 2 and the network depth is 4. The leaky factor, the dropout factor, the initial learning rate, the momentum parameters and regularization factor are determined as 0.1355, 0.6114, 0.0007, 0.8659 and 0.0081, respectively.
Since BOA is very time-consuming, it has only been implemented once using snapshots of all three views. For the CNNs with different snapshots and views, the architecture and parameters given as above are fixed.

B. CLASSIFICATION BY COLORFUL SNAPSHOT
The accuracy and loss of the training and validation processes in CNN-CD model (colorful snapshots, the dorsal view) are shown in Figure 3(a) as an example. It can be seen that the training loss of decreases continuously and reaches about 0.2 after 2200 iterations. Meanwhile, the training accuracy increases gradually and reaches above 0.98 after 2200 iterations. With the iteration, the accuracy of validation dataset  can increase to above 90% and the loss of validation dataset can drop to the range of 0.2-0.3. No big difference between the training and validation datasets is observed for both the accuracy and loss. Therefore, the overfitting has been alleviated within an acceptable range.
Using the ventral, isometric and dorsal view of colorful snapshots as inputs, the accuracy of CNN-CV, CNN-CI and CNN-CD achieves 86.8%, 86.7% and 87.5% (Table 3). CNN-CD outperforms CNN-CV and CNN-CI, indicating that snapshots from the dorsal view might present more differences of airway tree appearances between COPD and HC. After majority voting, the accuracy slightly rises to 88.2%.
The ROC and AUC of CNN-CV, CNN-CI and CNN-CD are given in Figure 4(a). As the accuracy, AUC of CNN-CD is the highest one, reaching 0.925. For CNN-CV and CNN-CI, it is 0.925 and 0.919. The confusion matrix of the final voting classifier is given in Figure 5(a).It is found that the sensitivity is as high as 0.963 (183/190). However, the false discovery rate is also as high as 0.289, i.e., 26 HCs among 90 are wrongly predicted as COPD.
For the colorful snapshot, three input channels (Red, Green and Blue) are utilized. Each folder of training will take about 41 minutes.

C. CLASSIFICATION BY GRAY SNAPSHOTS
For Scheme II. Gray snapshots, the accuracy and loss of the training and validation processes of CNN-GD model (the dorsal view) are shown in Figure 3 (Table 3). It is noted that, while using gray snapshot, the ACC of the final voting model increases to 88.6%, indicating that the color rendering adds confusing information for COPD identification. The AUC of CNN-GV, CNN-GI and CNN-GD is 0.906, 0.913 and 0.916 (Figure 4(b)), slightly smaller that of colorful snapshot. According to the confusion matrix ( Figure 5(b)), there are 27 false positives (HC is wrongly predicted as COPD) and 5 false negative (COPD is wrongly predicted as HC). Moreover, compared with Scheme I Colorful snapshots, the training time using gray snapshots is shorter (one folder, about 25 minutes) for there is only one input channel.

D. CLASSIFICATION BY BINARY SNAPSHOTS
In Scheme III. Binary snapshots, the ACC of CNN-BV, CNN-BI and CNN-BD is 86.4%, 83.2% and 87.1% (Table 3), VOLUME 8, 2020  smaller than that of the corresponding models with colorful and gray snapshots. After majority voting, the ACC reaches to 86.4%. It is suggested that abnormal appearances of airway tree in COPD are the major determinant of the identification and the gray contrast distributions characterizing the smoothness of airway tree have certain but limited contributions.
Similarly, the AUC of CNN-BV, CNN-BI and CNN-BD is 0.896, 0.893 and 0.894 (Figure 4(c)), smaller than that of models driven by colorful and gray snapshots. There are 32 false positives and 6 false negatives ( Figure 5(c)).
No significant difference is observed among the accuracy of models using the colorful, gray and binary snapshots. However, a significant AUC is observed for these models (p < 0.05, the colorful snapshots > the gray snapshots > the binary snapshots). The accuracy of model using dorsal view of snapshot is significantly higher than that of ventral view (p < 0.05), but no significant difference is observed in the other two comparisons: dorsal vs. isometric; ventral vs. isometric.

E. RATIO OF FULFILL AREA (RFA)
There is no significant difference of RFA between COPD and HC for any view (the ventral, dorsal or isometric view) ( Figure 6). The p value is 0.058, 0.153 and 0.515. These results indicate that the classification of COPD is independent on RFA. The appearances of airway tree take crucial role for distinguish COPD from HC, rather than the volume or surface area of airway tree.

F. CLASSIFICATION BY OTHER COMPARATIVE CNNs
Using the gray snapshots of airway tree from dorsal view, AlexNet, VGG-16 and Inception-V1 achieve an ACC of 77.9%, 74.6% and 76.79%, respectively, lower than our proposed CNN-GD (87.9%). Our specially designed and BOA optimized CNN models outperform the typical off-the-shelf CNN models. Moreover, our current models also outperform or show comparable performance with the CNN model using the original CT slices (ACC of 77.3%) [30] and other state-of-the-art models using transfer learning (AUC of 0.79-0.90) [31] and geometric tree kernels (ACC of 68.2%) [32].
While adopting dorsal view of gray snapshots, the final accuracy is 83.93% and 85.71% for pre-trained AlexNet and VGG-16, respectively. The accuracy of pre-trained GoogLeNet is 78.57% for the dorsal view of gray snapshots. The possible reasons lie in the current small dataset and the significant difference between snapshots of airway tree and natural images used in the pre-trained CNNs. The designed 3D CNN model has an accuracy of 78.6% for the classification of COPD and HC using the 3D airway trees extracted from CT images. Due to the small dataset, we have augmented the data by 16 times. For the deeper CNNs, no convergence can be obtained.

G. CLASSIFICATION INTERPRETATION AND GRAD-CAM
In form of snapshot of airway tree from dorsal view, some examples of classification are shown in Figure 7. It seems that the correctly predicted COPD subjects (the first row) have thicker trachea and main bronchi and less number of terminal bronchi, compared with HC (the second row). It is noted that CNN models may identify and utilize the pattern of abnormalities of airway tree in COPD in the agnostic way. For the false positives and negatives, it is very difficult to identify and represent the appearance differences even from perspective of computer vision. Figure 8 presents the overlap of airways and the classdiscriminative localization maps for HC and COPD groups. Different views of snapshots correspond to different patterns of CDL maps for both groups. Moreover, CDL maps in COPD show more class-discriminative regions and these regions are mainly located at the central airways including trachea, main bronchi and segmental bronchi. Comparatively, the class-discriminative regions in HC group are scattering and include some edges of peripheral airways. In general, the highlighted class-discriminative regions are the edges or textons with low abstract level because the last convolutional layer in our CNN still has the demission of 100 × 100.

IV. DISCUSSIONS
In the present study, deep CNN models have been successfully built to identify COPD using multi-view snapshots of 3D lung airway tree extracted from CT images. It is indicated that, first, the airway tree in COPD extracted from CT images presents distinct and abnormal morphological appearances from that of HCs. Second, these abnormalities can be represented efficiently by deep CNN from 2D multi-view snapshots of 3D airway tree, rather than the 3D airway tree, and therefore the classification accuracy between COPD and HC can reach an ACC of 88.6%. Third, abnormal appearances of airway tree in COPD are the major determinant of the identification and the local smoothness has certain but limited contributions.

A. PATHOPHYSIOLOGICAL RATIONALE OF AIRWAY ABNORMALITIES IN COPD
Rationale of using CT quantitatively measure airway morphology in COPD root in both aspects of airway innate anatomy and remodeling [43]. McDonough et al. [44] found that small conducting airways (<2.0 mm) in COPD become narrow and disappear before the onset of emphysema. Previous reports have asserted that these abnormalities in small airways will be reflected in large-and intermediate-sized airways [45]- [47]. A reduced tracheal section, total airway count on CT and airway branch variation have been proposed as the independent and innate risk factor of COPD [17], [19], [20], [43]. Our study has provided further evidence to the above rationale. VOLUME 8, 2020

B. CNNs EFFICIENTLY REPRESENT AIRWAY ABNORMALITIES IN COPD
Morphologic abnormality of large and small airways in COPD presents various different patterns and severities. To efficiently represent these abnormalities of the complex 3D structure of airways is a tough task. The difficulty originates from three aspects. First, abnormalities of airways in COPD are multiple dimensions. It is unclear that what and how many measurements would be representative of the entire airways tree [48]. Second, abnormalities of COPD airways vary with anatomy locations. To find the right locations to measure airways has been noted and the spatially matched method has been proposed [49], [50]. Third, quantitative measurements of COPD airway are not standardized. Different CT scanners and acquisition protocols of hospitals further increase the uncertainties and difficulties [8], [11], [51].
Our deep CNN method avoids the above three difficulties and results have shown that it can efficiently represent airway abnormalities in COPD. By the way of computer or machine vision, our deep CNNs lean a kind of effective, highly abstract and agnostic representation of COPD airway abnormalities from the data, instead of extracting ''explicit'' or ''semantic'' handcraft features [23]. González et al. [30] have shown the power of deep CNN in detection of COPD in the COPDGene testing cohort (ACC of 0.773). However, they directly used four canonical CT slices at predefined anatomic landmarks and 16 filters enhanced different structures in chest including the lungs, the chest wall, or the bone structures. No previous study has been done by using deep CNN and CT based airway trees for identification of COPD.

C. MULTI-VIEW SNAPSHOTS VS. SINGLE-VIEW SNAPSHOT
Inspired by the face recognition, we utilized multi-view (ventral, dorsal and isometric) 2D snapshots of 3D airway trees, instead of 3D airways themselves. Multi-view deep leaning has been widely used in face recognition and medical imaging analysis [52]- [54]. Recently, Keceli have utilized the pretrained deep CNNs to extract features from five 2D views of 3D volume for the action recognition [55]. The common ways to transfer 3D to 2D include the combination of axial, coronal and sagittal images, several axial slices, several canonical CT slices [30], [56], [57]. Given our dataset is relatively small, only includes 280 subjects, we adopted 2D snapshots of 3D airway trees from ventral, dorsal and isometric views. From Figure 8, one can see that different view corresponds to different class-discriminative localization map for COPD, indicating that the airway information extracted from snapshots of different views is complementary.
Our fusion strategy of three views is a kind of late or decision-level fusion, which is simpler than the early or data-level fusion and intermediate or feature-level fusion [58]. CNNs with multi-stream architectures (intermediate fusion) have been explored for multi-scale image analysis and multi-view (or 2.5D) classification [25], [54]. The method of fusion in the current study is the majority voting and other fusion methods may include average or even meta-learners (e.g., ELM) [59]. Our results have shown that the fusion can increase the identification accuracy certainly but slightly.

D. COLORFUL, GRAY AND BINARY SNAPSHOTS OF AIRWAY TREE
Colorful snapshots of airway are actually pseudo-color. Pseudo-color generated in the visualization rendering, as the name suggests, does not only contain any useful information, but also adds confusing information for COPD identification. Lesions of various severities have observed in the bronchi of patients with moderate and severe COPD and considered to be related with airway remodeling [60]. We suppose that these lesions make the inner surface of airway irregular and the smoothness is reduced compared with that of HC. The gray contrast distributions on the airway tree surface characterize the scales of smoothness. More importantly, lesions in bronchi and bronchioles will make both their lumen area and counts of CT total airway smaller [20], [61]. Therefore, the abnormal appearances of COPD airway will be mainly characterized by the irregular contours and being with fewer branches in binary snapshots. That explains why even CNN models using binary snapshots achieve rather high ACC. The reduced smoothness of COPD airway, presented as the alternations of gray contrast distributions in gray snapshot, only has limited contribution to COPD identification.

E. SPECIALLY DESIGNED CNNs VS. TYPICAL OFF-THE-SHELF CNNs
There are several ways to obtain effective CNN models for medical imaging analysis: the specially designed CNNs, the typical off-the-shelf CNNs trained from scratch [53], [62]. We have tried these strategies in the current studies. It is found that the specially designed CNNs outperform the typical offthe-shelf CNNs trained from scratch. The reasons might be two fold. First, we adopted an effective method of BOA to find the optimized architectures and super-parameters. Second, the off-the-shelf CNNs are usually developed for natural images and trained with millions of images. However, our study only includes 280 subjects, which makes the data insufficient to train CNNs from scratch.

F. LIMITATIONS AND FUTURE WORK
Our study has limitations. First, the sample size is relatively small and all subjects are from one center. Hence, it is unknown for the generalizability of our developed models across centers and populations. The independent test dataset and experiments are urgently required as done by González et al. [30]. In the next step, we will enroll more patients from different centers and update the dataset. Second, our CNN models utilize multi-view 2D snapshots, belonging to the 2D CNN or 2.5D CNN. Fully 3D CNNs are also available, but it is unfeasible because of the limited number of images and computational devices, at least at high resolutions [53]. More advanced methodological strategies of utilizing CNN can be adopted to precisely distinguish COPD from HC using CT images [63]. Third, we only include COPD and HC in the study. However, a real diagnosis model should be able to distinguish COPD from other airway diseases (e.g., asthma) and differentiate the severity of COPD (e.g., Global Initiative for Chronic Obstructive Lung Disease (GOLD) stage 2-4). Fourth, we only utilize the information of airway, but the lung parenchyma, patterns of LAA, air trapping, pulmonary blood vessels and medical record (e.g., demographic information, biomarkers, and some known medical diseases) are not included in the deep CNN model [4], [64]. Inclusion of the above information will likely improve the prediction performance. The influence of airway tree segmentation algorithm on the identification performance is also not clear [36]. In summary, these limitations are potentially solved and will be addressed in the future investigation.

V. CONCLUSION
It is feasible to identify COPD using snapshots of 3D lung airway tree extracted from CT images via deep CNN. The capability of identification mainly comes from the abnormal appearances of airway tree in COPD and the deep CNN with the ability of representing these abnormalities well. Fusion of multi-view CNN models by majority voting leads to better identification performance than CNN models from single view. Abnormal appearances of airway tree in COPD characterized by the irregular contours and being with fewer branches in binary snapshots are the major determinant of the identification and the local smoothness has certain but limited contributions. These methods and findings may help detect early COPD, reduce the rate of misdiagnosis and missed diagnosis, elucidate underlying disease pathogenesis, and leading to proper management of COPD. RAN   WEI QIAN has many years' experience on the biomedical imaging for breast cancer and lung cancer detection and diagnosis while he worked at the H. Lee Moffitt Cancer Center and Research Institute. In the past years, he has been conducting original studies in medical imaging, biomedical imaging, and molecular imaging, particularly in computer-aided cancer detection and diagnosis or remote detection and diagnosis for digital mammography, X-ray lung images, lung CT images, cellular images and molecular images, integrating intelligent techniques including the use of neural networks, fuzzy logic, genetic algorithms and evolution algorithms in simulation, modeling, and the design of high performance and robust systems. He is currently a Professor with The University of Texas at El Paso. He has published widely in these areas with more than 100 articles in peer-reviewed journals and premier international conferences (with H-Index 29). He has been awarded six U.S. patents, of which three patents have been licensed and commercialized by the Kodak Company (now Carestream Health, Inc). This is an important index for the translational research (the successful rate is 50%).