Bullous Parametric Response Map for Functional Localization of COPD

Advanced bronchoscopic lung volume reduction treatment (BLVR) is now a routine care option for treating patients with severe emphysema. Patterns of low attenuation clusters indicating emphysema and functional small airway disease (fSAD) on paired CT, which may provide additional insights to the target selection of the segmental or subsegmental lobe of the treatments, require further investigation. The low attenuation clusters (LACS) were segmented to identify the scalar and spatial distribution of the lung destructions, in terms of 10 fractions scales of low attenuation density (LAD) located in upper lobes and lower lobes. The LACs of functional small airway disease (fSAD) were delineated by applying the technique of parametric response map (PRM) on the co-registered CT image data. Both emphysematous LACs of inspiratory CT and fSAD LACs on expiratory CT were used to derive the coefficients of the predictive model for estimating the airflow limitation. The voxel-wise severity is then predicted using the regional LACs on the co-registered CT to indicate the functional localization, namely, the bullous parametric response map (BPRM). A total of 100 subjects, 88 patients with mild to very severe COPD and 12 control participants with normal lung functions (FEV1/FVC % > 70%), were evaluated. Pearson’s correlations between FEV1/FVC% and LAV%HU-950 of severe emphysema are − 0.55 comparing to − 0.67 and − 0.62 of LAV%HU-856 of air-trapping and LAV%fSAD respectively. Pearson’s correlation between FEV1/FVC% and FEV1/FVC% predicted by the proposed model using LAD% of HU-950 and fSAD on BPRM is 0.82 (p < 0.01). The result of the Bullous Parametric Response Map (BPRM) is capable of identifying the less functional area of the lung, where the BLVR treatment is aimed at removing from a hyperinflated area of emphysematous regions.


Introduction
COPD is a progressive, irreversible disease characterized by alveolar destructions due to small airways disease and loss of supporting tissue due to emphysema. Currently, the treatment options for patients suffering from emphysema are limited. Various techniques of bronchoscopic lung volume reduction treatment (BLVR) are recommended by the international expert panel as a minimally invasive procedure for improving pulmonary function, exercise capacity, and life quality in patients with COPD. These techniques include valves, coils, vapor thermal ablation, and sealant. Each method requires careful evaluation due to its limitation in treatment. In particular, endobronchial valves (EBVs) improve lung volumes and allow for the expansion of lung tissue. However, this lobar exclusion cannot enroll patients who have extensive interlobar collaterals in the emphysematous lung. On the other hand, endobronchial coil therapy which is independent of collateral ventilation improves the elastic recoil by the compression and redistribution of airflow toward healthier segments. Bronchoscopic thermal vapor ablation (BTVA, Uptake Medical Corporation, Seattle, WA, USA) reduces lung volume by the instillation of heated water in the most diseased segments to provoke irreversible parenchymal fibrosis and scarring of emphysematous tissue [1,2]. The challenge for the treatment success is, nevertheless, to define the specific signatures of the "most diseased segments or subsegments of the lobe" [3].
Quantitative computed tomography (QCT) can sensitively assess the size and spatial distribution of low attenuations clusters (LACs) as the emphysematous characteristics in patients with chronic obstructive pulmonary disease (COPD). The radiographic assessment has been developed for patient inclusion using quantitative computed tomography (QCT) as the predictors of treatment outcome, such as fissure integrity for evaluating the collateral ventilation and low attenuation clusters for evaluating emphysema heterogeneity. In the current paradigm of BLVR-specific quantitative CT imaging, emphysema lesion is segmented by the threshold of -950 Hounsfield units (HU) to extract low attenuation clusters (LACs) in the inspiratory CT data. Emphysema severity is then computed from the low attenuation volume percentage (LAV%) as the ratio of LACs to whole lung volume. However, there is no consensus definition for the emphysema heterogeneity score (HS), and many studies have proposed methods using interlobar LAV% [3,4]. Recently, Lor et al. have developed a LAC-based representation of emphysematous lesions grouped by four proportional scales in upper and lower lobes. The compositions of categorized LACs are used to build the predictive modeling of voxel-wise airflow limitation. As a result, instead of most "emphysema-destroyed lobe," the most "functionally affected" tissue can be identified for the treatment [5].
Although previous studies mainly take emphysema heterogeneity on inspiratory CT into consideration, BLVR is expected to treat the lobe with hyperinflation, which results from gas being trapped during the expiratory phase of the breathing. As suggested in the findings of recent ex vivo studies using high-resolution microcomputed tomography (microCT), the destruction of the terminal bronchiole is hypothetically the primary site in the lungs of patients with mild and moderate COPD, but with no sign of emphysema [6][7][8]. On the other hand, other studies have shown cases of emphysema without obstruction on spirometry [9]. The discordance between the results of the pulmonary function test (PFT) and structural alterations in emphysema and functional small airway disease only highlights partial aspects of the disease complexity. While the PFT provides a global assessment to diagnose COPD, CT imaging quantifies local structural and functional abnormalities, which induce the heterogeneity of COPD phenotypes, and further increase the disease complexity. As a result, both emphysema and small airway disease should be assessed for the COPD heterogeneity in the target lobe selection.
Emphysema can be obtained from LACs on inspiratory CT, whereas the morphology of the small airway (< 2 mm in diameter) is below the resolution of CT image data. Wellaccepted quantification of emphysema severity is the percentage of LAV segmented by the threshold of -900 ~ -950 Hounsfield units (HU) in total lung volume on inspiratory CT, whereas LAV% segmented by a threshold of -856 HU on expiratory CT has been used as a measure of air-trapping due to the functional small airway disease (fSAD). While both emphysema and fSAD are the main components of the heterogeneous COPD process, the surrogate measure of air-trapping is partially overlapped with emphysematous destruction. To quantify the COPD phenotypes, Galbán et al. co-registered inspiratory and expiratory CT scans to distinguish the relative contributions of fSAD and emphysema for a more accurate diagnosis [10]. Here, we co-registered the paired inspiratory and expiratory CT to identify the emphysema LACs as well as the fSAD LACs. After deriving the prediction model for the airflow limitation, the bullous parametric response map (PRM) (BPRM) is then constructed by accumulating the functional contributions of the regional low attenuation densities (LADs).

