1. Introduction

The diagnosis and treatment of skin disease is largely based on visual diagnosis by dermatologists who rely on extensive training to build their expertise. There are, however, many conditions that are visually ambiguous. Often, especially for cancer screening, biopsy followed by histological analysis is used to resolve these ambiguities. There has long been a hope that novel imaging methods could be used for a non-invasive or optical biopsy. This ambition has been one of the drivers for the development of techniques such as in-vivo confocal microscopy [1], in-vivo Raman imaging [2], in-vivo two-photon microscopy [3] and optical coherence tomography (OCT) [4,5]. Most of these techniques are limited to imaging the topmost 100–200 micrometers due to a host of issues including multiple scattering, numerical aperture and photon absorption. OCT, in contrast, can image as deep as one millimeter into the skin using commercially available high-quality instruments [6]. The images, examples of which are shown in Figure 1, clearly reveal skin structures including sweat ducts, a change in contrast around the junction between the dermis and the epidermis (the DEJ) and an impression of the layers in the various dermal plexi. Others have described seeing other important features, such as cancerous lesions and other skin diseases. Detailed examination of the images, however, is dominated by the random texture of speckle. Seeing through speckle to definitively identify structures and boundaries is notoriously challenging. Despite significant work to improve the technique to minimize speckle [7], the problem remains. While it seems obvious that image analysis (IA) should be part of the approach, it has turned out to be quite challenging. Specifically, the problem of identifying a contour is made difficult by the need to see it through the modulation of a random speckle field.

 

Fig. 1. OCT images of skin. (a) Volar forearm, showing skin layers. b) Thumb, with thick stratum corneum and visible sweat glands.

Download Full Size | PDF

Identification and monitoring of the epidermis is important because thickness of the epidermis and the local shape of the DEJ are important factors in monitoring skin health, aging, and photodamage [8]. Despite this, little published work addresses the measurement of epidermal thickness from OCT images. The most common approach, embodied, for example, in the software available with the commercial ISIS skin imaging system involves simply adding columns of the image together to produce the average depth profile of the OCT signal [9]. While this seems intuitively correct, it glosses over the fact that neither the skin nor the DEJ are flat. In skin, this approach, known as A-scan analysis, has not to our knowledge been subject to any systematic validation. Another approach is to use the boundary detection system of the human visual system; that is to train operators to find the DEJ and measure epidermal thickness using a ruler or potentially using software to convert the hand-drawn lines to an average thickness. We have found this technique to be challenged by two problems. First, it is time consuming. Second, no two lines are drawn exactly the same and variations exist when comparing different lines, even when drawn by the same observer. We report here an IA method that is sufficiently successful that we believe it justifies more detailed clinical comparisons between OCT and physical biopsy. Overall, it seems probable that the ability to automatically identify features in healthy tissue should enable the extension of this methodology to help determine the boundaries of both benign and cancerous lesions.

2. Materials and methods

The data considered in this report were all collected using a commercial OCT instrument (ISIS Skindex 300). This device has been designed specifically for imaging skin, and illuminates the skin using a set of 8 near- infrared light emitting diodes (LED) at 1300 nm. The 8 LEDs are used to collect data simultaneously from 8 channels. Mismatches in brightness and gain between the different LED channels cause the appearance of vertical bands in the displayed images. The Skindex 300 collects data from as deep as 900 µm. The data is collected and displayed as a vertical slice in the plane perpendicular to the skin surface. All data were collected from forearms or lower legs of twelve subjects as part of previous studies. To prepare the data for analysis, we imported the raw data into our analysis package (Matlab) using routines provided by the instrument manufacturer. By working with the raw data, we were able to access the full 12-bit dynamic range of the instrument. Close examination of the images will reveal that the pixel size is considerably smaller than a single speckle. This is, of course, a characteristic of any imaging technique based on coherent illumination [7].

 

Fig. 2. A-Scan analysis of an OCT image of skin.

Download Full Size | PDF

