Smart Region-Growing: a novel algorithm for the segmentation of 3D clarified confocal image stacks

Motivation Accurately mapping the brain at the micro-scale is still a challenge in cellular neuroscience. While notable success has been reached in the field of tissue clarification and confocal imaging to obtain high-fidelity maps of three-dimensional neuron organization, neuron segmentation is still far away of ground-truth and manual segmentation performed by experts may be needed. The need of an expert is in part related to the limited success of the algorithms and tools performing single-neuron segmentation from 3D microscopic image data available in the State of Art, in part to the non-complete information given by these methods, which typically perform neuron tracing and thus limit the interpret-ability of results. Results In this paper, a novel algorithm for segmenting single neurons in their own arrangement within the brain is presented. The algorithm performs a region growing procedure with local thresholds based on the pixel intensity statistics typical of confocal acquisitions of biological samples and described by a mixture model. The algorithm is developed and tested on 3D confocal datasets obtained from clarified tissues. We compare the result of our algorithm with those obtained by manual segmentation performed by 6 different experts in terms of neuron surface area, volume and Sholl profiles. Statistical analysis performed on morphologic features extracted from the segmented structures confirms the feasibility of our approach. Availability The Smart Region Growing (SmRG) algorithm used in the analysis along with test confocal image stacks is available on request to the authors. Contact alejandrocallara@gmail.com Supplementary information Supplementary data are available on request to the authors.


