Machine learning approaches to analyze histological images of tissues from radical prostatectomies.

2. Materials Radical prostatectomy specimens from 20 patients with a diagnosis of G3 or G4 prostate cancer according to the contemporary grading criteria [ 2 , 3 ] were retrieved from archives in the Pathology Department at our institution under an Institutional Review Board approval no. Pro00029960. Slides were digitized by a high resolution whole slide scanner SCN400F (Leica Biosystems, Buffalo Grove, IL) dedicated for pathology research. The scanning objective was set to 20× and the focusing was automatically adjusted by the scanner. The output was a color RGB image with the pixel size of 0.5 μm × 0.5 μm and 8 bit intensity depth for each color channel. We utilized freely available libraries from OpenSlide.org [ 20 ] to import Leica (.scn) images, select histopathologically important fields of view (FOV) that were converted to TIFFs for methods development. Furthermore, the FOVs were split into tiles of 1200 × 1200 pixels for image analysis. Of the total 5000 tiles, a subset of 210 images was selected by pathologists (SM, SB) who identified ST, BN glands, G3 cancer and G4 cancer containing cribriform and non-cribriform growth patterns. Depending on their content, tiles were categorized into groups ( Table 1 ), and then annotated manually using a custom graphical user interface (GUI) which we specifically developed for this task ( Fig. 1 ). The GUI facilitated: (a) import–export of images in common formats (tiff, jpg, png, etc), (b) free-hand contouring with the capability of an intuitive contour closing and space filling of the contour, (c) tissue labeling by colors: high-grade tumor (red), low-grade tumor (yellow), benign/normal glands (blue) and stroma (cyan), (d) easy correction of wrongly delineated areas, and (e) delineation of stroma by using a semi-automatic image flood-filling procedure after all glandular annotations are finished. At the end of this procedure all annotated image tiles were returned to pathologists for cross-evaluation.

3. Methods 3.1.

Overview of the image analysis strategy

