2.1.

Study Population

This was a retrospective study that included 50 solitary pulmonary nodule cases (22 malignant and 28 benign) from previous NIH funded studies: the NLST and the Chronic Obstructive Pulmonary Disease Genetic Epidemiology (COPDGene) study. Our study was approved by the University of Iowa institutional review board, and informed consent of participants was obtained. CT and demographic data were collected as a part of the parent studies. For the NLST study, lung cancer diagnoses were tracked as a part of the primary study outcomes; however, COPDGene subjects were separately contacted and consented in order to collect details related to nodules detected on their CT studies that were acquired to study chronic obstructive pulmonary disease (COPD).

Our study consisted of 14 cases from the Iowa cohort of the NLST; 6 nodules were malignant and 8 were benign, all diagnoses confirmed through pathology. The remaining 36 cases—16 malignant and 20 benign—were from the COPDGene study. For 12 of the malignant cases from the COPDGene study, the diagnoses were confirmed using the COPDGene Lung Cancer Registry, an ongoing ancillary study to collect diagnostic and treatment information from participants with both COPD and lung cancer. Additionally, 11 of the benign nodules remained benign after at least 2 years of follow-up. The diagnoses of the remaining four malignant and nine benign nodules from the COPDGene study were determined by comparing the self-reported presence/absence of lung cancer by the participants and the radiologists’ qualitative description of the nodules in their written reports. These 13 nodules will be reassessed upon the return of the participants for stage II of the COPDGene study.

2.2.

Computed Tomography Data

NLST and COPDGene acquired CT data with very different goals in mind. NLST, targeting annual screening for lung cancer, developed a low-dose protocol to minimize participant exposure to ionizing radiation. The low-dose CT scans from the NLST cohort were acquired from a single center site, the University of Iowa. All 14 NLST scans were collected with a tube voltage of 120 kVp and a tube current–time product of 40 to 80 mAs on either a Philips Mx8000, Siemens Sensation 64, or Siemens Sensation 16 CT machine. While the NLST protocol specified a slice reconstruction interval of 1.0 to 2.5 mm, additional high resolution reconstructions were made at our institution and used in this study.³ The three Philips scans had a slice thickness of 1.3 mm and were reconstructed with a B filter. The Siemens scans had slice thicknesses of 0.75 mm (four scans) or 0.6 mm (seven scans) and were reconstructed with a B30f kernel.

The imaging protocol for COPDGene was developed for high resolution lung imaging to best characterize emphysema, air trapping, and airway wall remodeling from the COPD disease process.²³ The 36 higher dose CT scans from the COPDGene subcohort were collected from Brigham and Women’s Hospital, Columbia University Medical Center, Morehouse School of Medicine, National Jewish Health, Reliant Medical Group, Temple University, the University of California San Diego, and the University of Iowa. Scans were collected with a tube voltage of 120 kVp and a tube current–time product of 200 mAs. Reconstructions were performed with B, B31f, or Standard kernel for Philips, Siemens, and GE respectively, with slice thicknesses between 0.6 and 0.9 mm.

For each of the 50 subjects with an identified solitary pulmonary nodule, the ROI was defined [Fig. 1(a)]. The ROI was selected within the CT data to include the solitary pulmonary nodule and immediate surrounding parenchyma. The amount of parenchyma included in the ROI was roughly proportional to the size of the nodule. Each ROI was manually segmented by a single user (J.U.) to label the individual voxels of the ROI as nodule, surrounding parenchyma, or invalid tissue. In addition, to explore the impact of different segmentations on feature extraction, a subset of 12 cases was independently segmented by three separate users (A.J., S.D., and K.K.). Voxels marked as invalid tissue included chest wall and blood vessels; these voxels were not used in feature extraction.

Fig. 1

