1.

Introduction

Skin cancer is among the most common forms of cancer in the USA, Australia, and Europe and is increasing in incidence in other parts of the world.¹ Visual inspection leading to biopsy followed by histology is the conventional method for clinical screening and diagnosis. Reflectance confocal microscopy (RCM) has been under development² for noninvasive imaging of skin for cancer screening and diagnosis.^{3, 4, 5, 6, 7} RCM enables imaging and visualization of the epidermis and superficial dermis layers below the surface of the skin. Maximum imaging depth is limited to the papillary dermis or superficial reticular dermis, depending on the state of the overlying epidermis and the dermis/epidermis junction. Nuclear and cellular detail is imaged with nominal (instrumental) optical sectioning of 1 to 3 μm and lateral resolution of 0.5 to 1.0 μm, which is comparable to that of conventional pathology. RCM is advancing toward clinical application in dermatology for detection of malignancies such as melanoma and basal cell carcinomas with high sensitivity and specificity.^{3, 4, 5, 6, 8} In these studies, sensitivity for detecting melanomas was reported to be 91% and specificity was 69%. Dermoscopy provides sensitivity of 80 to 100% and specificity of 70 to 100% but mainly for pigmented lesions.^{9, 10} Lightly pigmented, nonpigmented (amelanotic), and pink colored lesions are difficult to differentiate with dermoscopy.¹¹ When examining all such types of lesions, in a more generalized setting, dermoscopy performs with specificity of 32 to 39%, compared to 68 to 84% with reflectance confocal microscopy.⁵ For melanoma, while the sensitivity of RCM was comparable to that of the current standard of visual dermoscopy, the specificity was two times higher. Lentigo maligna melanomas were detected with sensitivity of 85% and specificity of 76%.⁶ For detecting basal cell carcinomas, the sensitivity was 94% and specificity 78%.³ These translational advances represent significant advances for RCM technology toward clinical utility. However, unlike histological sections, RCM images are oriented en face and are grayscale (unstained), and they do not visually resemble conventional stained pathology sections. Thus, differentiation of certain cellular features remains challenging. One such example of a clinically important feature is the dermis/epidermis junction.

Moreover, routine clinical use will require substantial training for clinicians. To accelerate such routine use and increase clinical acceptance, computer-automated methods to assist detection of morphologic features for screening and diagnosis may significantly expand the utility of RCM. Such methods may be either fully automated or semiautomated with user-assistance. However, as we describe below, this is a challenging application for automation with either method. One attempt toward automation was reported by Koller ¹² In their work, they applied an automated classification algorithm to identify RCM image regions as either benign nevi lesions or melanocytic skin tumors. Their algorithm classified 46.71 ± 19.97% of the tumor images in benign melanocytic skin lesions as “malignant,” in contrast to 55.68 ± 14.58% in malignant melanocytic skin lesions.

In this paper, as an initial target, we present a method for semiautomatic localization of the irregular three-dimensional dermis/epidermis junction (DEJ) in RCM images of human skin. The DEJ separates the superficial epidermis from the underlying deeper dermis. We chose to study the localization of the DEJ for three major reasons. First, the DEJ is a clinically meaningful target because cancers such as melanoma and basal cell carcinoma both originate and later spread from this junction. Thus, clinicians (when examining patients) and pathologists (when examining biopsy sections) need to accurately and repeatably evaluate the DEJ for screening, diagnosis, and staging of skin cancers and lesions. Secondly, the DEJ is also a clinically significant structure in many other types of epithelial tissue. Finally, we believe that localization of the DEJ is a useful surrogate for other desired forms of clinically relevant automation of RCM image analysis. Our findings may then be expected to inform and guide future attempts to detect and classify pathologies in skin and also in other tissues.

In dark skin or strongly pigmented skin types, melanin pigment in the basal cell layer creates a strong contrast. The basal cell layer lies directly on the DEJ. Thus, localization of the pigmented basal layer may offer a useful surrogate for the DEJ. By comparison, localization of the DEJ is difficult in RCM images of fair or lightly pigmented skin types due to lack of melanin pigmentation and lack of contrast at the basal layer. Indeed, this is a key motivation for this work: in this initial effort (of an anticipated long-term study), our choices were to address either strongly or lightly pigmented skin types. We chose to address the more difficult cases of lightly pigmented skin. This choice was carefully made based on advice from our clinical colleagues in dermatology.

