In this section, the skin explant assay used to generate the skin samples is described, followed by the image acquisition procedure and the methodology used to pre-process, threshold and fine-tune the segmentation of the epidermis tissue. Figure 2 summarises the main stages of the segmentation algorithm. The algorithm is based on the differences in the colour and intensity of staining in the different tissue types, the texture within each tissue and the overall shape and size of tissue regions. A colour normalisation step is also included to handle variations in the chromacity of the stained tissue caused by differences in sample thickness, staining procedure, and lighting during sample preparation and image acquisition. The algorithm was implemented using the Image Processing Toolbox™ in MATLAB®, Version 7.11, R2010b (The MathWorks, Inc.).
Data source: skin explant assay
The skin samples forming the test set were generated in a skin explant assay [31] developed to predict and investigate the immunobiology of graft versus host disease occurring post hematopoietic stem cell transplantation in HLA-matched siblings. Punch skin biopsies of 4 mm diameter were taken from the back below the waist-line of transplant patients after informed consent. Each 4 mm2 skin biopsy was equally dissected into 4 small sections and co-cultured with pre-primed alloreactive T cells. Skin sections cultured in medium alone were used as controls. After 72 hr incubation at 37°C, skin sections were fixed in formalin, sectioned and stained with H&E. The manual histopathological grading was assessed and confirmed independently by two experts. The histopathological grading system was as follows: grade I, mild vacuolisation of epidermal basal cells; grade II, diffuse vacuolisation of basal cells with scattered dyskeratotic bodies; grade III, subepidermal cleft formation; grade IV, complete epidermal separation. The assay is described fully in Sviland and Dickinson [31]. In addition to the vacuolisation, cleft formation and the presence of dyskeratotic bodies, some of the images also included regions of necrotic tissue. Ordinarily, samples with necrotic tissue are not manually scored and biopsies with such artefacts are not included in the standard assay readout, however these images were included here to enable the software to identify and ignore artefacts or necrotic regions. In the 40 sample data set of skin explant slides, there were eleven grade I examples, twelve grade II, nine grade III and eight grade IV.
Ethics statement
All skin samples were obtained from patients after informed consent. All human material was collected and stored using protocols approved by Local Research Ethics Committees (Newcastle and North Tyneside Hospitals NHS Trust).
Image acquisition
Digital images of the sectioned and stained skin samples were created using a Zeiss AxioImager II system with MOSAIC tiling facility. The image tiling and stitching facilities were used to create single x10 magnification RGB images of each skin sample. The image dimensions and aspect ratios varied, with heights of 1047 to 4819 pixels, and widths from 2676 to 5254 pixels. Figure 3 shows examples of stained skin sections with damage grades I-IV (Figure 3a-d) with the epidermis tissue outlined in green, and an example of a typical whole slide image generated by the Zeiss AxioImager system (Figure 3e).
The remainder of the methods section is used to describe the eight individual steps in the epidermis segmentation algorithm in detail.
Sample segmentation and image cropping
The first stage in the algorithm is to segment the pixels in the image representing the skin sample. This first segmentation increases the efficiency of the algorithm by limiting the number of pixels being processed during subsequent steps. While segmentation of the skin sample could be achieved by locating either the background or the sample pixels, the background pixels are used as they have lower inter and intra-image variance. The background pixels are located by creating a composite image, K, by summing the red, green and blue (R, G and B) intensities for each pixel in the RGB image, then taking the most frequently occurring value in K to approximate the background colour and storing this value as the background threshold, bg
thresh
. Black pixels present at the image edges due to the image tiling procedure are excluded by determining the mode of non-zero elements in K.
The black pixels at the image edge (as shown in Figure 3e) are replaced with the background threshold intensity, bg
thresh
, to create a consistent background. To reduce memory requirements in the implementation of this algorithm, excess outer rows and columns of background pixels which do not intersect the sample are cropped at this stage. For an m x n size image, this is done by cropping any rows where the sum of composite pixels in the row is less than bg
thresh
*n (an average value for a background row), and cropping any columns where the sum of composite pixels in the column is less than bg
thresh
*m. The process of cropping is shown in Figure 4 for a particularly challenging image showing grade IV damage. The original image is shown in Figure 4a, and the in Figure 4b all excess background has been cropped without cropping more than a few pixels of tissue at the sample edges.
Prior to thresholding the colour image is smoothed using an averaging filter which replaces each image pixel with the mean value of its pixel neighbours in each colour channel. This has the effect of reducing variation within the background pixel and the sample pixel sets and facilitating the choice of threshold when creating a binary sample mask. The operation is performed using convolution with a kernel filter to represent the pixel neighbourhood (K
smoothed
= K * kernel). Initially a filter size of 19×19 (where each element is 1/(19*19)) was chosen based on the typical size of vacuoles and clefts in the image which needed to be smoothed.
A new binary image sampleMask is created by thresholding the smoothed image, K
smoothed
. The value of the threshold was based on bg
thresh
, however since bg
thresh
is a measure of central tendency, the actual threshold used must be lower to ensure the majority of background pixels are below it. To calculate how much lower the threshold needed to be, the standard deviation of the background pixels in all 40 smoothed composite images was found. The standard deviation ranged between 1.2 and 3.2, subtracting three standard deviations of the image with the highest variance (3*3.2) from bg
thresh
will mean the threshold will fall below the intensity of the majority of background pixels.
The effect of the averaging filter and its size on the thresholding operation is shown in Figure 5. A variety of kernel filter sizes were tested to determine their effect on the thresholding step. The post thresholding binary mask created after pre-processing with no smoothing is shown in Figure 5a, with a 9×9 mean filter kernel in Figure 5b, a 29×29 filter in Figure 5c, and a 49×49 filter in Figure 5d. The dimensions of the mean filter must be large enough to smooth the coarse texture in the lower parts of the dermis and clefts at the dermal-epidermal junction so these features are included in the sampleMask, while minimising loss of accuracy at the perimeter of the sample. Without any smoothing the thresholding results in a mask which is very detailed but does not include any vacuoles or clefts. The smoothing creates a simpler mask which includes gradually more of clefts, vacuoles and white regions within the dermis as the filter size is increased. While mean filters of sizes between 9 and 49 can be used successfully, an intermediate value of 29 was chosen for two main reasons. Firstly, at this size the thresholded mask tends to have fewer separate objects than when smaller filters are used meaning a hole filling operation can be used subsequently to create a simple mask including clefts and vacuoles. Secondly, the 29×29 filter size was chosen in preference to a larger filter as the larger filter could result in a loss of segmentation accuracy at the sample boundary. The original RGB image and the image after smoothing with the 29×29 filter are shown in Figure 5e and f; there is a strong blurring effect which removes detail from the internal parts of the tissue.
In some images the dermal epidermal clefts are very large the smoothing operation is not enough to include them in the sample mask as foreground objects in the sampleMask (as seen in Figure 5). Mathematical morphology is used to in-fill these “holes” in the binary sample mask and also to remove small objects such as dust or tissue debris on the slide which have been captured during thresholding, but which are not useful for subsequent analysis. The sequence of operations used to refine the segmentation is as follows.
Step 1: Fill holes - Fills internal regions of background pixels within foreground objects in the binary image using the MATLAB, imfill function, an implementation of morphological reconstruction described in [32].
Step 2: Remove small objects - Removes foreground objects that consist of less than 25,000 connected pixels. This value was chosen so that the smallest fragments of tissue found in the image set used in this research were not excluded, but the value was high enough to exclude smaller objects such as dust or other tissue debris.
Colour normalisation
The initial optimisation of the method and parameters was carried out without colour normalisation, however the addition of this step improves the performance in terms of sensitivity, specificity and overall accuracy. The relative improvement is described in the results section. Staining inconsistencies in the input images are addressed by mapping the histogram for each individual colour channel of the RGB image to those of a target image, I
ref
, identified as well stained by an expert histopathologist. Only the sample pixels identified in the previous sample segmentation step are included in this colour normalisation.
The best approximation of the target image is obtained through the application of a grayscale transformation, T, to the sample pixel intensities, k, in the input image. A transform is selected for each colour plane so as to minimise the difference between the cumulative histogram of the transformed input image intensities, c
input
and the cumulative histogram of the well stained target histogram, c
ref
. The function to be minimised is
which can be implemented in MATLAB using the function histeq[26].
The effect of the colour normalisation on two images with significant differences in staining and lighting is shown in Figure 6. The original images are shown in Figure 6a and in Figure 6b the non-sample pixels have been changed to white and the sample pixels have been normalised. The final two images have a similar contrast between the epidermis and dermis, and a similar range of colour hues and saturation to each other.
Colourspace conversion
The next part of the segmentation procedure is a coarse segmentation of the epidermis based on the thresholding of a high contrast composite image. The composite image is created using a weighted linear combination of contrast stretched images in colour spaces which enhance the contrast between the regions of interest and the rest of the image.
A number of colour spaces were investigated to identify a representation that would maximise the contrast between the epidermis and the rest of the skin tissue. Those tested included CMYK (cyan, magenta, yellow and black) which is based on subtractive colour mixing, HSV (hue, saturation and value), YCbCr (luminance, blue chrominance and red chrominance) and the L*a*b* (lightness, red/green, yellow/blue) colour space [27].
The Cb channel in the YCbCr colour space and the b* channel in the L*a*b* colourspace were observed to provide good contrast between the epidermis and the rest of the tissue. The Cb (blue chrominance) and the b* (yellow/blue) colour channels both highlight the blue staining of haematoxylin which stains the nuclei in the cells of the epidermis. Although there are nuclei-containing cells present in the dermis, they are few in number.
When the two colourspaces were tested on a set of 16 images chosen to include different damage levels and staining variation, with optimised contrast enhancement, the b* image enhanced the contrast of the tissue types more consistently than the Cb plane. The b* plane of the L*a*b* colourspace was therefore chosen for use in the algorithm. The CIE L*a*b* specified by the International Commission of Illumination in 1976, is a perceptually uniform colourspace based around human perception of colour. In the CIE L*a*b*, colour is represented on three axes: (1) lightness, with the maximum value a perfect reflecting diffuser and the minimum representing black, (2) a red to green axis, with a positive a* value showing red, and a negative a* showing green, and (3) a yellow to blue axis, with a positive b* value showing yellow, and a negative b* showing blue.
It was noted that during optimisation of the subsequent contrast enhancement and thresholding stages that when the RGB image was converted to grayscale and contrast enhanced, good contrast for the epidermis was observed in some images where the b* colour channel was displaying poor contrast. The images showing poor contrast of the epidermis in the b* colour channel tended to have weak nuclear staining by haemotoxylin, which appears as a strong blue/purple colour and therefore stands out in this yellow/ blue colour channel. The complement image of the grayscale representation highlights more intensely stained areas, despite not being as specific to particular colours as the b* plane and tends to highlight both the cytoplasm and nuclei in the epidermis which are usually stained more intensely than the dermis tissue. When tested on a set of 25 images, the image variation meant that in different images the epidermis was highlighted best in either the b* or the grayscale image, it was therefore decided that a combination of the data in the grayscale and b* images could be used to enhance the robustness of the following steps.
Contrast enhancement
Contrast enhancement by stretching transforms a low contrast image, b*, to a high contrast image, b’, by remapping the gray levels to cover a wider dynamic range. A linear transformation maps the lowest gray level in b*, GL
min
, to a new minimum gray level GL’
min
, and the highest gray level in b*, GL
max
, to a new maximum gray level GL’
max
. The linear transformation that preserves the intensity histogram shape is given by:
where INT returns an integer value. GL
min
and GL
max
can be replaced with the penetration points P
min
and P
max
, to remap a band of intensities. Penetration points can be determined using a cumulative percentage histogram; for example, selecting the intensities to exclude the top and bottom 1% of the data. This is a useful technique if the intensity band of interest is at a mid-gray level rather than at an extreme gray level.
Only sample pixels were included in the contrast enhancement process, as the aim was to maximise the contrast between the dermis and epidermis and background pixels are not relevant to the rest of the process. Remapping a narrow, more specific band of intensities was investigated to try and improve the contrast. When tested manually using a variety of absolute intensity levels as penetration points, the optimal intensity band for remapping to enhance contrast of the epidermis varied significantly for different images. This issue was addressed by determining penetration points based on the cumulative percentage histogram so that a set percentage of low and high intensity pixels were saturated in the final image. This method is more able to cope with any remaining staining variation and pixel intensity outliers than choosing absolute intensity levels. This approach was used on both the b* and grayscale image planes to create new contrast enhanced images, b’ and G’.
Following the contrast enhancement, the two images G’ and b’ were smoothed using an averaging mean filter, as described for sample segmentation. The operation is performed using convolution with a kernel filter to represent the pixel neighbourhood (K
smoothed
= K * kernel). This has the effect of reducing variation within the sample pixels and smoothing minor variations within the epidermis and dermis regions and leaving the main variations due to differences in the two tissue types. This will facilitate the choice of threshold in method step 6.
As the contrast enhancement and smoothing are critical steps in the process, the optimal values for the upper and lower penetration points for the grayscale and b* images were determined along with the optimal sizes of smoothing mean filter and the structuring element used in the morphological processing (method step 7). This optimisation varied the 6 key parameters within defined ranges using a Design of Experiments approach. A quadratic model was then built and used to select parameter values which maximised sensitivity and specificity. The process and results are described in detail in the results section.
Based on the optimisation the 29.1% of pixels with lowest intensity in the grayscale image G’ were saturated and the remaining pixels remapped to utilise the full dynamic range. For the red-green b’ image, the 36.7% of pixels with the lowest intensity were saturated and the remaining pixels remapped to utilise the full dynamic range. A mean filter size of 41×41 (where each element is 1/(41*41)) was chosen during the optimisation.
Linear combination
An equal weighted linear combination is then applied to the two enhanced and smoothed images G’ and b’ to create a new image, Gb:
Combining the two images captures both staining intensity and staining colour information in a single grayscale image.
Thresholding
Following colour normalisation and contrast enhancement, Otsu’s automated method [23] was applied to determine the optimal threshold based on the intensity distribution of sample pixels in the combined, enhanced image, GB’. The Otsu method uses discriminant analysis to determine a threshold, Level, which maximises the separability of two pixel classes by minimising the intra-class variance. In this process, the aim is for this method to maximise separability of dermis and epidermis pixels. GB’ is converted into a binary image, BW, using the threshold, Level. Any non-sample pixels are changed to black (as background).
Figure 7 shows a histogram of a typical GB’ image. The method assumes an approximately bimodal distribution. The Otsu threshold, Level, generated using Otsu method is labelled in the figure at the intersection of the upper and lower intensity components.
Morphological processing
Morphological processing is used to further process the binary image, by removing small misclassified objects, merging multiple objects, in-filling holes and closing gaps. These operations are applied to the whole image, however once the operations are completed, any non-sample pixels which may have been affected are changed back to black. The sequence of operations summarised below is used to refine the segmentation. The choice of structuring element size (radius = 8) for the morphological closing and opening steps was optimised based on the final sensitivity and specificity of the algorithm as described in the results Section. A disk shaped structuring element was used as this shape reflects biological structures more accurately than sharp angles or linear shapes.
Step 1: Morphological closing - Morphological closing (dilation then erosion) enlarges the boundaries of foreground (bright) objects in the image and closes gaps between them, and shrinks background-coloured holes in the foreground objects. A disk shaped structuring element with radius = 8 pixels is employed.
Step 2: Morphological opening - Morphological opening (erosion then dilation) removes some of the foreground (bright) pixels from the edges of foreground objects, which will break fine bridges between objects while preserving the object size. A disk shaped structuring element with radius = 8 pixels is employed.
Steps 1 and 2 combine to smooth and simplify the objects edges without changing the size of objects. Simple object perimeters are useful in the image classification algorithm for which this segmentation algorithm is being developed; however these steps could be omitted if this is not an important output of the segmentation.
Step 3: Remove small objects - Remove foreground objects that consist of less than 4000 connected pixels. The threshold of 4000 pixels was decided based on a size assessment of any regions of epidermis identified as epidermis objects at this stage of the process. The majority of epidermis regions at this stage consist of more than 4000 pixels.
Step 4: Fill holes - Fills internal regions of background pixels within foreground objects in the binary image that consist of less than 7000 connected pixels. A threshold is required as in some images there are regions of dermis tissue surrounded by epidermis tissue (due to tissue slicing technique) and if these regions are filled the specificity of the final algorithm is compromised. Again the value was chosen based on the typical size of enclosed dermis regions and other structures within the epidermis which should not be included.
Object classification
Following morphological processing, the binary mask includes objects that are not part of the epidermis. These include collections of cells within the dermis, which have a similar appearance to epidermis tissue, and parts of the dead surface layer, the stratum corneum, which may be segmented with the epidermal tissue in cases where it is highly stained. For each object, Z, the object area, Z
Area
, and the area of the object’s bounding box, Z
BoundBox
, are determined. The ratio of Z
Area
to Z
BoundBox
gives the extent, Z
Extent
, of the object, which is a useful measure for differentiating between the long thin objects of the epidermis and the more compact, circular clusters of cells within the dermis:
The Z
Area
and Z
Extent
can be used to classify the remaining objects as either epidermis or non-epidermis. The thresholded objects are either retained or removed based on their area and their extent and so the effect of changing the area and extent thresholds used in this rule on algorithm sensitivity and specificity was tested. Only objects with an area greater than the area threshold and an extent less than the extent threshold were retained. The values used in this classification were optimised together once all the other critical factors had been set at optimal values, the optimisation is described in the results section. Based on the optimisation the following classification rules were used to classify each object pixel, z, in the binary mask.
where z are pixel elements in the object Z.
User interaction
The object classification step can be used to tune the specificity and sensitivity of the final algorithm, however the algorithm also includes the option for the user to interactively include or remove objects for the final epidermis mask. Epidermis segmentation is critical to the performance of the subsequent steps in the skin damage classification that is the ultimate aim of this research, and so this step is available to improve the performance of the algorithm. It is relatively straightforward for a user to determine whether a given object is part of the epidermis when shown next to an image of the RGB image. The user has the option to (1) approve the object selection, (2) click on objects to remove them, or (3) select additional objects, which were removed during object classification.