2.1.

HSI Acquisition System

A custom acquisition system was designed, with a 2D spectroradiometer at its core. The system was developed and calibrated specifically for applications in skin gross pathology.²³ An HSI camera (SR-5000, TOPCON Corporation, Japan; 1.4 Mpixel CCD image sensor, lens 32 mm, $1376\times 1024\text{pixel}$, 380 to 780 nm) was positioned above the capture base, and an LED source (SERIC Solax-iO LE-9ND55, SERIC, Japan) was arranged at 45 deg from the capture base. Two polarizers (TS High Contrast Plastic Polarizing Plate, Edmund Optics, Japan) were positioned in front of the lens and the light source, in a Crossed Nicols configuration. The system was enclosed in a dark box. The HSI data cube had spectral range (420, 750) nm with step 1 nm, in accordance with the range limits of the polarizers. The spatial dimensions varied per target sample, thus the largest data cube had dimensions $500\times 800\times 311$, corresponding to an 8 cm long tissue sample. The HSI system was integrated into the medical practice and operated by pathology experts (TI and MH). The operator could choose an appropriate region of interest to be imaged and the capture duration was adjusted accordingly. Data capture and postprocessing lasted 2 to 3 min. Additionally, a white reference surface (X-Rite ColorChecker White Balance sheet) and dark signal (lights off state) were imaged with the same system after the capture of the tissue. This resulted in three HSI data cubes per sample, the tissue image $I_t$, the white image $I_w$, and the dark image $I_d$. A black porous material was used as the capture surface, to avoid leaks and contamination of the camera system.

2.2.

Experimental Dataset

An original HSI dataset was prepared using the custom acquisition system. The experiments in this study were approved by the Ethical Committee of Saitama Medical University (977). All participants gave informed consent for the scientific use of their data. Pigmented skin lesions (PSL) were selected as target pathologies for this study. Lesions that fit the conditions were imaged immediately after surgical excision, before any kind of processing. Tissue samples were imaged with the surface facing up on the capture board. While this process constitutes ex-vivo imaging, the fast acquisition ($\le 30\min$ after surgery) is analogous to in-vivo imaging. Each sample was formalin-fixed and cross-sectioned two days later. The cross-sections were captured both with the HSI system and conventional RGB cameras. The sectioning direction and positioning of cross samples were recorded on the RGB image. Afterward, samples were pressed into paraffin blocks, shaved, and converted into tissue slides, for histopathological evaluation. Next, histological margins were traced on the images of sectioned tissue by expert pathologists (TI and MH). These margins were then traced on the HSI of the raw tissue surface. This resulted in 2D masks of tumor margin that coincided with the spatial dimensions of the HSI. Additionally, a full diagnosis was provided, including disease name and malignancy status (malignant or benign).

The data collection process resulted in the HSI of 32 tissue samples of PSL. Four of the samples were not used for the following reasons. The sample of Paget’s disease had a very different structure compared to the rest of the PSL. Two samples of melanocytic nevus were discarded due to heavy staining with blue dye. A basal cell carcinoma (BCC) sample was discarded due to difficulty in creating the label mask. This resulted in 28 tissue samples from 21 patients. An HSI data cube, a tumor segment mask, disease name, and malignancy status were recovered for each sample. Table 1 shows a detailed breakdown of the dataset. A total of 384,014 spectral signatures were collected for validation, 33% of which were within the tumor margin detected by histology. Another set of 283,206 spectral signatures was collected in the same manner and used for testing. Melanocytic Nevus is closely associated with Malignant Melanoma in terms of structure. These tumors differ from cancer tumors like BCC, due to the absence of keratin in the epithelial layer. Nonetheless, both types of tumor can potentially evolve into malignancy.

Table 1

Breakdown of skin tumors included in the dataset.

2.3.

Data Preprocessing

HSI data cubes were preprocessed individually as follows. Normalized reflectance $r(\lambda )$ for a tissue sample was recovered from each pixel as in

