2.1.

Hyperspectral Imaging Instrumentation

An in vivo small animal imaging system called Maestro (PerkinElmer Inc., Waltham, Massachusetts) was used for acquiring the hyperspectral dataset. This system mainly consists of the light source, wavelength dispersion device, and area detector. A Cermax-type, 300-W xenon light source is used for white-light excitation, which spans the electromagnetic spectrum from 450 to 800 nm. The interior near-IR light source can be used to extend the spectral range of excitation up to 950 nm. Four fiber optic, adjustable illuminator arms yield an even light distribution to the subject. The light from the excitation source illuminates the subject. The reflected light from the surface is split into a series of narrow spectral bands by liquid crystal tunable filters with a bandwidth of 20 nm and collected by a 12-bit, high-resolution charge-coupled device. The hyperspectral camera can simultaneously acquire a full dataset from as many as three mice, with a spatial resolution of $25\mu \mathrm{m}/\text{pixel}$, in only a few seconds. The acquisition wavelength region ranges from 450 to 950 nm, with varying step sizes such as 2, 5, or 10 nm.

2.2.

Animal Imaging and Tumor Surgical Experiments

A head and neck tumor xenograft model using head and neck squamous cell carcinoma (HNSCC) cell line M4E with green fluorescence protein (GFP)²² was adopted in the experiment. The HNSCC M4E cells were maintained as a monolayer culture in Dulbecco’s modified Eagle’s medium/F12 medium (1:1) supplemented with 10% fetal bovine serum.¹⁴ M4E-GFP cells, which are generated by transfection of the pLVTHM vector into M4E cells, were maintained in the same condition as M4E cells. The animal experiment was approved by the Animal Care and Use Committee (IACUC) of Emory University.

In this experiment, eight female nude mice aged 4 to 6 weeks were injected with $2\times {10}^6$ M4E cells with GFP on the lower back. Surgery was performed $\sim 1$ month after tumor cell injection. Before surgery, mice were anesthetized with a continuous supply of 2% isoflurane in oxygen. After the anesthesia administration, the skin covering the tumor was removed to expose the tumor to simulate a surgical situation.

2.3.

Evaluation Methods for Cancer Detection

In this study, GFP fluorescence images provide the in vivo gold standard for tumor margin delineation. Although the current gold standard for cancer diagnosis remains histological assessment of H&E stained tissue, the ex vivo tissue specimen undergoes deformations, including shrinkage, tearing, and distortion, which makes it difficult to align the ex vivo gold standard with in vivo tumor tissue. However, the in vivo GFP fluorescence images provided a much better alignment with hyperspectral reflectance images in the animal model, since they were acquired in vivo immediately after the acquisition of reflectance images for each mouse while the tumor position and shape in the reflectance images are exactly the same as those in the fluorescence images. The tumor and surrounding normal tissue exhibited high contrast in the GFP images. Since histopathology is the gold standard for cancer diagnosis in clinics, we also acquired the H&E stained histological image as the ex vivo gold standard.

In this paper, tumor regions were identified manually on the GFP images, and the classification results were then compared with the manual maps. Since human tissue does not contain GFP naturally, registration methods are desirable to align the in vivo hyperspectral images with ex vivo histological images as discussed in our previous publication²³ in order to move forward for future human studies.

3.1.

Preprocessing Methods for Hyperspectral Images

The preprocessing steps of the hypercube consist of data normalization, image registration, glare detection and removal, and curvature correction.

3.1.3.

Glare removal method for hyperspectral images

Glare, also called specular reflection, is the mirror-like reflection of incident light from a moist surface. Optical images such as those from an endoscope and colposcopy images acquired during surgery are often strongly affected by glare spots in the images, which presents a major problem for surgical image analysis.²⁴ In HSI, glare alters the intensities of the pixels in each image band and consequently changes the spectral fingerprint, which could introduce artifacts in feature extraction and hence deteriorate classification results. Therefore, glare pixels need to be detected and removed before feature extraction and classification.