The DEJ is the boundary in skin between the dermis and the epidermis above it, and marks a difference in tissue function and morphology. The epidermis contains the melanin pigment responsible for skin coloration and epidermal thickness is an indicator of skin aging and health. We tested the existing methods for measuring epidermal thickness in OCT data by comparing the results of line-tracing by three trained observers with the A-scan method. The observers traced the DEJ using a simple computer program that recorded the line they drew and then calculated the thickness in each column of pixels, returning an average. All of the observers were research scientists specializing in skin and were trained by a process of repeat observation, discussion and repetition of the DEJ tracing process on at least 100 images. The DEJ was identified by comparing OCT images to traditional histological images of skin and to results of previous investigations of skin structure using OCT [10]. The variability between different lines drawn in repetition was tested by having each observer draw five lines marking the DEJ on each of a series of test images. Each observer was fairly consistent with RMS error between different attempts to draw the DEJ falling between 6.7 and 7.5 mm. Also, inter-operator agreement was good for most files, with an average RMS error of less than 6 micrometers for 4 out of 5 files. The A-scan method was performed using the ISIS OCT Evaluation software (Fig. 2). The DEJ was identified in the A-scan as the significant local minimum with the best correspondence to the observed depth in the OCT image.

The IA presented here is organized into four stages (Fig. 3). First, the images are first preprocessed to slightly smooth over speckle. Then a shapelet decomposition algorithm is used to produce intermediate images corresponding to changes in the primary image at three different lengthscales. Third, contours are traced. The shortest length scale shapelet is used to find the bottom of the stratum corneum (SC), while the medium and long lengthscale shapelets help identify the DEJ and calculate the confidence in the accuracy of the line identified as the DEJ. The final step consists of calculating an average epidermal thickness, using those only columns in which confidence estimates are sufficiently large. The resulting lines are displayed as superimposed features on the original image.

 

Fig. 3. Shapelet decomposition of OCT skin images. a) Raw image. b) Short lengthscale shapelet highlights stratum corneum. c) Mid and long lengthscale shapelets find DEJ. d) Reconstructed image with shapelet results.

Download Full Size | PDF

It is helpful to pre-process OCT images as they are speckled and contain other artifacts that will interfere with image analysis. The pre-processing step removes some of the noise and prepares the data for shapelet analysis, which adequately handles most subtle noise features. The preprocessing algorithm is divided in four steps. First, we reduce the image size by half to speed processing. Second, the eight acquisition channels are rescaled to the global channel intensity. Random single-pixel noise is then removed by factor analysis denoising (FAD) [11], which is implemented by treating the image columns as vectors and then performing an eigenvalue decomposition. The 35 most significant factors are preserved in the reconstructed image. This step significantly reduces random pixel noise. Finally, after the FAD step, the speckles are softened by a least-median-of-squares (LMS) adaptive spatial filter that preserves the larger scale edge features. These preprocessed images are suitable inputs for shapelet decomposition.

Shapelet decomposition is well known in the astrophysics literature as a method for modeling deep-space image features [12], and has also been successfully applied in the computer vision literature to the reconstruction of 3D spatial features from surface normals [13]. The method we have chosen for our work is an adaptation of this latter methodology. Shapelet decomposition can be thought of as a two-dimensional analog of wavelet decomposition in which an image containing a particular “shape” is used as the basis kernel rather than an orthogonal or bi-orthogonal 1D or 2D waveform. In our application, it decomposes an OCT image using a set of unique basis images such that the gray level information is transformed into a set of surface contours. This is achieved through convolution of the original image with specially chosen shapelets representing different features at different lengthscales. This type of shapelet analysis is useful for examination of OCT images of skin because it helps find boundaries between the different skin layers that are otherwise difficult to identify by traditional IA techniques due to speckle and other sources of noise.

Shapelet analysis uses one or more basis image “shapes” as the kernel in a mathematical correlation with an image. Correlations are generated using a convolution filter that applies the kernel to the entire image. In this implementation of shapelet analysis, kernels are chosen that best represent contours corresponding to ridge peaks and troughs in OCT images related to the SC and DEJ features in the original image. Different lengthscale kernels are employed to distinguish between these two OCT features and to exclude interfering features such as random noise and speckle.

The shapelet we chose for our analysis was based on the depth profile of the two features of interest. A simple gaussian shapelet was good match for the stratum corneum at short lengthscales but was too sharp to identify the contour of the DEJ. Better performance was obtained with a damped gaussian-like shapelet, s, based on a Butterworth function: s(x,y)=1/(1+((x ²+y²)^1/2/L_i)²ⁿ). Lengthscales (L_i) are specified in ranges and are allowed to vary within the range in order to provide the best possible fit, and the subscript i refers to a particular lengthscale. Four different ranges of length-scales were used to perform feature identification within the skin OCT images. The sharp boundary of the stratum corneum was found using shapelets with the smallest L_i values (3.6 µm, 5.4 µm, and 7.2 µm). Five shapelet functions of longer lengthscales were used together to identify the DEJ. These functions had L_i values of 14.4 µm, 18.0 µm, 21.6 µm, 28.8 µm, and 32.4 µm.