This difficulty is compounded by high intra/intersubject, and even intralayer variability, resulting from the natural biological heterogeneity of skin tissue. In particular, the epidermis layer has a changing depth-dependent layered structure composed of cells with different morphology at each layer, while the dermis is mostly collagen fibers and blood vessels and is inconsistent in appearance. Due to heterogeneity of skin, at some sites dermis region including collagen fibers appear very bright and at other sites, dermis region appears dark, and at some sites dermis regions looks very much like epidermis. Moreover, since the DEJ is highly corrugated with a hills-and-valleys topography (Fig. 1), its location is not a simple function of depth from the skin surface. This topography of the DEJ surface also causes a lack of consistent appearance and a lack of discriminative features of interest for important local regions, especially regions where the DEJ surface has a large slope. Experience with translational and clinical studies with RCM has shown that clinicians often cannot visually find a single well-defined and consistent junction with accuracy and repeatability without a significant amount of training; even with training, experts may find it difficult to reliably locate the DEJ.

Fig. 1

Left figure shows the DEJ in a vertical histology cross-section image and the middle and right figures show lateral slices from a RCM stack with the epidermis/dermis boundary marked. The DEJ is a thin membrane, shown with a blue solid line, that separates the epidermis from the dermis. a single layer of basal cells lies directly on the DEJ. The basal cell layer is typically at average depth of 100 μm below the surface in normal skin and 10 to 15 μm in thickness (Ref. 1). (Color online only.)

Indeed, after our steps in initial algorithm development¹³ and subsequent feedback from clinical collaborators, we determined that localization of a single DEJ boundary (i.e., a single DEJ) is not consistently supported by the data. Instead, the goal adopted here is to locate two 3D surfaces that together bound a transition region in which the DEJ is located. The region above the first surface is determined to be epidermis with high confidence, and the region below the second surface is determined to be the dermis, also with high confidence. We reported some preliminary results in Kurugol ¹⁴

In this paper, we report further advances to solve this DEJ localization problem. To attack this problem, we found it necessary to incorporate ideas from three fields: texture segmentation, pattern classification, and sequence segmentation (SS). One component of our algorithm grew out of a class of standard approaches to segment regions based on their different textural features by using a classifier¹⁵ (see Ref. 16 for a general review of texture segmentation algorithms). Beyond such standard 2D texture segmentation, we incorporated ideas from dynamic texture segmentation and texture change detection. Dynamic texture segmentation¹⁷ is used to segment a sequence of 2D images with identical spatial statistics but dynamic changes along a third dimension. This dimension can be time, as in image sequence segmentation; here, the third dimension is depth through the skin, along the optical axis. Such problems can be solved by sequence segmentation algorithms¹⁸. For example, a change detection solution for image sequences is reported in Ref. 19.

Our algorithm works on small regions or “tiles.” The algorithm partitions each vertical stack of tiles (see Fig. 2 for an example) into a variable number of homogeneous segments in z (depth). We define homogeneity in terms of a multivariate dynamic model on a chosen set of image texture features. Based on the resulting set of break points between segments, we then identify two of them as the transition region boundaries using a multiclass classifier in the z-direction that assigns segments to epidermis, transition, or dermis. This process starts from the upper and lower images in each stack and moves toward the middle until the classifiers locate the respective region boundary at a segment end point.

Fig. 2

An example stack (sequence) of 60 tiles is shown, with increasing depth indicated by increasing slice number in the figure. For this stack, an expert evaluator (see Sec. 3 for details) located the epidermis boundary at slice 19 and the dermis boundary at slice 29.

We evaluated two different training scenarios for our algorithm. The first scenario was based on training and testing the algorithm on the same stack. In this scenario, we labeled a small training set for a given stack and used it to train the classifiers, and then we applied them to the rest of that stack (same stack scenario). In the second scenario, we labeled and trained on an entire stack and applied the resulting classifiers to different stacks from subjects with similar skin types (cross stack scenario). Results from both scenarios showed a reasonable performance considering the difficulty of the problem. Our methods and results are reported below.

2.

Methods