Currently, no glare detection technique exists for intraoperative hyperspectral images. To better understand the difference of glare pixels, we compared the spectral fingerprint of glare pixels with nonglare pixels and identified the characteristics of glare pixels. It can be seen from Fig. 2 that glare pixels not only had higher intensities but also showed larger variations along the spectral bands. Therefore, the first-order derivatives of the spectral curves of glare pixels were much larger than those of nonglare pixels, as shown in Fig. 2(d).

Fig. 2

Rationale for the proposed glare detection method. (a) Image band at 758 nm. (b) Enlarged image of a selected glare region in the image band in (a). Glare pixels are much brighter than nonglare pixels. (c) Normalized reflectance curve of glare pixels (G1 to G3) and nonglare pixels (NG1 to NG3). Spectral curve of glare pixels varies significantly in many wavelengths. (d) First-order derivative curves corresponding to the spectral curves shown in (c).

Based on the above observations, we proposed a three-step glare detection method: (1) estimate the first-order derivatives of spectral curves with a forward difference method. (2) Calculate the standard deviation (std) of each derivative curve and generate an std image for each hypercube. Glare pixels show higher std than normal pixels. (3) Compute the intensity histogram of each std image, fit the histogram with a “log-logistic” distribution (MATLAB function), and experimentally identify the threshold that separates glare and nonglare pixels.

3.2.

Feature Extraction from Hyperspectral Data

Feature extraction and representation is a crucial step for image classification tasks. Efficient feature extraction could lead to improved classification performance. However, HSI has only recently been applied to medical applications, and it is not well understood what features are the most effective and efficient to differentiate cancer from normal tissue in medical hyperspectral images. Currently, one of the most frequently used features is the normalized reflectance curve of each pixel. With only the reflectance features, classification performance is still far from optimal. Therefore, there is still space for adding new features to improve the distinction of cancer from normal tissue.

In this study, we explored the utility of several spectral features, which are derived from the reflectance curve of each pixel. Spatial features are not explored in this study but will be included in our future work. The extracted spectral features included (1) first-order derivatives of each spectral curve, which reflect the variations of spectral information across the wavelength range; (2) second-order derivatives of each spectral curve, which reflect the concavity of the spectral curve; (3) mean, std, and total reflectance at each pixel, which summarize the statistical characteristics of the spectral fingerprint; and (4) Fourier coefficients (FCs), which were initially found to be effective for target detection in the remote sensing field²⁶ and later were adopted for breast cancer margin classification from ex-vivo breast cancer hyperspectral images.²⁷

FC feature extraction involves transforming the original spectral $f(n)$, $n=\mathrm{0,2},\dots N-1$ into the Fourier domain as $F(k)=\sum _{n=0}^{N-1}s(n)e^{-\frac{2\pi i}{N}kn}$ and then combining the selected real and imaginary components of FCs $F(k)$ using the following formula:

Different features may have very different numerical ranges, so each feature was standardized into its z-score (MATLAB function) by subtracting the mean from each feature and then dividing by its std.

3.3.

Feature Selection Method

After feature extraction, the spectral feature dimension was increased from 226 to 904. Such a high spectral dimension poses significant challenges to the analysis of the hypercube. High dimensionality can significantly increase the computational burden and storage space, leading to increased data processing time, which is against the requirement of real-time tumor detection during surgery. Depending upon the wavelength range of imaging systems, different studies may have different reflectance features. It is not clear which wavelength is more relevant in characterizing cancerous tissue and would provide a better contrast between cancer and normal tissue. In addition, with a narrow wavelength increment, there is likely spectral redundancy between adjacent bands. Last but not least, increasing the feature dimensionality without increasing the number of training samples may lead to a decrease in classification performance due to the curse of dimensionality, that is, the Hughes phenomenon.²⁸ Therefore, it is desirable to analyze the spectral redundancy in the high dimensional data and select the most characterizing compact feature set.

The goal of the feature selection is to find a feature set $S$ with $n$ wavelengths $\{{\lambda }_i\}$, which optimally characterize the difference between cancerous and normal tissue. To achieve the optimal condition, we used the maximal relevance and minimal redundancy (mRMR)²⁹ framework to maximize the dependency of each spectral feature on the target class labels (tumor or normal) and minimize the redundancy among individual features simultaneously. Relevance is characterized by mutual information $\mathrm{I}(x;y)$, which measures the level of similarity between two random variables $x$ and $y$:

We represent each pixel with M features $\mathrm{\Lambda }=\{{\lambda }_{\mathrm{i}},i=1,\dots M\}$, $M=904$, and the class label (tumor or normal) with c. Then the mR condition is to search features satisfying Eq. (3) which maximize the mean value of all mutual information values between individual features ${\lambda }_i$ and class $c$:

The features selected by the mR condition are likely to have redundancy, which means that the dependency among these features could be large. When two features highly depend on each other, the respective class-discriminative power would not change much if one of them were removed. So, the minimal redundancy condition can be added to select mutually exclusive features:

The simple combination [Eqs. (5) and (6)] of these two conditions forms the criterion mRMR, which can optimize $D$ and $R$ simultaneously:

In practice, incremental search methods can be used to find the near-optimal features defined by $\mathrm{\Phi }(·)$. Suppose we already identified a feature set $S_{m-1}$ with $m-1$ features. The task is to select the $m$th feature from the set $\{\mathrm{\Lambda }-S_{m-1}\}$. This is done by selecting the feature that maximizes $\mathrm{\Phi }(·)$. The respective incremental algorithm optimizes the following condition:

3.4.

Hyperspectral Image Classification

To determine the optimal feature set, a K-nearest neighbor (KNN) classifier was employed to evaluate the effectiveness of the selected features with leave-one-out cross validation. The workflow of the feature selection and classification was shown in Fig. 3. First, eight hypercubes were partitioned into a training dataset with seven hypercubes and a testing dataset with the rest of the hypercube. Next, 904 features were extracted from the training data, and mRMR was used to select the optimal feature set on the training dataset as the feature number varied from 1 to 904. The selected feature set was applied to the testing dataset correspondingly. Next, the KNN classifier was used to train the training data and to predict the labels of the testing dataset with feature sets of different sizes. Finally, the feature set that gave the best classification performance was chosen as the optimal feature set for differentiating cancer from normal tissue.

Fig. 3

Flowchart for feature selection and classification.

Accuracy, sensitivity, specificity, and F-score, as defined in Refs. 30 and 31, were chosen as the metrics to evaluate the classification performance of the features. Accuracy represents the percentage of the correctly detected tumor and normal pixels relative to the total number of tumor and normal pixels in an image, respectively. Sensitivity represents the percentage of correctly detected tumor pixels relative to the total number of tumor pixels in an image. Specificity is the percentage of correctly identified nontumor pixels relative to the total number of nontumor pixels in an image. F-score is the harmonic mean of precision and sensitivity.

4.1.

Results on Glare Detection and Removal

Figure 4 shows one example of glare detection results. The bright area in the std image represents the location of glare pixels, which were characterized by large spectral variations and distributed mostly along the long tail of the histogram. The key issue was to identify an appropriate threshold on the histogram that could separate the glare pixels from nonglare pixels in the std images. We compared the traditional thresholding method such as the Ostu method³² and entropy method³³ with the proposed loglogistic curve fitting method. The Ostu and entropy methods were most suitable for histograms with a bimodal shape, but the histogram here contained a very long tail, which was not in the typical bimodal shape. The curve fitting methods produced much a lower threshold than the Ostu and entropy methods, which enabled the detection of the relatively bright glare margins as well as the isolated glare pixels. The threshold was set to be the intensity value, which yielded a certain percentage ($ϵ$) of the peak value in the fitted loglogistic distribution curve. The value of $ϵ$ was experimentally set to 5% through trial and error, which was found sufficient to detect most of the glare pixels. After the glare masks were generated, glare pixels along all the spectral bands were removed from the training and testing datasets, since they did not contain useful diagnostic information.

Fig. 4