2.4.

Use Case Scenarios

Two frameworks were proposed, one for single-pixel processing using machine learning and one for patch-wise processing using deep learning. The reason for the two frameworks is that they are expected to be used in different scenarios. A single-pixel framework can be used when a small number of images/large number of spectral signatures is available and fast system development is required. It can also be used for the initial stages of system development for exploratory analysis. Moreover, a model based on well-defined feature extractors is more easily interpretable and the decision process can be associated with pathophysiological characteristics. On the other hand, previous experiments showed that there is a benefit in using the joint spatio-spectral information of HSI data cubes. Moreover, deep learning classifiers that perform abstract feature extraction have shown exceptional results in applications in medicine. For this reason, the second framework has deep learning at its core. Due to the large size of HSI data cubes and to fit memory and time constraints, a patch-based approach is suggested. This framework is suitable for cases where a large dataset with accurate labels is available. The two frameworks are presented in Fig. 1. The available models are presented in Fig. 2. Each framework is described in detail in the following subsections.

Fig. 1

The proposed frameworks for skin tumor segmentation. (a) The pixel-wise framework and (b) patch-based framework have different models at their core. Processes in blocks with vivid blue color are the same for both frameworks.

Fig. 2

The models that were implemented in this study. Both machine learning and deep learning models are used, each with the respective framework.

Fig. 3

The database of end members used for the calculation of abundance maps.

Fig. 4

The database of reference signatures used for SAM calculation.

2.5.

Pixel-wise Framework

The pixel-wise framework in Fig. 1 receives as input the HSI data cube. Pixel-wise normalization is applied and the foreground signatures are extracted based on the foreground mask, after background removal. The signatures from multiple HSI data cubes are reshaped into a 2D matrix, where each row represents the foreground pixels and columns represent the spectral dimension. Nontissue pixels were ignored because their inclusion affects model training. Then, feature extraction is performed. Although there is a wide variety of feature extractors, a simple approach was assessed. Two sets of features were prepared (a) reflectance signatures and (b) abundance maps.²⁷ The former refers to the spectral signatures after normalization. Abundance maps were calculated using eight endmembers from a basal cell carcinoma sample. This sample was selected because it contained a variety of colors, textures, and underlying structures. In essence, endmembers represent “pure pixels” or individual substances in the spectral profile of an image target. Endmembers were calculated by the N-Findr algorithm and eight abundance maps were produced, respectively.²⁸ The endmembers are shown in Fig. 3. The signature vector was 311-point long, whereas the abundance vector was eight-point long. Then, the feature vectors are fed into a support vector machine (SVM) classifier. This resulted in two schemes: $\text{Signature}+\mathrm{SVM}$ and $\text{Abundance}+\mathrm{SVM}$. On datasets where the two classes have a small amount of overlap and are a bit unbalanced, like this skin HSI dataset, SVM tends to perform well. The SVM model produces a prediction of tumor or normal tissue. The predicted values are reshaped back to the original spatial dimensions of the input images. For comparative purposes, an unsupervised method based on image processing is also implemented. For $\mathrm{Kmeans}+\mathrm{SAM}$, the HSI signatures were grouped into seven clusters with Kmeans clustering and then grouped based on SAM similarity.⁸ Six signatures (three of normal tissue and three of PSL tissue) were used as references for SAM calculation, as shown in Fig. 4. MATLAB^® 2020b (The MathWorks Inc., Natick, Massachusetts) was used for implementation and experimentation.

2.6.

Patch-based Framework

For the patch-based framework in Fig. 1, the initial stages of data handling are the same as before. However, in this case, the normalized HSI data cube is split into $32\times 32\times 311$ patches. Deep models need a lot of memory to work, and HSI data only exacerbate this problem. By splitting data into patches, implementation in reasonable memory levels can be achieved. The patches are reshaped into a 4D matrix of batches, where the dimensions represent (batch number, height, width, and spectrum). The data are fed to a deep-learning model that performs both abstract feature extraction and segmentation. The model produces segmentation results for each patch. The predictions are stitched back together to form the original images.

