Multi-Pass Adaptive Voting for Nuclei Detection in Histopathological Images

Nuclei detection is often a critical initial step in the development of computer aided diagnosis and prognosis schemes in the context of digital pathology images. While over the last few years, a number of nuclei detection methods have been proposed, most of these approaches make idealistic assumptions about the staining quality of the tissue. In this paper, we present a new Multi-Pass Adaptive Voting (MPAV) for nuclei detection which is specifically geared towards images with poor quality staining and noise on account of tissue preparation artifacts. The MPAV utilizes the symmetric property of nuclear boundary and adaptively selects gradient from edge fragments to perform voting for a potential nucleus location. The MPAV was evaluated in three cohorts with different staining methods: Hematoxylin & Eosin, CD31 & Hematoxylin, and Ki-67 and where most of the nuclei were unevenly and imprecisely stained. Across a total of 47 images and nearly 17,700 manually labeled nuclei serving as the ground truth, MPAV was able to achieve a superior performance, with an area under the precision-recall curve (AUC) of 0.73. Additionally, MPAV also outperformed three state-of-the-art nuclei detection methods, a single pass voting method, a multi-pass voting method, and a deep learning based method.

Nuclei detection is often a critical initial step in the development of computer aided diagnosis and automated tissue grading schemes in the context of digital pathology images [1][2][3][4][5][6][7][8][9][10][11][12] . However, nuclei detection is a challenging task in images with poor staining and noise. In breast cancer diagnosis, the Nottingham Histologic Score system is highly correlated to the appearance of cancer nuclei 13 . In melanomas, the quantitative assessment of melanocytes within the epidermis is an important cue for disease presence 14,15 . It is clear therefore that features pertaining to nuclear shape and spatial distribution of cell nuclei have important diagnostic value 12 . However, manually identifying the location and extent of melanocyte invasion or breast cancer nuclei can be subjective and time consuming. With recent interest in developing automated and computerized cancer diagnosis and grading schemes 7,[16][17][18] , there is a need for improved methods for nuclei detection and counting.
With recent advent of digital pathology, a number of groups have recently proposed approaches for nuclei detection 1,14,15,[19][20][21][22][23][24][25][26] . Most of these approaches make assumptions about the staining quality of the tissue and the resulting image quality and consequently may yield sub-optimal detection performance if the tissue preparation and staining is less than ideal. Due to the uneven absorption of the staining dyes by the tissue, variations in staining procedures adopted in various labs, variations in the exposure time for stain absorption by the tissue, differences in the quality of biochemical tissue staining and resulting appearance of the slides can be substantial. This can adversely affect the appearance of the tissue and the associated histologic primitives such as nuclei. For instance, over-staining of hematoxylin can cause adjacent nuclei to appear to be clumped together, over-staining with eosin can cause the nuclei to blend in with the stroma and result in under-emphasis of the nuclear boundaries. Uneven and imprecise staining can also cause hollow cores within the centers of individual nuclei. An example is shown in Fig. 1(b). Unlike the image shown in Fig. 1(a) where all nuclei are uniformly stained and have a consistent appearance, most nuclei in Fig. 1(b) appear to have an ill-defined nuclear contour, with a hollow interior. While a number of gradient based approaches for nuclear segmentation have been proposed, e.g., active contour 19 and level set schemes 21 , these approaches typically require model initialization. However, because gradient based approaches rely on the presence of strong edges for the model to latch onto, they will invariably yield worse performance on poorly stained images where the nuclei are either not uniformly stained or there are cavities within the nuclei. While deep learning (DL) based approaches 13 have the ability to learn discriminative features, the heterogeneous appearance of nuclei in different cohorts and often the availability of only a small number of representative training samples may hinder these supervised learning based approaches. Thus, it is necessary to develop improved, efficient methods to identify individual nuclei on pathology images.
Apart from the gradient based approaches, another class of nuclear detection approaches is the voting family of methods. By utilizing the inherent symmetric property of individual nuclei as a cue, voting based methods aggregate the votes provided by each pixel on the nuclear contour. Since the nuclear contour is generally symmetric, the center region is expected to have a higher value than other regions. While the hollow interior and ill-defined nuclear contour may affect other methods, the voting based methods 22,24,27 are typically able to handle the limitations of other families of nuclei detection methods due to their ability to ignore the heterogeneous appearance of nuclei and infer the center of nuclei from broken edge fragments.
In this paper, we present a Multi-Pass Adaptive Voting (MPAV) method that will adaptively select/modify the gradient information on nuclear contours. Unlike existing voting based methods, pixels that will join in the voting procedure will be adaptively selected and refined. This could potentially help alleviate the issue of incorrect voting, resulting in more accurate detection results. In MPAV, we assume that the shape of the objects of interest, i.e., nuclei, is convex. This is generally speaking a valid assumption since nuclei typically have elliptical shapes, which are convex. Table 1 enumerates some recent techniques for nuclei detection. Most approaches typically tend to use image derived cues, such as color/intensity 25,28-31 , edges 19,21,24,[32][33][34] , texture 35 , self learned features 13,36 , and symmetry 22,24,27,37 .

