2.1.

Acquisition System

The system employed in this study is depicted in Fig. 1(a). It can automatically capture large WSHSIs. A more detailed description of the system is shown in our previous publication.³² The system consisted of an upright microscope (Olympus, Tokyo, Japan) and a customized hyperspectral snapscan camera that captures images by moving a sensor tiled with Fabry-Perot interferometers.³³ The hyperspectral camera can image from 467 to 920 nm in 150 spectra, with a full width at half maximum value ranging from 3.1 nm (487 nm) to 11 nm (900 nm). The maximum image resolution was $2048\times 2048\text{pixels}$. The stage consisted of an $X\text{-}Y$ motorized stage (Thorlabs, Newton, New Jersey, United States) and a $Z$-direction stepper motor (Prior, Rockland, Massachusetts, United States). The stage could be controlled in three spatial directions using the computer. A C++ program was developed to integrate all components together to scan hyperspectral images using the flowchart shown in Fig. 1(b).³² First, the stage and the camera identified the location of maximum focus using the Brenner filter.³⁴ Then, the program moved the $X\text{-}Y$ stage in a predetermined pattern and captured hyperspectral images.

Fig. 1

Overview of the acquisition system. (a) Photograph of the acquisition system with its components. (b) Flowchart of the automated acquisition software.

2.2.

Histological Specimen

The histological specimens consisted of 65 slides of thyroid cancer, which included three common types of thyroid carcinoma: papillary thyroid carcinoma (24 slides), follicular thyroid carcinoma (eight slides), and medullary thyroid carcinoma (three slides). The remaining 30 slides were normal non-carcinoma thyroid tissue. The specimen collection was described in our previous publication.⁶ The gross specimens were fixed, embedded, sectioned, and stained with H&E.

2.3.

Hyperspectral Dataset

Using the instruments described in Sec. 2.1, we produced a dataset of 2599 images of size $2000\times 2000\times 84\text{pixels}$, with each image being 2000 pixels in length and width and each image containing 84 spectra from 467 to 721 nm. The size of the image was predetermined by the camera manufacturer’s settings. Each image was divided into non-overlapping sub-images of size $250\times 250\times 84\text{pixels}$ (Fig. 2) and then resized into $224\times 224\times 84\text{pixels}$ to fit into pretrained networks. The size of the sub-images was chosen to best balance between enough coverage of tissue structure and good margin accuracy. An anti-aliasing resizing technique was used to not introduce resizing artifacts. A 10× magnification objective was used, so each image had a spatial resolution of $0.556\mu \mathrm{m}/\text{pixel}$.

Fig. 2

Example of hyperspectral images. Each image is $250\times 250\text{pixels}$ and corresponds to a spatial field of size $139\mu \mathrm{m}\times 139\mu \mathrm{m}$.

Before the scanning of each slide, a white reference hyperspectral image was obtained at a blank area on the slide, and a dark current image was acquired automatically by the camera with the camera shutter closed. Afterward, each raw hyperspectral image was calibrated with the white reference and dark current images to obtain the transmittance of the tissue as

2.4.

Transformer Architecture

We used the “TimeSformer” transformer architecture that was described in our previous work.³⁵ Transformers are multi-layer neural networks that use the mechanism of self-attention. In the self-attention mechanism, the input data (image, word sequence) is divided into small tokens. Each token has its own key-value pairs and its own query, so it can be matched with other tokens. The transformer network trains itself through gradient descent, so appropriate keys, values, and queries emerge. Queries can be multiplied with keys to produce a matrix of relevant weights, the result of which is retrieved as “values.” Within a transformer layer, they are represented using the attention equation, given as

Fig. 3

TimeSformer architecture modified for HSI. (a) Spatial attention compared the patch with other patches within the same wavelength band, whereas spectral attention compared the patch with other patches at the same spatial location but with different wavelengths. (b) The overall architecture for the TimeSformer network modified for HSI.

2.5.

Augmentation of Spectral Data