Patient and Public Involvement
The volumetric CT scans were taken at full inspiration and full expiration for each subject. However, the interpretation and conclusions contained in this study are those of the authors alone. A total of 88 subjects with symptoms of chronic obstructive pulmonary disease (GOLD stage 1 ~ 4) and 12 normal subjects was included. Subject demographics are summarized in Table 1. All subjects are evaluated in the correlational studies, and 80 of randomly selected subjects are trained in the predictive modeling for testing on the remaining 20 subjects as the setup for repeated fivefold cross-validations.

CT Image Analysis and PRM
The co-registration of expiratory and inspiratory was performed in three main steps: lung mask segmentation, threedimensional non-rigid point set registration, and image deformation. CT Automated airway and lung segmentation were performed using in-house software. In particular, the method of lung segmentation affects the accuracy of LAV%. The lung mask was segmented by applying the adaptive region growing with the upper bound value of HU 200 and using the airway mask as the starting point. After all the dark regions (all voxels with HU less than 200) were segmented and marked as the initial lung mask, the enclosed vascular regions were included by morphological operations (such as closing and filling by dilation and erosion). The separation of left and right lobes was performed by erosion to detach the anterior borders of the lobes. Then, the completed left and right lobes were obtained by the dilation separately.
To minimize the contribution of the airway lumen, the airway was removed from the lung mask. Once the lung masks were extracted from the paired CT dataset, the registration of expiratory CT (floating) aligning with inspiratory CT (reference) takes place in the left and right lobes separately. The evenly distributed sparse points of the lung mask were used as the input of the registration algorithm. We applied the coherent point drift (CPD) to transform the floating-point set to the reference point set. The outstanding performance of CPD was evaluated thoroughly in the work of Wang et al. [11] using COPDgene dataset of DIR-lab (http:// www. dir-lab. com/ index. html) [12]. While preserving the shape context, CPD also transforms the spatial distribution of the interior point set within the lung mask. The transformed interior-point set can be used as the landmarks for the image deformation, where the co-registered volume data is the result of interpolation within landmarks.
To furtherly investigate the impact of the regional emphysema on lung functions, not only the lobes are labeled as left and right lobes, but also the upper and lower lobes. The upper lobes include the right upper lobe and the left upper lobe, while the lower lobes are the right middle lobe, the right lower lobe, and the left lower lobe. The method of segmenting the upper and lower lobes is derived from the fissures as the results of thin plate structures of the hessian filter. The sample points of the fissure are the control points of thin plate spline, TPS, which can be used to delineate the borders between lobes.

Measurements of Low Attenuation Volume Percentage
The low attenuation volume is the targeting voxel of the binarized image using a fixed HU threshold value on CT. LAV% is the percentage of targeting voxels in the whole lung. In this study, we include four types of targeting voxels: mild emphysema on inspiratory CT using the threshold of HU-920, severe emphysema on inspiratory CT using the threshold of HU-950, air-trapping on expiratory CT using the threshold of HU-856, and fSAD of PRM which is a result of air-trapping on co-registered expiratory excluding overlapped emphysema voxels on inspiratory. The notations are LAV% Emph920 , LAV% Emph950 , LAV% AirT , and LAV% fSAD , respectively.

Modeling of the Low Attenuation Cluster using the Local Maxima of Bulla Voxels
This study utilized an algorithm that was developed in the previous work to find the LACs in the binarized image using a fixed HU threshold value on inspiratory CT and expiratory CT. By applying iterative erosion to the binarized image data, each LAC will accumulate the number of eroded voxels from previous steps of erosion. The method of obtaining LACs is described in the works of Lor et al. [5]. The fraction density of LACs or low attenuation density (LAD) is the total number of bulla voxels divided by the total number of parenchymal voxels. The summation of LADs is the same as the low attenuation volume percentage (LAV%). The corresponding LAD notations for each type of targeting voxel are LAD% Emph920 , LAD% Emph950 , LAD% AirT , and LAD% fSAD . While LACs of emphysema and air-trapping are segmented directly from binarized image using single-value thresholding, LACs of fSAD are not based on the fSAD of PRM. The excluding emphysematous voxels (HU < -950 on both inspiratory CT and co-registered expiratory CT) destroy the completeness of bullous structures and create more fragments than the number of visually observed clusters. The shortcoming of extracting LACs in the PRM approach can be registration algorithm dependent and thus lowers the reproducibility. More accountable approaches are referring to the classification based on the predominant ratio between emphysematous voxels (HU-950) and air-trapping voxels (HU-856). The emphysematous predominant LACs are excluded from co-registered expiratory CT, whereas the LAD% fSAD of LAC fSAD is the percentage of fSAD voxels determined by PRM.

Statistical Analysis
Pearson's test was performed to report the correlations between LAV% and FEV 1 % predicted and FEV1/FVC% on inspiratory CT and expiratory CT. The cluster analysis of LADs is performed by ckmeans.1d.dp method (which is based on dynamic programming for optimal one-dimension K-means clustering) to cluster into 10 size ranges [13]. A total of 5.3 million LACs from all the datasets on inspiratory CT and 2.8 million LACs from expiratory CT were collected to illustrate their relationship with lung functions.
The initial correlation study is conducted to FEV 1 /FVC% as the predicting target. Although the GOLD classification proposes a COPD grading system, the categorical staging result has a weak correlation with LAV%. In this study, the full dataset was evaluated for the correlation study of the disease, and the predictive model derived from the training dataset will be used to predict the outcome of the testing dataset. The final model was tested with repeated fivefold cross-validation. The full model of multivariable linear regression associated with FEV 1 /FVC% was tested for CT emphysema and air-trapping using 40 possible predictor variables. We then used backward selection to reduce the collinearity and use only 16 predictors in the final model. The final model has p-value of less than 0.05 in each predictor. All statistical analyses were performed in R statistical software (version 3.6.1; R Foundation for Statistical Computing, Vienna, Austria).

Subject Characteristics
The study included 100 participants with inspiratory and expiratory CT taken at the same time before the treatment. The dataset has shown the normal distribution of lung function. The regions of the whole lung are partitioned into upper and lower lobes. In contrast to the conventional emphysema threshold of HU-950, the proposed method uses single HU-930 to segment the mild emphysematous regions [14]. The baseline characteristics of participants are shown according to GOLD stages 0 through 4 in Table 2. Both the mean values of LAV% Emph and LAV% AirT are increased relative to the COPD severity graded by GOLD as expected.
There is no significant difference in LAV% between upper lobes and lower lobes at each stage. While FEV 1 % predicted is the GOLD criteria for staging severity (stage 1, mild: FEV 1 % predicted ≥ 80%; 2, moderate: 50 ~ 70%; 3, severe: 30 ~ 49%; 4, very severe: < 30%), the mean values of LAV% at each stage increase accordingly in the categorical subjects as shown in Fig. 1. BMI declines gradually as the COPD gets more severe.

Correlation Study Between Lung Functions and LAV%
In the past decades, Pearson's correlations between interlobar LAV% and PFT have been studied extensively, and Spearman's rank analysis was used to assess the categorical correlation between interlobar LAV% and GOLD stage classification [15]. In contrast to previous studies, we sought to characterize the heterogeneity of COPD phenotypes in terms of contribution at different stages of development. Figure 2 shows correlation coefficients between LAV% and PFT at the different GOLD stages. It can be seen that although the total dataset (n = 100) has significant correlation with PFT (r >|0.5|, p < 0.01), only patients at GOLD stage 3 (n = 22) have contributions with significant correlation (r >|0.5|, p < 0.05) with FEV 1 % pred. Table 3 shows correlation coefficients between various LAV% and PFT at the different GOLD stages. In contrast to the results of FEV 1 % pred, it can be seen that while severe emphysema (GOLD 3) has the most contribution to FEV 1 /FVC% in lung function decline, air-trapping and other composite measurements using both inspiratory and expiratory CT have a better linear association with earlier stages of COPD. Interestingly, the weaker correlation at stage 2 (r <|0.3|), as comparing to stage 1 and 3 (r >|0.5| and r >|0.8|), suggested a non-linear relationship between LAV% and FEV1/FVC% in patients with the mild syndrome. The resulting correlation comparison of Table 3 is shown in Fig. 2. While the ideal Pearson's correlation coefficients are expected to be equally strong in all GOLD stages, FEV 1 % predicted only has a moderate correlation with severe COPD, and FEV 1 /FVC% was shown stronger in all stages as compared to FEV 1 % precited.
Other than a simple thresholding method for quantitative measurement of air-trapping, some studies aimed at differentiating air-trapping from those remaining in emphysematous regions by the commonly used ratio between the mean lung density at expiration and inspiration (E/I) ratio [16] and more sophisticated method such as parametric response mapping (PRM) [10]. With the fact of LAV% on expiratory CT is the mixture of air-trapping and fraction of emphysema in mind, we studied the correlation between PFT and the mixing effect of LAV% Emph and LAV% AirT by simple addition. The result is shown in the second last category of Table 3. Despite the possible introduction of collinearity, it can be seen that the linearity is strengthened, achieving higher correlation coefficients in the total dataset (n = 100). Furthermore, we studied the correlation between PFT and the composite effect of emphysema and functional small airway disease, two main components of COPD. In the last category, LAV% Emph950+fSAD has shown the strengthened correlation analysis as expected for eliminating the overlapped area of emphysema from expiratory CT. These findings affirm the concept that the heterogeneity of airflow limitation (FEV 1 % predicted and FEV 1 /FVC% in this study) is the result of the non-linear contribution from the low attenuation lesions revealed in parenchyma on inspiratory and expiratory CT. Another interesting result showed that Emph% and AirT% have a noticeable stronger correlation with FEV 1 /FVC than FEV 1 % predicted (in agreement with the findings of previous studies [7,9,15,17,18]). Hence, we use FEV 1 /FVC % as the predicting target in this study.

Cluster Analysis of LAD
In the previous study [5], we have developed predictive modeling of FEV1/FVC using LAD distribution on inspiratory CT to predict the functional severities of subjects, differentiated in spirometry but equivalent in LAV%, as well as those equivalents in spirometry but differentiated in LAV%. In this study, we extend the concept to apply the cluster analysis on expiratory CT. Unlike the 4 empirical scales used in the previous study, the LADs were clustered into 10 scales by univariate K-means clustering method. A noticeable similarity in boxplots of 10 scalar clusters is illustrated in Fig. 3 A (inspiratory CT) and B (expiratory CT). Figure 3C, D are the boxplots showing the total FEV1/FVC distribution of each scalar clusters for all the data (n = 100). The total number of clusters obtained on inspiratory and expiratory CT is more than 8 million. Although the studied subject can be composed of LADs varied in sizes within lobes, the distribution shown in boxplots in Fig. 3 can be seen as the increased probability of having severe airflow limitation when the heterogenous lung destruction is composed of more LADs in cluster sizes 6 to 8.

Linear Model Fitting and Derivation of the Predictive Model
It is well understood that the mechanism and transition of how low attenuation relates to disease progression are not yet clear. We hypothesized that the severity of airflow limitation can be partially attributed to the effects of LAD% at a different level according to the size and location of the LAD. To evaluate the effects of the LADs on the severity of lung function, we formulate the full model of multiple linear regression using 10 LAD scales from upper lobes and 10 LAD scales from lower lobes on  Table 3. Plot A compares the average LAV% with FEV1% pred, whereas plot B compares with FEV1/FVC% inspiratory and expiratory CT. A total of 40 LAD scales formulate the predictors used in the full model. Substantial collinearity is expected because of the overlapped effect due to the partial LAC Emph inevitably existing within the LAC fSAD on expiratory CT. Backward selection is applied to improve the statistical significance of the model.

Prediction Performance Evaluation
In order to evaluate the robustness of the final model, we initially derived the estimated coefficients of linear models on the complete dataset using fivefold cross-validation repeated for 30 times. The Pearson's correlation between FEV1/

Discussion
In this study, the correlation between LAV% Emph920 and PFTs is only statistically significant in the lungs with a severe GOLD stage. The result showed that LAV% Emph920 does not have a linear association with PFTs, and LAV% AirT correlates better than LAV% Emph920 with PFTs at each GOLD stage. LAV% AirT encloses areas of air-trapping as well as emphysematous regions, which can be decomposed into LAV% Emph950 and LAV% fSAD as shown in Fig. 5. Accessing the functional distribution of low attenuation clusters (LACs) in different sizes, we found that giant bullous emphysema (cluster-10) does not necessarily attribute more impact on the lung function when comparing to smaller bullae. Weaker correlation with PFTs is therefore expected in the lungs consisting of giant bullous emphysema. Clinical studies have shown the heterogeneous functional outcome among the subtypes of emphysema. Larger bullae do not necessarily attribute to severe COPD, and smaller, evenly distributed bullae might be more common in the lobe with advanced destructive emphysema.
From the pathophysiological perspective, pursuing a linear relationship between PFTs and LAV% would be more challenging. The centrilobular emphysema (CLE), panlobular emphysema (PLE), and paraseptal emphysema (PSE) are three main phenotypes of emphysematous destruction seen in COPD. PLE refers to diffuse emphysema across the lobule, whereas CLE is considered the primary lesion of destruction at respiratory bronchioles. Not only that each subtype of emphysema differs in sizes, distribution, and location, but they may also have different pathophysiologic contributions to the increase in small airway resistance in COPD. For example, PLE is considered less associated with small airway obstruction than CLE, because the alveolar destruction in PLE is milder than that in CLE, which at the later stage spreads to fuse destroyed lobules and becomes large bullous lesions [19,20]. Other studies reported that PLE although represents more advanced emphysema has a weaker association with FEV 1 % as compared to CLE [21]. Moreover, previous studies reported the confounding functional effect of bullous emphysema (advanced CLE) in patients with diffuse emphysema (PLE) [22].
Recent studies have shown that the presence of fSAD might be predictive of spirometry decline [18,23]. From the functional perspective, the extent of hyperinflation can be quantified by the decline in FEV 1 % as the result of increasing bronchitis/bronchiolitis in small airway disease. Detection of LAC fSAD might not only provide detection of COPD at its early stage but also differentiate the functional severity among heterogenous COPD phenotype by the assessment of regional association with PFTs. Another study reported that current smoker with regional air trapping tends to have less emphysema and better lung function than those without [16]. More studies have hypothesized that small airway disease might be a precursor to emphysema [7,18,24]. As a result, we extend our previous work of emphysema on inspiratory CT to include fSAD on expiratory CT, in order to provide better functional correlation. Our findings of the better correlation between combined LAV% Emph950+fSAD and FEV1/FVC% support our hypothesis that emphysema and air-trapping are complementary to each other in association with airflow limitation.
The current target lobe selection in BLVR is based on densitometric data on inspiratory CT scans. In the work of Kloth et al. [25,26], the cluster analysis of LACs representing connected voxels of HU values under -950 was conducted using Pulmo 3D software (Siemens Healthcare). Unlike the fraction of density integrated into this study, Kloth et al. categorized the LACs of inspiratory CT in four empirical volume clusters to visually differentiate the homogeneous and heterogeneous of the emphysema phenotypes. Studying the treatment response of BLVR's endobronchial coiling, they found that the emphysema phenotypes of the target lobe had no significant impact on Fig. 4 First row: left is inspiratory CT, the center is the original expiratory CT, right is the co-registered result of expiratory CT. Second row: left is the vascular tree rendering of inspiratory CT, and center is the expiratory CT. The visual inspection shows extensive pruning in lower left lobes and relatively mild pruning in lower right lobes. Third row: Red color marks the emphysema regions (HU-950), green color marks the normal regions, yellow marks the air-trapping of expiratory CT, distinguished from emphysema of inspiratory CT using technique of parametric response map (PRM). The fourth row is the result of second-row using low attenuation clusters. Fifth row: left is the visualization result of PRM in a three-dimensional model, the center is the mapping of the PRM on HU values of inspiratory CT, right is the result of regional functional prediction mark on the grids, which has blue color as the most severe candidates for the target lobes the outcome [26]. However, it has to be considered that determining the emphysema heterogeneity in our study is more than qualitative analysis. Lacking the objective image descriptor for emphysema phenotypes with pathophysiological meaning, we pursue the regional COPD severity with the functional prediction from the compositions of the LADs. The proposed predicting model for airflow limitation using low attenuation data on respiration CT has shown a remarkably strong correlation with statistical significance. Although the emphysematous destruction might be caused by more than one phenotype, making it more difficult to identify the disease trajectory, the longitudinal changes using the proposed model should be studied to provide clinical insight and utilization.
We are aware of the limitations of this retrospective study. We note that the major limitation is the number of subjects included in the study, particularly those at GOLD stage 4. Secondly, the study does not include the measurement of airway wall thickness, although it has been shown to have an association with lung function only in CLE [20]. Another limitation, as well as the future work of this study, is the co-registration of inspiratory and expiratory CT. The transition of LACs in respiration has been studied previously to assess the expansion and collapse of the LACs [10,27,28]. Future studies will extend the utility of LACs localization and qualification to distinguish emphysema subtypes.

Conclusions
In conclusion, the pattern of lung destruction revealed in the low attenuation clusters (LACs) of co-registered CT images has a direct impact on the lung functions. The result of analyzing PFT and LAV% correlation, which is also supported by other studies, shows that LAV% AirT on expiratory CT has a better association than LAV% Emph on inspiratory with spirometry measurements. We have extended the analysis and showed the appreciated correlation contributed by combing both fSAD on expiratory CT and emphysema on inspiratory CT. On top of the PFT correlational analysis, in the categorical analysis of low attenuation volume percentage (LAV%) at a different GOLD stage, we have shown that only patients suffering from severe COPD have a moderate correlation with LAV% on inspiratory CT. The heterogeneity of the disease phenotype is thus evidently affecting the functional outcome measured by PFT. After collecting more than 8 million LACs from 100 subjects, the cluster analysis divides the LACs into 10 scales according to the occupying fractions. The percentage occupying the parenchyma by each LAC can be seen as a fraction of total low attenuation density (LAD). The distribution of lung function in each LAD has shown that the size of LACs has a non-linear relationship with lung function. Particularly, the subjects with mediumsized LACs tend to have weaker lung function than those with diffusive lesions which form into much larger LACs. Furthermore, these LACs are distributed across the entire lung, and their spatial information also has a functional contribution to the disease. Together with the size and the location of LACs, the heterogeneity of emphysema and fSAD can be generalized by the total volume of each LAD category in upper and lower lobes. Utilizing the derived predictive modeling of airflow limitation, we are able to achieve a much stronger correlation between CT image predictors and PFT measurements. The low attenuation density distribution was furtherly incorporated into the BPRM to show a promising result in identifying regional COPD severity for the treatment planning of BLVR.

Declarations
Ethics Approval and Consent to Participate This retrospective study was based on the CT dataset and PFT measurements obtained in the clinical routine and approved by the institutional review board (IRB) of the National Taiwan University Hospital (NTUH) and Cardinal Tien Hospital (CTH). The need for individual informed consent was waived by the Institutional Review Board because of the retrospective nature of the study, and the medical data accessed were anonymized and deidentified prior to analysis. The volumetric CT scans were taken at full inspiration and full expiration for each subject. Disclaimer The interpretation and conclusions contained in this study are those of the authors alone. A total of 88 subjects with symptoms of chronic obstructive pulmonary disease (GOLD stage 1 ~ 4) and 12 normal subjects was included; 55 cases of NTUH, and 45 cases of CTH.

Competing Interests
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.