Previous Related Work and Novel Contributions
The color and texture-based methods require consistent color/texture appearance for the individual nuclei in order to work optimally. The method presented in ref. 31 applied the Laplacian of Gaussian (LoG) filter to detect the initial seed points representing nuclei. However, due to the uneven distribution of nuclear stain, the response of LoG filter may not reflect the true nuclear center. Filipczuk et al. applied circular Hough transform to detect the nuclear center 34 , however the circular Hough transform assumes that the shape of the underlying region of interest can be represented by a parametric function, i.e., circle or ellipse. In poorly stained tissue images, the circular Hough transform is likely to fail due to the great variations in appearance of nuclear edges and the presence of clusters of edge fragments.
Recently, there has been substantial interest in developing and employing DL based methods for nuclei detection in histology images 13,36 . The DL methods are supervised classification methods that typically employ multiple layers of neural networks for object detection and recognition. They can be easily extended and employed for multiple different classification tasks. Recently a number of DL based approaches have been proposed for image analysis and classification applications in digital pathology 13,36 . For instance, Xu et al. proposed a stacked sparse autoencoder (SSAE) to detect nuclei in breast cancer tissue images. They showed that the DL scheme was able to outperform hand-crafted features on multi-site/stain histology images. However, DL methods required a large number of dedicated training samples since the learning process requires a large number of parameters to be learned. These approaches therefore tend to be heavily biased and sensitive to the choice of the training set.
The key idea behind voting based techniques is to cluster circular symmetries along the radial line/inverse gradient direction on an object's contour in order to infer the center of the object of interest. An illustrative example is shown in Fig. 2(a,b). Figure 2(a) shows a synthetic phantom nucleus with foreground color as grey, and the background color in white. A few sample pixels/points on the nuclei contour with their inverse gradient  Figure 2(b) illustrates the voting procedure with three selected pixels on the contour. Note that for each pixel, a dotted triangle is used to represent an active voting area. The region where three voting areas converge can be thought of as a region with a high likelihood of containing a nuclear center.  Several effective symmetric voting-based techniques have been developed employing variants of the same principal. Parvin et al. 27 proposed a multi-pass voting (MPV) method to calculate the centroid of overlapping nuclei. Qi et al. 22 proposed a single pass voting (SPV) technique followed by a mean-shift procedure to calculate the seed points of overlapping nuclei. In order to further improve the efficiency of the approach, Xu et al. 24 proposed a technique based on an elliptic descriptor and improved single pass voting for nuclei via a seed point detection scheme. This initial nuclear detection step was followed by a marker-controlled watershed algorithm to segment nuclei in H&E stained histology images. In practice, the MPV procedure tends to yield more accurate results compared to the SPV procedure in terms of nuclei detection. The SPV procedure may help improve overall efficiency of nuclear detection 24 , however, it needs an additional mean-shift clustering step to identify the local maxima in the voting map. This additional clustering step requires estimating additional parameters and increases overall model complexity.
Since existing voting-based techniques typically utilize edge features, nuclei with hollow interiors could result in incorrect voting and hence in generation of a spurious detection result. One example is shown in Fig. 2(c), where we can see a color image, its corresponding edge map and one of the nuclei, denoted as A. Nucleus A has a hollow interior so that it has two contours, an inner and an outer contour, which results in two edge fragments in the edge map (see second row of Fig. 2(c)). For the outer nuclear contour, the inverse gradients are pointing inwards, whereas for the inner nuclear contour, the inverse gradients are pointing outwards. As one may expect, the inverse gradient obtained from the inner contour minimally contributes towards identifying the nuclear centroid (because the active voting area appears to be outside the nucleus, while the nuclear center should be within the nucleus). Another synthetic example of a nucleus with a hollow interior is shown in Fig. 2(c), and a few inverse gradient directions are drawn on the inner contour. In most cases, those inverse gradients from the inner contour will lead to a spurious result in regions of clustered nuclei. In Fig. 2(e), three synthetic nuclei with hollow regions are shown. It is clear that due to the vicinity of these three nuclei, the highlighted red circle region has received a large number of votes and thus could lead to a potential false positive detection. In section, we will show that in real histopathologic images, existing voting-based techniques tend to generate many false positive detection results.
In this paper, we present a Multi-Pass Adaptive Voting (MPAV) method. The MPAV is a voting based technique which adaptively selects and refines the gradient information from the image to infer the location of nuclear centroids. The schematic for the MPAV is illustrated in Fig. 3. The MPAV consists of three modules: gradient field generation, refinement of the gradient field, and multi-pass voting. Given a color image, a gradient field is generated by using image smoothing and edge detection. In the second module, the gradient field is refined, gradients whose direction leads away from the center of the nuclei are removed or corrected. The refined gradient field is then utilized in a multi-pass voting module to guide each edge pixels for generating the nuclear voting map. Finally, a global threshold is applied on the voting map to obtain candidate nuclear centroids. The details of each module are discussed in the next section and the notations and symbols used in this paper are summarized in Table 2.