Glare detection results: (a) Standard deviation (std) image of the first order derivative for a hypercube. (b) Binary glare map generated by the classical Otsu method. (c) Binary glare map generated by the entropy method. (d)–(f) Glare map generated by the proposed method with ratios 0.01, 0.05, and 0.1. (g) Histogram of the std image with blue color and loglogistic fitting curve with red color. The five vertical lines represent the five thresholds generated by the five methods in (b)–(f).

4.2.

Results on the Comparison Between Green Fluorescent Protein and Non-Green Fluorescent Protein Images

Figures 5 and 6 show the exposed tumors under the white excitation with and without GFP, respectively. GFP emission peaks under blue excitation occur at the wavelengths of 508 and 510 nm. It was found that under white excitation, the spectral images at these two bands did not exhibit enhanced contrast between tumor and normal tissue compared to other spectral bands. This is consistent with the observation in our previous study.³⁰ Therefore, we did not remove GFP spectral bands at 508 and 510 nm in the preprocessing steps.

Fig. 5

Visualization of a tumor with green fluorescence protein (GFP). The upper part is the image of the tumor in a mirror. (a) RGB composite image of the hypercube. (b) Preprocessed spectral images at wavelengths 450 nm, 508 nm, 510 nm, 542 nm, 554 nm, 576 nm, 600 nm, and 650 nm. (c) Spectral curve of cancerous and healthy tissue.

Fig. 6

Visualization of a tumor without GFP. (a) RGB composite image of the hypercube. (b) Preprocessed spectral images at wavelengths 450 nm, 508 nm, 510 nm, 542 nm, 554 nm, 576 nm, 600 nm, and 650 nm. (c) Spectral curve of cancerous and healthy tissue.

Both figures demonstrated that hyperspectral imaging can probe vessels at different depths below the visual surface of the tumor. Light at a shorter wavelength region of the visible range is more sensitive to superficial vascular information due to limited light penetration into the tissue. As the wavelength becomes longer, information from deeper tissue can be acquired. Hence, changing the illumination wavelength may enhance vascular contrast and allow visualization of angiogenesis at the tumor region.

The RGB color image shows the highly vascularized tumors. Figures 5 and 6(c) illustrate the characteristic spectrum of hemoglobin at 542 and 577 nm. These characteristics may contribute to the distinction between cancerous and normal tissue by HSI.

4.3.

Results on Feature Extraction and Visualization

The most commonly utilized feature for cancer detection with HSI is the normalized reflectance spectra, which reflects the physiological and pathophysiological state of tissue at each pixel. However, this feature alone may not have enough discriminative information to minimize the classification error. It would be very interesting to explore the usefulness of other features besides reflectance. We derived a series of features based on the spectral curve of each pixel and boosted the original feature dimension from 226 to 904.

Figure 7 shows an RGB composite color image of an example hypercube. The tumor exhibited a white necrotic appearance, which was confirmed by the histological image in the interface between viable and necrotic tumor tissue. The reflectance spectrum of a tumor with necrosis differs from that of a tumor without necrosis. The image also shows that the pre-existing blood vessels in normal tissue reached the tumor region and grew new blood vessels into cancerous tissue. The highly vascularized tumor region grew well due to adequate oxygen and nutrient supply and due to removal of waste through blood vessels, while the nonvascularized tumor region became necrotic due to the lack of vessels. These observations were consistent with their spectra curves (Fig. 7). The viable cancer tissue is heterogeneous with higher spectra variations compared to the necrotic tissue. The average reflectance spectra of the viable tumor were lower than those of the normal tissue, which reflected the higher amount of hemoglobin in the highly vascularized cancer tissue.

Fig. 7

Reflectance spectral curve of a tumor with necrosis. (a) RGB composite image of hypercube. The white region looks necrotic, and the other part of the tumor contains many vessels; (b) Histological image of the rectangular tissue region in (a). The upper part is the necrotic tissue without nuclei, and the lower part is the viable cancerous tissue. (c) The average reflectance spectra of the tumor, necrosis, and normal tissue with std. The red solid line represents the average spectra of cancerous tissue, and the blue dotted line represents the average spectra of the normal tissue. The green dashed line represents the average spectra of the necrotic tissue. The error bars are the std at a certain wavelength of the three curves.