Tissue classification based on H&E images is a challenging and computationally expensive process [ 13 , 19 ]. Our approach involves the stratification of prostate tissue in two sequential steps: (1) separating stroma (ST) from the epithelium (EP) and (2) differentiating of benign/normal glands (BN) from prostate cancer (PCa) ( Fig. 2 ). In each step only two tissue classes are analyzed at once. In step 1, a mask covering epithelial areas is generated to facilitate recognition of BN and PCa tissues in the next step. This strategy is significantly different from other published approaches [ 8 – 10 , 13 ], because it does not involve the per se glandular segmentation and extraction of glandular features. Many gland segmentation approaches known to date [ 9 , 21 – 24 ] are based a priori on the segmentation of nuclei, and use the glandular lumen to recognize the presence of glands. However, nuclear segmentations often fail in dense, overstained or understained tissue areas and glandular lumens can be occluded by corpora amylacia or secretions. Furthermore, empty spaces caused by tissue retraction artifacts or blood vessels can falsely be recognized as glandular lumens. In areas of high-grade cancer lumens are absent or difficult to detect. To overcome the limitations of existing methods, we formulate the problem as pixel-wise semantic segmentation based on intensity and texture descriptors that are extracted from images after color deconvolution. We apply novel and robust descriptors to distinguish BN from PCa. Due to the complexity of the cellular architecture in prostate tissue, different classifiers and descriptors are utilized at each step ( Fig. 2 ). Since non-epithelial tissues such as ST are morphologically less heterogeneous than BN glands and PCa, the computational expense can be decreased by removing ST before analyzing glandular structures. Our methodology was developed based on the following premises: (a) the density of epithelial nuclei is higher than the density of stromal nuclei, (b) the intensity of eosin is different in ST versus BN glands and areas of PCa and (c) the texture of hematoxylin and eosin in glands differs from the texture in ST. Hence, we utilize the difference in pixel intensity characteristics of eosin and hematoxylin to classify areas of prostate tissue. The classifiers are able to distinguish the three tissue classes based solely on the descriptors (intensity and/or texture histograms) that are extracted from hematoxylin (H) and eosin (E) image stains after color deconvolution [ 25 ]. The color deconvolution algorithm developed by Ruifrok and Johnston, was used to process RGB color images for digital separation of immunohistochemical dyes in a tissue image and performs well even if the dyes have overlapping absorption spectra or co-localize in same cellular compartment. Using optical density vectors of the pure dyes and an ortho-normal transformation of image intensities the algorithm deconvolves optical density of dyes for all image pixels. As a result, blue hematoxylin and pink eosin images are obtained as single channel intensity images. For details regarding color deconvolution the reader is referred to [ 25 ]. Histogram-based features are extracted by sliding a window W over H and E images with the interval of W /4 in horizontal and vertical directions. Descriptors are calculated for W and assigned to its central pixel. Descriptors for the remaining pixels are obtained via bicubic interpolation and then classified. Descriptors obtained this way form a three dimensional descriptor matrix with x and y dimensions that are same as in the input image and the z -dimension determined by the descriptor length. W slides from the top left corner to the bottom right corner vertically and horizontally throughout the image and stops systematically at W/4 intervals (defined by x and y ) to calculate the descriptor (histogram). Subsequently, the descriptor is inserted to the three dimensional matrix and at x and y positions defined by W /4 intervals. Descriptors at image borders are calculated after mirroring pixels that are up to W /2 away from the image border. The mirroring provides a sufficient number of pixels to calculate a descriptor at pixels for which W positioned at the upper left corner and along the borders of the original image. A bicubic interpolation is then applied to the three dimensional descriptor matrix to obtain values at positions not visited by W . It is applied to each z-plane separately using a third degree polynomial and values from a 4 × 4 grid of descriptor components at which W stopped. Prior to color deconvolution, images are color-normalized [ 26 ] to account for variability in staining intensity. Bright areas (background and lumens) are removed by thresholding of gray-level images obtained from color-normalized H&E images. This intensity threshold was manually optimized and set to the value of 210. The resulting binary mask of the background is cleaned up by removing isolated blobs with areas smaller than 50 pixels followed by a morphological closing with a flat disk-like structuring element with the radius of 7. Pixels under the background mask are excluded from analysis. 3.1.1. Intensity and texture histograms of hematoxylin and eosin stained images Prediction of tissue content requires a set of descriptors–numerical features extracted from images that have the capacity to capture differences in tissue morphology. Tissue descriptors utilized in our study are divided into two groups: intensity histograms and joint distributions of the uniform rotation-invariant local binary patterns and rotation invariant local variance. The first type of histograms represents a distribution of pixel intensities, whereas the second histogram represents the spatial relationships (mostly differences) between intensities of pixels, which have a close proximity. Thus, the second type is an image texture descriptor. Our approach utilizes texture and intensity based features that are extracted from hematoxylin and eosin images after color deconvolution ( Fig. 3 ). For each image pixel the descriptors: (a) intensity histograms: Hist( H ) and Hist( E ) and (b) joint distributions of the uniform rotation-invariant local binary patterns (LBP) [ 27 ] and rotation invariant local variance (VAR) [ 27 ] denoted as LBP P , R riu 2 × VAR P , R ( H ) ( H ) and LBP P , R riu 2 × VAR P , R ( E ) ( E ) respectively for hematoxylin and eosin. P and R denote neighborhood size and radius for LBP and VAR. riu2 stands for the uniform rotation-invariant type of LBP. To obtain the LBP P , R riu 2 × VAR P , R descriptor the intensity and contrast insensitive histogram LBP P , R riu 2 is combined with a single value representing image contrast VAR P,R . The resulting LBP P , R riu 2 × VAR P , R descriptor is a one-dimensional histogram with length k × n , where: k is the number of quantizing levels of VAR P,R and n is the length of the LBP P , R riu 2 histogram. Since each image pixel can contribute one LBP P , R riu 2 and VAR P,R it makes sense to calculate LBP P , R riu 2 × VAR P , R for pixels in a window W . For instance, for a 10 bin long LBP P = 8 , R = 1 riu 2 histogram and VAR P=8,R=1 quantized as 1 the bins of L B P P = 8 , R = 1 r i u 2 histogram are added to the first 10 bins of LBP P = 8 , R = 1 riu 2 × VAR P = 8 , R = 1 . If VAR P=8,R=1 is quantized as 3, then the bins of LBP P = 8 , R = 1 riu 2 are added to the third tenth of bins in LBP P = 8 , R = 1 riu 2 × VAR P = 8 , R = 1 . We quantized VAR P,R into 8 bins. After calculating all LBP P = 8 , R = 1 riu 2 histograms in W and sorting them according to values of VAR P=8,R=1 , the joint LBP P = 8 , R = 1 riu 2 × VAR P = 8 , R = 1 histogram will be 8 × 10 = 80 bin long. Source codes for calculating LBP and VAR can be found here [ 28 ]. The Hist( H ), Hist( E ) and LBP P , R riu 2 × VAR P , R descriptors are derived from a N × N pixel window W centered at each image pixel. Intensity histograms Hist(.) are sorted into 18, and each LBP P , R riu 2 × VAR P , R into ( P + 2) × 8 equally spaced bins. Paired values of P = {8,24} and R = {1,3}, are tested for best classification performances. For deconvoluted H and E images the histograms: Hist( H ), Hist( E ), LBP P , R riu 2 × VAR P , R ( H ) and LBP P , R riu 2 × VAR P , R ( E ) are concatenated, and then used as descriptors for tissue classification. The window size N depends on scanner resolution and needs to be chosen to provide enough data points to derive meaningful descriptors. According to [ 17 ] several G3 glands fit W with N = 64 and may be sufficient to obtain spatial precision in predicting ST and EP. Since its impact on epithelial stratification is unknown, we tested performance of descriptors obtained for N = {64, 128, 256}. The prediction of ST and EP regions is conducted using full size images (highest resolution). Prior to extracting the descriptors for the epithelial stratification ( Fig. 2 ) the image was downsized by 50% with bicubic interpolation to reduce computation and storage costs.