Algorithm 1 shows the pseudo-code implementation of RandAugment for HSI. For each image, RandAugment successively applied three random transformations from the set of transformations.²⁹ The strengths of the transformations were uniformly random values between the ranges shown in Table 1. In addition to the list of transformations seen in Faryna et al.²⁷ (identity, rotation, shearing, translation, brightness, and contrast), we proposed the following tailored transformations: spectral noise, shifting, zeroing, and NMF estimation. We detail the specifics of each transformation. Rotation rotated the image by its center clockwise anywhere from 0 to $\pi /2$. Shearing sheared the image either horizontally or vertically from 0% to 50% of the image length. Translation translated the image either horizontally or vertically from 0% to 50% of the image length. Brightness increased the intensity of the image from 0% to 100% of the image’s maximum intensity. Contrast applied the contrast adjustment equation: $I^{\prime }=1.0157(C+1)(I-0.5)/(1.0157-C)+0.5$, where $I$, $I^{\prime }$, and $C$ are the image before contrast enhancement, the image after contrast enhancement, and the contrast score, which ranges from 0 to 2.0, respectively.³⁶ Sharpen applied the unsharp masking matrix, which can be scaled to increase or decrease the sharpness.³⁷ Spectral noise applied a uniformly random shift to each image spectra that ranges from 0% to 50% of the maximum image intensity. Shifting introduced a linear deviation of the form $I^{\prime }(\lambda )=I(\lambda )\mathrm{A}(\lambda -600\mathrm{nm})$, where $I(\lambda )$ and $I^{\prime }(\lambda )$ are the spectra before and after modification that are dependent on wavelength $\lambda$, respectively, and $A$ is the strength of the modification that ranges from 0 to 2. Zeroing randomly selected up to three channels and replaced them with black image bands. This augmentation method mimicked data occlusion schemes such as Cutout.²⁵ NMF estimated the content of different tissue stains and modified them to simulate different stain concentrations.³⁶ Using NMF, we decomposed the hypercube into three components, which corresponded to hematoxylin, eosin, and hemoglobin in the stained slide. The decomposed components were modified using gamma correction $I^{\prime }=I^\gamma$, where $I$ and $I^{\prime }$ are the image before and after modification, respectively, and $\gamma$ is the gamma score that ranges from 0.5 to 1.5. We then reconstructed the hypercube from the decomposed cube. Table 1 shows the range of modifications applied for each transformation. We chose the values from a combination of experimental values and from values obtained by Faryna et al.²⁷ Figure 4 shows the effect of augmentations in the spatial domain. Figure 5 shows the effect of augmentations in the spectral domain.

Algorithm 1

Hyperspectral RandAugment.

Table 1

Ranges of different types of data augmentation used. The intensity of an image is measured on a scale from 0 to 1, where 0 is the total black and 1 is the total white.

Fig. 4

Effect of augmentations on the spatial domain of hyperspectral images. For visualization purposes, we applied transformation from HSI into RGB images using the method from Ma et al.¹⁴

Fig. 5

Effect of augmentations on the spectral domain of hyperspectral images. The spectra were generated by taking the spectrum average of the image patch. Blue dashed lines mark the original average spectra. Black solid lines mark the transformed average spectra.

2.6.

Training

To avoid data leakage, we used stratified sampling to split the data into training (28 slides), validation (seven slides), and test (26 slides) sets. Data was separated by patients, so no patients contributed to two splits. Table 2 shows the number of $250\times 250\times 84$ images used for each dataset. In the training dataset, there were 47,630 images (56% of the training data) labeled normal and 37,627 images (44% of the training data) labeled cancer. In the validation dataset, there were 10,273 images (47% of the validation data) labeled normal and 11,379 images (53% of the validation data) labeled cancerous. In the test dataset, there were 26,539 images (56% of the test data) labeled normal and 20,638 images (44% of the test data) labeled cancer.

Table 2

Number of images used for the training, validation, and test datasets.

