Histological classification of non-small cell lung cancer using Gabor filtering Minkowski functionals and Neighborhood Component Analysis on CT images.

Information regarding the histological types of non-small cell lung cancer is essential to determine the treatment strategy. Although several radiomics studies using almost similar feature variables were reported, a considerable variation in the performances has been observed. In this study, as novel radiomic features, 2D Gabor filtering Minkowski functionals were used. They were calculated in rotational invariant and both scale and rotational invariant ways using circular shift operations of Gabor filters on nonenhanced computed tomographic images. Eighty-six patients (47 adenocarcinomas, 39 squamous cell carcinomas) were analyzed. Two independent observers manually delineated a single slice segmentation of a tumor. Feature selection was made by neighborhood component analysis. Among various classifiers, 1-nearest neighbor gave a promising performance. The observer-averaged accuracy of rotational invariant analysis was 86.28% and that of both scale and rotational invariant one was 88.27%. However, there was no common feature among the ten top-ranked features of each observer with the identical Gabor filtering type. Hence further study of the robustness is necessary to create a more reliable model. This study aims to develop a novel radiomics method for the histological classification of NSCLC on nonenhanced CT images. 2D Gabor filtering (GF) Minkowski functionals (MFs) were used as quantitative imaging features. In the previous studies, a vast range of radiomics features in which MFs were not included have been investigated. MFs are fundamental quantities of Integral-geometry morphological analysis 10 . In 2D, area, perimeter, and Euler number characterize the geometrical and topological pattern formed by black-and-white (objects and background) binary images. MFs have proven to be very useful for statistical analysis of microstructures in material science 11 . They have also been used to detect tumor response heterogeneity to treatment and differentiate pseudoprogression from glioblastoma progression on MR images 12-13 . In the previous studies 2-9 , the majority of essential features were wavelet decomposed features. It implied that a single resolution approach was not sufficient and a multiresolution approach was useful. GF is one of the most commonly used methods to extract features at different scales and directions 14 . Unlike the discrete wavelet decomposition, GF can select the filtering direction arbitrarily and perform multi-directional analysis. Moreover, GF circular shift operation enables the analysis of object’s rotational invariant (RI) and both scale and rotational invariant (SRI) properties 15-16 . . features from high-dimensional data. Several nearest neighbor (NN)-based feature selection methods have been successfully developed. et al. showed that reliefF, a popular NN-based method, was optimal among 24 feature selection methods for the histological classification of NSCLC. As a novel NN-based feature selection method, neighborhood component analysis (NCA) was proposed 18-19 . The better performance of NCA than that of the state-of-the-art NN-based methods was also reported 19 . In the present study, NCA was used as the feature selection method and the performances of various machine learning models were evaluated. only ground-grass opacity (GGO)


Introduction
Accurate histologic classification of non-small cell lung cancer (NSCLC) between adenocarcinoma (AD) and squamous cell carcinoma (SCC) is mandatory in selecting the best systemic therapeutic option for patients. The classification is not always feasible technically because of limited tissue samples of lung cancer 1 . A computed tomographic (CT) scan is a commonly exploited medical imaging tool for lung cancer patients. However, on the CT images, the visual differentiation between lung cancer histological subtypes is difficult.
This study aims to develop a novel radiomics method for the histological classification of NSCLC on nonenhanced CT images. 2D Gabor filtering (GF) Minkowski functionals (MFs) were used as quantitative imaging features. In the previous studies, a vast range of radiomics features in which MFs were not included have been investigated. MFs are fundamental quantities of Integral-geometry morphological analysis 10 . In 2D, area, perimeter, and Euler number characterize the geometrical and topological pattern formed by black-and-white (objects and background) binary images. MFs have proven to be very useful for statistical analysis of microstructures in material science 11 . They have also been used to detect tumor response heterogeneity to treatment and differentiate pseudoprogression from glioblastoma progression on MR images [12][13] . In the previous studies 2-9 , the majority of essential features were wavelet decomposed features. It implied that a single resolution approach was not sufficient and a multiresolution approach was useful. GF is one of the most commonly used methods to extract features at different scales and directions 14 . Unlike the discrete wavelet decomposition, GF can select the filtering direction arbitrarily and perform multi-directional analysis. Moreover, GF circular shift operation enables the analysis of object's rotational invariant (RI) and both scale and rotational invariant (SRI) properties 15-16 . In radiomics study, a proper feature selection method improves the performance of machine learning 17 . Feature selection is a technique eliminating irrelevant and redundant features from high-dimensional data. Several nearest neighbor (NN)-based feature selection methods have been successfully developed. Wu et al. 4 showed that reliefF, a popular NNbased method, was optimal among 24 feature selection methods for the histological classification of NSCLC. As a novel NN-based feature selection method, neighborhood component analysis (NCA) was proposed [18][19] . The better performance of NCA than that of the state-of-the-art NN-based methods was also reported 19 . In the present study, NCA was used as the feature selection method and the performances of various machine learning models were evaluated.

