1. Training samples

For training in the classification methods, 27 histological slides of 4 or 5 µm thickness were digitized at 20X magnification (0.46 µm/pixel) using the NanoZoomer HT Scan System (Hamamatsu Photonics, Japan). This scanner uses the line-scanning method (i.e., the slide is scanned by lines of 4096 pixels) and time delay integration (i.e., each line is scanned 64 times and averaged) to improve the signal to noise ratio. As mentioned in the introduction, the image sharpness is evaluated on high-resolution subregions of the VS, also denoted as tiles. We thus collected a total of 48,000 tiles of 200×200 pixels sampled from the 27 VSs taken at 20X magnification, which represented typical histological samples. The tile size enabled us to target specific VS regions, which presented different tissue and staining patterns and displayed either sharp or blurred areas.

As detailed in Table 1 and illustrated in Figure 2, tiles were sampled from VSs to account for two morphologic patterns: structured morphologies (SM, i.e., presenting well-defined histological structures such as colonic glands) and unstructured morphologies (UM, i.e., lacking clear histological structures such as brain tissue). These morphologic characteristics were combined with three staining patterns: HE, nuclear IHC (NS) and cytoplasmic IHC (CS), presenting variable expression levels of the targeted markers. We also included membranous staining patterns in the CS category (see Figure 2H), considering that (trans)membranous antigens, such as receptors, are also often expressed in the cytoplasm (due to internalization and protein trafficking).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Tile origins.

https://doi.org/10.1371/journal.pone.0082710.t001

The combination of morphologic and staining patterns resulted in six categories (see Table 1). For each category, we selected 4000 blurred and 4000 sharp tiles from a panel of VSs showing the patterns of interest (see Table 1). These different tile categories were included in the supervised dataset to best represent the large diversity of patterns that are encountered in HE and IHC slides.

2. Validation samples

To validate the method using independent slides, we selected another 97 VSs from our daily routine and imaged from 4- or 5-µm-thick tissue samples. For each slide, an expert randomly selected up to 20 blurred tiles and up to 20 sharp tiles. This new dataset of 3438 tiles (1462 from 39 HE slides and 1976 from 58 IHC slides) was then used as a validation set to estimate the performance of the classification algorithms (see section 3). This validation set included tiles that covered at least the six pattern categories included in the training data. The tiles were sampled from various tissue origins (e.g., skin, bladder, colon, pancreas, brain tumor, liver, and thyroid) and various IHC markers (e.g., KI67, glucagon, insulin, IGF1R, E-cadherin, CD79A, FHL2, GFAP, EGFR, HER2, etc.), including samples different from those encountered in the training set (see Table 1).

In addition to quantitative validation, we illustrate an application of this method to complete VSs (included in the independent set of 97 VSs), as detailed in section 3.3.

1. Tile blur characterization

A set of 16 different blur features was determined per tile to locally evaluate the image sharpness of histological VSs.

In [15], the authors presented two image analysis metrics: sharpness (SH) and noise (NO). SH was computed for tiles that had been converted to grayscale (i.e., the Y component of the CIE XYZ color space). The gradient direction was determined for each pixel whose gradient was larger than a threshold value (determined by Otsu’s method [14]). The edge width () was defined as the distance (in pixels) between the local minimum and local maximum of the gradient. As described in [15], the edge width is modified as follows:

with θ being the edge direction sampled at . The edge width was then averaged across the tile to yield SH:

where is the edge width of the i-th edge.

NO was computed by first subtracting a Gaussian-filtered (σ = 2) tile from the original tile. Using the resulting tile, NO was defined as the mean square of the minimum difference between each pixel and its eight nearest neighbors (in a 3×3 window). In contrast to SH, NO is the average of the NO value determined for the individual red, green and blue channels.

where I(i, j, c) and G₂(i, j, c) are the intensity values of channel c at row i and column j of the original and Gaussian-filtered tile (with σ = 2), respectively. R and C are, respectively, the numbers of rows and columns in the tile (i.e., R  =  C  =  200 in our application).

We added three other features. The first two used the difference image, (I - G₂), generated during NO computation and computed the mean value, or mean blur difference (MBD), and standard deviation, or standard blur difference (SDBD), across the difference image, respectively. The third feature is the mean value of the gradient, or the mean gradient magnitude (MGM), of a Gaussian-filtered tile (previously converted to gray-scale) with σ = 0.5 (G_0.5).