TimeSformer networks were trained for 15 epochs using a batch size of 16 and an initial learning rate of $1\times {10}^{-3}$. We initialized the model with the pretrained weights from Bertasius et al.³⁰ The optimizer used was stochastic gradient descent with a Nesterov momentum value of 0.9 and a learning rate decay of 0.5 every five epochs. The epoch that achieved the maximum validation $F_1$ score was selected for testing. The training was implemented in PyTorch using graphic processing units on high-performance computers. We trained two TimeSformer networks, one with and one without the augmentation methods described in Sec. 2.5.

We trained three additional networks for comparison with TimeSformer: a ConvNext trained on HSI, a ConvNext trained on RGB images, and a ViT trained on RGB images. We modified the “patchify” block of ConvNext³⁸ to expand the input kernel into 84 input channels. This technique of increasing the number of input layers in pretrained CNN to fit HSI is shown to be effective in classifying H&E-stained hyperspectral data.^11,14 To train ConvNext on HSI, a batch size of 16 and an initial learning rate of $1\times {10}^{-3}$ were used. All training data were augmented using the method described in Sec. 2.5. The RGB images were synthesized from HSI using the method by Ma et al.¹⁴ To train ConvNext and ViT on RGB images, a batch size of 64 and an initial learning rate of $1\times {10}^{-3}$ were used. In all networks, we used stochastic gradient descent with a Nesterov momentum value of 0.9 and a learning rate decay of 0.5 every five epochs. For augmentation, we used the default RGB RandAugment transformations in the PyTorch library.

2.7.

Evaluation Metrics

We compared the performances of each network on the task of identifying thyroid cancer per image. The networks classified an image as either containing or not containing cancer cells based on labels provided by a clinically trained pathologist. For each model, the accuracy, the $F_1$ score, the AU-ROC, and the sensitivity score (also known as recall) were calculated. We then tested the performance of each network in identifying the cancer margin. The test dataset consisted of nine WSHSIs. To generate the cancer margin for each whole-slide image, we performed morphological opening (erosion followed by dilation) on each prediction. The metrics used to evaluate margins were the Jaccard index (intersection over union) and the Hausdorff distance. The Jaccard index $J(A,B)=|A\cap B|/|A\cup B|$ calculates the intersection between two sets $A$ and $B$ over their union. The Hausdorff distance measures the largest distances between two morphological segmentations:

2.8.

Attention Unrolling

We visualized the weights of the trained attention modules using the attention unrolling technique.¹⁹ This method aggregates all attention weights across all layers into one. At each layer $L$, a single head produced an attention map $A_{h,L}$, which was the SoftMax multiplication of the keys $K_{h,L}$ and queries $Q_{h,L}$. The attention maps were averaged across all heads and added with the identity matrix $I$ to produce the layer attention matrix ${\hat{A}}_L$. The rollout attention was calculated as the dot product of all attention maps at all layers.

3.1.

Classification Results

TimeSformer trained on the HSI dataset achieved an accuracy score of 90.87%, an $F_1$ score of 89.79%, a sensitivity score of 91.50%, and an AU-ROC score of 97.04%. Without data augmentation, TimeSformer achieved an accuracy score of 88.23%, an $F_1$ score of 86.46%, a sensitivity score of 85.53%, and an AU-ROC score of 94.94%. ConvNext trained on the HSI dataset achieved an accuracy score of 90.38%, an $F_1$ score of 88.74%, a sensitivity score of 86.27%, and an AU-ROC score of 96.23%. ViT trained on a synthesized RGB dataset achieved an accuracy score of 89.98%, an $F_1$ score of 88.14%, a sensitivity score of 84.77%, and an AU-ROC score of 96.17%. ConvNext trained on a synthesized RGB dataset achieved an accuracy score of 89.31%, an $F_1$ score of 87.84%, a sensitivity score of 87.91%, and an AU-ROC score of 95.63%. Table 3 shows that TimeSformer achieved the highest accuracy, $F_1$ score, sensitivity, and AU-ROC score of all networks on the test dataset. The results achieved by TimeSformer on the thyroid dataset are greater than those achieved by Halicek et al.⁶ on the same set of H&E-stained thyroid slides, which achieved an accuracy of 89.4% and an AU-ROC score of 95.4%.