Introduction
Digitalizing a complete and high-fidelity map of the neurons populating a mammalian brain is a central endeavour for neuroscience research, since it represents the first step for the delineation of the full Connectome (Alivisatos et al., 2012). Single-neuron reconstructions from empirical data can aid generating models to make predictions about higher-level brain organization, studying the normal development of dendritic and axonal arbours and documenting neuro-pathophysiological changes (Budd et al., 2015). Doubtless, confocal and two-photon microscopy are the cardinal tools for studies of the cellular morphology and brain cyto-architecture (Wilt et al., 2009;Ntziachristos, 2010). These cutting-edge technologies, coupled with clarification methods, that act on the penetration limits imposed by light scattering (see the review of (Richardson and Lichtman, 2015) for more details) making the biological tissues transparent, have enabled the inspection of the global arrangement of large brain neuron populations (i.e. up to 500 µm in thickness) at cellular (i.e. sub-micrometeric) resolution. A fundamental limitation to unravel the intricate brain micro-architecture is the limited success of the algorithms and the tools available in literature. In fact, despite the efforts reported by the scientific community, as witnessed by the DIADEM (DIgital reconstructions of Axonal and DEndrite Morphology) challenge in -2010(Gillette et al., 2011 and the BigNeuron project in 2005 (Peng et al., 2015), a robust algorithm or tool performing single-neuron segmentation from datasets representing dense-packed neurons in their native arrangement within the brain still lacks.
Moreover, most of these approaches just perform three-dimensional neuron tracing, giving scarce information about neuron morphology (i.e. surface and volume of the whole cell, dendrite diameters).

A. L. Callara et al.
The neuron tracing methods available can be categorized, as suggested by (Acciai et al., 2016), into global processing, local processing and metaalgorithm approaches. While the global approaches process the whole datasets, the local processing considers local features around relevant structures and meta-algorithms just enhance particular aspects of existing methods to manage large-scale images.
Due to their ability in managing signal variability, local tracing methods represent the ideal candidate for single-neuron segmentation from confocal or two-photon microscopy datasets. In fact, even after the clarification of the biological samples, deeper layers in the specimen are imaged with lower photon intensity due to light scattering, refraction and absorption, as well as photobleaching and spherical aberration (Diaspro, 2001). All these phenomena, along with the biological fluorescence variability (i.e. the non-uniform distribution of the fluorophores through the sample), lead to both on-plane and intra-plane pixel intensity heterogeneities within the dataset acquired.
The core idea of local methods is to segment objects based on local hard (Li et al., 2008) or soft (Wang et al., 2009) clustering such as local binary fitting (LBF) or gaussian fitting (LGF). This procedure is generally applied without specific spatial constrains, even though some efforts have been done to integrate hard clustering in a region growing approach (RG) (Mostafa et al., 2016). Region growing is a pixel intensity and seed generationbased image segmentation method (Brice and Fennema, 1970). First, a seed point belonging to a region of interest (ROI) is chosen, manually or automatically. Then, the neighboring pixels of the seed are iteratively examined to determine whether they should be added to the ROI or not based on some predefined rule (i.e. a homogeneity predicate). The performance of this procedure may be influenced by the seed selection and by the choice of the rule used by the RG scheme (Baswaraj et al., 2012), which typically depends on the intrinsic properties of the image.
Image properties derive principally from the image formation process. Confocal microscopy is based on a sampling procedure of successive points in a focal plane that creates a spatial representation of the distribution of fluorescent probes within a sample. Hence, each pixel contains a discrete measure of the detected fluorescence within a sample period, represented by a photon count, and certain amount of noise, deriving from different sources (Pawley and Pawley, 2006;Calapez and Rosa, 2010). Given the nature of the measure, statistical methods represent a well suited way for describing the image and to date different models have been proposed to depict the confocal image properties such as the Poisson model (Pawley and Pawley, 2006), the Gaussian model (Calapez et al., 2002) and mixture models (MM, Calapez and Rosa, 2010). A quantitative comparison between these models has been reported elsewhere, defining MM as the best descriptor of the sharp peaks and the long tails characteristic of background and low fluorescence distributions (Calapez and Rosa, 2010). In this light, we developed a new Smart Region Growing (SmRG) algorithm to segment single neurons from 3D confocal image stacks obtained from clarified samples. SmRG is an automatic (semi-automatic) framework that combines a RG procedure and MM describing the confocal microscopy signal statistics. We describe the SmRG algorithm workflow for single-neuron segmentation. To assess the algorithm accuracy, we evaluate the SmRG performance on segmented Purkinje cells (PCs) from 3D image stacks of clarified mouse cerebella acquired using a confocal microscope.
The results are compared with those obtained from manual segmentation (Magliaro et al., 2017).

An outline of SmRG
The SmRG is a region growing algorithm based on the local features of confocal images. Precisely, it exploits the statistics of the background and signal distribution coming from confocal acquisitions and a linear mixture model to determine the probability at which a given pixel (voxel) can be considered as part of the specimen of interest or the background. The homogeneity predicate to grow regions is then designed depending upon these probabilities.
The SmRG algorithm is developed in Matlab (The Mathworks-Inc, USA) and exploits the negative binomial parametrization as in https://github.com/mscipio/MATLABtools___nbin_mu.

The mixture model
In its original version (Calapez and Rosa, 2010), the model is supposed to describe K different fluorescence levels or classes; the k-th-class is described by the linear mixture model where , and αk denote respectively the pixel intensity level, the system offset and the mixture parameter. is the distribution for the background pixels and is modeled according to a discrete normal distribution, with variance and mean , and is the intensity distribution of the k-th class pixels, described by a negative-binomial distribution with parameters = (2) and = (3)

SmRG
For region growing purposes however, it is reasonable to assume the presence of a single class k of pixels, at least locally. In this case, the complete model for a pixel is described by the 5-parameter distribution where all the parameters are real valued except for which is an integer and α ∈ [0,1]. The model fitting is done by means of an Expectation-Maximization (EM) algorithm in which: 1. and are obtained by the method of moments 2. and are given by the maximization of the log-likelihood 3. is given by the posterior density The region growing procedure The SmRG region growing procedure begins by selecting a seed, manually or automatically. In the first case, the user is asked to give the seed position by clicking on the image, while in the latter two different routines are performed. Initially, the algorithm searches for spherical objects within the stack and, if at least one sphere is found, the seed is chosen as the center of the sphere. Otherwise, the MM described in (equation) is fitted globally and seeds are determined as those pixels that certainly belong to the signal ( =0). After seed detection, the homogeneity predicate is derived locally on a 32x32x3 sub-volume (crop) centered on the seed. The smart behavior of SmRG performs at this step: a Hartigan's dip test (Hartigan and Hartigan, 1986) is performed on the crop intensity distribution to test for unimodality against multimodality. In the latter case, the segmentation proceeds with Otsu's method (Otsu, 1979), a well-known thresholding technique for multimodal distributions (Guo et al., 2012). Otherwise, the MM is fitted on the sub-volume pixel intensity distribution and the homogeneity predicate is defined as those pixels with < 0.01. Every pixel that satisfies the homogeneity predicate and is spatially connected to the seed is added to the ROI. New seeds are chosen from the segmented ones. For each segmented plane the regional maxima of the distance transform (Maurer and Raghavan, 2003) are taken as new seeds. The algorithm iterates for each new seed and the process stops when there are no more pixels to add. Moreover, anytime the MM is fitted, the mean absolute error (MAE) between the fitted model and the empirical distribution is calculated. If the MAE is bigger than a settable threshold, the cropped volume in which the fitting fails is saved for visual inspection.
The result of the SmRG is a segmented 3D structure that allows for subsequent morphological operations such as the calculation of neuron surface and volume and dendrite diameters and lengths (Figure 1).

PCs Segmentation evaluation
For evaluating the performance of the algorithm in terms of segmentation precision and accuracy, we used confocal datasets representing neurons in their native arrangement within the brain. The same n=3 PCs from 3D image stacks of murine cerebellum slices clarified as in (Magliaro et al., 2016) and acquired using a confocal microscope (Nikon A1) were segmented using both the SmRG algorithm and ManSegTool (Magliaro et al., 2017), a tool for manual segmentation of 3D stacks (each neuron was segmented manually by 6 experts). To quantify the goodness of neuron segmentation we evaluated the differences between the SmRG segmentation and a manual segmented gold standard by comparing features extracted from both segmentations. As relevant morphological features for our algorithm we considered the surface area, the volume and the Sholl analysis (Sholl, 1953) of segmented structures.
To compare Sholl profiles we calculated the total area under the curve (AUC) using the trapezoidal rule: this allows to have a single measure for each profile (Binley et al., 2014). Statistical differences between the features in the gold standard and those coming from automatic segmentation were evaluated by means of the Friedman's test with replicates. Specifically, surface area, volume and AUC of Sholl profiles represented blocking factors (rows, nuisance effects) with replicates (neurons #1, #2 and #3) and users and SmRG represented treatments (column). Figure 2 shows an example of the same neuron segmented by an expert and by the SmRG algorithm: SmRG can follow the neurite arborization, without missing primary branches, as well as smaller structures; moreover, the structure obtained with SmRG is characterized by a smooth three-dimensional structure that allows the extraction of different morphological information.

Results
The morphological features extracted from the segmented neurons (i.e. surface area, volume and AUC of Sholl profiles) are shown in Figure 3, while a detailed Sholl analysis is reported in Figure 4. The Friedman's test showed no statistically significant differences between automatic (SmRG) and manual (ManSegTool) segmentation in terms of surface area, volume and Sholl profiles of the segmented structure; a detailed ANOVA table of the Friedman's test is reported in Figure 5.

Discussion
The reconstruction of single-neuron morphology from three-dimensional image stacks represents a primary step in neuro-research, aiming at revealing the complex relationships between brain structure and function. Despite remarkable attempts to address this issue, the problem is still far from being solved. On one hand, automatic general-purpose robust methods to deal with the large variability of neuro-image datasets are still lacking; on the other, while considerable efforts have been made in the field of neuron tracing, little attention has been paid to the morphology of segmented structures as if of secondary importance.
SmRG is an open-source Matlab-based algorithm purposely developed to segment complex structures in a three-dimensional environment represented by confocal image stacks, while preserving important and useful information about the morphology of the segmented objects. The algorithm is based on a region growing scheme that allows to follow the signal variability locally and on mixture models aiming at defining a homogeneity predicate for the segmentation. Our results ( Figure 5) suggest that SmRG performance in terms of extracted features from segmented structures is comparable to those obtained from a manual segmentation gold-standard. SmRG was developed and tested on confocal image stacks of Purkinje cells from clarified mouse cerebella. The higher signal-to-noise ratio (SNR) and contrast-to-noise ratio (CNR) of clarified images (Magliaro et al., 2016), along with the results of this work suggest that this procedure should always be performed before segmentation when possible. Unlike other tools and algorithms, the SmRG output allows the extraction of additional features about the segmented neurons (i.e. neurite thickness, neuron surface and volume). Preserving this information may be helpful for studies that rely on neuron morphology, such as simulating the neuron electrophysiological behavior based on empirical microscopic data and cell characterization.
It remains an open question whether a general-purpose algorithm should be possible, considering the great variability of different cell-imaging methods. However, it is a common view that segmentation procedures depend on the goal of the segmentation itself (Zhang, 1996). In this view, SmRG is not supposed to fill this gap, but to represent the last step of an imaging workflow that comprises tissue clarification and confocal imaging.

Conclusions
We integrated a region growing procedure and mixture models to isolate single neurons from confocal image stacks obtained after tissue clarification. To our knowledge, although statistical methods represent a natural candidate for describing confocal image intensity distribution, their integration in spatial constrained segmentation procedures such as region growing are limited. We believe that such integration may help considerably developing new approaches to neuron segmentation. Moreover, SmRG output allows for subsequent 3D analysis of segmented structures in terms of surface area and volume of segmented structures. It is our view that this kind of information, which is often neglected in the state of art, may lead to a more accurate neuron classification.
In (Magliaro et al., 2016) a quantitative analysis of the effects of clarification on the SNR and CNR of confocal images is reported. It may be of particular interest evaluating the effects of clarification on neuron segmentation, since algorithms for this kind of data seems to lack in the state of art. From this point of view, our algorithm can represent a first step in this direction. From a biological perspective, one of the problems that could benefit from digital neuron reconstruction is the characterization of cell type by morphology, including in wild-type and model animals, as well as in healthy and diseased human individuals (Acciai et al, 2016).