Flowchart of the computer-aided diagnosis (CAD) tool. (a) Segmentation: The region of interest containing the nodule and the surrounding parenchyma is labeled as nodule or parenchyma. (b) Feature extraction: Over 300 features are extracted from these two areas, and an additional 54 global features are extracted based on characteristics of the lung. (c) Classification: Finally, the feature set is reduced and the selected features are used to train the linear discriminant analysis (LDA) and artificial neural network (ANN) classifiers, labeling a nodule as likely malignant or likely benign.

2.3.

Feature Extraction

Voxels labeled as parenchyma and nodule were used in the extraction of four classes of features [Fig. 1(b)]: intensity, shape, border, and texture. Intensity characteristics included histogram features, such as mean, minimum, and maximum intensity (in Hounsfield Units), as well as measures of heterogeneity, such as entropy and kurtosis; these were extracted from the nodule and parenchyma voxels separately. Shape characteristics such as sphericity and features comparing the surface of the nodule to a sphere of equivalent volume were extracted using the nodule mask. The size of the nodule was also extracted using the RECIST maximum diameter and the effective radius, calculated from the volume of the nodule.

The texture of both the nodule and the parenchyma were each quantified using a three-dimensional (3-D) extension of Laws texture energy measures (TEM).²⁴ These measures interrogate the grey-scale images by looking for patterns in different gradients. In the two-dimensional (2-D) application, five 5-element vectors that describe levels, edges, spots, waves, and ripples are convolved with each other to form 25 $5\times 5\text{matrices}$ called kernels. These kernels locate a different texture type in each dimension, such as edges in the horizontal and spots in the vertical. By convolving the five-element vectors with the 25 2-D kernels, 125 3-D kernels in the form of $5\times 5\times 5\text{matrices}$ were generated. These 125 kernels were convolved with each ROI, generating 125 texture volumes that, when combined and normalized to account for the directional texture information, simplified into 34 rotationally invariant 3-D TEMs. Using the same method for intensity feature extraction, the mean, variance, kurtosis, and skewness of the textures for the nodule and parenchyma were extracted for each 3-D TEM. This resulted in 272 texture feature values.

A lung CAD tool developed by Way et al.²² examined the nodule border through the use of a rubber band straightening transform, originally implemented by Sahiner et al.²⁵ This process straightens the nodule based on the perimeter identified during segmentation. While the identification of the boundary can vary greatly between users and segmentation methods, the centroid is minimally affected, as illustrated in Fig. 2. Quantitatively, this was tested in the subset of 12 nodules segmented by three different users, yielding a mean centroid difference of $0.45\pm 0.41\mathrm{mm}$. By shifting the straightening from the border of the nodule to the nodule’s centroid, the dependence on segmentation is reduced. To accomplish this task [illustrated in Fig. 3(a)], the axial slice of the ROI containing the centroid was selected. For this 2-D slice, the distance between the centroid and the closest bounding box edge was computed; this length was determined to be the ray length. Rays were cast from the centroid toward the edge of the bounding box radially, resulting in 360 rays of uniform length, with each representing a single row in the straightened image. Once the straightened image was created, the border was analyzed [Fig. 3(b)]. As benign nodules tend to have well-defined borders,²⁶ we expect the average rate of border change to be steeper, indicating a quicker transition to parenchyma. However, for malignant nodules, which tend to be spiculated, it is predicted that the intensity would change more slowly and be more varied across the rows of the straightened image (higher values in both of the variance border features).

Fig. 2

Comparison of segmentation of a nodule compared to placement of centroid. The nodule (a) was manually segmented by two users: (b) A.J. and (c) S.D. The grayscale image of slice 20 is shown in the top image, while the two segmentations with centroid indicated are shown in the lower images. Though the segmentations differ, the centroid placement remained constant at pixel location (29, 33, 20).

Fig. 3