Similar to SH and MGM, the following features were computed for tiles converted to gray-scale. As suggested by [10], we computed a set of Haralick features, namely, the correlation and entropy features, computed on five gray-level co-occurrence matrices (GLCM). For a given tile I of size n x m, the GLCM computed for the offset (Δx, Δy) is defined as follows:

The normalized GLCM, labeled P, is obtained by dividing each element of C by the sum of all the elements of C:

The contrast feature is then calculated as follows:

with L being the maximum possible gray value (for our images, L  =  255).

Similarly, the entropy is calculated as follows:

For each tile in our database, we computed the contrast and entropy values corresponding to the GLCMs for the set of offsets , resulting in ten Haralick features.

For thoroughness, we also used the Tenengrad function (TG_t), as described in [12]. TG_t sums the square of the gradient values above a threshold t [13]. Because all of the tiles are the same size, the difference between TG₀ (i.e., TG_t with t  =  0) and MGM is that the gradient is computed for the original and Gaussian-filtered tiles, respectively. To prevent a strong correlation between the TG_t and MGM values, we altered the threshold of each tile by applying Otsu’s method to the squared norm of the gradient of the tile.

All of the features had larger values for sharp images compared to blurred images, with the exception of SH, which had larger values for blurred images [15]. All image features were computed in Python using SciPy [20] and the scikit-image image processing toolkit [21]. For performance purposes, we optimized the computation of the SH, NO, MBD, SDBD, MGM, HE_Δx,Δy and TG_t features by compiling the corresponding code to LLVM bytecode using the Numba library [22].

To compare the discriminatory abilities of the 16 features between sharp and blurred images, we used the θ measure, a non-parametric measure of effect size used to characterize the degree of separation of two distributions. This measure is easily estimated by U/mn, i.e., the Mann–Whitney U statistic divided by the product of the two sample sizes [23]. This normalized statistic ranges between 0 and 1, with values near 0.5 indicating similar distributions and values near 0 and 1 indicating strong separation. To simplify the comparison, we ranked the features by means of their discriminatory ability (or separation degree) evaluated by max(U/mn, 1-U/mn), which ranges between 0.5 (complete distribution overlap) and 1 (no overlap).

2. Tile classification and model selection

We used two different approaches of supervised classification: the decision tree (DT) approach, for which we tested three different split criteria, and the SVM approach which has been shown to be very efficient in cell classification [16] and for which we tested different kernels and hyper-parameter values. As detailed below, the DT approach enabled us to include feature selection in the classification process. Another advantage of this approach is that it generates relatively simple classification rules as a series of if-then statements, which are easy to implement in a decision process and are especially fast to apply in production (a particularly beneficial property). We also tested the SVM, which is a powerful, state-of-the-art algorithm with strong theoretical foundations that is particularly well adapted for binary classification in a numerical space [24]. This latter classifier enabled us to approximate the best results we could expect for our classification problems. As detailed below, for each classifier, we used a nested (5-fold x 5-fold) cross-validation to choose the hyper-parameters (inner 5-fold cross-validation) and to estimate the performance of the resulting model (outer 5-fold cross-validation), as is usually recommended to avoid overfitting [25].

The DT approach executes recursive partitioning of the feature space until it identifies subspaces in which the training cases predominantly belong to the same class (i.e., the majority class, which is allocated to the corresponding subspace). The DT was implemented using the classification tree module of the Statistica package (StatSoft, Tulsa, OK), except when the complete training set (i.e., grouping HE and IHC tiles) was considered. For this very large set (48,000 tiles), we used the QUEST Classification Tree software (version 1.10) [26]. These software packages allowed for the testing of different split criteria (for recursive partitioning), including univariate and multivariate split options. The univariate options are based either on an exhaustive search evaluated by the Gini index of node impurity [28] or on discriminant-based split criteria (based on ANOVA and k-means in the case of quantitative predictors, see [29]). The univariate split options resulted in a selection of the most discriminating features. The multivariate option determined a discriminant-based linear combination of all of the quantitative predictors provided for the model. In addition, minimal cost-complexity cross-validation pruning [28] was used to generate the “right-sized” tree for each split criterion. The complexity of a DT model is usually determined by its maximum depth, i.e., by the maximum number of splits required to determine the majority class in any region of the feature space. Finally, outer cross-validation (with inner cross-validation pruning performed on each training set) was used to evaluate the predictive performance of each algorithm.