In this section, we describe our algorithm for the semiautomatic localization of the DEJ in RCM images. The stages of the algorithm are summarized in the flow chart (Fig. 3). We will describe each of these stages in turn. In particular, in Sec. 2.1, we describe the acquisition and preprocessing of the RCM data stack. Section 2.2 describes the set of texture features extracted from the training data and the feature selection algorithm used to select the most relevant, least redundant features. We explain sequence segmentation and classification algorithms employed in the two parallel stages in Secs. 2.3, 2.4 respectively. In Sec. 2.5, we explain how the algorithm combined the results of the sequence segmentation and classification stages to identify which segment boundaries correspond to skin layer boundaries. Finally, we describe some post-processing applied to smooth the two resulting boundaries in Sec. 2.6.

Fig. 3

Flow chart of the algorithm.

2.3.

Sequence Segmentation in z-Direction

We now explain the two parallel stages used to extract boundaries and classifiers, which are subsequently combined to determine boundaries. We describe the sequence segmentation algorithm, which uses the features selected for that algorithm as described above. To start, we denote that set of d features for all n tiles in a given stack (or sequence) of tiles as F = {f ₁, …, f _n}, of length n, where each point is d-dimensional. A segmentation partitions F into k contiguous segments such that each segment is as homogeneous as possible as in the sense below.

We applied the sequence segmentation algorithm to each sequence of tiles in the 3D RCM stack.

In detail, for each k in a fixed range (here, k takes values 3, 4, …, 8), using an optimal dynamic programming (DP) algorithm, we partitioned the multidimensional feature vector sequence into k homogeneous segments. We modeled each dimension of our d-dimensional feature vector with a piecewise affine model for each segment. The cost function used by the DP algorithm is the l-2 norm of the error between that model and the features:

Eq. 2

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} {\rm Cost}=\sum _{j=1}^k \sum _{i\in s_j} ||\mathbf {f_i}-(\mathbf {a_j} i+\mathbf {b_j})||_2. \end{equation} \end{document} $$\mathrm{Cost}=\sum _{j=1}^k\sum _{i\in s_j}\left|\right|{\mathbf{f}}_{\mathbf{i}}-({\mathbf{a}}_{\mathbf{j}}i+{\mathbf{b}}_{\mathbf{j}}){\left|\right|}_2.$$

Here, given the desired number of segments k, the depth indices of the tiles, denoted by i, are partitioned into a set of k contiguous subsets, denoted s _j, one for each segment. For each i’th tile within a segment, the d-dimensional feature [TeX:] \documentclass[12pt]{minimal}\begin{document}$\mathbf {f_i}$\end{document} ${\mathbf{f}}_{\mathbf{i}}$ ( [TeX:] \documentclass[12pt]{minimal}\begin{document}$\mathbf {f_i}=[f_{1i},f_{2i}, \dots, f_{di}]$\end{document} ${\mathbf{f}}_{\mathbf{i}}=[f_{1i},f_{2i},\cdots ,f_{di}]$ ) is approximated by the affine combination a _j i + b _j, where the d-dimensional vectors of affine function parameters [TeX:] \documentclass[12pt]{minimal}\begin{document}$\mathbf {a_j}$\end{document} ${\mathbf{a}}_{\mathbf{j}}$ ( [TeX:] \documentclass[12pt]{minimal}\begin{document}$\mathbf {a_j}=[a_{1j}, a_{2j} \dots a_{dj}]$\end{document} ${\mathbf{a}}_{\mathbf{j}}=[a_{1j},a_{2j}\cdots a_{dj}]$ ) and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\mathbf {b_j}$\end{document} ${\mathbf{b}}_{\mathbf{j}}$ ( [TeX:] \documentclass[12pt]{minimal}\begin{document}$\mathbf {b_j}=[b_{1j}, b_{2j} \dots b_{dj}]$\end{document} ${\mathbf{b}}_{\mathbf{j}}=[b_{1j},b_{2j}\cdots b_{dj}]$ ) are held constant throughout the j’th segment but the affine “variable” i varies with tile index. We collect the approximation errors across the k segments as our cost. We then optimized for the model parameters for each segment along with the segment boundaries (definitions of s _j), for a given k, by jointly minimizing the cost across all variables.