Process for quantifying the border from the straightened image. Nodule straightening: the slice containing the centroid of the nodule mask is selected (a). The minimum distance from the centroid to the bounding box edge is computed [circle in (b)] and rays cast radially for 360 deg [arrow head in (b)]. The intensity values are interpolated along each ray for the straightened image [right panel in (b)]. Each row of the straightened image contains the same number of pixels as the uniform-length cast ray, with the bold line in (b) corresponding to the horizontal line in the row of the straightened image. The same process is repeated for the nodule mask. The straightened mask is combined with the straightened image to remove the invalid tissue (c). Feature extraction: for each row in the straightened image (c), a plot of intensity versus column position is generated (e). Shown by the thick vertical lines, the minimum and maximum column of the border is found from the straightened mask (d). A first-order polynomial was fit for each intensity plot (e) over the column range in (d). The mean and variance of the slopes as well as the width in the border was computed for features (f).

As an additional area of interest, global features were included to determine whether imaging data captured within the same lobe, third of the lung, and/or lung contain meaningful information related to nodule diagnosis. Measures including region-specific (lobe, third, or lung) mean and standard deviation intensity, tissue volume, air volume, and percent emphysema (percent of voxels $<-950\mathrm{H}\mathrm{U}$ on inspiratory imaging) were extracted through the commercially available lung analysis software, Apollo Software (VIDA Diagnostics, Inc., Coralville, Iowa). About 54 global measures were included in feature extraction, leading to a total of 165 nodule, 139 parenchymal, and 54 global features.

2.4.

Feature Selection and Classification

The 358 features collected from the ROIs were then reduced through statistical analysis and feature selection [Fig. 1(c)]. A statistical analysis module was created to determine which of the features were statistically different between the malignant and benign cases. Each feature was tested for normality using the Jarque–Bera normality test. Additionally, the variances for each feature were tested for equality. If the feature followed a normal distribution and the variances for the malignant and benign groups were equal, a two-sample $t$-test for equal variances was performed. If the feature followed a normal distribution but did not have equal variances between the two groups, the two-sample unequal variances $t$-test was used. Finally, if the normality test showed the feature distribution was not normal, the nonparametric Wilcoxon rank sum test was used. A feature with a $p$-value $<0.05$ was determined to be statistically significant. This statistical testing was used to reduce the number of features used in feature selection in order to decrease computation time while ensuring meaningful features were included in analysis.

To determine which of the significant features should be used for classification, stepwise forward selection was used. A framework was developed for optimizing feature selection when a small number of ROIs were available in an effort to reduce both performance and selection biases (Fig. 4). The performance of the selected features at each step was estimated using leave-one-case-out cross-validation on 49 cases, with the 50th case completely left out of training and testing. This was repeated so every case was completely removed from the process. Majority voting was then used to select the best feature set at each stage of stepwise forward selection. Features were added until the best combination of features was found, with the total number of features equivalent to 1 feature for every 10 ROIs to prevent overfitting;²⁷ with 50 ROIs, features were added until five features were selected or until performance was not improved with the addition of more features. An artificial neural network (ANN), the structure of which is illustrated in Fig. 1(c), was used for classification. Performance was estimated using a leave-one-out cross-validation method where 49 cases were used for training and the 50th used for testing. In the hidden layer, two nodes were used in addition to the bias node. Due to the random initialization of weights during training of the ANN, this was repeated 10 times to better approximate the classifier’s performance. Classification using the features selected was also performed using a linear discriminant analysis (LDA) classifier [Fig. 1(c)]. Due to the deterministic nature of the LDA classifier, the leave-one-out cross-validation method was performed once.

Fig. 4

Process for selecting the best set of features at each step along stepwise forward selection. For each iteration, one case is left completely out of the internal loop, where accuracy approximations are computed for each subset of features being tested. After acquiring accuracy approximations for all possible subsets for each set of 49 ($N-1$) subcohorts, the feature set with the highest approximate accuracy across all subcohorts is selected.

2.5.

Classification Evaluation