A deep learning model based on 3D convolutions is proposed. 3D Convolution has been used in other segmentation problems in medicine such as MRI diagnosis. The advantage of 3D convolution compared to 2D convolution is the direction of the kernel, which can move along the spectral direction with independent steps. Thus, the kernel for convolution can have different values in the spatial directions and the spectral direction. This way, the information in the spectral direction is preserved and processed individually. The 3D Xception model is based on Xception architecture²⁹ with 3D Convolutions instead of 2D. Other deep models, such as Inception need a lot of data to be trained effectively, but Xception can reduce the parameters of the model. Xception architecture outperforms other state-of-the-art (SOTA) segmentation models in problems with natural images. Its feature is that it uses depthwise convolutions instead of common convolutions. Depthwise convolution is performed in two stages to decouple the process. The structure of the proposed 3D Xception is shown in Fig. 5. It preserves the 3D structure in the entire part of the decoder, which enables independent adjustment of the kernel size for convolutions in the spatial and spectral dimensions. A too short kernel size in the spectral dimension could be insufficient to detect patterns in the spectral signatures. At the same time, a too wide size could obscure such patterns. The spectral dimension is dropped at the bottleneck layer. Two different approaches to dropping the 3D at the bottleneck layer were proposed: applying (a) the max function at the spectral dimension or (b) the mean function at the spectral dimension. Dropout layers with a 40% rate were introduced at the end of each Xception block. An additional convolutional neural network based on 3D convolutions (3D CNN) was designed for baseline comparison. The model was structured with a U-Net backbone as shown in Fig. 6, which is commonly used in the medical field for segmentation problems.¹³

Fig. 5

The architecture of the proposed 3D Xception network.

Fig. 6

The architecture of the 3D CNN network.

The proposed models were compared with various SOTA models. Visual Geometry Group (VGG) and Inception networks showed poor performance in initial experiments, so only the residual network (Resnet) architecture was considered moving forward. A precompiled Resnet34 model^30,31 was used in two different configurations, as shown in Fig. 7. The first configuration is pretrained Resnet. Because Resnet is designed to accept three-channel inputs, an additional convolutional layer is added before the input layer to reduce the 311 spectral channels into 3. Only the input convolution weights are trained, whereas the rest of the parameters are frozen to weights trained on the ImageNet dataset. Alternatively, a second configuration is prepared, referred to as custom Resnet. The input dimensions can be passed as a parameter to the model’s input layer, the spectral dimension is reduced with 2D convolutions, and weight training for all layers is performed from scratch. The models were built and tested in an Anaconda (Anaconda, Inc., Austin, Texas) environment, using Python 3.9, Tensorflow 2.7, and a system with NVIDIA RTX6000 GPU.

Fig. 7

The architecture of the two Resnet-based networks.

2.7.

Training and Evaluation

Three different data splits were used for experimentation, according to the flowchart in Fig. 8. The input dataset was split into validation and testing datasets. Validation data were split in two different ways. The first data split (data split 1) consisted of 16 training samples and 3 testing samples as described in Table 2. Data split 1 was used for the initial system building and model optimization. Validation data were separated into folds, with each fold consisting of the data derived from a single tissue sample. This structured dataset, data split 2, was used for leave-one-out cross-validation (LOOCV) with 19 folds in total. The training was performed on eighteen folds and performance was evaluated on the last. The training/testing process was repeated until each fold was used for testing once. Average performance was derived from each iteration’s test performance. Finally, the testing data and data split 3 were used for final testing and visual confirmation of the produced segmentation results.

Fig. 8

The flowchart of system building, optimization, cross validation, and testing.

Table 2

Data splits.