Figure 8 shows the visualization of all the spectral features explored in this study, which include the RGB composite image of hypercube, and the mean, std, and sum of the selected tumor and normal tissue regions of interest (ROI). The texture of the tissue could be clearly visualized from the std image. The mean and sum image only differ by a constant. Due to the large dynamic range of the spectra along the wavelength, the average of the hypercube in the ROI did not reflect the correlation with the RGB image. All these features captured different aspects of the differences between tumor and normal tissue, which are visualized in these figures.

Fig. 8

Feature extraction and visualization. (a) RGB composite image of hypercube, mean, std, and sum of the selected tumor and normal tissue ROI. (b)–(e) Average normalized reflectance curve, first derivative, second derivative, and the difference of FCS between normal and tumor tissue in the selected ROI in (a).

4.4.

Feature Selection and Classification Results

The objective of the feature selection is to identify the features which are critical to minimize the classification error. The mRMR feature selection method selected a compact feature set with the mR to the target class and the MR within the feature set. Figure 9 shows the mutual information of the extracted spectral features with the class labels (tumor or normal), which reflects the relevance of each feature with respect to the class labels. The highest mutual information was achieved by FCs. Normalized reflectance above 850 nm showed higher relevance than other wavelengths. Derivatives, mean, std, and sum generally exhibited lower mutual information with class labels.

Fig. 9

Mutual information between features and class labels. The $x$ axis represents the feature number, and the $y$ axis represents the mutual information.

Figure 10 displays visualization of the mutual information among individual features, which represents the redundancy within the 904 features. It can be seen clearly that the lower left square of size $226\times 226$ shows relatively higher mutual information than other regions, which demonstrates that wavelength features were highly correlated with each other. Therefore, the normalized intensity values across the wavelength range contain complementary information as well as redundancies.

Fig. 10

Mutual information between individual features. Color bar on the right shows the color map corresponds to the value of mutual information. Higher mutual information indicates more redundancy between features.

Figure 11 shows the evaluation of feature sets of different sizes by supervised classification. The metrics initially increased with the feature number, reached a maximum, and then decreased as the feature set went to its maximum size of 904. Classification with the full feature set of dimension 904 was not as good as the feature set of dimension 20. We found that a feature set of size 20 gave the best performance, with an average accuracy, sensitivity, and specificity of 67.2%, 77.5%, and 54.0%, respectively. In this experiment, classification was used only for evaluating and comparing the effectiveness of feature sets of different sizes. Since the focus of this study was not on classification, we will try to improve classification in the future.

Fig. 11

Feature selection and classification.

As shown in Fig. 11, we have run mRMR with cross validation to select the optimal feature set $\mathrm{F}={\{f_i\}}_{i=\mathrm{1,2},\dots m}$ of size $m$ ($m=1,3,5,\dots ,750,904$) from 904 features. We have generated a series of feature rankings, $\{f_1\},\{f_1,f_2,f_3\},\dots ,\{f_1,f_2,\dots ,f_{750}\},\{f_1,f_2,\dots ,f_{904}\}$, where $f_i$ represents the feature that has been ranked in the $i$th position of each optimal feature set. It should be noted that $f_i$ could be different on different cross validation folds or different among optimal feature sets of different sizes. Therefore, we can analyze the composition of the $i$th ranked feature $f_i$ and count the ranking frequency of individual features. Figure 12 summarizes the ranking frequency of different feature types. Each bar represents the normalized frequency of different feature types being selected as the $i$th ranked feature $(i=1,2,\dots ,20)$. FC is the only selected top ranking feature type, normalized reflectance is the only feature type that ranked the second, and the mean spectrum is the third-ranked feature type. Spectral derivatives ranked fourth and fifth even though their mutual information with corresponding labels was low. This could be explained by the fact that mRMR selects not only the features with high mutual information but also with low redundancy within the feature set.

Fig. 12

Feature ranking. The $x$ axis is the ranking from 1 to 20, and the $y$ axis is the percentage of the selection frequency for different features on each rank. Each bar represents the normalized frequency of different feature types being selected as the $i$th ranked feature.