The accuracy, sensitivity, specificity, and the standard deviations of each measure were computed by averaging the leave-one-case-out cross-validation results from 10 trials with each ANN classifier. Performance was evaluated through ROC analysis in which the area under the ROC curve (AUC) was obtained. Additional performance measures of positive predictive value (PPV) and negative predictive value (NPV) were also computed. Classifier ROC curves were compared using the Delong method to assess for differences in performance.²⁸ Predicted classification probabilities were averaged over the 10 trials for comparisons of ANN performance. All statistical comparisons were two-sided and assessed for significance at the 5% level.

3.1.

Statistical Analysis

Of the 358 features (165 nodule, 139 parenchymal, 54 global), 52 were found to be significant predictors of malignancy ($p<0.05$). In summary, the significant predictors incorporated features from both the nodule and the surrounding parenchyma including the following 47 features:

Comparatively, when only the nodule features were considered, only nine features returned as significant:

There were five significant global features:

The effects of this manual segmentation were analyzed by comparing the feature values across three different segmentations (A.J., S.D., and K.K.). Based on the results of the one-way analysis of variance test, none of the features varied significantly from the others when segmented by a different user. The intensity histogram ($p=0.16$ to 1.00) and shape ($p=0.16$ to 0.78) features were most affected by manual segmentation, as expected due to reliance on the nodule mask. However, none of the features were statistically significantly affected.

3.2.

Classification

From the subset of 47 statistically significant nodule and parenchyma features, the features chosen by stepwise forward selection for the ANN classifier can be seen in Table 1, along with the performance of the classifier using leave-one-case-out cross-validation. Likewise, the features selected for classification from the nine-feature subset for the nodule only are also shown in Table 1. Both classifiers selected the nodule minimum intensity as a feature. While the nodule and parenchyma classifier selected intensity measures of both the nodule and the parenchyma and a shape feature, the nodule-only classifier selected size and border features, in addition to intensity measures. For the nodule-only classifier, only four features were selected when up to five could be used to maintain the 1 feature per 10 cases rule of thumb. In this instance, performance was not improved with the addition of a fifth feature. The ANN built with nodule and parenchyma outperformed that built with features from the nodule alone, correctly classifying 46 of the 50 nodules on average with an average AUC of 0.938, compared to the nodule-only performance of $44.9/50$ and average AUC of 0.918. Both classifiers had similar sensitivities (90.9% to 93.6%) or $20/22$ to $20.6/22$ correctly classified as malignant.

Table 1

Summary of performance measures for the artificial neural network (ANN) classifications. Due to the random initialization of the neural networks during training, performance is reported in percent as the mean ± the standard deviation over the 10 trials. The positive predictive value (PPV), negative predictive value (NPV), and the area under the ROC curve (AUC) are also reported.

As mentioned previously, the classifier built with features only from the nodule selected the size feature: the RECIST maximum diameter. Due to a size bias existing within the cohort, as indicated by the statistical significance of this feature ($p=0.01$), the size feature was eliminated from feature selection. Without this feature, the ANN relied on the border feature and nodule intensities and only selected three features before performance was not improved. This ANN had a resulting accuracy of 87.8% ($43.9/50$) and an average AUC of 0.872 and suffered a reduction in sensitivity and NPV compared to the ANN built with nodule features including size. Overall, the ANNs built with nodule-only features had worsened outcomes, including decreased accuracy, specificity, PPV, and AUC.

A final trial looking at the impact of nodule, surrounding parenchyma, and global features resulted in a similar set of features being selected as the nodule and parenchyma CAD tool (Table 1). Instead of parenchyma variance being selected as the fifth feature, the kurtosis of the third of the lung containing the nodule was selected. The ANN incorporating all three types of features (nodule, surrounding, and global) had a performance accuracy of 96.0% ($48/50$), improved sensitivity and specificity compared to the other classifiers, as well as improved PPV and NPV. It had a similar but reduced AUC compared to that of the nodule and parenchyma ANN, with an AUC of 0.932.