Multi-Pass Adaptive Voting
For a pixel located at (x, y), an active voting area A(x, y; r min , r max , α, δ) is generated, where r min and r max indicate the radial range of the active voting area. r min is usually set to 1 pixel without loss of reliability (an active voting area can now be represented as A(x, y; r max , α, δ)). Let α represents the inverse gradient direction at pixel (x, y); δ indicates the angular range. An active voting area is a fan-like area illustrated in Fig. 4. Within the voting area, a kernel, normally a 2-D Gaussian Kernel K(u, v; σ, A) 22,24,27 , is placed at (u, v) to diffuse the votes. That is, pixels close to (u, v) will have received a higher number of votes and vice versa.

Gradient Field
Generation. An unsupervised sparse non-negative matrix factorization (SNMF) based color deconvolution technique 38 is first applied to separate different stains. We denote the single stain image as I and its 2-D image coordinates as (x, y). In order to reduce image noise, a smooth image I s is obtained by convolving I with a Gaussian kernel ) The inverse gradient map  is generated using the following equations: s represent the gradients along the X and Y axes, respectively. In this approach, an edge detection method based on Canny operator 39 is applied to I s to obtain nuclear edges. An example of edge detection is shown in Fig. 5(c), where white regions represent pixels that belong to nuclear edges. From here on, we only consider pixels that belong to nuclear edges. We denote the pixels set that belongs to a nuclear edge as E. An updated gradient map ′ is then generated as follows: Refinement of the Gradient Field. In the last step, we obtained an edge map that reflects which pixels will join in the voting procedure (shown in Fig. 5(c)). Since in the scenario that there are staining inconsistencies (see Fig. 2(c)), the inverse gradients on the inner edge fragment may not contribute to the voting of nuclear centroid and may introduce spurious votes in the voting map, we refine the inverse gradient of pixels that will join in the voting procedure. In other words, we do not blindly use all pixels and their corresponding gradient information during the voting step. Instead, we adaptively select the valid pixels and only use their associated gradient information during the voting stage. The steps corresponding to the gradient field refinement are described below.
Convex Hull Center Computation. For an edge fragment, E i,i∈{1,…,Q} , where Q is the total number of edge fragments in the edge map, we compute the corresponding convex hull, denoted as H i , and identify the center of the convex hull, denoted as C i . The convex hull region corresponding to the edge fragment in Fig. 2(c) is shown in Fig. 6(a). The red circle reflects the center of C i . The white region represents the convex hull region for the outer edge shown in Fig. 2(c), whereas the original pixels belonging to the edge fragment, i.e., pixels p(x, y) ∈ E i , are marked with rectangles and stars.
Angle Difference Computation. For a pixel in the current edge fragment, p(x, y) ∈ E i , we compute a vector, → C p i , formed by the current pixel, p, and the center of C i . We then compute the angle difference θ(x, y) between the vector → C p i and inverse gradient direction as follows: Two illustrative examples of the procedure for angle difference computation are shown in Fig. 7(a,b), respectively, where a hollow circle represents a pixel p(x, y) ∈ E i and a solid red circle represents the center of C i . The angle difference θ(x, y) between the vector → C p i and inverse gradient D(x, y) (shown via the solid blue arrow) is also shown in Fig. 7(a,b). Given a color image, the gradient filed generation module generates a gradient field by using of image smoothing and edge detection. In the second module, the gradient field is refined, gradients whose direction are pointing outwards the outerior region of nuclei are removed or corrected. The refined gradient field is then utilized in the multi-pass voting module to guide each edge pixels for generating the nuclear voting map. Finally, a global threshold is applied on the voting map to obtain candidate nuclear centroids.
Identifying Valid Voting Pixels. We employ the following heuristic to determine whether the gradient of a pixel is appropriate for participating in the voting strategy, where L p(x,y) denotes the label for pixel p(x, y). The rationale behind Eq. 4 is that if the inverse gradient points towards the hull centroid, it has a positive contribution towards the symmetry voting. On the other hand, if the inverse gradient points outwards from the hull centroid, it contributes negatively and may result in a false positive detection result. Illustrative examples are shown in Fig. 7. Figure 7(a) shows an example where the pixel is determined to be an invalid pixel, while Fig. 7(b) shows an example where the pixel is determined to be a valid pixel. In Fig. 6(a), all the invalid pixels are marked with black rectangles whereas all valid pixels are marked with red stars.
Invalid Voting Gradients Refinement. We present two gradients refinement strategies, Strategy 1 (GS1) and Strategy 2 (GS2), to refine the inverse gradient filed (we denote the MPAV with GS1 as M GS MPAV 1 , the MPAV with GS2 as M GS MPAV 2 , and the original MPV method as M MPV ). GS1 simply ignores the pixels which are labeled as 'invalid' . Mathematically, the refined inverse gradient map ′ x y ( , ) 1  is then calculated as follows: GS2 involves setting the inverse gradient directions of 'invalid' pixels to the reverse direction. Mathematically, the refined inverse gradient map ′ x y ( , ) 2 is calculated as follows:       Figure 6(b,c) show the refinement result calculated by using Eqs 5 and 6, respectively. It is clear that for GS1 shown in Eq. 5, the inverse gradients of all invalid pixels are set to 0. For the GS2 shown in Eq. 6, the the directions of inverse gradients for all invalid pixels are reversed. GS1 ignores pixels with 'invalid' gradients and reduces the number of pixels that used in the voting procedure, whereas GS2 tries to refine the gradient direction of pixels with 'invalid' gradient and attempts to infer a more accurate result. We do not suggest that either of GS1 or GS2 is necessarily better than the another, rather they represent two alternative strategies which the user might want to invoke for a specific application.
Multi-Pass Voting. The steps of MPV are summarized as follows.
Initialize voting parameters. Radial maximum range r max , maximum angular range δ max , and iteration number N are the voting parameters that need to be initialized. The angular ranges set is computed as follows, The initial voting directions α for all pixels are based on the inverse gradient map G′ (note that we have two different refinement strategies which correspond to two different refined inverse gradient maps ′ G 1 and ′ G 2 , respectively). That is, Update the vote image. For each pixel (x, y) ∈ E and G′ (x, y) ≠ 0, compute the vote image V i , for the ith iteration, using Eq. 8. In Eq. 8, K(u, v; σ, A(x, y; r max , α i , δ i )) represents the kernel placed at location (u, v) within an active voting region A(x, y; r max , α i , δ i ). Gradient map G′ can then be refined gradient map ′ G 1 or ′ G 2 .
Update voting direction. For each pixel (x, y) ∈ E and G′ (x, y) ≠ 0, we employ the following equation to update the voting direction, Intuitively, Eq. 9 involves updating the voting direction for pixel (x, y) with the pixel located at (u*, v*), which is the pixel that has the maximum vote value within the active voting region A, calculated in Eq. 10. In this way, we can perform the next pass/iteration of voting with a more precise voting direction, since the voting from the first iteration may be affected by the presence of noisy edge fragments that do not belong to the nuclei.
Iteration parameter adjustment. Using a smaller angular value, i.e., δ i+1 , set i = i + 1, and repeat steps 2.3.2 to 2.3.4, until i = N. The reason to reduce the voting area angle is to reduce the size of voting area. The idea is that by doing this we can identify a more precise nuclear centroid candidate region.
After N iterations, a final vote image V N is obtained. The final vote image V N contains the voting information from all pixels that join in the voting procedure. Figure  have high values near the ground truth nuclear centers (the manually labeled nuclei centers are shown in Fig. 9(d)). However, it should be noted that in the voting image appear to yield more specific and accurate detection results. Finally, a global threshold T D is applied to the voting image V N to obtain nuclear centers. Specifically, we use the automated Ostu's threshold method 40 to determine a global threshold T D . We then apply this threshold on V N to obtain a binary image. The centroids of the connected components are then identified as the final nuclear centers. Figure 9 on three different histopathologic image datasets. Table 3 describes the properties of these datasets. Dataset A is a publicly available histopathologic image dataset (H&E-stained digital images of breast cancer tissue slides) 41 . Dataset B consists of 10 breast cancer images stained with CD31 antibody and hematoxylin 13 . Dataset C consists of 5 skin cancer images stained with Ki-67 42 . Dataset D consists of 21 breast cancer tissue images stained with H&E. Note that the dataset A and D contain images with noisy background and nuclei with heterogeneous appearance due to the uneven staining. The nuclei locations were manually labeled by two experts, one a breast cancer pathologist with over 10 years of experience and the second an image computing scientist with over 3 years of experience in working with breast pathology. Figure 10 illustrates example images from A, B, C, and D, respectively. Evaluation Metrics. The main objective of the evaluation was to determine if the nuclear centers detected by the MPAV technique are concordant with the manually labeled nuclear centers. We calculate the center of each segmented region obtained by the automated technique. A nucleus was identified as having been correctly detected if its center is within a range of 15 pixels, i.e., about 3.6825 μm, of the manually labeled nuclei center location.
Denoting  ML as the total number of manually labeled nuclei centers,  DO as the total number of detected nuclei centers, TP  as the number of true-positives, (i.e., correctly detected objects compared to the manually labeled ones), FP  as the number of false-positives. (i.e., falsely detected objects compared to the manually labeled ones).
The performance is evaluated with respect to recall ( REC  ), precision ( PRE ), and F-measure ( F ) which are defined as follows:   13 . The M SSAE is a 4 layers deep neural network. The input patch size used in M SSAE was 34 × 34 pixels, therefore, the input layer receives 34 2 × 3 input units. The first and second hidden layers have 400 and 255 hidden units, respectively. The output layer is a soft-max function which classifies the input patches as nucleus or non-nucleus. The M SSAE was trained with a different set, other than datasets A, B, and C, of 37 H&E images of 2200 × 2200 pixels at 40x magnification 13 .