We used a heuristic approach to find the number of segments to use: We first ran the algorithm for k from 3 to 8. Then the cost ratio between using k − 1 and k segments was calculated for k’s starting from 4 and repeating until this ratio fell below a threshold t (set to 1.3 experimentally).

In Fig. 4, we show an example multivariate z-sequence of features. For illustration purposes, only four features were plotted. Each feature is represented by a different color. The z-sequence is of length 60, i.e., 60 slices. The segment boundaries of the eight segments found by the sequence segmentation algorithm are shown with solid blue vertical lines. The dashed vertical red lines show the epidermis and dermis boundaries located by the expert. We note that one can visually observe a qualitative difference between each adjacent pair of regions, which presumably is what determined where the sequence segmentation algorithm placed the boundaries automatically.

Fig. 4

An example multivariate z-sequence of features. For illustration purposes, only four features are shown. The segment boundaries of the eight segments found by the sequence segmentation algorithm are shown with solid blue vertical lines. The dashed vertical red lines show the epidermis and dermis boundaries located by the expert. (Color online only.)

2.5.

Combined Sequential Decision

The basic approach we adopted to combine these two sets of results is to use the sequence segmentation endpoints as the candidate set of final layer boundaries for each tile. We used the classifiers to determine which two endpoints to choose as those boundary locations. In particular, we employed the following procedure once from top to bottom for the epi-trans boundary and again from bottom to top for the trans-der boundary.

Starting from the top (bottom) segment of tiles in the stack we identified the epi-trans (trans-der) boundary as the segment end point beyond which the average probability of belonging to epidermis (dermis) from the epidermis versus rest (dermis versus rest) classifier was below 0.4. This number was selected to be less than 0.5 intentionally to account for the noise and the effect of averaging the probabilities, i.e., if probabilities of most of the tiles within a segment are around 0.5 but a few tiles in that segment are below, we intend to include that segment to that class. Figure 5 illustrates the procedures of sequence segmentation and boundary location decision. For the sequence in Fig. 4, the algorithm located the epidermis boundary and dermis boundary both at the second segment boundary (with slice number 37).

Fig. 5

Left panel shows the tile sequence and an example output of the sequential segmentation algorithm. Right panel shows the resulting epidermis and dermis boundaries (yellow longer horizontal lines) of the combined sequential+classification decision algorithm. (Color online only.)

3.

Experiments

As explained in Sec. 1, we used two different scenarios for training. In the same-stack-training scenario, out of 41 tile sequences in the stack, 8 tile sequences were used, 5 for validation and 3 for training, to train our classifier and tune the classifier parameters. We performed standard N-fold (N = 6) cross validation experiments to tune our classifier parameters. In these experiments, a “brute force” search in the parameter space was performed, followed by choosing the parameter combination with the best classification performance over the validation set. In order to provide the LS-SVM algorithm with the neighboring samples it uses for local smoothing, tiles were chosen in 2 × 2 adjacent blocks, of size 100 × 100 pixels, which were then subdivided into 50 × 50 pixel tiles. The reason for choosing this small neighborhood size of 2 × 2 is due to a constraint of the expert markup in scenario 1. In particular, in scenario 1 the expert only marked a small number of training samples. To ensure the availability of adjacent tiles for training, we had the experts mark regions of size 100 × 100 pixels (giving us 2 × 2 neighborhoods of 50 × 50 tiles). The expert marked the two boundary depths in these selected tile stacks. The rest of the processing was automated once these labeled tiles were provided. In scenario 2, the cross-stack-training scenario, we assumed that data stacks from similar skin types have similar features. Thus, a classifier pretrained on one completely labeled RCM stack from a particular skin type was applied to a new RCM stack from a similar skin type without any need to label tiles in the new stack.

4.

Results

We report results from experiments on four RCM test data sets from four different subjects with fair skin types (Ref. 29)—skin type I and II—which as noted above are the most difficult cases due to the small amount of melanin pigment and, hence, contrast. Results are reported for both training scenarios (scenario 1: same stack and scenario 2: cross stack). Note that the expert fully labeled both stacks. In both scenarios, the expert markings not made available to the algorithm were used to evaluate the performance.