Across the four trials, AUCs for the CAD tools using ANN classifiers ranged from 0.872 to 0.938. Their differences were not statistically significant, although the modest sample size in this study limits power to detect small differences. AUCs for LDA classifiers ranged from 0.710 to 0.857 and were significantly lower than those for the ANN classifiers ($p<0.05$). A comparison of the ROC curves for these classifiers can be found in Fig. 5. The measures of performances for the LDA classifiers, including the accuracy, sensitivity, specificity, PPV, NPV, and AUC, are shown in Table 2. As seen with the ANN classifiers, the nodule and parenchyma LDA ($\mathrm{AUC}=0.847$) and the LDA incorporating nodule, parenchyma, and global features ($\mathrm{AUC}=0.857$) performed better than the LDAs that included only nodule features ($\mathrm{AUC}=0.710$ to 0.713), providing additional support for the inclusion of parenchyma features in a CAD tool.

Table 2

Summary of performance measures for the linear discriminant analysis (LDA) classifications. The accuracy, sensitivity, specificity, PPV, NPV, and the AUC are reported.

Fig. 5

ROC curves for the (a) ANN and (b) LDA classifiers. The ROC curves generated from the ANN classifiers show better performance compared to the LDA classifiers, which is supported by the area under the ROC curve values seen in Table 1 and Table 2.

4.1.

Limitations

The ANN classification approach incorporated a stepwise forward selection method utilizing majority voting. This approach was designed to take advantage of a smaller training dataset while still allowing the maximum number of features to be selected without leading to over-fitting or over-training. While we reported a trend of increased CAD performance through incorporation of parenchymal features, the modest sample size in this study did not allow us to demonstrate statistically significant improvements in AUC. We are currently collecting a retrospective clinical cohort to which this CAD approach will be applied for further testing and validation in a larger cohort.

The nodule cases in this study were manually segmented, which is a time-consuming process. However, we present stability of our features despite separate segmentations by three independent users, as support for the adaptability of our system to accommodate automated or semiautomated segmentation approaches in the future.

Ensuring accurate truth in the diagnoses of the cases is always a major concern in CAD development. While COPDGene’s primary goal is to study COPD, the identification and determination of lung cancer is a secondary purpose that, until recently, relied heavily on participant reporting. However, recent support has been acquired for the collection of clinical diagnosis, classification, and treatment approach for participants with lung cancer. This ongoing project has led to confirmation of lung cancer diagnoses for 18 of the 22 malignant nodules for the cohort. Additionally, the COPDGene study is currently re-recruiting participants for a second phase of testing (including CT) scheduled for 5-year postphase I testing to assess progression of COPD within the cohort. This has led to 19 of the 28 benign nodules being confirmed as benign with at least 2 years of follow-up. The diagnoses of the remaining four malignant and nine benign nodules will be verified upon the participant’s return for phase II of COPDGene. Additionally, as COPDGene participants return for phase II, we will continue to expand the number of nodules in the cohort to better assess the statistical significance of the differing classifiers' performance.

4.2.

Conclusions

In this study, we report the development of an expanded feature extraction method and accompanying ANN CAD approach for lung nodules detected with CT; incorporating QCT features from not only the solid nodule but also the surrounding lung tissue. The resulting ANN classifier was based on five features—two nodule intensity features, one size feature, and two parenchyma intensity features that also quantify texture and resulted in excellent lung nodule classification (AUC of 0.94, PPV of 0.91, and NPV of 0.93). While the AUC was not improved when global features were included, the PPV and NPV of the nodule, surrounding, and global ANN were improved (AUC of 0.93, PPV of 0.96, and NPV of 0.96). The ANN tools that incorporated only nodule features achieved lower performance values (AUC of 0.92, PPV of 0.85, and NPV of 0.95). This study was underpowered to show statistically significant improvement in the resulting AUC values of the systems, due to a modest sample size. However, the trend toward increased performance when the parenchyma is included, coupled with the large number of significant parenchymal features indicate the value of the developed, expanded QCT feature set.