Automated Neuron Reconstruction from 3D Fluorescence Microscopy Images Using Sequential Monte Carlo Estimation

Microscopic images of neuronal cells provide essential structural information about the key constituents of the brain and form the basis of many neuroscientific studies. Computational analyses of the morphological properties of the captured neurons require first converting the structural information into digital tree-like reconstructions. Many dedicated computational methods and corresponding software tools have been and are continuously being developed with the aim to automate this step while achieving human-comparable reconstruction accuracy. This pursuit is hampered by the immense diversity and intricacy of neuronal morphologies as well as the often low quality and ambiguity of the images. Here we present a novel method we developed in an effort to improve the robustness of digital reconstruction against these complicating factors. The method is based on probabilistic filtering by sequential Monte Carlo estimation and uses prediction and update models designed specifically for tracing neuronal branches in microscopic image stacks. Moreover, it uses multiple probabilistic traces to arrive at a more robust, ensemble reconstruction. The proposed method was evaluated on fluorescence microscopy image stacks of single neurons and dense neuronal networks with expert manual annotations serving as the gold standard, as well as on synthetic images with known ground truth. The results indicate that our method performs well under varying experimental conditions and compares favorably to state-of-the-art alternative methods. Electronic supplementary material The online version of this article (10.1007/s12021-018-9407-8) contains supplementary material, which is available to authorized users.


Introduction
The brain is regarded as one of the most complex and enigmatic biological structures. Composed of an intricate network of tree-shaped neuronal cells (Ascoli 2015), together forming a powerful information processing unit, it performs a myriad of functions that are essential to living organisms  (Kandel et al. 2012). Obtaining a blue print of the architecture of this network, including the morphologies and interconnectivities of the neurons in various subunits, helps to understand how the brain works (Ascoli 2002;Donohue and Ascoli 2008;Cuntz et al. 2010), including how neurodegenerative disease processes alter its function. A key instrument in this endeavor is microscopic imaging, as it allows detailed visualization of neuronal cells in isolation and in tissue, thus providing the means to study their structural properties quantitatively (Senft 2011).
Quantitative measurement and statistical analysis of neuronal cell and network properties from microscopic data rely on the ability to obtain accurate digital reconstructions of the branching structures (Halavi et al. 2012) in the form of a directional tree of connected nodes (Ascoli et al. 2007). The ever increasing amount of available image data calls for automated computational methods and software tools for this purpose, as manual delineation of neurons is extremely cumbersome even in single image stacks, and is downright infeasible in processing large numbers of images (Svoboda 2011;Senft 2011). Automating neuron reconstruction requires solving fundamental computer vision problems such as detecting and segmenting tree-like image structures (Meijering 2010;Donohue and Ascoli 2011;Acciai et al. 2016). This is complicated by the large diversity of neuron types, imperfections in cell staining, optical distortions, inevitable image noise, and other causes of ambiguity in the image data. Consequently, with the current state-of-theart, manual proof-editing of automatically obtained digital reconstructions is often necessary (Peng et al. 2011b). Recent international initiatives such as the DIADEM challenge (Gillette et al. 2011) and the BigNeuron project (Peng et al. 2015a, b) have catalyzed research in automated neuron reconstruction but have also clearly revealed that further improvement is still very much needed before computers can fully replace manual labor in performing this task.
With this paper we aim to contribute to the developments in the field by proposing a novel fully automated neuron reconstruction method based on probabilistic filtering techniques. Starting from seed points that have a high probability of being centered at neuronal branches, our method recursively traces these branches by sequential Monte Carlo estimation, using state transition and measurement models designed specifically for this purpose. This results in a series of possibly overlapping but probabilistically independent estimates of the branches, which are subsequently combined into a refined estimate of the actual branch centerlines using mean-shifting. We presented early versions of the method at conferences (Radojević et al. 2015;Radojević and Meijering 2017b) and donated one implementation of it (named Advantra) for inclusion in the BigNeuron benchmarking study (Peng et al. 2015a, b). Since then we have improved the method and its software implementation and have significantly extended its experimental evaluation. Here we provide a detailed description of the method, its implementation, and the experimental results, and show that it performs favorably compared to several state-of-the-art neuron reconstruction methods from the BigNeuron project as well as an alternative probabilistic method (Radojević and Meijering 2017a). The source code of our software implementation will be released along with this paper.