The models in this study were minimally optimized because this is an initial proof of work. The parameters of each SVM (kernel, box constraint, and kernel scale) were optimized with Bayesian optimization. Kernel functions were limited to linear and radial base function (RBF) kernels. Only 10% of available signatures (sampled at equal intervals) were used for training, to avoid a long training duration of the SVM. For the per-patch framework, data augmentation resulted in 1780 HSI cubes for the validation dataset. The deep learning models were optimized using grid search for the following hyperparameters: Optimizer among Adam, RMSProp, learning rate among ${10}^{-3}$, ${10}^{-4}$, ${10}^{-5}$, ${10}^{-6}$, exponential decay among none, ${10}^{-5}$, ${10}^{-6}$, loss function among binary cross entropy (BCE), $\text{binary cross entropy}+\text{Jaccard coefficient}$ loss ($\mathrm{BCE}+\mathrm{JC}$). BCE is the most common loss for segmentation problems. However, in cases where the classes are unbalanced, model training converges at a good accuracy for the prevalent class, instead of focusing on both classes. To alleviate this problem, JC loss has been proposed.^17,32 The losses are defined using the predicted value PR and the ground truth value GT as follows:

Accuracy, sensitivity, and specificity metrics were used for evaluation. These metrics are commonly used in classification problems but are lacking when evaluating performance concurrently in the two spatial dimensions. Accuracy expresses the ratio of correct prediction. Sensitivity expresses the ratio of positive cases being detected. Specificity expresses the ratio of a negative prediction being correct. Sensitivity and specificity are commonly compromised with each other because in a medical setting no positive case should be missed, but high false positives induce an unnecessary burden to the patient. Additionally, Jaccard coefficient (JC), also known as intersection-over-union was used as the main performance metric. JC is more appropriate for segmentation tasks because its calculation is unaffected by a small number of true positives or true negatives. The performance metrics are calculated as:

3.1.

System Building

Due to the high-dimensionality of the HSI data space, the validation dataset is visualized using t-distributed stochastic neighbor embedding (t-SNE). This technique uses a distance metric to reduce each spectral signature to a 2D or 3D embedding, and it is well suited for high-dimensional spaces. Two-dimensional embeddings and euclidean distance were used in this case. The embeddings of the labeled dataset are displayed in Fig. 9. The two-dimensional embedding shows a great overlap among pathologies, tumor tissue and non-neoplastic tissue. This fact highlights the difficulty of the problem of tumor segmentation, which requires a nontrivial solution.

Fig. 9

t-SNE visualizations of the dataset (a) according to tumor status and (b) according to PSL type.

The average signature for each pathology in the validation dataset is presented in Fig. 10. Average signatures show a clear difference among pathologies, both in curve shape and curve amplitude. Melanocytic nevus and Melanoma show a straight appearance and a sudden upward slope around 600 and 650 nm, respectively. Basal cell carcinoma shows a distinct double valley at 540 and 570 nm, but at a much lower reflectance compared with benign pathologies. Additionally, Basal cell carcinoma and melanocytic nevus show a much higher standard deviation compared to other pathologies. However, this fact could be attributed to the large number of available tissue samples and the signatures for those two categories.

Fig. 10

Average signatures for different pathologies. The error bars represent 20% of the standard deviation across signatures for each pathology.

Each model was optimized using the train/test splits of data split 1. The SVM model was optimized using Bayesian optimization at 200 evaluations. This resulted in observed parameters 44.184 for box constraint and 4.0136 for the kernel scale of the RBF kernel. For $\mathrm{Kmeans}+\mathrm{SAM}$, values proposed in the literature were used. Specifically, signatures were split into seven clusters and the sum of SAM similarity per pixel in the cluster was used for segmentation.

Four alternative structures were investigated for the 3D Xception network. The balance of sensitivity and specificity in terms of AUC for the different structures is shown in Fig. 11. Networks with max function at the bottleneck layer display higher AUC compared with the alternative. Initially, four Xception blocks show improved AUC performance compared with three blocks. The convergence of the training process at 200 epochs is also affected. High loss performance during training and testing for the four-block network shows that it suffers from overfitting. The three-block network has lower training and testing loss, as well as smoother convergence. Therefore, a network with three Xception blocks is used in subsequent experiments.