Since in scenario 1 the expert only marked a small number of training samples, we made sure that marked regions included at least one neighborhood; to do so we had the experts mark regions of size 100×100 pixels (giving us 2 × 2 neighboring 50×50 tiles). Thus, we used a small neighborhood of size four (three adjacent tiles). Note that scenario 2 does not require this restriction, and as the reviewer pointed out, it would make more sense to use 4 or 8. To be consistent for both training scenarios, we utilized the neighborhood of size four (three adjacent tiles) consistently in the results reported here.

The results from the automatic supervised feature selection algorithm indicate that log-Gabor and wavelet features were the most commonly selected features. These features are consistent with the visual features used by an expert to discriminate epidermis from dermis, which are typically the blurriness of dermis versus the cellular texture pattern of epidermis. We speculate that this information can be captured well by obtaining frequency information localized in space provided by the features calculated from log-Gabor filtered and wavelet packet decomposed tiles. The energy and variance features calculated from the wavelet transform at different scales provided texture information at those scales. Log-Gabor filters, which are a product of a Gaussian (with frequency represented in log-scale) and a sinusoid, are also used for texture discrimination at various scales. For example, for structured cellular epidermis regions, we obtained higher values from our energy feature at high frequency bands, while for blurry collagen fibers within dermis, we obtained higher values from our energy feature at lower frequency bands. This kind of texture exploration at multiple scales provided us the information needed to discriminate epidermis from dermis. However given the variability of the data and the stochastic nature of the feature selection process, we cannot make definitive statements either about which features are selected (or even the number of features selected) nor a physical interpretation of this outcome.

To quantify performance we calculated the distance in micrometers between the expert labeled boundaries and classification boundaries on each test set. For each experiment, we calculated the number of z-sequences of tiles for which this distance was smaller than both 10 and 15 μm. We report this number as well as the mean and standard deviation of the error distances for both stacks in Table 2 for scenarios 1 and 2, respectively. We also report the results as confusion matrices. The confusion matrices for the first and second scenarios are shown in Tables 3, 4, respectively. The (j, k)'th entry of the confusion matrix indicates the number of tiles found to be class k [1-epidermis, 2-transition (transition) or 3-dermis] by the algorithm given that they belong to class j according to expert markings. For both data sets we also plot the resultant 3D epidermis and dermis boundary surfaces for scenario 2 in comparison to the expert marked boundary surfaces in Fig. 6. Single-color surfaces [Figs. 6a and 6c] show the expert labeled epidermis (dermis) boundary while the multicolored surfaces [Figs. 6b and 6d] indicate the boundaries found by the algorithm. The color maps on the multicolored surfaces show the distance from the expert labeled boundary. In Fig. 7, for both RCM stacks and scenarios 1 and 2, we show the epidermis and dermis boundaries located by the algorithm in comparison to the expert located boundaries for all 164 ¹ z-sequence of tiles. In Fig. 8, a comparison of expert markings with the smoothed algorithm results are shown for scenario 1 in two orthogonal vertical views at the locations from the RCM stack 1 indicated by the solid lines drawn in the axial views (on the left).

Fig. 6

Surface plot of the epidermis boundary and the dermis boundary in 3D in comparison to the expert labeled boundaries of RCM stack 1 and 2 for scenario 2 (cross training). Top blue (bottom red) surfaces show the expert labeled epidermis (dermis) boundary for (a) RCM stack 1 and (c) RCM stack 2. The colored surfaces indicate the resultant boundaries of the algorithm for (b) RCM stack 1 and (d) RCM stack 2. The color maps indicate the distance from the expert labeled boundary. The z-axis is in micrometers. x and y axes are in pixels, where the pixel spacing is 0.5μm. Flat regions are the masked out wrinkles. For the smooth visualization purpose, the boundaries are plotted after interpolating them twice in 2D with spline interpolation. (Color online only.)

Fig. 7

For scenario 1 and 2 and RCM stacks 1 and 2, the figure show the epidermis and dermis boundaries located by the algorithm in comparision to the expert located boundaries for all of the 164 tile-sequences that were processed by the algorithm. The boundaries shown are 2D Gaussian filtered for smoothness as explained in the post-processing step in Sec 2. The dotted vertical lines in (c) indicate the location of the vertical slice shown in Fig. 8b. (Color online only.)

Fig. 8