OCT images lack the necessary graylevel contrast to allow direct shapelet analysis of the raw image. We chose instead the tilt (τ) and slant (σ) image representations that are commonly used in similar analyses in the computer vision literature [13]. Tilt and slant images are computed directly from the 2D gradient components of the image, G_x and G_y. τ is the solid angle the gradient makes with the horizontal plane. σ is the slope, or the arctan of gradient image. The slant image is related to the magnitude of the image gradient at each pixel. The tilt image is used to provide feature direction information with respect to the image plane.

Shapelet images are generated from image gradient components by computing gradient correlations (C _∇) between the magnitude of image gradient (|∇_G|) and the magnitude of the shapelet kernel gradient (|∇_si |).

Gradient correlation values by themselves contain no information regarding the orientation of the gradient, and the tilt correlation information must be used to determine whether a feature is pointing in or out of the image plane. The tilt correlation between a shapelet kernel at scale i and the image is formed by summing the cosine of the tilt differences between points on the image and the shapelet. Thus the gradient slant correlation must be multiplied by the tilt correlation measure that varies between 1, when tilts are aligned, to -1 when the tilts are in opposite directions. Cosine of the tilt angle difference is used to compute this relationship:

where τ_G and τ_si denote the tilts of the image and shapelet at scale i respectively. The overall correlation measure between image and shapelet kernel at scale i is obtained by taking the point-wise product of the gradient correlation and tilt correlation:

where · denotes point-wise multiplication and C _i is the ith shapelet image. The individual C _i’s are the images for different lengthscales and are the images used for quantification of the stratum SC and DEJ boundaries.

The boundaries of the epidermis are found by using a column-wise search algorithm that looks for local minima in the shapelet images corresponding to the different lengthscales. A piece-wise line search algorithm compares the median value of a sliding 3×3 window to the previous global median to find the first minimum and the first maximum in each column of image data. The mid-point between the two locations determines the SC surface position for each of the three short length-scale images. Then the SC curves are drawn across the image by connecting adjacent columns. The median of the three curves is used as the final SC line. A similar computation is performed on the eight long length-scale images to determine the DEJ curve. However, the least median of squares (LMS) statistic is used to determine the final DEJ curve [14], because of the large variation in position values between curves. The LMS is the most robust estimator of a curve under these conditions. The epidermal thickness calculation is then performed using the final SC and DEJ lines.

One advantage of using an IA method that includes information from several lengthscales is the ability to include an internal test of the quality of the results. The data at this point was examined for reproducibility at each point. We observed that individual point values relative to the LMS value was quite large for some points and not for others. On further examination of the data it became clear that these points of large variance in general, produced unreliable values. We therefore decide to add a variance test for each point. The standard deviation between the LMS point value and each of the eight individual point values was determined for each point on the lines. All points with a standard deviation between lines greater than 3 were excluded from the epidermal thickness calculation. By including results from only those columns that have good agreement in the location of the DEJ for all of the lengthscales, the shapelet analysis automatically excludes regions in the OCT image of high uncertainty and bolsters the robustness of the method.

3. Results

Results from shapelet analysis of OCT images of skin are displayed in Figure 4 and Figure 5. Figure 4 shows an OCT image where the DEJ is clearly visible and the line drawn by shapelet analysis matches well in many respects with those drawn by the human observers and produce thicknesses similar to values expected by A-scan analysis. The OCT image in Figure 5 is more difficult to analyze. The skin surface is not flat, and bubbles in the ultrasound gel used to reduce surface reflections degrade the signal in areas directly below. A-scan analysis is difficult in such an image because a single peak corresponding to the DEJ is not clearly identified. The shapelet analysis does not follow all of the contours of the dermal papillae as closely as the manual traces and instead charts a more slowly varying course across the image. These problems do not compromise the shapelet result because the varying skin surface is taken into account and columns with questionable results are excluded and marked with red circles.

Table 1. Comparison of epidermal thickness measured from 47 OCT images using A-scan, shapelet analysis, and manual line-drawing by three different people. Shapelet analysis performs similarly to human operators, and has a much lower error rate as compared to the A-scans.

View Table