Fig. 11

The influence of the structure of the 3D Xception network on the AUC performance. (a) In-full view, and (b) zoomed-in at the knee of the curve.

The influence of kernel size in the spectral dimension on testing performance is shown in Fig. 12. While JC increases as the kernel size increases, this comes at a cost of precision at 30 nm. After further experimentation, the kernel (3, 3, 10) was used for the first block and (3, 3, 20) for the last two blocks of the 3D Xception network. Dropout layers with a dropout rate 40% were also applied at the end of each block, to avoid overfitting during training. Furthermore, the choice of training parameters greatly affects model performance, as expected. Training with $\mathrm{BCE}+\mathrm{JC}$ loss converged more smoothly and at a lower value compared to training only with BCE loss. The combination of RMSProp with $\mathrm{BCE}+\mathrm{JC}$ loss achieved the best performance for all deep learning models. The learning rate was optimized at 0.0001 for ResNet, at 0.001 for 3D CNN, and at 0.00001 for 3D Xception.

Fig. 12

The influence of kernel size at the spectral dimension in testing JC and precision.

3.2.

Validation Performance

LOOCV was performed twice with the same model parameters for both the preprocessed and the denoised data. This way the effectiveness of spatial denoising with the NGMeet algorithm can be evaluated. Performance metrics without denoising and with denoising are provided in Tables 3 and 4, respectively. Accuracy and specificity metrics paint a misleading picture about the overall performance of segmentation, due to the high prevalence of the negative class. For this reason, more emphasis should be given to JC, AUC, and sensitivity metrics. Out of all the metrics, JC performance is deemed as the most important, due to the balance it provides between false positives and false negatives. Denoising did not improve and even slightly reduced the performance for pixel-wise models. JC performance of preprocessed data without denoising came with increased sensitivity and slightly increased specificity for $\text{Abundance}+\mathrm{SVM}$ and $\text{Signature}+\mathrm{SVM}$. On the other hand, spatial denoising improved JC performance for patch-based models. However, JC performance improved at a cost of specificity.

Table 3

Cross validated performance with LOOCV (without denoising).

Table 4

Cross validated performance with LOOCV (with denoising).

The fluctuation of JC across folds is presented in Fig. 13. For pixel-wise models, preprocessed data is used, whereas patch-based models were trained and tested with denoised data. A tissue sample that was diagnosed as Melanophage aggregate (benign neoplasia), two BCC, and a melanocytic nevus had the worst prediction in terms of JC. It needs to be emphasized that the ground truth labels are approximate because they refer to the histological evaluation of specific cross-sections on the tissue. For the Melanophage Aggregate sample, since the true positive area is just a small growth of melanocytic cells, it is difficult to detect. SVM models detected both the true positive and another false positive area on the left. Pretrained-Resnet detected only the false positive. This false positive area was not evaluated by a cross-section. Therefore, there is a chance that the false positive was a true positive that remained undetected due to the benignity of the diagnosis. The BCC test sample at fold 6 has an unusual appearance where the epidermal layer of the tumor tissue was shed and the deeper cells are prominent. Almost all models missed its detection apart from pretrained-Resnet, which detected almost the entire tissue sample as tumorous. For the other BCC sample (fold 9) where models failed, the mapping of histological ground truth was not successful. The same holds for the Melanocytic Nevus at fold 18. This may affect the metric calculation and should be assigned lower importance during model evaluation.

Fig. 13

Performance in terms of Jaccard coefficient per fold.