Results and Discussions
Performance comparison with existing voting based techniques. Qualitative results. Figure 11    are able to generate more accurate results with fewer false positives (automatically generated nuclear centers which are not concordant with the manually labeled nuclear centers).
Quantitative results. The detection performance in terms of PRC for datasets A, B, C, and D are shown in Fig. 12. The first row shows the performance of M SPV 24 compared to the M GS SPV 1 and M GS SPV 2 techniques (Fig. 12(a-c)). The second row shows the performance of M MPV 27 compared to M GS MPAV 1 and M GS MPAV 2 ( Fig. 12(d-f)). Figure 12 achieve better performance in such noisy images since they adaptively refine the gradient information for voting. In Fig. 12(b,f), it may be observed that the improvements obtained by the gradient refinement strategies  for dataset B are less substantial. In Fig. 12(c,g), the detection performances are almost identical. This is because the images in datasets B and C contain less noise that may affect the voting procedure. The refinement procedure thus has little effect on the voting result. In dataset A and D, there are many clustered nuclei, and hence, M GS for datasets A, B, and C, comprising a total of 11844 annotated nuclei. It is worth nothing that for M GS MPAV 2 , where the gradients direction are reversed (in Eq. 6), provides the best performance.
In Table 4 in terms of the best  F and its associated  PRE and  REC (fifth, thrid and fourth column, respectively). The area under the PRC curve for all plots, shown in Fig. 12, are also calculated and presented in Table 4. are shown in the first three columns (all nuclear centers are marked with blue crosses). An example image from dataset A is shown in the first row, M SSAE appears to produce a number of false positive detection results for this poorly stained image. Additionally M SSAE is prone to treat the boundary point as the nuclei center due to the heterogeneous intensity distribution within the nuclei and nuclei boundaries. In the second row, one may observe that the MPAV appears to produce more false positive errors for nuclei that have a thin and long shape (see Fig. 13(f,g)).
Quantitative results. The detection performance in terms of PRC for all 4 datasets are shown in Fig. 14. As shown in Fig. 14(a), the nuclei detection result for M SSAE on dataset A is poor, since the M SSAE is learned from dataset with well stained nuclei (and we also trained M SSAE with dataset A and then tested on dataset A, the result is even worse, possibly on account of too few training instances). In Fig. 14(c), M SSAE marginally outperforms are unsupervised methods that do not require training, whereas the M SSAE requires a large number of samples to learn the underlying patterns, and the heterogeneous appearance of nuclei in different stain conditions may hence affect the learned features and final detection performance.  Fig. 14 are presented in the last column of Table 5.
Statistics evaluation of detection results. We calculated the AUC for every image in the datasets and plotted the AUCs as a point for each image, as shown in Fig. 15. Figure 15(a-d) show the AUC cloud plots for dataset A-D respectively. A higher AUC corresponds to a better performance. From these cloud plots, one may observe the detection performance for each image within a dataset. The two-sample t-test p-values between each technique are calculated based on the AUCs obtain from different techniques and presented in Table 6. In this work, we used 0.05 as the significance level for the p-valudes to determine statistical significance. One may observe that in dataset A and D, the M GS  Fig. 15, we generated the PRC and calculated the AUC for each image whithin the dataset spearately.