Table 3

Accuracy, F1 score, sensitivity, and AU-ROC scores for the test dataset.

3.2.

Cancer Margin Detection in Whole-Slide Images

Out of the four networks tested, TimeSformer achieved the overall best average Hausdorff distance of $0.76\mathrm{mm}\pm 0.33\mathrm{mm}$, with a maximum of 1.47 mm and a minimum of 0.44 mm. In comparison, ConvNext trained on HSI achieved an average Hausdorff distance of $1.07\pm 0.73\mathrm{mm}$, ViT trained on RGB achieved an average of $0.99\pm 0.64\mathrm{mm}$, and ConvNext trained on RGB achieved an average of $1.08\pm 0.52\mathrm{mm}$. The results showed the ViT performs better for hyperspectral images. Table 4 shows the Hausdorff distance achieved on each slide. The table shows that typically slides that achieved good performances in TimeSformer also achieved good performances on other networks, and vice versa. The table also shows that TimeSformer produced good margins more consistently compared with other networks, as evidenced by the lower standard deviation value.

Table 4

Performance in terms of the Hausdorff distance on each slide. Data in bold refer to the network with the best performance.

Figure 6 shows the carcinoma margins in four representative slides. Figure 6(a) depicts slide 7, which contains papillary thyroid carcinoma. In this slide, the network achieved a Hausdorff distance of 0.67 mm and a Jaccard score (intersection over union) of 0.95. Figure 6(b) depicts slide 1, which contains follicular thyroid carcinoma. In this slide, the network achieved a Hausdorff distance of 0.65 mm and a Jaccard score of 0.88. Figure 6(c) depicts slide 2, which is a follicular carcinoma slide and the slide with the worst performance. Examination of the slide reveals that the classifier associated some fibrous structures with normal tissues, which might not be accurate for this particular instance. Figure 6(d) depicts slide 5, which is a papillary carcinoma slide and the slide with the best performance. Examination of the slide reveals that other networks also achieved good results on slide 5 and that slide 5 has a very clear delineation between carcinoma and normal tissues.

Fig. 6

Cancer margins predicted on four thyroid whole slides. The black border shows the predictions made by the pathologist, and the red border shows the predictions made by TimeSformer. Whole-slide images were captured by a whole-slide histology scanner. (a) Whole-slide image of thyroid diagnosed with papillary carcinoma. (b) Whole-slide image of thyroid diagnosed with follicular carcinoma. (c) Whole-slide image of thyroid diagnosed with follicular carcinoma. (d) Whole-slide image of thyroid diagnosed with papillary carcinoma. The RGB images were taken using a whole-slide histology scanner.

3.3.

Visualization of Attention in Transformer

Figure 7 shows the spatial attention map of the trained network on an image of papillary thyroid carcinoma, for which the network produced an accurate prediction with a confidence of 96.51%. It could be seen that spatial attention was highest among the clusters of nuclei in the middle section; this attention was consistent across all wavelengths and highest in the blue-green wavelengths (470 to 540 nm). We attribute this to the fact that hematoxylin stains the nucleus blue, so nuclei were highlighted for the ViT. Our findings were consistent with the findings of Halicek et al.,⁶ which showed that deep neural networks focus on nuclei to predict the probability of cancer.

Fig. 7

Attention map of an image that showcased papillary thyroid carcinoma. The RGB image was taken using a whole-slide histology scanner.

Figure 8 shows the attention map of the trained network on a patch of healthy and normal thyroid cells. The network produced a correct prediction with 99.53% confidence. It could be seen that spatial attention was highest in the regions of follicular colloid, which are the fluid-filled sections of the thyroid follicle. Less attention was paid to the surrounding follicular thyroid cells. Spatial attention was consistent across all wavelengths and highest in the blue-green wavelength bands (470 to 540 nm).

Fig. 8

Attention map of an image that showcased regular healthy thyroid. The RGB image was taken using a whole-slide histology scanner.