Related Work
Early methods and tools for digital neuron reconstruction were semi-automatic and required extensive manual intervention for their initialization and operation or the curation of faulty results (Glaser and Van der Loos 1965;Capowski and Sedivec 1981;Glaser and Glaser 1990;Masseroli et al. 1993). With the increasing capabilities of computers it became possible to store and process 3D images of neurons (Cohen et al. 1994;Belichenko and Dahlström 1995). More recently, the state-of-the-art in the field has moved towards full automation of neuron reconstruction, and various freely available software tools are now available for this purpose (Peng et al. 2010;Longair et al. 2011;Peng et al. 2014a, b), though the need for flexible editing tools has remained unabated (Luisi et al. 2011;Dercksen et al. 2014).
Neuron reconstruction methods typically have a modular design where each module or stage of the processing pipeline deals with different structural objects. Depending on the subproblems being solved, modules can operate independently, or work together for example to combine local and global processing, possibly requiring multiple iterations. Several subproblems that can be identified in the literature include image prefiltering and segmentation (Zhou et al. 2015;Türetken et al. 2011;Sironi et al. 2016;Mukherjee and Acton 2013), soma (cell body) detection and segmentation (Quan et al. 2013), landmark points extraction (Al-Kofahi et al. 2008;Wang et al. 2011;Choromanska et al. 2012;Baboiu and Hamarneh 2012;Su et al. 2012;Radojević et al. 2016), neuron arbor tracing (Zhao et al. 2011;Liu et al. 2016;Leandro et al. 2009;Radojević and Meijering 2017a;Xiao and Peng 2013), and assembling the final tree-like graph structure Türetken et al. 2011;Yuan et al. 2009). In the remainder of this section we briefly review techniques for solving each of these subproblems. Since our primary goal in this paper is to present a new method, the review is not meant to be exhaustive, but to put our method into context.
The pool of neuron reconstruction methods is very diverse (Meijering 2010;Donohue and Ascoli 2011;Acciai et al. 2016;Peng et al. 2015a) but there are also many commonalities. For example, image prefiltering to enhance tubular structures is typically carried out using Hessian or Jacobian based processing (Xiong et al. 2006;Al-Kofahi et al. 2008;Yuan et al. 2009;Wang et al. 2011). And to cope with uneven staining, adaptive thresholding (Zhou et al. 2015), perceptual grouping , and vector field convolution (Mukherjee et al. 2015) have been used. For image segmentation (separating foreground from background), a wide variety of methods has been proposed, including the use of feature-based classifiers (Türetken et al. 2011;Chen et al. 2015;Jiménez et al. 2015), tubularity based supervised regression (Sironi et al. 2016), and even deep learning (Li et al. 2017). The general difficulty of supervised methods, however, is their need for extensive manual annotation for training to arrive at usable segmentation models. In our proposed method we have chosen to avoid this by using carefully designed explicit models.
For the detection and segmentation of the neuronal somas, which typically have a much larger diameter than the dendritic and axonal branches, a simple and efficient solution is to apply morphological closing and adaptive thresholding ). An alternative is to use shape fitting approaches (Quan et al. 2013). Next, to initialize and/or guide the segmentation of the arbor, landmark points are often extracted using image filters that specifically enhance tubular structures Türetken et al. 2011;Choromanska et al. 2012;Su et al. 2012;Radojević et al. 2016), a popular one being the so-called "vesselness filter" (Frangi et al. 1998). In our proposed method we have adopted classical approaches for soma and seed point detection as detailed in the next section.
Segmentation or tracing of all branches of the dendritic and axonal trees is the main challenge of the reconstruction problem. A widely used approach to overcome the difficulties caused by imperfect staining and image noise is to use techniques that find globally optimal paths between seed points by minimizing a predefined cost function (Meijering et al. 2004;Peng et al. 2011a;Longair et al. 2011;Quan et al. 2016). But many other concepts have been proposed as well, including model fitting (Schmitt et al. 2004;Zhao et al. 2011), contour extraction (Leandro et al. 2009), active contour segmentation Luo et al. 2015), level-set or fast-marching approaches (Xiao and Peng 2013;Basu and Racoceanu 2014), path-pruning from oversegmentation (Peng et al. 2011a), distance field tracing , marching rayburst sampling (Ming et al. 2013), marked point processing (Basu et al. 2016), iterative back-tracking , and learning based approaches (Chen et al. 2015;Gala et al. 2014;Santamaría-Pang et al. 2015). In recent works we have shown the great potential of probabilistic approaches to neuron tracing (Radojević et al. 2015;Radojević and Meijering 2017a, b) which formed the basis for the new fully automated neuron reconstruction method presented and evaluated in the next sections.
The final aspect of neuron reconstruction is the assembling of the complete neuronal tree structure from possibly many partial or overlapping traces and putting it into a format that is both representative and suitable for further automated analysis. This is typically solved by graph optimization strategies such as the minimum spanning tree (MST), the alternative K-MST (Türetken et al. 2011;González et al. 2010), or integer programming (Türetken et al. 2013). To deal with very large data sets it has also been proposed to assemble the 3D graph representation through tracing in 2D projections and applying reverse mapping . However, with the advent of sophisticated assemblers such as UltraTracer , it is possible to extend any base tracing algorithm to deal with arbitrarily large volumes of neuronal image data . Therefore, in our proposed method, we do not use projections but perform the tracing in the original image (sub)volumes. And to obtain the graph representation we propose a new approach to refining and grouping the individual traces.