The SVM approach aims to identify the boundary that represents the largest separation, or margin, between two classes. Initially developed for linear separation, the method was extended to solve non-linear problems by transforming the initial feature space (using kernels) into a new space in which the classes are (nearly) linearly separable [24]. Training a SVM classifier requires the selection of a kernel type (and its parameters) and a soft-margin parameter C, i.e., a penalty parameter for allowing some classification errors inside the margin (for nearly linearly separable classes in the transformed space). It can be shown that the SVM classification rule is a function of a subset of the training data that lie on the margin, i.e., the support vectors that determine the complexity of the class boundary implemented by the model [24]. The SVM was implemented using the SVC (Support Vector Classification) scikit-learn Python module [30]. A grid search was performed to identify the optimal hyper-parameters set by cross-validation. The grid search space was defined as

- kernels: { linear, rbf }

- C: {1, 2, 5, 10, 50, 100, 500, 1,000}

- : {10⁻⁵, 10⁻⁴, 10⁻³, 10⁻², 10⁻¹},

where is the coefficient of the radial basis function (rbf) used as the kernel [24].

Outer cross-validation was used to evaluate the predictive performance of the optimal SVM classifier (defined by the set of hyper-parameter values that gave the best predictive performance in the inner cross-validation step).

As detailed in the results, the different classification algorithms were trained with different feature sets as inputs. We thus implemented a model selection process based on a trade-off between the predictive performance (evaluated by means of nested cross-validation, see above) and model complexity to favor both accurate and quick decision-making in the application phase. Candidate models were first selected on the basis of a tolerance of 0.5% applied to the best accuracy observed among all the generated models. Among these successful candidates, we kept the model(s) that required the lowest number of blur features for computation (which is the most time-consuming task during the application). If several models remained in the resulting selection, we then chose the model with the lowest complexity, evaluated by the maximum depth for the DT and by the number of support vectors for the SVM. This complexity criterion also prevented overfitting in the model selection [25].

3. Quality map production and VS annotation

The map image was initialized by thresholding the image at 0.1X magnification of the entire slide. In the image, each pixel corresponds to a tile of 200×200 pixels at 20X magnification. To avoid computing the features and classifying tiles outside the tissue area, we first performed an Otsu thresholding procedure on the map image. Pixels corresponding to tiles outside versus within the tissue are were set to 0 versus 1. Tiles whose corresponding pixels were set to 1 were extracted at 20X magnification to compute the blur features and were then classified as being sharp or blurred. As illustrated in Figure 3, this classification resulted in a raw map showing black pixels (out-of-tissue area), grey pixels (sharp tiles) and white pixels (blurred tiles). This raw map may be difficult to interpret, especially in the case of scattered misclassification errors. Therefore, a post-processing tool, known as an alternate sequential filter [27] and consisting of a combination of a gray-scale morphological closing followed by a gray-scale morphological opening (with a 3×3 square as the structuring element), was applied to smooth the raw map and to ease the work of the operator. The contours of the blurred regions on the smoothed blur map were converted to an annotation compatible with NDP.View, the Hamamatsu annotation viewer. This conversion enables one to outline the blurred regions on the navigable VSs (Figure 3).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Application to complete HE and IHC VSs.

The first row shows the raw classification results obtained by classifying all of the tiles constituting an HE and an IHC VS, by means of the HE-specialized DT and the IHC-specialized DT (see results for details), respectively. Gray or white pixels correspond to tiles classified as being sharp or blurred, respectively. The second row shows the final blur map obtained after grayscale morphological closing followed by an opening. In the third row, the borders of the blurred (gray) regions identified in the map are transferred as annotations on the VS view provided by our scanner viewer. Magnified details for the small rectangles (see arrows) are provided in the last line. In these two magnified fields (scale bar  =  200 µm), the red lines indicate the boundary between a sharp and a blurred region. The blurred region is located in the left-hand region for HE and at the bottom for IHC.

https://doi.org/10.1371/journal.pone.0082710.g003

1. Feature analysis and selection