Comparison of expert markings with the algorithm results shown in vertical views y-z (top) and x-z (bottom). The solid line on the left of both (a) and (b) indicate the vertical slice location. Transition region is located by the algorithm in between epidermis algorithm (green) and dermis algorithm (purple) curves. The green (purple) curve is the epidermis (dermis) boundary found by the algorithm. The blue (red) curve is the dermis (epidermis) boundary marked by the expert. If there is no epidermis expert (blue), the expert found no transition region and the upper and lower boundaries coincide. For visualization purposes, algorithm boundaries computed for each tile are linearly interpolated to the same grid (pixel grid) that the expert used in their mark-up. (Color online only.)

Table 2

Results of Scenario 1 and Scenario 2 are given for both RCM stack 1 (column 2) and RCM stack 2 (column 3). Rows labeled N give the number (and ratio) of test tiles (out of the total tested) for which the detected boundary was within 10 μm (N10) or 15 μm (N15) of the expert marked boundary. Rows labeled m ± std give the mean and standard deviation of the error in micrometers between detected and expert marked boundaries across the test set.

Table 3

Scenario 1: Confusion matrices for RCM stack 1 RCM stack 2 as test set. Confusion matrices show the results of algorithm (columns) given the expert results (rows), as explained in the text.

Table 4

Scenario 2: Confusion matrices for RCM stack 1 RCM stack 2 as test set. Confusion matrices show the results of algorithm (columns) given the expert results (rows), as explained in the text.

We also show average results for experiments performed on stacks 1 to 4. In those experiments, a classifier is trained on one stack and applied to the remaining three stacks. Therefore, three different classifiers were applied to each stack and the results were calculated. The results for each stack averaged over all three classifiers applied to that stack are reported in Tables 5, 6. Table 5 shows average confusion matrices for RCM stacks 1 to 4. The average was calculated over three classifiers applied on a stack. For each stack, the classifier trained on that same stack was not used in the testing. Table 6 also reports average results calculated similarly. Rows labeled N show the mean ratio of test tiles (out of the total tested) for which the detected boundary was within 10 μm (N10) or 15 μm (N15) of the expert marked boundary. Rows labeled m ± std give the average mean and standard deviation of the error in micrometers between detected and expert marked boundaries across the test set.

Table 5

Average confusion matrices for RCM stack 1 to 4 as test set. For each stack, the average was calculated over three classifiers trained on the other stacks and applied on the remaining stack. For each stack, the classifier trained on that same stack was not used in the testing. Confusion matrices show the results of algorithm (columns) given the expert results (rows), as explained in the text.

Table 6

Average results for stack 1 to stack 4 as test set. For each stack, the average was calculated over three classifiers trained on the other stacks and applied on the remaining stack. For each stack, the classifier trained on that same stack was not used in the testing. Rows labeled N show the mean ratio of test tiles (out of the total tested) for which the detected boundary was within 10 μm (N10) or 15 μm (N15) of the expert marked boundary. Rows labeled m ± std give the average mean and standard deviation of the error in μm between detected and expert marked boundaries across the test set.

The results from the same experiments are also reported by averaging the performance of each classifier over the three stacks that it was tested on. For each classifier, the stack on which the classifier was trained on was not used in the testing. Table 7 shows average confusion matrices for classifiers 1 to 4 and Table 5 shows average distance results calculated similarly for each classifier. Thus with these tables we can examine both the performance of each classifier across all other stacks and the performance of all classifiers on each stack.

Table 7

Average confusion matrices for classifiers 1 to 4. The average was calculated over three out of four data stacks, on which the same classifier was applied. For each classifier, the stack on which the classifier was trained on was not used in the testing. Confusion matrices show the results of algorithm (columns) given the expert results (rows), as explained in the text.

The video, which shows the classification results of scenario 2 applied on RCM data stack 1, is provided as a multimedia file and a snapshot from the video is shown in Fig. 9. The left panel shows epidermis and dermis surfaces and a cutting data slice that moves from the top of the stack to the bottom. The right panel shows the original data slice (bottom) and the same slice with the overlayed algorithm results (top). The video starts from a superior slice of the stack, where all regions were either classified as epidermis (red shaded) or were masked out (dark gray shaded) in the preprocessing stage. Then the cutting plane proceeds to deeper slices. Moving deeper in the stack, first the epidermis regions shrink and the transition regions (light gray shaded) start. Then the transition regions shrink and the dermis regions (blue shaded with solid boundary) start. The deepest slices in the stack include only dermis regions.