Proposed Method
The pipeline of our proposed method consists of six steps ( Fig. 1) described in detail in the following subsections. We assume that image stacks contain a single neuron (one soma) or just an arbor (no soma) as in the DIADEM  and BigNeuron data (Peng et al. 2015a). In short, we first extract the soma and a set of seeds, which serve to initialize our probabilistic branch tracing scheme. The resulting traces are iteratively refined and their corresponding nodes spatially grouped into a representative node set that is traversed to form the final reconstruction.

Soma Extraction
The soma typically has a considerably larger diameter than the individual branches of the neuronal arbor (Fig. 1a). Thus it can be easily extracted using morphological filtering operations . Specifically, in our method, we apply grayscale erosion to remove all branches and leave only the (eroded) soma. To this end, the radius r s of the structuring element needs to be larger than the largest expected branch radius in a given data set, and smaller than the expected soma radius. The resulting image is then smoothed using a Gaussian filter with standard deviation equal to r s and segmented using max-entropy thresholding (Radojević et al. 2016) to obtain a blob corresponding to the a b c d e f Fig. 1 Schematic overview of the six main steps of the proposed method: a soma extraction, b seed extraction, c branch tracing, d trace refinement, e node grouping, f tree construction soma. For computational efficiency both the erosion and the Gaussian smoothing operation are carried out by separable filtering. In this paper we model the soma in the final graph representation of the neuron as a single spherical node with position equal to the centroid of the segmented blob and radius equal to the average distance of the blob voxels to the centroid. Alternatively, we could model the soma with a set of nodes that together represent the blob as accurately as we like, but in our applications this is not needed.

Seed Extraction
To initialize the branch tracing we extract a set of seed points (Fig. 1b). These seeds are points with very high likelihood of being centered on a branch. In our method we estimate this likelihood using a Hessian-based multiscale tubularity filter (Frangi et al. 1998). 1 It computes for every voxel in an image I the Gaussian-smoothed second-order derivatives: where * denotes convolution, G σ the Gaussian filter at scale σ , and H σ the resulting local Hessian matrix at that scale. The eigenvalues |λ 1 | ≤ |λ 2 | ≤ |λ 3 | of H σ are indicative of the geometry of the local image structure and are used to quantify its tubularity as (Frangi et al. 1998): with the free parameters typically set to α = β = 0.5 and c to half the maximum Hessian norm and where For each voxel location p = [x, y, z] the spatial scale of the local image structure is taken to be the Gaussian σ at which the filter yields the highest tubularity value υ.
And the orientation of the structure is then taken to be the eigenvector v = [v x , v y , v z ] corresponding to the smallest absolute eigenvalue λ 1 of the Hessian matrix H σ . From the resulting tubularity map, initial seed points s i = [p i , v i , σ i ] are selected whose tubularity value is the highest in a cylindrical neighborhood with radius 3σ i and length σ i , centered at p i , and oriented along v i . For this 1 http://bitbucket.org/miroslavradojevic/vess purpose we use a find-maxima function ported from ImageJ, which applies a noise tolerance τ to prune insignificant local maxima (Ferreira and Rasband 2012). The final set of seeds is subsequently obtained by excluding the maxima where the correlation of the image with a cylindrical template model is too low, using exactly the same criterion as for termination of branch tracing, described next.