Table 2 shows the ranking of the features (in decreasing order) based on the discriminatory ability between blurred and sharp tiles. The discriminatory indices were evaluated by max(U/mn, (1–U)/mn) and ranged between 0.5 and 1 (see section 3.1). These indices and their rankings were computed for the complete training set as well as for the HE and IHC groups, taken separately. The results in Table 2 demonstrate high discriminatory ability for the HC_Δx,Δy family, followed by the gradient-based features (MGM and TG), NO and the HE_Δx,Δy family. SDBD, SH and MBD exhibit the weakest discriminatory ability. Interestingly, the feature ranking carried out for the complete training set is equivalent to the ranking average for the HE and IHC groups. This group distinction highlights some slight variations. For example, TG is more discriminatory in the HE tile group, whereas MGM is slightly more discriminatory in the IHC group.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Discriminatory ability of the blur features.

https://doi.org/10.1371/journal.pone.0082710.t002

Figure 4 details the data distributions of four discriminatory features across the six pattern categories constituting the training set. It should be noted that the level of discrimination between sharp and blurred images allowed by these features can vary from one category to another for a given staining group. These data demonstrate the benefit of including these different patterns in the dataset.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Distributions of HC₁₁, MGM, TG and NO values.

The feature values are computed for the training set. The results distinguish blurred and sharp tiles and are displayed by morphological patterns (structured vs. unstructured, i.e., SM vs. UM) in rows and staining types (HE, CS and NS) in columns. The data are shown by the following markers: median (black squares), 25–75% percentiles (boxes), non-outlier minimum and maximum (whiskers), and outliers (open dots).

https://doi.org/10.1371/journal.pone.0082710.g004

We also observed that the HC_Δx,Δy features obtained for the different GLCMs were highly correlated among themselves (Spearman’s rank correlation coefficients r_S > 0.95 for all training data). Even higher correlations were observed for the HE_Δx,Δy features (r_S > 0.98 for all training data). We thus selected the most discriminant feature in each family (i.e., HC₁₁ and HE₁₀, cf. Table 2) for further analysis.

We ultimately kept a set of eight features (HC₁₁, MGM, TG, NO, HE₁₀, SDBD, SH and MBD) to train the classification models. Because some classification algorithms do not include a feature selection process (see section 3.2), we trained the models with different feature sets provided as input. To generate these feature sets, we implemented a backward elimination procedure, which is based on the discriminatory ability shown in Table 2 (i.e., beginning with the complete set and then removing the features, one by one, from the least to the most discriminatory feature).

2. Supervised classification models

While detailed results are provided in the supplementary data (Table S1), the top of Table 3 summarizes the characteristics of the DT classifiers selected on the basis of their performance (evaluated by nested cross-validation) and complexity (see section 3.2). This table describes the results obtained with and without consideration of the two staining groups (i.e., HE and IHC vs. ALL), together with the results pooled for the two staining groups (see row “Pooled”).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Characteristics of the DT models selected on the basis of both their performance with the training set and their complexity and results obtained on the validation set.

https://doi.org/10.1371/journal.pone.0082710.t003

The accuracies (i.e., the percentage of correctly classified tiles) of the HE and IHC classifiers were 98.11% and 96.38%, respectively. Training only one classifier on all the data resulted in a global drop of approximately 1% in accuracy and sensitivity (i.e., the percentage of blurred tiles correctly classified as blurred by the classifier) (see lines “ALL” and “HE + IHC Pooled” in Table 3). All of these classifiers (selected for providing a good trade-off between accuracy and complexity) were trained with the Gini (univariate) split criterion. The difference between the HE and IHC classifier accuracies can be explained by the higher heterogeneity encountered in IHC slides. In addition, HE staining helps to introduce contrast in VSs because hematoxylin colors all cell nuclei in blue, whereas eosin stains most of the cell cytoplasm (and other extracellular eosinophilic structures) in various shades of pink (or even red). In contrast, IHC VSs mix IHC staining (with variable locations) with hematoxylin counterstaining, making the cell cytoplasm nearly invisible in areas in which the IHC biomarker is not expressed (see Figure 2). This difference may explain the increased difficulty for distinguishing between blurred and sharp tiles in the IHC group.

We also tested the impact of training more specialized classifiers for IHC VSs, one for nuclear staining and another for cytoplasmic and/or membranous staining (denoted NS and CS categories in Table 1, respectively). The accuracies of the selected classifiers (see Table S1) were 96.68% (NS) and 97.34% (CS). Their pooled performance (97.01%) was slightly higher than that of the IHC classifier (see Table 3). However, the specialized NS classifier required a computation of the eight blur features for the tiles (see Table S1), instead of only three for the generic IHC classifier (see Table 3). In addition, this generic classifier requires no distinction between staining locations. This latter property is interesting for future applications because subcellular antigen locations may vary in a single VS (e.g., due to possible protein translocation between the cell cytoplasm and the cell nucleus, as we observed in a previous study [31]).