Fig. 9

A snapshot from the video file which shows the classification results of scenario 2 applied on RCM data stack 1. The left panel shows epidermis and dermis boundary surfaces and a cutting data slice that moves from the top of the stack to the bottom. The right panel shows the original data slice (bottom) and the same slice with the overlayed algorithm results (top). The video starts from a superior slice of the stack, where all regions were either classified as epidermis (red shaded) or were masked out (dark gray shaded) in the preprocessing stage. Then the cutting plane proceeds to deeper slices. Moving deeper in the stack, first the epidermis regions shrink, and the transition regions (light gray shaded) start. Then the transition regions shrink and the dermis regions (blue shaded) start. The deepest slices in the stack include only dermis regions. (MPEG, 21.1MB)

5.

Discussion

In this work, we developed an algorithm to locate the DEJ in RCM image stacks for lightly pigmented skin types. The proposed hybrid algorithm locates the epidermis and dermis boundary surfaces with a decision combining sequential segmentation and classification stages. The results show that the algorithm performed reasonably well with epidermis/dermis misclassification rates smaller than 10% and average distance from the expert labeled boundaries around 7 to 12 μm. When considering the fact that the standard epidermis—dermis classification approaches perform poorly due to a number of reasons including the heterogeneity of skin layers, the proposed combined decision algorithm achieves a significant improvement by making use of the spatial structure and dynamics inherent in the tissue.

Due to this heterogeneity of skin structure, for less than 10% of all tile stacks ², the classifier effectively could not distinguish between the different layers. Tiles in these tile-stacks were found upon visual inspection to have an ambiguous appearance which was dissimilar to the large majority of the tile stacks. Hence, the decision function that uses the LS-SVM classifiers failed to discriminate either epidermis from the rest or dermis from the rest due to lack of discriminative features between layers for these tile stacks and the concomitant lack of dynamics in the z-direction.

The results from all tile stacks are shown in comparison to the ground truth as determined by expert markups in Table 2. The presentation of our results implies that we accept an accuracy of within 10 to 15 μm as useful. The rationale for using these thresholds is that because the size of a basal layer cell is about 10 to 15 μm, we are within about one cell distance of the boundary. We note that when clinicians visually evaluate a data stack, they typically evaluate with 5 μm slice separation, which is considered a reasonable distance for a visually detectable change to occur from slice to slice.³⁰ In that case, 10 to 15 μm would correspond to 2 to 3 slices.