The segmented masks at each test sample of 19 folds of the LOOCV are provided in Figs. 14 and 15. The preprocessed dataset without denoising was used for $\mathrm{Kmeans}+\mathrm{SAM}$, $\text{Abundance}+\mathrm{SVM}$ and $\text{Signature}+\mathrm{SVM}$. The denoised dataset was used to train the patch-based models. All pixel-wise models suffered from false positive detection at tissue boundaries. Patch-based models did not suffer from such a problem. This effect was even more pronounced for $\mathrm{Kmeans}+\mathrm{SAM}$. Moreover, $\mathrm{Kmeans}+\mathrm{SAM}$ failed to detect any tumor at all in four cases. $\mathrm{Kmeans}+\mathrm{SAM}$ suffers from detecting illumination artifacts and areas with high contrast from the surroundings, instead of detecting tumor growths. This is especially evident in folds 1, 8, 10, 13, 14, and 15. $\text{Abundance}+\mathrm{SVM}$ detected more homogeneous tumor masks compared to $\text{Signature}+\mathrm{SVM}$. However, it suffered from increased false positives, for example in the cases of BCC (fold 10) and melanoma (fold 11).

Fig. 14

LOOCV segmentation results for folds 1 to 10. From top to bottom: Mucinous Carcinoma, Melanophage aggregate, BCC, BCC, Bowen’s disease, BCC, Melanocytic Nevus, BCC, Melanocytic Nevus, BCC.

Fig. 15

LOOCV segmentation results for folds 11 to 19. From top to bottom: Melanoma, Melanocytic Nevus, BCC, Dermatofibroma, Melanocytic Nevus, Melanocytic Nevus, Melanocytic Nevus, Melanocytic Nevus, Melanocytic Nevus, BCC.

Patch-based models provided more compact and dense masks compared to pixel-wise models. Dense masks are following the biological structure of the tissue. However, they all suffered from false positives and overestimation of the tumor area. Pretrained-Resnet performed better in terms of JC compared to custom-Resnet. Pretrained-Resnet was the only model that detected tumor at the center of the BCC sample at fold 6. However, both it produced much larger tumor masks compared with custom-Resnet, showing a tendency of high false positives. Additionally, pretrained-Resnet tended to label an entire patch as only tumor or only non-neoplasia, resulting in square-shaped segments. This fact is discordant with the shape of tumor borders, which tends to be circular or wavy. 3D Xception produced more curvy tumor segments, although not as dense. False positives were more controlled compared to Pretrained-Resnet. However, the 3D Xception model falsely detected a non-neoplastic edge at the melanoma sample, with more emphasis compared with the Resnet models. This could be attributed to the influence of staining, which was reduced but not completely eradicated by removing the more heavily stained patches. 3D CNN fails to detect the tumor area and mostly detects tissue pixels against the background.

3.3.

Testing Performance

After the validation stage, the models are trained on the entire validation dataset and tested on an unknown testing dataset, data split 3. Segmentation results for each sample are demonstrated in Fig. 16. 3D CNN results are not included during testing, because it was falsely detecting skin instead of the tumor during validation. It should be noted that the test dataset contained one blue Nevus sample, a pathology that was not present in the validation and training dataset. Performance fluctuates a lot depending on the pathology of the test sample. Performance is exceptionally good for melanocytic PSL. The issue of SVM detecting tissue edges as false positives is still present. In terms of technical evaluation, deep learning models provide clearer tumor masks for BCC samples compared to SVM models, with fewer false positives in non-neoplastic areas. However, in terms of medical evaluation, SVM-based pixel-wise masks are more appropriate in cases where tumor detection was successful. Although margins are not delineated well and false positives appear in areas with non-neoplastic tissue, all models apart from $\mathrm{Kmeans}+\mathrm{SAM}$ detect the general area of the tumor.

Fig. 16

Segmentation results for test samples 1 to 9. From top to bottom: BCC, BCC, BCC, Bowen’s disease, Melanocytic Nevus, BCC, Blue Nevus, Melanocytic Nevus, BCC.

LOOCV and testing results were examined by two expert pathologists (TI and MH). The consensus was that the tumor masks were still not sufficient for computer-assisted diagnosis and further improvements are needed. SVM-based models were deemed as closer to the expected ground truth, despite the false positives around the tissue border.