The six pattern categories constituting the training set enabled us to identify which pattern(s) contributed the most to the classification errors. When using the HE-specialized DT to classify the HE tile group, the SM (structured morphology) pattern exhibited the largest (but acceptable) misclassification rate (slightly more than 2%, with a similar impact on the sensitivity), whereas the tiles presenting an UM (unstructured morphology) pattern were nearly always correctly classified (less than 1% of error) with a sensitivity near 100% for the blurred tiles. It should also be noted that the unspecialized DT strongly deteriorated the results for the SM-HE tiles (classification accuracy of approximately 93% and sensitivity of approximately 88%), whereas the results for the UM-HE tiles remained nearly similar (with a slight decrease in specificity). In the IHC group, the IHC-specialized classifier had the greatest difficulty in correctly classifying tiles presenting a combination of the UM pattern and the nuclear staining (NS) location (approximately 6% of errors relatively equidistributed between sharp and blurred tiles). In contrast, the UM-CS combination presented the best classification rate (near 99%). The two other categories (SM-NS and SM-CS) presented intermediate error contributions. In contrast to the HE group, the unspecialized classifier provided similar performances for all categories in the IHC group, without particular deterioration.

For completeness, we compared the DT approach to the SVM, a state-of-the-art classification approach. The SVM showed slightly improved performances compared to the DT approach, at the cost of a higher model complexity. Indeed, the best SVM classifiers, which were separately trained on the HE and IHC groups (with all of the features as input), exhibited classification accuracies of 99.06% and 97.23%, respectively. However, these two SVM models (kernel: rbf, C = 1000 and  = 0.0001) required the computation of dot products with 510 and 2216 support vectors, respectively, to classify a tile. Comparatively, the most accurate DT models, which were generated with all of the features as input, provided similar performances (i.e., 98.56% and 96.63% classification accuracy, respectively) but with maximal DT depths of only 6 and 8, respectively (see Table S1). These results imply that DT algorithms provide very simple classification rules, which allow for accurate and quick decision-making in tile classification.

In view of all these results, we therefore chose to use the DT approach for blur detection application and selected the two generic DT classifiers, one for each staining type (HE vs. IHC), as described in Table 3. The resulting models were tested on an independent set of VSs, as described in section 4.3.

3. Validation on an independent tile set and application to entire VSs

We first verified that the ability to discriminate blur features was conserved in the validation set of 3438 tiles. We obtained a ranking similar to that previously observed for the training set (see Table 2), with a decrease in the discriminatory indices, i.e., between 0.96109 for HC₁₁ and 0.65573 for MBD (evaluated on a much smaller dataset than the training set).

Table 3 (at the bottom) details the results obtained with the selected HE and IHC DTs on this completely independent dataset extracted from our daily routine in WSI and including samples vey different from those used for classifier training (see section 2.2). These results show an accuracy of near 90% for the two DTs (with a better sensitivity in the case of the HE DT). As previously observed on the training set, these performances decreased when the unspecialized (ALL) DT was used (in particular, the sensitivity for the HE tiles decreased to 82.54%).

In the final application (see Figure 3), all of the tiles constituting a VS should be classified. A post-processing step that considers the spatial consistency of the resulting tile labels (i.e., sharp or blurred) is then applied, withdrawing scattered classification errors. We found that this complete process takes approximately 15 minutes per image (averaged for the 97 slides of the validation dataset). After transferring blurry region borders as annotations in the VS view, the operator then only needs to verify consistently blurry regions to identify those that require additional focusing points to improve the VS quality (as shown in Figure 1C). This information may thus strongly reduce the manual reviewing time by focusing the operator’s attention on regions identified as blurry and for which the addition of new focusing points could be beneficial for locally improving the VS sharpness. If the operator confirms that the regions outlined in the VS are blurry, he/she can redistribute the focusing planes and the focusing points around the blurry regions and can then queue the slide for a rescan (see Figure 1C).

Code Availability

The source code for the Python tool presented in this work can be found at:

https://bitbucket.org/diapath/sharpa-lite.