Absolute comparison of the epidermal thickness as measured by different techniques is not straightforward due to the lack of an absolute standard for such thickness measurements. Without a gold-standard, we look for general agreement between techniques, or parsimony, to compare the different methodologies used for measuring epidermal thickness- A-scan, manual line drawing, and shapelet analysis. The concept is to identify which techniques most often disagree with the majority. As our effort is to match a feature apparent to a trained human observer, this is a reasonable test. The epidermal thickness was measured by each technique for a set of 47 OCT images. Manual line drawing was performed by three researchers. For each image, the median epidermal thickness was calculated from the five data points, and the absolute difference between the measurement and the median value was recorded for each technique. Table 1 shows the average deviation from the median for each technique, and examination of the results indicates that shapelet analysis performs similarly to the manual results, and that the A-scan method produces the largest deviation from the other techniques. The values of epidermal thickness as measured by all techniques falls within the range 60–110 µm with am average thickness of 88.2 µm and a standard deviation of 11.9 µm, which corresponds well with typical values reported in literature.

 

Fig. 4. Comparison of DEJ location as determined by trained operators (cyan, red, and magenta lines) and the shapelet analysis program (yellow line). The stratum corneum position is shown by the top green line, as determined by shapelet analysis.

Download Full Size | PDF

4. Discussion

The deviation of the A-scan results from the other methods is not surprising, given the inability of the A-scan to take into account the effect of features such as lines and wrinkles in the skin structure that influence only parts of an image. The similarity between results from shapelet analysis and manual processing, and the deviation from A-scan values is highlighted by performing a principal components analysis on the results of epidermal thickness measurements from each technique (Fig. 6). The automated measurement provides average results that are more reliable than A-scan results and comparable to those obtained by human observers. Furthermore, this approach increases reproducibility compared to manual traces and requires no human intervention.

 

Fig. 5. Identification of the DEJ in a non-ideal OCT image. The software correctly traces the DEJ even with a curved surface and bubbles in the ultrasound gel.

Download Full Size | PDF

Most skin examined in this paper has a relatively flat DEJ that does not display the characteristic egg-carton pattern associated with the dermal papillae. This flattening of the DEJ, normally associated with photodamage and aging, is typical of thin skin found on the face and forearm. Increased sun exposure also means that facial skin is more likely to develop cancerous lesions and therefore is the type of skin which needs to be studied to reach the ultimate goal, the identification of tumor margins for optical biopsy and to assist in Mohs surgery. The algorithm described here, however, fails to closely follow all of the contours of a wavy DEJ because the shapelet kernel is only processed in one direction - normal to the skin surface. OCT images that show a clearly defined wavy DEJ could be analyzed more accurately by allowing the shapelet function to rotate throughout all angles and highlight edges in all directions. This approach would require more computer processing time, however, and is not expected to significantly alter the epidermal thickness measurements for facial skin. We do plan to incorporate such functionality in future versions of this analysis.

The ultimate goal of this exercise is not the measurement of epidermal thickness, but rather to enable automatic and accurate feature identification in the complex OCT images of skin. The IA presented in this letter generally finds contours in these images despite the complicating presence of a speckle mask. It not only has potential application in the analysis of OCT of other tissues including the esophagus and the retina but also potential to analyze ultrasound images. Epidermal thickness is not only a parameter of skin health, but also a feature overlaying the boundary of basal and squamous cell tumors. Any robust method designed to perform optical biopsy or identify tumor margins will certainly need to identify basic structural features such as the DEJ. The work presented in this paper is a first step towards the fully automated analysis of skin OCT images.

5. Conclusion

The results presented here show that OCT images of human skin can be successfully interpreted by automated image analysis. The use of shapelet decomposition techniques to find soft boundaries within noisy images represents a novel approach that allows automated segmentation of the image and measurement of important parameters such as epidermal thickness. While there are many other features that are important targets for automated identification within these images, this work represents an important step towards fully automated image analysis of OCT images of skin. The existence of a successful and easily automated IA approach gives us confidence that there is further clinical utility of OCT for human skin. Advances such as phase contrast and polarization analysis [15] have shown potential in the laboratory but have yet to be incorporated in commercial equipment. Further improvements in both hardware and analysis will bolster use of the technique.

 

Fig. 6. Results of principal components (PC) analysis on epidermal thickness measurements for A-scan, shapelet analysis, and manual line drawing. Clustering of manual and shapelet results indicate good agreement in the majority of cases.

Download Full Size | PDF