Concluding Remarks
In this paper, a Multi-Pass Adaptive Voting technique was presented for automated nuclei detection on histopathological images. While most previous nuclear detection approaches tend to work well on well stained images, these approaches tend to fail on unevenly and imprecisely stained images. In practice, due to the variations of staining and slide preparation methods, a number of images tend to stain poorly. The MPAV is able to adaptively ignore or refine and identify the pixels that may lead to inaccurate nuclei centers in the voting procedure. To show the effectiveness of the MPAV, we compared it with existing voting based methods as well as a state-of-the-art deep learning method. The main contribution of the presented MPAV is to adaptively utilize the gradient information of pixels on edge fragments to generate more accurate detection results by exploiting the symmetry of   nuclei. Both qualitative and quantitative evaluation results show that the MPAV appears to be able to address many of the limitations of existing voting based techniques for nuclei detection in unevenly and imprecisely stained histology images. While compared to the DL-based method, MPAV provides a consistently superior detection performance whereas the DL method requires more training samples for poorly stained images. Note that the MPAV approach aims solely to detect nuclei and not to explicitly segment them, though the result of detection could serve as the initialization for other (e.g., watershed 24 ) nuclear segmentation approaches. However, our approach does have its limitations. Firstly, due to the noise and extraneous tissue components in the image, there may exist some noisy edge pixels outside the true nuclei on the edge segments. In these cases, the estimated center of the convex hull may fall outside the true nuclei, which leads to an inaccurate determination of which pixels are either valid or invalid for the purposes of voting. In our experiments, the number of such cases was small. Secondly, while our new gradient refinement strategy is able to improve the detection performance in noisy histopathological images, for relatively good quality images with well stained nuclei, like datasets B and C shown in Fig. 10(b,c), there is marginal or no improvement (detection result shown in Fig. 12(b,e,c,f)). However, the MPAV technique only has a marginal increment in terms of computational cost compared to the original MPV method (A Matlab implementation of MPAV algorithm is provided in the supplementary material C). In future work, we aim to integrate the MPAV with segmentation methods such as level sets and active contour schemes to explicitly extract nuclear boundaries.