Methods
The radiomics analysis was performed using MATLAB (Release 2020a, The MathWorks, Natick, MA).

Patients
Ethics Committee, University of Toyama approved this study and waived informed consents.
This study was performed in accordance with the relevant guidelines and regulations.
Radiation oncologist searched the clinical and radiological database of the institution. The data of NSCLC patients who received radiation therapy at the University of Toyama Hospital from April 2014 to May 2020 was collected. The inclusion criteria were as follows: (1) pathologically confirmed AD and SCC except for favor AD and SCC; (2) available nonenhanced CT images with a voxel size 0.39 ×0.39 × 1 mm 3 before chemotherapy or radiotherapy. The patients with only ground-grass opacity (GGO) lesions were excluded.
Definitive diagnosis was confirmed by cytology, biopsy, or surgical resection. Voxel size resampling was not used in this study, since the resampling was not sufficient for some intensity histogram and texture features 20 . The numerical values of those features were highly correlated with a voxel count of a tumor and feature normalization using the voxel count was proposed 21 .
Eighty-six patients (47 AD and 39 SCC ) were included. Table 1 summarizes the characteristics of patients. AD had a significantly higher female rate ( = 0.003, 2 test). In T2 and T4 stage, there were significant differences between AD and SCC ( =0.0006, and = 0.002, 2 test). SCC had a significantly higher rate of T3 or more tumor than AD ( = 0.0004, 2 test), while AD had a significantly higher rate of T2 or less tumor than SCC ( = 0.0009, 2 test).

Region of interest (ROI) segmentation
A tumor was segmented manually using image segmentation software (ITK-SNAP, version 3.6.0). Two independent observers did the segmentation as radiation oncologist (Obs. 1) and radiologist (Obs. 2). We delineated a single slice ROI of a tumor instead of 3D whole tumor segmentation because SCC had a significantly higher rate of advanced T-stage tumor than AD. Obs. 1 selected the slice so that the difference between voxel counts of AD and SCC becomes small. Obs. 2 selected the segmented slice arbitrarily. The following regions were not included in ROI: (1) separated or narrowly connected region; (2) GGO region; (3) invasive region out of lung field area; (4) non-tumor region like vessels, cavity and calcification to the extent possible. As seen in Supplementary Fig. S1, each observer selected a different slice for the same tumor and excluded the separated or narrowly connected regions.

Gabor filter
A two dimensional Gabor function ( , ) can be written as: ,where ′ = cos + sin , ′ = − sin + cos . and are the standard deviations of the Gaussian factor of the Gabor function and preferred wavelength, respectively.
specifies the orientation normal to the parallel stripes of a Gabor function.
Aspect ratio specifies the ellipticity of the Gabor function's support. Theγvalue of 0.5 was used. The half-response spatial frequency bandwidth (in octaves) of a Gabor filter was set to 1, in which case, and are connected as = 0.56 . The wavelength was logarithmically spaced as: where minimum wavelength = 2 and the scaling factor = 2 , respectively.
Four orientations were adopted. Sixteen Gabor filters , ( , ) were applied to the original CT image ( , ).
The convolution gives GF images:

Minkowski functionals
To remove the highly noisy data effects, 1% of the data were saturated at low and high intensities. The intensities of segmented tumor images were linearly remapped onto uniform interval 0 to 1. The simplest method of reducing gray-scale images to two-valued images (objects and background) is to use a threshold that is any number between 0 and 1. Integralgeometry morphological image analysis consists of the calculation of MFs as a function of the threshold. Larkin et al. 12 reported that thresholding steps of more than ten did not improve the classification performance. We also adopted ten thresholding steps from 0 to 1 with 0.1 intervals ( = 0,1, … ,10) to obtain 11 binary images. In Supplementary Fig. S2 The schematic illustration for the process of MFs computation is shown in Fig. 1.