Branch Tracing
For each seed s i , our method traces the local image structure in two directions, +v i and −v i , producing a pair of local traces (Fig. 1c). A trace is considered to consist of a sequence of hidden states, x 0:L = (x 0 , . . . , x L ), where x 0 is the initial state extrapolated from the seed s i , and x L is the last state of the trace. Similar to the seeds, the states and the scale σ i of the underlying neuron branch. The states are estimated sequentially in a probabilistic fashion using Bayes' rule: where p(x i |z 0:i ) is the posterior probability distribution of the state x i given measurements z 0:i from the first to the current iteration, p(x i |x i−1 ) is the state transition prior, and p(z i |x i ) is the likelihood of measuring z i given state x i . It is assumed that the state transition is a Markovian process and the measurements are independent. To allow for nonlinearities in the process, we solve the estimation problem (4) using sequential Monte Carlo (SMC) filtering (Doucet et al. 2001), also known as particle filtering (Arulampalam et al. 2002). Here the posterior is approximated using a set of N samples x k i with corresponding weights w k i as: where the weights are normalized so that k w k i = 1 and Each iteration in SMC filtering consists of a prediction step and an update step. In the prediction step, given the samples x k i−1 from the previous iteration, N new samples x k i are drawn using the state transition prior. The importance sampling distribution that we use for this is ( Fig. 2a): Functions used in the prediction and update steps of the SMC filtering: a the prediction importance sampling distribution (for ease of visualization a 2D example is given) and b the measurement likelihood function for different values of K. Reprinted with permission from Radojević and Meijering (2017b) x where I 0 denotes the zero-order Bessel function of the first kind, κ is the circular variance parameter, η is a normalization factor that makes the prediction over all N samples integrate to unity, d i = ||p i − p k i−1 || is the Euclidean distance between the predicted position and the sample position in the previous iteration, d is the tracing step size, and ζ the scale variance parameter. Each predicted state is assigned a unit direction v represents the difference in scales, which contributes to the importance sampling function by a Gaussian component, giving a higher value to state samples that retain the scale.
In the update step, the newly drawn samples are updated using the following likelihood function (Fig. 2b): where K determines the sensitivity to the normalized crosscorrelation c x ∈ [−1, 1], which quantifies the similarity of the underlying image structure for x = p, v, σ to a cylindrical template model with Gaussian profile (Fig. 3): where (k, l, m) are the template coordinates, which transform to p in image coordinates since the template is centered at p and is oriented in the direction v and has scale σ of x, and by definition u ⊥ v, w ⊥ v, and u ⊥ w. The summation is limited to −3σ ≤ k, l ≤ 3σ and σ ≤ m ≤ σ which corresponds to the spatial extent of the template.Ī andḠ denote the mean of the image intensities and of the template intensities, respectively, within the mentioned limits. The value of c x is independent of intensity scalings and offsets and thus provides us with a robust measure of structural resemblance, which may range from −1 (inverse correlation), to 0 (no correlation), to +1 (full correlation). The weights of the samples are updated accordingly as: and renormalized so that k w k i = 1. To avoid weight deterioration, systematic resampling (Kitagawa 1996) is performed each time the effective sample size N eff (Kong et al. 1994) falls below 80% of N. The final state estimate after each iteration i, which constitutes a node of the trace, is computed from the weighted samples as the centroid: Filtering is terminated if the average correlation value k c x k i /N drops below the threshold c min , indicating the end of the underlying neuron branch in the image, or if the iteration limit L is reached. Since the filtering is done for each seed, and in both (opposite) directions, the same neuron branch may be traced many times over, but in a probabilistically independent way, providing accumulating evidence about the presence and location of the branches. However, to avoid excessive over-tracing and to reduce the computation time, we also monitor the node density D n per image volume unit of n voxels and terminate the tracing if the density in the current position exceeds the limit δ n . In principle n can be any number, but in our work we typically consider n = 1 (single voxel), n = 5 (4-connected in 3D), and n = 9 (8-connected in 3D), which is sufficient given that the image stacks often have a large voxel anisotropy.