For the dermis boundary, as shown in Table 2, this acceptable distance from the expert boundary is achieved for more than 82% of all tile stacks. Similar results can be observed from the plots in Fig. 7 (i.e., the algorithm's dermis boundaries are almost touching the expert dermis boundaries for most of the tile stacks). The algorithm performance is reported in terms of the confusion matrices in Tables 3, 4. Regions labeled as dermis both by the algorithm and the expert are given in the (3,3) entries of the confusion matrices, and again these values are around 80%. The epidermis/dermis misclassification rates for dermis are shown in (3,1) entries of the confusion matrices and are less than 10%. For experiments performed on all four RCM stacks, similar results averaged over a stack and averaged over a classifier were also reported in Tables 5, 6, 7, 8. The average classifier performance on each stack reported in terms of confusion matrices indicate a comparable performance over all stacks even though some classifiers might work better on some stacks and worse on other stacks due to variability from subject to subject. Mean distances between expert and algorithm boundaries were around 7 μm for dermis and 9 μm for epidermis and the acceptable distance of 15 μm was achieved for around 90% of tiles for dermis and 85% of tiles for epidermis.

Table 8

Average results for classifier 1 to classifier 4. The average was calculated over three out of four data stacks, on which the same classifier was applied. For each classifier, the stack on which the classifier was trained on was not used in the testing. Rows labeled N give the mean ratio of test tiles (out of the total tested) for which the detected boundary was within 10 μm (N10) or 15 μm (N15) of the expert marked boundary. Rows labeled m ± std give the average mean and standard deviation of the error in micrometers between detected and expert marked boundaries across the test set.

On the other hand, the epidermis boundaries of the algorithm are generally chosen more conservatively than the expert labeled epidermis boundaries. The algorithm was designed to label the regions as epidermis only if it finds high probability

that they are epidermis; otherwise it labels the regions as transition. Hence, the epidermis boundary of the algorithm is on average about 20 μm away from the expert labeled epidermis boundary. Consistent epidermis results are shown in Tables 3, 4, where confusion matrix entries (1,1) indicate epidermis classification accuracies. These accuracies for epidermis are around 45%, while the (1,2) entries of around 40% indicate conservative detection of epidermis boundaries as regions labeled as epidermis by the expert are labeled as transition by the algorithm. Although the epidermis classification accuracy is not as high as for dermis due to this phenomenon of epidermis being classified as transition in many cases, the epidermis/dermis misclassification for epidermis, as shown in the (1,3) entries of the matrices, is very low, with values less than 10%, which suggests successful epidermis/dermis classification. If we decide, based on this analysis, to include the transition regions into the epidermis class—that is place our detected epidermis boundary below the transition region instead of above it—then the epidermis classification accuracy of the algorithm is increased to around 90% in all cases and dermis accuracy stays around 76%. In that case, the detected DEJ boundary would be a single surface, the dermis boundary found by the algorithm. Similar analysis for experiments performed on all four RCM stacks show that including transition regions into the epidermis results in accuracy values larger than 75% for dermis for all stacks and 90% for epidermis for all stacks except stack 3.

The transition region labeled by the algorithm generally includes some deeper epidermis regions lacking strong epidermis features and some superior dermis regions lacking strong dermis features. For such regions, the expert labeled transition region is likely somewhat subjective, i.e, the height of the transition region in z may change from expert to expert. However, experts generally try to mark the transition region as narrowly as possible. In particular, we speculate that experts may, in effect, add in more complex perceptual features to their decision process, or use information from a large number of neighboring tiles to make decisions in uncertain regions, even if the sufficient low-level textural clues are not present. This suggests that the algorithm could be refined in the future to, similarly, identify regions of larger uncertainty, and in those cases employ more complex features, rely more on spatial context, and/or use hierarchical classification approaches.

In the current version of our algorithm, the LS-SVM classifier leverages the spatial correlations between tiles. However, during the sequential segmentation of each tile stack based on dynamics, the spatial correlation between neighboring tile stacks is not used in the segmentation. Moreover, in the combined decision step, where we locate the epidermis and dermis boundaries as one of the segment end points, we again do not consider the variation of tile stack boundaries across stacks within a neighborhood; we only enforce smoothness in our post-processing step. These limitations will also be addressed in future development of our approach.

One significant limitation of the algorithm is parameter tuning and sensitivity. The SVM parameters are tuned automatically with a grid search over a labeled validation set. This is the method commonly used in the SVM literature.³¹ There are also two thresholds in the feature selection algorithm that are also tuned by using the validation set. As future work, a more robust feature selection algorithm such as the one in Ref. 32 may replace the current fast filter approach.

We have begun to test our algorithm on a larger database, since the proof of concept reported in this paper indicated that training on one stack and testing on another stack from a similar skin type (cross training scenario) will likely be a useful method. This will be especially important if we wish to process a large number of RCM stacks without the need to mark up a training set in each individual stack.

In this work, we only treated data stacks from subjects with fair skin types, which have only a small amount of melanin pigment and hence almost no contrast at the DEJ location. These were, by far, the most challenging cases and we felt it was important to establish whether this approach could work under these low contrast conditions. For darker skin types, the availability of strong contrast at the basal layer provides a clearly detectable surrogate feature to localize the DEJ. In such skin, an algorithm based on peak detection in z-profiles of mean intensities of each tile stack may be adequate. A peak in the intensity profile is likely to be localized in the basal layer (located right above DEJ). In preliminary studies on a few cases it appeared that the first strong jump in that z-profile corresponded to the superficial stratum corneum layer, provided that the stack started from the skin surface, and the second strong jump corresponded to the basal layer location superficial to the DEJ. We note that a similar idea was used in Ref. 33 to detect pagetoid melanocytes.

In our future work, we plan to test the DEJ localization algorithm in skin lesions and cancers in vivo. After locating the DEJ, the location of the lesion can be identified with respect to DEJ location. A similar or extended set of texture features can be calculated from the lesions and physical and physiological interpretation of most relevant features can be examined.