Neighborhood Component Analysis
The weighted distance between two d-dimensional feature vectors and (the corresponding class labels and ) of two samples is defined as: where is a weight associated with t h feature. The probability of selects as its  (9) is differentiable with respect to w, and used to optimize the feature weights by gradient ascent method. In this study, the regularization parameter was tuned in the vicinity of 1/n (n=86) to minimize the mean squared error (MSE) of leave-one-out CV.
The features with weights of more than 0.5 were selected, although the threshold value was not optimized.

Performance evaluation
The classification performances in the original image,

Data availability
The datasets generated and analyzed during this study are available from the corresponding author on a reasonable request.

Segmentation
The voxel counts of each observer for AD and SCC are shown in Fig. 2. Obs. 2 delineated ROI for an arbitrarily selected slice of a tumor. It is natural that SCC has a significantly larger voxel count than AD in Obs. 2, because SCC has a higher rate of advanced T-stage tumor than AD. In contrast, Obs. 1 selected the slice so that the differences between voxel counts of AD and SCC become small. As a result, there was no significant difference between the voxel counts of AD and SCC. Figure 3 shows the area, perimeter, and Euler number calculated on the original image and 0,0 , respectively. GF MFs began to change at a smaller threshold than those of the original images. As the threshold increases, the increase of holes can increase the perimeter and decrease the Euler number. The increase in isolated regions increases the Euler number. At higher thresholds, the area and perimeter decrease as the number of remaining pixels decreases.

Neighborhood Component Analysis
Before NCA, each feature's values were standardized to have a zero mean and a standard deviation of one to avoid bias caused by the difference of the value range of each feature.

Performance Evaluation
Among various machine learning models in MATLAB classification learner app, 1-NN with Manhattan distance measure (MD) and Euclidean distance measure (ED) performed most favorably. The classification performance of 1-NN are summarized in Table 3. Both types of GF improved the classification performance substantially and gave approximately the same observer-averaged performance. However, the performance should be evaluated using the independent test dataset preferably.
Zhu et al. 7 reported the best performance as AUC 0.893 using the training datasets (n= 81) and the independent test datasets (n=48). Bashir et al. 9 investigated a radiomic-based and semantic feature-based model using a training data set (n=106) and an independent test data set (n=100). On training data, the radiomics-based and semantic feature-based model gave AUC=1 and 0.78, respectively, while on the test data, the semantic feature-based model outperformed the radiomics-based model (AUC=0.82 and 0.5). Thus, the reported results differ considerably; therefore reproducibility of the results and generalizability to new independent data must be verified to create a reliable radiomics-based model.
NN-based feature selection methods are unable to detect redundant features. Therefore, the correlation-based feature elimination before the feature selection stage was performed 4 .
As another approach, the two-stage feature selection method combining minimal redundancy-maximal-relevance criterion (mRMR) and other feature selection method was proposed [22][23] . In addition to NCA alone, we performed the two-stage method combining mRMR and NCA, but the performance weakened considerably in Obs. 2 (see Supplementary   Fig. S4, S5 and Table S1).
The stability of extracted features must be evaluated in radiomics study. Several studies estimated the interobserver stabilities of various features by intra-class correlation coefficient [24][25][26] . MFs were not included in those studies. In the brain tumor analysis using MFs on T2-weighted MR images, MFs showed good interobserver agreement in Bland-Altman analysis 13 . Although those methods are applied to normally distributed variables, the normality of many features in the present study was not confirmed by Jarque-Bera test. Hence, a nonparametric method was required. Aerts et al. 27 used Friedman test, which is a nonparametric repeated measurement test. Haga et al. 6 introduced an interobservervariability-based feature selection filter that selects features common to each observer. In this study, no common feature was found among ten top-ranked features of each observer. It might be because, in Obs. 2, the voxel count of SCC was significantly higher than that of AD.
To remove any feature's dependence on the voxel count of ROI, feature normalization for the area and perimeter was performed as equation (6), but it might not be enough to remove that effect. A dataset without significant difference in T-stage between histological types is required to investigate further about the robustness.

Conclusion
We proposed the novel radiomics method using 2D RI and SRI GF MFs as the radiomics features and NCA as the feature selection method for histological classification of NSCLC on nonenhanced CT images. The promising results (about 88% accuracy) was obtained, although a more detailed study about the robustness is required.