Trace Refinement
After the tracing step, each neuron branch may have multiple corresponding traces, and each trace node has bidirectional links to neighboring nodes (Figs. 1d and 4a) to allow trace traversal in any of the possible directions in but in the sequel we write the elements of N more generally as n k , k = 1, . . . , M, where M = T t=1 M t . Each node n k contains an estimate of the center position p = (x, y, z) and the cross-sectional radius (r) of the underlying branch structure, as well as the cross-correlation (c) with the cylindrical Gaussian template model, and a set (I) containing the indices in N of the neighboring nodes: where I k has either two elements (in the case of a body node) or just one (in the case of a terminal node).
The goal of the trace refinement step is to exploit the cumulative evidence provided by the over-tracing in the previous step to improve the individual node estimates. Specifically, we update each node n k to: by applying mean-shifting (Cheng 1995), resulting in an updated node setN . Mean-shifting iteratively moves each node element to the local mean of the nodes in its vicinity: This reduces the variance of the estimates but preserves the linking of the nodes:Ī = I. In practice, five iterations are sufficient to reach satisfactory radial trace alignment  (Fig. 4b). The kernel size used in the mean-shifting process is taken to be the initial radius of each node. In our implementation, prior to mean-shifting, we resample all traces with a step size of one voxel to get a more finegrained result.

Node Grouping
Although the previous step results in refined node estimates, it keeps the total number of nodes and corresponding multiple traces. The next step is to merge overlapping traces and obtain a single trace for each neuron branch. In our method this is accomplished by the process of node grouping (Figs. 1e and 4c) as detailed in Algorithm 1. It iteratively takes from the refined setN an as-yet ungrouped node with the highest cross-correlation value, finds all its neighboring nodes within the predefined Euclidean distance r g , and groups them by calculating the mean value of each element. In our implementation we use unweighted averaging for this. Alternatively, weighted averaging could be used, based on the cross-correlation scores. The node links within a group are accumulated and their indexes mapped to the group node index list. This results in a new setN = {n 1 , . . . ,n P }, P ≤ M, of group nodes: and any twon i andn j are connected if there exists a link between any of the refined nodes captured by these two, as revealed by the accumulated index setsÎ i andÎ j . Thus, all existing inter-node connectionsĪ are preserved, and are projected into the inter-group connectionsÎ (see Supplementary Fig. S13 for an example case).

Tree Construction
The final step of our method is the construction of a graph representing the complete neuronal arbor. This is facilitated by the bidirectional connectivity of the group nodes inN . However, similar to a real neuron, the final graph must be a tree, in which the nodes are unidirectionally linked (Figs. 1f and 4d), as also required by the SWC file format for storing digital neuron reconstructions (Stockley et al. 1993;Cannon et al. 1998). Starting from the soma node, or from the group node with the highest cross-correlation value if no soma was found in the image, the nodes inN are iteratively traversed using a breadth-first search (BFS) algorithm. In this process it is possible to discard any isolated branches and single-node terminal branches (false positives).

Implementation Details
Our method, which we call Probabilistic Neuron Reconstructor (PNR), was implemented in C++ as a plugin for the freely available and extendable bioimage visualization and analysis tool Vaa3D (Peng et al. 2010(Peng et al. , 2014a. 2 The source code of PNR is freely available for non-commercial use. 3 As mentioned in the preceding sections, the method has a number of free parameters, which are summarized in Table 1, where we also list default values.  Iteration limit δ n δ 9 = 4 δ 9 = 4 δ 9 = 4 δ 1 = 3, δ 9 = 4 δ 1 = 3, δ 9 = 4 Node density limit [count/volume] r g 2 2 2 2 2 G r ouping radius [voxels] The ordering is according to first mention in the main text

Experimental Results
The performance of our PNR method was evaluated using both synthetic and real fluorescence microscopy image stacks of single neurons and was compared to several alternative 3D neuron reconstruction methods that yielded favorable performance in the BigNeuron project (Peng et al. 2015a). These included the second all-path pruning method (APP2) (Xiao and Peng 2013), NeuroGPS-Tree (GPS) (Quan et al. 2016), BigNeuron's minimum spanning tree (MST) method, and we also added our recently published alternative probabilistic method based on probability hypothesis density filtering (PHD) (Radojević and Meijering 2017a).
To quantify performance we adopted the commonly used measures of distance and overlap of neuron reconstructions with respect to the ground truth (in the case of synthetic images) or the gold-standard reconstructions obtained by manual annotation (in the case of real images). The distance measures were the average minimal reciprocal spatial distance (SD) between nodes in the reconstructions being compared, the substantial spatial distance (SSD) using only the nodes with a spatial distance larger than a threshold S, and the percentage of these substantially distant nodes (%SSD), all computed after densely resampling each reconstruction to reduce the distance between its adjacent nodes to one voxel (see Peng et al. 2010 for details). The overlap measures were precision (P), recall (R), and the F score (Powers 2011), computed from the numbers of truepositive (TP), false-positive (FP) and false-negative (FN) nodes according to the spatial distance threshold S.
All experiments were performed on a MacBook Pro with 2.2 GHz Intel Core i7 processor and 16 GB RAM memory to test the practicality of the methods on a typical computer system. For each method we optimized the score for each performance measure by exploring a grid of possible parameter values around the default ones (see Table 1 for our method and the cited papers for the other methods). To keep the experiments feasible, we set the maximum allowed processing time per stack and method to 2 hours. In the sequel, to save space, we show only the F scores (higher is better) and SSD scores (lower is better), while the P, R, SD, and %SSD scores are given in the supplement. Our conclusions are based on the complete body of results.

Experiments on Synthetic Neuron Images
Prior to evaluating how well our method emulates expert manual reconstruction in real neuron images, we first performed a controlled experiment using synthetic neuron images, with known ground-truth reconstructions and predefined levels of signal-to-noise ratio (SNR) and inter-voxel correlation (COR). This allowed us to study the robustness of our method compared to the others as a function of these image quality factors. For this experiment we selected 10 neurons from the BigNeuron training data set (Peng et al. 2015a), representative of the range of morphological complexities in the data set, and for which node radius information (non-default) was available in the corresponding gold-standard reconstructions in SWC format.
We developed a plugin for ImageJ (Schneider et al. 2012) called SWC2IMG, 4 which takes any SWC file as input and simulates fluorescence microscopy imaging of all neuronal branches in the file at a specified SNR and COR level, producing an image stack whose true digital reconstruction is the very input. It assumes that in practice, because of the relatively large spatial extent of even a single neuron with its complete arbor, the combination of optical magnification factor and digital image matrix size in real neuron images is typically such that the voxel size is larger than the point spread function (PSF), implying that the partial-volume effect of digitization is more prominent than the optical blurring by the microscope. Based on this, the plugin simulates the imaging simply by estimating for each voxel which fraction of its volume is occupied by the neuron. Next, it simulates noise by using the Poisson noise model representative of optical imaging, which defines SNR as the image intensity inside the neuron above the background, divided by the standard deviation of the noise inside (Sheppard et al. 2006). And finally, to allow for correlated signal and noise, which we found to improve the visual realism of the simulated images, the plugin also offers the possibility to apply Gaussian smoothing at a specified scale, being the COR parameter, while preserving the SNR level. Generally, the lower the SNR and/or the higher the COR level, the more challenging the data and the reconstruction problem.
Using this plugin we created a synthetic data set containing image stacks for a range of SNR and COR values for each neuron (Fig. 5). Specifically, we considered SNR = 1, 2, 3, 4, 5, 10, 20, and COR = 0, 0.5, 1, 1.5, 2. Thus, our synthetic data set consisted of 10 (neurons) × 7 (SNR levels) × 5 (COR levels) = 350 image stacks, 5 which we attempted to reconstruct optimally using the five considered methods (APP2, GPS, MST, PHD, PNR) and a parameter grid-search approach. However, some of the images were very challenging, especially the ones with many branches and low SNR or high COR values, causing the methods to sometimes 5 https://bitbucket.org/miroslavradojevic/pnr require excessive computation times or even to get stuck altogether. Because of the mentioned time constraint, not all methods were able to complete all the reconstructions, and it turned out that only 7 out of the 10 neurons could be reconstructed by all the methods for all SNR and COR values. Therefore we present the results only for those.
From the average F and SSD scores of the methods as a function of SNR for a few sample values of COR and S (Figs. 6 and 7) we generally observe that, as expected, increasing the SNR improves the performance of all methods (increasing F and decreasing SSD scores). We also note that the two probabilistic methods (PHD and PNR) are more robust against noise (especially according to F) and that our proposed method (PNR) is often superior overall. The results also show that, as expected, increasing the value of COR (which yields more difficult images) has a strong negative impact on the performance of all methods (lower F and higher SSD scores for the same SNR). This is confirmed when looking more in-depth at the results as a function of COR (Figs. 8 and 9). Additionally, again as expected, in all cases we observe that increasing the value of S (meaning being more lenient in matching reconstructions to the ground truth) may also strongly affect the scores of all methods (meaning higher F scores, but in this case also higher SSD scores, as the latter includes only node distances larger than S). This is confirmed when looking explicitly at the performance of the methods as a function of S (Figs. 10 and 11). These results reveal that both the absolute and the relative performance of different methods being compared may depend on S. This is an important observation, since in all studies we are aware of, the somewhat arbitrary value of S = 2 is taken for granted in calculating performance and ranking the methods. Our results (Figs. 10 and 11) show that taking other values of S may yield a different ranking. Notwithstanding this finding, our results also show that under most experimental conditions (SNR, COR, S), the proposed method (PNR) repeatedly and in a statistically independently way, indeed yields more evidence about the underlying neuron branches and leads to better reconstructions. This also follows from a visual comparison of the reconstructions (Fig. 12).
Especially at low SNRs, pruning and fast-marching based methods tend to oversegment the images, while our probabilistic methods still perform relatively well regardless. Even at high SNRs, when most of the methods

Experiments on Real Neuron Images
In addition to synthetic data we also used three real neuron image data sets to evaluate the absolute and relative performance of our method. The first two are the olfactory projection fibers (OPF) data set (9 image stacks) and neocortical layer-1 axons (NCL1A) data set (16 image stacks) from the DIADEM challenge , and the third is part of the BigNeuron (BGN) training data set (76 image stacks) (Peng et al. 2015a), all imaged with fluorescence microscopy (confocal or two-photon) and manually annotated as described in detail in the cited works and corresponding resources. Being the smallest of the three, in terms of both neuronal volume and complexity, OPF is probably the most often used data set in the field. NCL1A is often used as it contains neuronal network-like structures and no clear somas. And BGN is the largest, most diverse, and thus most challenging data set for evaluating neuron reconstruction methods. Together, the 100+ image stacks in these data sets have a wide variety of image qualities and volumes (10 MB to 2 GB per stack) and portray a wide range of neuronal shapes and complexities  Fig. 17 Example neuron reconstructions of an image stack from the OPF data set. Shown are the original arbor (volume rendering on the left) and the reconstructions (overlaid surface renderings in red) of the different methods (indicated at the top) corresponding to the best F score (given below each reconstruction) for S = 2 with respect to the available manual reconstruction (Fig. 13), representative of many studies. For some stacks in the BGN data set, the voxel size was unknown, and in these cases we used a default x:y:z voxel aspect ratio of 1:1:2, reflecting the typically lower resolution in the depth dimension. Also, because of the mentioned processing time constraint, 3 of the 76 image stacks could not be reconstructed by all methods (see Supplementary Fig. S14 for these and other hard cases), so the presented results are based on the remaining 73.
From the results of the experiments on these three real data sets (Figs. 14, 15 and 16) we observe that, as in the experiments on synthetic data, the probabilistic methods PHD and PNR typically show superior performance in terms of both F and SSD score. Of these two methods, our proposed PNR method consistently shows the smallest performance spread, indicating it is more robust than our previously published PHD method. For the BGN data set, being the most diverse of the three, the performance spread (including outliers) of all methods is the largest, and the increase in performance as a function of S is the smallest, as expected. Nevertheless, the PNR method consistently shows the best overall performance especially for this data set. In other words, for any given data set similar to those considered in this study, PNR is the favorable method a priori. Obviously this does not necessarily mean that PNR will give the best reconstruction for each and every image stack in the data set, but simply that the chances are higher. This is confirmed when we look at a few example image stacks from the three data sets and the corresponding best reconstructions produced by the different methods by maximizing the F score in the parameter grid search (Figs. 17,18 and 19). As these examples show, although PNR often outperforms the other methods, in specific cases one of the other methods may give better results.
Finally we investigated the sensitivity of PNR with respect to two of its parameters that one might suspect to be rather critical. The first is the noise tolerance parameter τ used to prune insignificant local maxima in the seed extraction (Seed Extraction). To obtain the best possible results while keeping the computation times manageable, we considered different sets of values for this parameter, depending on the data set (Table 1). For example, in the case of the relatively small-sized OPF data set we considered values τ = 4, 6, 8, 10, 12, while for the larger NCL1A and BGN data sets and the synthetic images we examined τ = 6, 8, 10. The results (Supplementary Figs. S15-S18) show that τ is in fact not a very sensitive parameter of the method and that the suggested default value is a suitable choice. The second parameter is the node density limit δ n used to terminate the tracing (Branch Tracing). Here, too, we considered different sets of values depending on the data set, for n = 5 and n = 9. The results ( Supplementary Fig. S19) show that the method is also not very sensitive to this parameter and its default value is suitable. Notice that due to the Fig. 18 Example neuron reconstructions of an image stack from the NCL1A data set. Shown are the original arbor (volume rendering on the left) and the reconstructions (overlaid surface renderings in red) of the different methods (indicated at the top) corresponding to the best F score (given below each reconstruction) for S = 2 with respect to the available manual reconstruction Fig. 19 Example neuron reconstructions of four image stacks from the BGN data set. Shown are the original arbors (volume renderings on the left) and the reconstructions (overlaid surface renderings in red) of the different methods (indicated at the top) corresponding to the best F score (given below each reconstruction) for S = 2 with respect to the available manual reconstruction probabilistic nature of the method there is inherently some randomness in the results. But altogether we believe the results justify the conclusion that PNR is a robust method and a valuable addition to the neuron reconstruction toolbox.

Conclusions
We have presented a new fully automated probabilistic neuron reconstruction method (PNR) based on sequential Monte Carlo filtering. It traces individual neuron branches from automatically detected seed points repeatedly but statistically independently to acquire more evidence and to be more robust to noise and other artifacts. The traces are subsequently refined, merged, and put into a tree representation for further analysis. We evaluated the method on both synthetic and real neuron images and compared it to various other state-of-the-art neuron reconstruction methods (APP2, GPS, MST, PHD) using commonly used quantitative performance measures (we presented F and SSD scores). To obtain realistic synthetic data we developed a novel simulator (SWC2IMG) that can turn any given SWC file into an image stack of specified quality whose ground truth reconstruction is the input. For the evaluation on real data we used about 100 single-neuron fluorescence microscopy image stacks of widely varying quality and complexity, with corresponding manual reconstructions serving as the gold standard, from three different data sets used in the DIADEM and BigNeuron studies. The results show conclusively that the proposed method is generally favorable and also outperforms our own alternative neuron reconstruction method based on probability hypothesis density (PHD) filtering we presented recently. Nevertheless, there still remains much room for further improvement, as none of the quantitative scores were near perfect for any of the considered methods even for high SNR levels and very lenient distance thresholds. Possible directions for future work within the presented probabilistic framework would be to explore other state transition and measurement models. Alternatively, since no single method always performs best on all images of a given data set, and the results of different methods are likely complementary, another possible direction could be to combine multiple methods either during tracing or in a post-processing step. The latter approach is already being explored in the BigNeuron project. But regardless of the outcome of this effort we conclude that the method proposed in this paper may already prove to be of great use in many cases. Our software implementation of the method will be made freely available for non-commercial use upon publication.

Information Sharing Statement
The source code of the presented method is freely available for non-commercial use from https://bitbucket.org/miroslav radojevic/pnr.