# Hidden Markov Modeling for Maximum Probability Neuron Reconstruction

Thomas L. Athey  Center for Imaging Science & Kavli Neuroscience Discovery
Institute
& Department of Biomedical Engineering
The Johns Hopkins University<EMAIL_ADDRESS>,  Daniel J. Tward  Department
of Computational Medicine & Department of Neurology
University of California at Los Angeles<EMAIL_ADDRESS>,  Ulrich
Mueller  Department of Neuroscience
The Johns Hopkins University<EMAIL_ADDRESS>,  Joshua T. Vogelstein
Center for Imaging Science & Kavli Neuroscience Discovery Institute
& Department of Biomedical Engineering
The Johns Hopkins University
301 Clark Hall
Baltimore, MD 21218<EMAIL_ADDRESS>and  Michael I. Miller  Center for Imaging
Science & Kavli Neuroscience Discovery Institute
& Department of Biomedical Engineering
The Johns Hopkins University
301 Clark Hall
Baltimore, MD 21218<EMAIL_ADDRESS>
###### Abstract.

Recent advances in brain clearing and imaging have made it possible to image
entire mammalian brains at sub-micron resolution. These images offer the
potential to assemble brain-wide atlases of neuron morphology, but manual
neuron reconstruction remains a bottleneck. Several automatic reconstruction
algorithms exist, but most focus on single neuron images. In this paper, we
present a probabilistic reconstruction method, ViterBrain, which combines a
hidden Markov state process that encodes neuron geometry with a random field
appearance model of neuron fluorescence. Our method utilizes dynamic
programming to compute the global maximizers of what we call the “most
probable” neuron path. Our most probable estimation method models the task of
reconstructing neuronal processes in the presence of other neurons, and thus
is applicable in images with several neurons. Our method operates on image
segmentations in order to leverage cutting edge computer vision technology. We
applied our algorithm to imperfect image segmentations where false negatives
severed neuronal processes, and showed that it can follow axons in the
presence of noise or nearby neurons. Additionally, it creates a framework
where users can intervene to, for example, fit start and endpoints. The code
used in this work is available in our open-source Python package brainlit.

##  1\. Introduction

Neuron morphology has been a central topic in neuroscience for over a century,
as it is the substrate for neural connectivity, and serves as a useful basis
for neuron classification. Technological advances in brain clearing and
imaging have allowed scientists to probe neurons that extend throughout the
brain, and branch hundreds of times (Winnubst et al., 2019). It is becoming
feasible to assemble a brainwide atlas of cell types in the mammalian brain
which would serve as a foundation to understanding how the brain operates as
an integrated circuit, or how it fails in neurological disease. One of the
main bottlenecks in assembling such an atlas is the manual labor involved in
neuron reconstruction.

In an effort to accelerate reconstruction, many automated reconstruction
algorithms have been proposed, especially over the last decade. In 2010, the
DIADEM project brought multiple institutions together to consolidate existing
algorithms, and stimulate further progress by generating open access image
datasets, and organizing a contest for reconstruction algorithms (Peng et al.,
2015). Several years later, the BigNeuron project continued the legacy of
DIADEM, this time establishing a common software platform, Vaa3D, on which
many of the state of the art algorithms were implemented (Peng et al., 2010b).
Acciai et al. (2016) offers a review of notable reconstruction algorithms up
to, and through, the BigNeuron project.

Previous approaches to automated neuron reconstruction have used shortest
path/geodesic computation (Peng et al., 2010a; Wang et al., 2011; Turetken et
al., 2013), minimum spanning trees (Yang et al., 2019), Bayesian estimation
(Radojević and Meijering, 2017), tracking (Choromanska et al., 2012; Dai et
al., 2019), and deep learning (Friedmann et al., 2020; Zhou et al., 2018; Li
et al., 2017). Methods have also been developed to enhance, or extend existing
reconstruction algorithms (Wang et al., 2019, 2017; Peng et al., 2017; Chen et
al., 2015). Also, some works focus on the subproblem of resolving different
neuronal processes that pass by closely to each other (Xiao and Peng, 2013; Li
et al., 2019; Quan et al., 2016).

Both DIADEM and BigNeuron initiatives focused on the task of single neuron
reconstruction so most associated algorithms fail when applied to images with
several neurons. However, robust reconstruction in multiple neuron images is
essential in order to assemble brainwide atlases of neuron morphology.

We propose a probabilistic model-based algorithm, ViterBrain, that operates on
imperfect image segmentations to efficiently reconstruct neuronal processes.
Our estimation method does not assume that the image outside the
reconstruction is background, and thus allows for the existence of other
neurons. Our approach draws upon two major subfields in Computer Vision,
appearance modeling and hidden Markov models, and generates globally optimal
solutions using dynamic programming. The states of our model are locally
connected segments. We score the state transitions using appearance models
such as exhibited by Kass and Cohen’s early works (Kass et al., 1988; Cohen,
1991) on active shape modeling and their subsequent application by Wang et al.
(Wang et al., 2011) for neuron reconstruction. For our own approach we exploit
foreground-background models of image intensity for the data likelihood term
in the hidden-Markov structure. We quantify the image data using KL-divergence
and intensity autocorrelation in order to validate our model assumptions.

Our probabilistic models are hidden Markov random fields, but we reduce the
computational structure to a hidden Markov model (HMM) since the latent axonal
structures have an absolute ordering. Hidden Markov modeling (HMM) involves
two sequences of variables, one is observed and one is hidden. A popular
application of HMM’s is in speech recognition where the observed sequence is
an audio signal, and the hidden variables a sequence of words (Rabiner and
Juang, 1986). In our setting, the observed data is the image, and the hidden
data is the contour representation of the axon or dendrite’s path.

The key advantages to using HMMs in this context are, first, that neuronal
geometry can be explicitly encoded in the state transition distribution. We
utilize the Frenet representation of curvature in our transition distribution,
which we have studied previously in Athey et al. (2021); Khaneja et al.
(1998). Secondly, globally optimal estimates can be computed efficiently using
dynamic programming in HMMs (Forney, 1973; Khaneja et al., 1998). The well-
known Viterbi algorithm computes the MAP estimate of the hidden sequence in an
HMM. Our approach, inspired by the Viterbi algorithm, also computes globally
optimal estimates. Thus, our optimization method is not susceptible to local
optima that exist in filtering methods, or gradient methods in active shape
modeling.

Figure 1. Left panel shows an image subvolume from the MouseLight project
containing a single neuron. The middle panel shows the same image overlaid
with a binary image mask in brown. This mask was generated by the random
forest based software Ilastik (Berg et al., 2019) and illustrates the typical
output of an image segmentation model. The right panels shows the same binary
image mask, with a different color for each connected component. The variety
of colors shows that the neuron has been severed into several pieces. All
panels are maximum intensity projections (MIPs).

In this work, we apply our hidden Markov modeling framework to the output of
low-level image segmentation models. Convolutional neural networks have shown
impressive results in image segmentation (Ronneberger et al., 2015), but it
only takes a few false negatives to sever neuronal processes that are often as
thin as one micron (Figure 1). Our method strings together the locally
connected components of the binary image masks into a reconstruction with a
global ordering. Thus, our method is modular enough to leverage state of the
art methods in machine learning for image segmentation.

We apply our method to data from the MouseLight project at Janelia Research
Campus (Winnubst et al., 2019), and focus on the endpoint control problem in a
single neuronal process i.e. start and end points are fixed. We introduce the
use of Frechet distance to quantify the precision of reconstructions and show
that our method has comparable precision to state of the art, when the
algorithms are successful.

##  2\. Results

###  2.1. Overview of ViterBrain

Viterbrain takes in an image, and associated neuron mask produced by some
image segmentation model. The algorithm starts by processing the mask into
neuron fragments and estimating fragment endpoints and orientation to generate
the states. Next, both the prior and likelihood terms of the transition
probabilities are computed to construct a directed graph reminiscent of the
the trellis from the HMM work in Forney (1973). Lastly, the shortest path
between states is computed with dynamic programming. We have proven that the
shortest path is the “most probable” estimate state sequence formed by a
neuronal process (Statement 1). An overview of our algorithm shown in Figure 2
and details are given in the Methods section. We validated our algorithm on
subvolumes of one of the MouseLight whole-brain images (Winnubst et al.,
2019). The image was acquired via serial two-photon tomography at a resolution
of $0.3\mu m\times 0.3\mu m\times 1\upmu$m per voxel. Viterbrain is available
in our open source Python package brainlit: http://brainlit.neurodata.io/.

Figure 2. Summary of our algorithm. The algorithm takes in an image and a
binary mask that might have severed, or fused neuronal processes. First, the
mask is processed into a set of fragments. For each fragment, the endpoints
and endpoint orientations are estimated and added to the state space. Next,
transition probabilities are computed from both the image and state data to
generate a directed graph reminiscent of the trellis graph in classic hidden
Markov modeling. Finally, a shortest path algorithm is applied to compute the
maximally probable state sequence connecting the start to the end state.

###  2.2. Modeling Image Intensity

Figure 3. Characterization of voxel intensity distributions in three different
subvolumes of one of the Mouselight whole-brain images. a) Correlation of
intensities between voxels at varying distances from each other. The curves
show that intensities are only weakly correlated ($\rho<0.4$) at a distance of
$>10$ microns for foreground voxels, or $>2$ for background voxels. Error bars
represent a single standard deviation of the Fisher z-transformation of the
correlation coefficient. Each curve was generated from all pairs of 5000
randomly sampled voxels. b) Kernel density estimates (KDEs) of foreground and
background intensity distributions. A subset of the voxels in each subvolume
was manually labeled, then used to train an Ilastik model to classify the
remaining voxels. Each KDE was generated from $5000$ voxels, according to the
Ilastik classifications. KDEs were computed using scipy’s Gaussian KDE
function with default parameters (Virtanen et al., 2020).

Figure 3a shows the correlations of image intensities between voxels at
varying distances of separation. As is typical for natural images, voxels that
are close by each other have positively correlated intensities, and those
farther away are uncorrelated. In the case of foreground voxels, correlations
become weak beyond a distance of about $10$ microns, with background voxel
correlation decaying rapidly. This lends support to our assumption that voxel
intensities are conditionally independent processes, conditioned on the
foreground/background model (Eq. 3b). This assumption is one of the central
features of our model because it provides for computational tractability.

Figure 3b shows kernel density estimates (KDEs) of the foreground and
background image intensity distributions. The distributions vary greatly
between the three image subvolumes, implying that modeling the image process
as homogeneous throughout the whole brain would be inappropriate.
Additionally, the distributions do not appear to be either Gaussian or
Poisson. Indeed, Kolmogorov-Smirnov tests rejected the null hypothesis for
both Gaussian and Poisson goodness of fit in all cases, with all p-values
below $10^{-16}$. For that reason, we exploit the independent increments
properties of Poisson emission conditioned on the underlying intensity model,
but do not assume that the marginal probabilities are Poisson (or Gaussian),
instead, we estimate the intensity distributions from the data itself using
KDEs (denoted $\alpha_{0}(\cdot),\alpha_{1}(\cdot)$ in Section 4.1).

###  2.3. Maximally Probable Axon Reconstructions

Figure 4 demonstrates the reconstruction method on both a satellite image of
part of the Great Wall of China, and part of an axon. Different image
segmentation models were used to generate the fragments in the two cases, but
the process of joining fragments into a reconstruction was the same. The
algorithm reconstructs the one dimensional structure in both cases.

Figure 4. Demonstration of maximally probable reconstruction on isolated
linear structures from a) a satellite image of part of the Great Wall of China
and b) a neuronal process from the MouseLight dataset (MIP). Left panels show
the original images. Middle panels shows the space of fragments,
$\mathcal{F}$, pictured in color. The green and red arrows indicate the start
and end states of the reconstruction task, respectively. The right panels show
the most probable fragment sequences, where the fragments are colored and
overlaid with a blue line connecting the endpoints of the fragments.

Since the state transition probabilities are modeled with a Gibbs distribution
indexed by the state (hence giving the Hidden Markov property),
reconstructions are driven by relative energies of different transitions.
Thus, our method can still be successful in the presence of luminance or
fragment dropout, as long as the neuronal process is relatively isolated
(Figure 5). Our geometric prior has two hyperparameters, $\alpha_{d}$ and
$\alpha_{\kappa}$, which determine the influence of distance and curvature,
respectively, on the probability of connection between two neuronal fragments.

Figure 5. ViterBrain is robust to image intensity and fragment dropout when
axons are relatively isolated. In this image, we censored the image intensity
and the fragment space periodically along an axon path (red arrows in panels c
and d show censored sections). Nonetheless, our algorithm was able to jump
over the censored regions to reconstruct this axon (panel e). All images are
MIPs

Figure 6a shows various examples of maximally probable reconstructions. The
algorithm was run with the same hyperparameters in all cases: $\alpha_{d}=10$
and $\alpha_{\kappa}=1000$.

Figure 6. a) Successful axon reconstructions; the ViterBrain reconstructions
are shown by the blue line; the manual reconstructions are shown by the red
line. The algorithm was run with the same hyperparameters in each case:
$\alpha_{d}=10$ and $\alpha_{\kappa}=1000$. b) Different hyperparameter values
lead to different results. Panel i) shows the neuron of interest. Panels ii-
iv) are close-up views of reconstructions with different hyperparameter values
that weigh transition distance ($\alpha_{d}$) and transition curvature
($\alpha_{\kappa}$). The red circle in Panel ii) indicates where the
reconstruction deviated from the true path by jumping $\sim 10$ $\upmu$m to
connect the gray fragment to the light blue fragment. Panel iii) shows how a
higher $\alpha_{d}$ value avoids the jump in panel ii), but takes a sharp turn
to deviate from the true path (red circle). Finally, in panel iv), the
reconstruction avoids both the jump from panel ii) and the sharp turn from
panel iii) and follows the true path of the axon back to the cell body. All
images are MIPs.

In some cases, reconstruction accuracy is sensitive to hyperparameter values
(Figure 6b). Higher values of $\alpha_{d}$ penalize transitions between
fragment states with large gaps; higher values of $\alpha_{\kappa}$ penalize
state transitions with sharp angles as measured by their discrete curvature.
It is important to note that the transition distributions depend on the exact
values of $\alpha_{d}$ and $\alpha_{\kappa}$, not just, for example, the ratio
between them. We found $\alpha_{d}=10$ and $\alpha_{\kappa}=1000$ to be
effective for the reconstructions shown throughout the paper, but these values
should be adjusted according to the quality of fragments, and the geometry of
the neurons being reconstructed.

###  2.4. Signal Dropout Causes Censored Fragment States

Though our method is robust to fragment dropout in certain contexts, it is
less robust when there is dropout near a parallel neuronal process (Figure
S2). Since the reconstruction is ultimately a sequence of fragments, missing
fragments forces the algorithm to consider alternative paths thereby trading
off the curvature of jumping to an adjacent state (nearby alternative) or
continuing on course on what might appear to be the correct neuronal
trajectory.

The most obvious source of fragment dropout is low image luminance, leading to
false negatives in the initial image segmentation. However, the reason for
fragment dropout depends on the underlying segmentation model. Empirically, we
found that the reconstruction algorithm fails when when the fragment
generation process neglects portions that are greater than $\sim 10\>\upmu$m
in length.

###  2.5. Comparison to State of the Art

We examined the accuracy of ViterBrain compared to state of the art
reconstruction algorithms. We identified four algorithms that have
accompanying publications, and freely available open-source implementations.
The first method is APP2 which starts with an oversegmentation of the neuron
using a shortest path algorithm, then prunes spurious connections (Xiao and
Peng, 2013). The second is Snake which is based on active contour modeling
(Wang et al., 2011). The third is called Advantra, based on the particle
filtering approach by Radojevic and Meijering (2019). The final reconstruction
software we use is GTree (Zhou et al., 2021), which uses the algorithm
outlined in Quan et al. (2016). This algorithm is similar to APP2 in that it
starts with an initial reconstruction that spans several neurons, then
identifies false connections. APP2, Snake, and Advantra were both used with
their default settings and hyperparameters in Vaa3D 3.2 for Mac. GTree version
1.0.4 was used on Linux. For GTree, the binarization threshold was set to
$1.0$, which was qualitatively identified as a good threshold to capture the
neuron. The default soma radius, $3\;\upmu$m, was used in the soma
identification step.

Shown in Figure 7a are the results of the various methods on a dataset of 35
subvolumes of a MouseLight whole brain image. Each subvolume contains a cell
body, and the initial part of its axon that is covered by the first ten points
of the Janelia reconstruction. So, the subvolumes vary in size but usually
encompass around $10^{6}$ cubic microns. The algorithms are evaluated on how
well they can trace the axon between fixed endpoints (cell body and tenth axon
reconstruction point). The outcomes were classified as either successful (if
the axon was fully traced), partially successful (if more than half of the
axon was reconstructed as evaluated visually), or failures. According to two
proportion z-tests the success rate of ViterBrain ($11/35$) was higher than
all other methods at $\alpha=0.05$. Also, APP2 had a higher success rate
($4/35$) than Advantra at $\alpha=0.05$. The success rates for several of the
algorithms are discouragingly low, so we discuss the possible reasons for this
in the Discussion, and demonstrate the typical failure modes in Figures S3, S4
and S5.

For the successful reconstructions, we examined the precision of the
reconstructions using spatial distance (SD) and Frechet distance. We manually
upsampled the Janelia reconstructions to provide the most precise ground
truth. All reconstructions were sampled at every $1\;\upmu$m before distances
were measured. Spatial distance represents the average distance from a point
on one reconstruction too the closest point on the other reconstruction (Peng
et al., 2010b). Frechet distance represents the maximum distance between two
reconstructions, and satisfies the criteria to be a mathematical metric and is
invariant to reparameterization (see section 4.5 for a derivation of Frechet
distance). The spatial distance between ViterBrain and manual reconstructions
are typically between 1 and 3 microns which is roughly the resolution of the
image. Successful reconstructions by the other algorithms are also in this
range. Frechet distances are larger, since they indicate maximal distance, not
average distance, between curves. We observe that the $\sim 5$ micron
deviations occur most often near the axon hillock where the axon broadens to
merge with the soma.

Figure 7. We applied several reconstruction algorithms to a dataset of 35
subvolumes of a MouseLight whole brain image (Snake was only applied to 10
subvolumes due to incoherent results and excessively slow runtimes, see Figure
S4). Each subvolume contained a soma and part of its axon. The task was to
reconstruct the portion of the axon that was contained in the image (no
branching). First, the algorithms were evaluated visually and classified as
successful, partially successful (over half, but not all, of the axon
reconstructed), or failed. The table in panel a) shows these results, along
with markers showing statistical significance in a two proportion z-test
comparing success rates of the algorithms. For each successful reconstruction,
we measured the Frechet distance and spatial distance from the manual ground
truth in order to evaluate the precision of the reconstructions. These
distances are shown as blue points in b), overlaid with standard box and
whisker plots (center line, median; box limits, upper and lower quartiles;
whiskers, 1.5x interquartile range; points, outliers).

###  2.6. Proof of Concept Graphical User Interface

Since our method is formulated as an optimal control problem conditioned on
the start and states it is most suited for semi-automatic neuron
reconstruction where a tracer could click on two different fragments along a
neuronal process, and the algorithm would fill in the gap. The user could
proceed to trace several neuronal processes until the full neuron is
reconstructed. We designed a proof of concept graphical user interface based
on this workflow and show an example set of reconstructions made using our GUI
in Figure 8. The GUI relies on the visualization software napari (napari
contributors, ).

Figure 8. We designed a proof of concept graphical user interface where the
user is presented with the image (a)) and the fragments (b)). The user can
then click on two fragments and generate the most probable curve between them.
c) shows three partial reconstructions (red, green and blue) of different
neurons using the GUI.

##  3\. Discussion

This paper presents a hidden Markov model based reconstruction algorithm that
connects fragments generated by appearance modeling. Our method converts an
image mask into a set of fragments and thus can be applied to the output of an
image segmentation model. We chose Ilastik to generate image masks because of
its convenient graphical user interface, and high performance on a small
number of samples (Berg et al., 2019). However, masks could also be generated
using a deep learning based model such as Liu et al. (2018); Li and Shen
(2019) or Wang et al. (2021).

These fragments are assembled based on the associated appearance model score
of the observed image and the discrete numerical curvature and distance of
adjacent fragments. In the Methods, we derive the Bayes posterior distribution
of the hidden state sequence encoding the axon reconstruction. We also show
that applying the polynomial time Viterbi algorithm is not possible for
maximum a-posteriori estimation since the state space has to grow to account
for possible cycles in the unordered image domain. We address the problem with
cycles with a modified procedure for defining the path probabilities making it
feasible to efficiently calculate the globally optimal neuron path. The
solution for efficiently generating the globally optimal path implies that the
local minimum associated to gradient and active appearance models solutions is
resolved in this setting.

We apply the algorithm to the fixed endpoint problem in two-photon images of
mouse neurons. In a dataset of 35 partial axons, our algorithm successfully
reconstructs more axons than existing state of the art algorithms (Figure 7).
We observed that the most common failure mode in this dataset was when there
are extended ($>\sim 10$ micron) stretches of the putative axonal path where
there is significant loss of luminance signals leading to highly censored
fragment generation. This implies that the maximally probable HMM procedure is
only effective if paired with effective voxel classification tools. The
algorithm can also fail in areas densely populated with neuronal processes. We
demonstrate that proper selection of the hyperparameters to reflect the
density of the fragments and the geometry of the underlying neurons can
resolve these issues. Our algorithm is specifically adapted for reconstructing
axons in projection neurons in datasets such as MouseLight in two ways. First,
the high image quality as indicated by large KL-divergence values (Figure 3a),
makes it straightforward to build an effective foreground-background
classifier which is an essential part of our HMM state generation. Secondly,
our algorithm encodes the geometric properties of axons such as curvature,
which allows our solutions to adapt to the occasional sharp turns in
projection axons.

The success rates of the other algorithms on our dataset is quite low,
considering the performances that they achieved in their accompanying
publications. There are two likely reasons for this, dataset differences, and
sub-optimal algorithm settings.

When an algorithm is validated on one type of data, the results do not
necessarily hold for datasets of a different type. Since none of the existing
algorithms had been designed for the MouseLight data, the unique details of
our dataset could lead to reduced performance. For example, several subvolumes
in our dataset contain multiple neurons, while two of the algorithms, APP2 and
Advantra, were explicitly designed for images containing single neurons. Snake
is designed to handle multiple neurons, but is largely validated on DIADEM, a
single neuron dataset (Peng et al., 2015). Other dataset differences include
different image resolutions (or levels of anisoptropy), different image
encodings (8 bit vs. 16 bit), and different signal to noise ratios. Lastly,
the existing algorithms are designed to reconstruct all dendrites and axons
simultaneously while our task of reconstructing a single section of the axon
does not give the algorithms credit for successfully reconstructing other
parts of the neuron.

It is also important to note that reconstruction algorithms can be sensitive
to hyperparameter settings. All algorithms had different hyperparameter
options except for Snake. While we tried various hyperparameter settings for
all algorithms, the only non-default setting that clearly improved
reconstruction performance was the binarization threshold setting in GTree,
all other settings we left to default. It is likely that these settings were
not optimal for our dataset, but it is quite time intensive for a typical user
to quantitatively determine the optimal settings. Detailed and accessible
software documentation makes the process of choosing effective algorithm
settings more efficient.

Figures showing the common failure modes for some of the algorithms are shown
in the Supplement S5. Advantra and Snake produced incoherent reconstructions
on the dataset, and the common failure modes for GTree were early termination
of reconstruction, or severing an axon into multiple components. Despite the
possible reasons for sub-optimal performance of the other algorithms, we
provide evidence that our algorithm is competitive with state of the art
methods for reconstructing neuronal processes. Future benchmark comparisons
could include reinforcement learning, or recurrent neural network approaches,
which have become prevalent in sequential decision processes. However, there
is not much scientific literature on these approaches to neuron reconstruction
with accompanying functional code.

We have scaled up the ViterBrain pipeline to process images of $3332\times
3332\times 1000$ voxels, representing one cubic millimeter of tissue (Figure
S6), and we will continue to improve the pipeline until it can be run on whole
brain images. The traces generated by our pipeline could also be paired with
tools like the one presented in Li et al. (2020) which can turn traces into
full neuron segmentations complete with axon and dendrite thickness
measurements.

The code used in this work is available in our open-source Python package
brainlit:
http://brainlit.neurodata.io/, and a tutorial on how to use the code is
located at:
http://brainlit.neurodata.io/link_stubs/hmm_reconstruction_readme_link.html.

##  4\. Methods

###  4.1. The Bayesian Appearance Imaging Model

Our Bayes model is comprised of a prior which models the axons as geometric
objects and a likelihood which models the image formation process.

We model the axons as simply connected curves in $\mathbb{R}^{3}$ written as a
function of arclength

$c(\ell),\ell\geq 0\ ,c(\ell)\in\mathbb{R}^{3}\ .$

We denote the entire axon curve in space as
$c(\cdot)\coloneqq\\{c(\ell),\ell\geq 0\\}$. To interface the geometric object
with the imaging volume we represent the underlying curve $c(\cdot)$ as a
delta-dirac impulse train in space. We view the imaging process as the
convolution of the delta-dirac impulse train with the point-spread kernel of
the imaging platform. We take the point-spread function of the system to be
roughly one micron in diameter, implying that the axons are well resolved. The
fluorescence process given the axon contour is taken as a relatively narrow
path through the imaging domain ($\sim 1\;\upmu$m diameter) with relatively
uniform luminance.

We take the image to be defined over the voxel lattice
$D=\cup_{i\in\mathcal{Z}^{m^{3}}}\Delta y_{i}\subset\mathbb{R}^{3}$ with
centers $y_{i}\in\Delta y_{i}$. We model the image as a random field
$\\{I_{y_{i}},\,\Delta y_{i}\in D\\}$ whose elements are independent when
conditioned on the underlying axon geometry, similar to an inhomogeneous
Poisson process as described in Snyder and Miller (2012). We denote the image
random field associated to any subset of sites $Y\subset D$, with joint
probability conditioned on the axon:

(1a) $\displaystyle I_{Y}=\\{I_{y_{i}}:\Delta y_{i}\in Y\\}\ ;$ (1b)
$\displaystyle P(I_{Y}|c(\ell),\ell\geq 0)=\prod_{y\in
Y}p(I_{y}|c(\ell),\ell\geq 0)\ .$

Because of the conditional independence, the marginal distribution determines
the global joint probabilities. Of course, while the conditional probabilities
factor and are conditionally independent, the axon geometry is unknown and the
measured image random field is completely connected if the latent axon process
is removed. We adopt a two hypothesis formulation $\\{f,b\\}$ corresponding to
a foreground-backround model for the images where the marginal probability of
a voxel intensity is:

(2a) $\displaystyle P(I_{y})=\alpha_{1}(I_{y}),y\in\text{foreground}$ (2b)
$\displaystyle P(I_{y})=\alpha_{0}(I_{y}),y\in\text{background}$

Our conditional independence assumption implies that the joint distribution of
group of foreground or group of backround voxels can be decomposed into the
corresponding marginal distributions over the foreground-backround models
(2a),(2b); defining the foreground and backround sets $Y=Y_{f}\cup Y_{b}$,
then

(3a) $\displaystyle P(I_{Y})=\prod\nolimits_{y\in
Y_{f}}\alpha_{1}(I_{y})\,\prod\nolimits_{y\in Y_{b}}\alpha_{0}(I_{y})\ .$ We
define the notations for the joint probability of the set in the foreground
for example (or background) (3b)
$\displaystyle\alpha_{1}(I_{Y_{f}})\coloneqq\prod\nolimits_{y\in
Y_{f}}\alpha_{1}(I_{y}),\ Y_{f}\subset\text{foreground}\ .$

Despite the Poisson nature of the image acquisition process, simple scaling or
shifting of the imaging data would mean the image intensities are no longer
Poisson. To accommodate this effect, we estimate the foreground-background
intensity distributions ($\alpha_{1}(\cdot)$ and $\alpha_{0}(\cdot)$
respectively) nonparametrically. The simplest nonparametric density estimation
technique is using histograms. However, it can be difficult to choose the
origin and bin width of histograms, so we opt for a kernel density estimate
(KDE) approach. We estimate $\alpha_{0}(\cdot)$, $\alpha_{1}(\cdot)$ by
labeling a subset of the data as foreground/background then fitting Gaussian
KDEs to the labeled data (see Figure 3b). We use the scipy implementation of
Gaussian KDEs (Virtanen et al., 2020), with Scott’s rule to determine the
bandwidth parameter (Scott, 2015). Under some assumptions on the derivatives
of the underlying density, our approach converges to the true density as the
number of samples increases (Theorem 6.1 in Scott (2015)). Further, the
Scott’s rule choice for bandwidth is (approximately) optimal with respect to
mean integrated square error.

Our approach allows us to empirically estimate the intensity distributions
while maintaining the independent increments property of a spatial Poisson
process, which is consistent with the autocorrelation curves depicted in
Figure 3a. This is the key property that allows us to factor joint
probabilities into products of marginals for sets of voxels.

We note that the foreground-background imaging model allows us to estimate the
error rate of classifying a voxel as either foreground or background. In the
Neyman-Pearson framework, foreground-background classification is a simple two
hypothesis testing problem and the most powerful test at a given type 1 error
rate is the log likelihood ratio test. The Kullback-Leibler (KL) divergence
between the foreground and background distributions gives the exponential rate
at which error rates converge to zero as the number of independent,
identically distributed samples increases (Cover and Thomas, 1991). In the
case of using Gaussians to model the foreground-background distributions, the
KL divergence reduces to the squared signal to noise ratio. In the absence of
normality, the KL Divergence is the general information theoretic measure of
image quality for arbitrary distributions. We propose KL Divergence as an
important statistic in evaluating the quality of fluorescent neuron images.

###  4.2. The Prior Distribution via Markov State Representation on the Axons
Fragments

Our representation of the observed image is as a hidden Markov random field
with the axon as the hidden latent structure. Given the complexity of sub-
micron resolution images, we build an intermediate data structure at the
micron scale that we call fragments $F\subset D$ defined as collections of
voxels without any assumed global ordering between them. Each fragment
represents a portion of a neuronal process, with a natural orientation given
by one end that is closer to the soma. Depicted in Figure 1 are the fragments
shown via different colors. The fragments are a coarser scale voxel and can be
viewed analogously as higher order features.

The axon reconstruction problem becomes the reassembly of the fragments along
with the imputation of the censored fragments. In the image examples presented
in this work, the complexity of the space of fragments is approximately
$|\mathcal{F}|=100$,$000$ for a cubic millimeter of projection neuron image
data.

We exploit the computational structures of hidden Markov models (HMMs) when
the underlying latent structure is absolutely ordered so that dynamic
programming can efficiently compute globally optimal state sequence estimates.
From the set of fragments we compute a set of states for the HMM. The states
are a simplified, abstract representation of the fragments that contain the
minimum information required to specify the HMM. Each state includes endpoints
$x^{0},x^{1}\in\mathbb{R}^{3}$ in order to compute "gap" or "censored"
probabilities, and unit length tangents $\tau^{0},\tau^{1}\in\mathbb{R}^{3}$
associated with the endpoionts in order to compute curvature. Each fragment
generates two states, one for each orientation. The two states are identical
except their endpoints are swapped, and their tangents are swapped and
reversed. We denote the natural mapping from state to fragment by
$F:s\in\mathcal{S}\mapsto F(s)\in\mathcal{F}$.

The collection of states $\mathcal{S}=\\{s\\}$ is the finite state space of
the HMM, and our goal is to estimate the state sequence $(s_{1},...,s_{n})$
that follows the neuronal process.

Our algorithms exploit two splitting properties, the Markov nature of the
state sequence and the splitting of the random field image conditioned on the
state sequence. We use the notation $s_{i:j}:=(s_{i},s_{i+1},\dots,s_{j})$ for
partial state sequences. We model the state sequence $s_{1},\dots,s_{n}$ as
Markov with splitting property:

$p(s_{i+1:n},s_{1:i-1}|s_{i})=p(s_{i+1:n}|s_{i})p(s_{1:i-1}|s_{i})\ ,\
i=2,\dots,n-1\ ,$

which implies the 1-order Markov property
$p(s_{i}|s_{i-1},s_{1:i-2})=p(s_{i}|s_{i-1})$.

We define the transition probabilities with a Boltzmann distribution with
energy $U$:

$\displaystyle p(s_{i}|s_{i-1})$
$\displaystyle=\frac{e^{-U(s_{i-1},s_{i})}}{Z(s_{i-1})}\ ,\ \text{with}\
Z(s_{i-1})=\sum_{s_{i}\in\mathcal{S}}e^{-U(s_{i-1},s_{i})}\ ,$ with energy
given by $\displaystyle U(s_{i-1},s_{i})$
$\displaystyle=\alpha_{d}|x_{i}^{0}-x_{i-1}^{1}|^{2}+\alpha_{\kappa}\kappa(s_{i-1},s_{i})^{2}\
,$

where $|\cdot|$ is the standard Euclidean norm. The two hyperparameters
$\alpha_{d}$ and $\alpha_{\kappa}$, determine the influence of distance and
curvature, respectively, on the probability of connection between two neuronal
fragments. The term $\kappa(s_{i-1},s_{i})$ approximates the curvature of the
path connecting $s_{i-1}$ to $s_{i}$ as follows:

Define
$\tau_{c}\coloneqq\frac{x_{i}^{0}-x_{i-1}^{1}}{||x_{i}^{0}-x_{i-1}^{1}||}$
which is the normalized vector connecting $s_{i-1}$ to $s_{i}$. Then we can
approximate the squared curvature at $x_{i-1}^{1}$ and $x_{i}^{0}$ with
$\kappa_{1}(s_{i-1},s_{i})^{2}=1-\tau_{i-1}^{1}\cdot\tau_{c}$ and
$\kappa_{2}(s_{i-1},s_{i})^{2}=1-\tau_{c}\cdot\left(-\tau_{i}^{0}\right)$
respectively. These formulas are derived in Supplement S2 and they approximate
curvature as modeled in Athey et al. (2021). Finally, we approximate average
squared curvature with the arithmetic mean:

$\displaystyle\kappa(s_{i-1},s_{i})^{2}=\frac{\kappa_{1}(s_{i-1},s_{i})^{2}+\kappa_{2}(s_{i-1},s_{i})^{2}}{2}$

We control the computational complexity associated to computing prior
probabilities for all $|\mathcal{S}|^{2}$ states by restricting the possible
state transitions. We set to 0 probability any transitions where the distance
between endpoints is greater than $15\,\upmu$m or the angle between states is
greater than $150^{\circ}$. We also set the probability that a state
transitions to itself to zero.

###  4.3. Global Maximally Probable Solutions via $O(n|\mathcal{S}|^{2})$
Calculation

The maximum a-posteriori (MAP) state sequence is defined as the maximizer of
the posterior probability (MAP) over the state sequences
$s_{1:n}\in{\mathcal{S}}^{n}$, with $|\mathcal{S}|$ finite:

(4) $\hat{s}_{1:n}:={\mathrm{argmax}}_{s_{1:n}\in{\mathcal{S}}^{n}}\log
p(s_{1:n}|I_{D})\ .$

The solution space has cardinality $|\mathcal{S}|^{n}$ so it is infeasible to
compute the global maximizer by exhaustive search. Our approach is to rewrite
the probability recursively in order to use the Viterbi algorithm and dynamic
programming with $O(n|\mathcal{S}|^{2})$ time complexity. We rewrite the MAP
estimator in terms of the joint probability:

$\hat{s}_{1:n}={\mathrm{argmax}}_{s_{1:n}\in\mathcal{S}^{n}}p(s_{1:n}|I_{D})={\mathrm{argmax}}_{s_{1:n}\in\mathcal{S}^{n}}p(s_{1:n},I_{D})\
.$

The image random field is split or conditionally independent conditioned on
the fragment states:

$p(I_{F(s_{i})},I_{D\setminus
F(s_{i})}|s_{i},)=p(I_{F(s_{i})}|s_{i})p(I_{D\setminus F(s_{i})}|s_{i})\ ,$

which implies $p(I_{F(s_{i})}|s_{i},I_{D\setminus
F(s_{i})})=p(I_{F(s_{i})}|s_{i})$.

Define the indicator function $\delta_{A}(x)=1$ for $x\in A$, $0$ otherwise.

######  Lemma 1.

Defining the shorthand notation identifying fragment sequences with the state
sequence:

$F_{1:n}\coloneqq(F(s_{1}),\dots,F(s_{n}))\ ,$

then for $n>1$ we have the joint probability:

(5) $\displaystyle p(s_{1:n},I_{D})$
$\displaystyle=\prod\nolimits_{i=2}^{n}\left(\frac{\alpha_{1}(I_{F_{i}})}{\alpha_{0}(I_{F_{i}})}\right)^{\delta_{D\setminus
F_{1:i-1}}(F_{i})}p(s_{i}|s_{i-1})\ p(s_{1},I_{D})\ .$

See Supplement S3 for the proof which rewrites the probability recursively
giving the factorization. The proof is similar to the classic HMM
decomposition in how it uses the two splitting properties, but there are two
differences. The first is that the probability needs to account for the full
image, including areas outside of the axon estimate, which explains the
presence of both $\alpha_{1}$ and $\alpha_{0}$ in Eq. 5. The second difference
is that if the state sequence $s_{1:n}$ contains repeated states, then the
corresponding image data should not be double counted in the probability. This
is enforced by the delta function $\delta(\cdot)$.

It is natural to take the negative logarithm of Eq. 5 to obtain a sum that
represents and the cost of a path through a directed trellis graph (Forney,
1973). Several algorithms exist that solve the shortest path problem in
$O(n|\mathcal{S}|^{2})$ complexity. However, we cannot use these algorithms
directly because the cost function is not “sequentially-additive” due to the
dependence of the indicator function on previous states in the sequence. In
the Supplement, section S4, we offer a example demonstrating that directly
applying the Viterbi algorithm to this problem does not generate the MAP
estimate.

We adjust our probabilistic representation on the $|\mathcal{S}|^{n}$ paths in
order to utilize shortest path algorithms such as Bellman-Ford or Dijkstra’s
(Bellman, 1958; Dijkstra, 1959). For this we note that the
$\frac{\alpha_{1}(I_{F_{i}})}{\alpha_{0}(I_{F_{i}})}$ term in Eq. (5) may
often be greater than 1. In the directed graph formulation (negative logarithm
transformation), this can lead to negative cycles in the graph of states. When
negative cycles exist, the shortest path problem is ill-posed. To avoid this
phenomena we remove the background component of the image from the joint
probability, which converts
$\frac{\alpha_{1}(I_{F_{i}})}{\alpha_{0}(I_{F_{i}})}$ to
$\alpha_{1}(I_{F_{i}})$, and converts our global posterior probability to our
path probability formulation.

######  Statement 1.

Define the most probable solution $s_{1:n}\in\mathcal{S}^{n}$ by the joint
probability
${\mathrm{argmax}}_{s_{1:n}\in\mathcal{S}^{n}}p(s_{1:n},I_{F_{1:n}})$. Then we
have

(6)
$\max_{s_{1:n}\in\mathcal{S}^{n}}p(s_{1:n},I_{F_{1:n}})=\max_{s_{1:n}\in\mathcal{S}^{n}}\prod\nolimits_{i=2}^{n}\left(\alpha_{1}(I_{F_{i}})\right)^{\delta_{D\setminus
F_{1:i-1}}(F_{i})}p(s_{i}|s_{i-1})\,p(s_{1},I_{F_{1}})\ .$

Further, if $\alpha_{1}(I_{y})\leq 1$ for all $y$, then the globally optimal
solution to the fixed start and end point problem is a nonrepeating state
sequence and can be obtained by computing the shortest path in a directed
graph where the vertices are the states, and the edge weight from state
$s_{i-1}$ to $s_{i}$ is given by:

(7) $\displaystyle e(s_{i-1},s_{i})=-\log\alpha_{1}(I_{F_{i}})-\log
p(s_{i}|s_{i-1})$

See Supplement S3 for proof. Our reconstruction problem has now become a
shortest path problem, and can be solved using one of the several dynamic
programming algorithms.

We note that since the path ($s_{1:n}$) defines the subset of the image in the
joint probability ($I_{F_{1:n}}$) we can define the probability as a function
of only the state sequence $\bar{p}(s_{1:n})\coloneqq p(s_{1:n},I_{F_{1:n}})$
emphasizing that we are solving the most probable path problem.

###  4.4. Implementation

####  4.4.1. Fragment generation:

Fragments are collections of voxels, or “supervoxels,” and can be viewed
analogously as higher order features such as edgelets or corners. As described
in section 4.2, identifying the subset of fragments that compose the axon,
then ordering them becomes equivalent to reconstructing the axon contour
model.

The first step of fragment generation is obtaining a foreground-background
mask, which could be obtained, for example, from a neural network, or by
simple thresholding. In this work, we use an Ilastik model that was trained on
three image subvolumes, each of which has three slices that were labeled (Berg
et al., 2019). During prediction, the probability predictions from Ilastik are
thresholded at 0.9, a conservative threshold that keeps the number of false
positives low.

The connected components of the thresholded image are split into fragments of
similar size by identifying the voxel $v$ with the largest predicted
foreground probability and placing a ball $B_{v}$ with radius $7\,\upmu$m on
that voxel. The voxels within $B_{v}$ are removed and the process is repeated
until the component is covered. The component is then split up into pieces by
assigning each voxel to the center point from the previous step, $v$, that is
closest to it. This procedure ensures that each fragment is no larger than a
ball with radius $7\,\upmu$m. At this size, it is reasonable to assume that
each fragment is associated with only one axon branch since no fragment is
large enough to extensively cover multiple branches.

Next, the endpoints $x^{0},x^{1}$ and tangents $\tau^{0},\tau^{1}$ are
computed as described in Supplement S1. Each fragment is simplified to the
line segment between its endpoints which is rasterized using the Bresenham
algorithm (Bresenham, 1965). Briefly, the Bresenham algorithm identifies the
image axis along which the line segment has the largest range and samples the
line once every voxel unit along that axis. Then, the other coordinates are
chosen to minimize the distance from the continuous representation line
segment.

####  4.4.2. Imputing Fragment Deletions

In practice the imaging data may exhibit significant dropouts leading to
significant fragment deletions. While computing the likelihood of the image
data, we augment the gaps between any pair of connected fragments in
$F_{1},F_{2},\dots$ by augmenting the sequence with imputed fragments
$F_{1},\bar{F}_{1},F_{2},\bar{F}_{2},\dots\ .$, with $\bar{F}_{i}\subset D$
the imputed line of voxels which forms the connection between the pair
$F_{i},F_{i+1}$. For this define the start and endpoint of each fragment as
$x^{0}(F)\in\mathbb{R}^{3},x^{1}(F)\in\mathbb{R}^{3}$ with line segment
connecting each pair:

$L_{i,i+1}=\\{y:y=ax^{1}(F_{i})+(1-a)x^{0}(F_{i+1}),\ a\in[0,1]\\}\ .$

The imputed fragment $\bar{F}_{i}\subset D$ for each pair $(F_{i},F_{i+1})$ is
computed by rasterizing $L_{i,i+1}$ with the Bresenham algorithm.

The likelihood of the sequence of fragments augmented by the imputations
becomes

$p(s_{1:n},I_{F_{1:n}})=\prod_{i=2}^{n}\alpha_{1}(I_{F_{i}})^{\delta_{D\setminus
F_{1:i-1}}(F_{i})}\alpha_{1}(I_{\bar{F}_{i}})p(s_{i}|s_{i-1})p(s_{1},I_{F_{1}})$

####  4.4.3. Initial and Endpoint Conditions

We take the initial conditions to represent

$p(s_{1},I_{F_{1}})=\pi(s_{1})p(I_{F_{1}}|s_{1})\ ,$

with $\pi$ the prior on initial state. For all of our axon reconstructions we
specify an axonal fragment as the start state $s_{start}$ and set
$\pi(s_{1})\coloneqq\delta_{s_{start}}(s_{1})$.

The endpoint conditions are defined via a user specified terminal state
$s_{term}$ where the path ends giving the maximization:

$\max_{s_{1:n}\in\mathcal{S}^{n}}p(s_{1:n},I_{F_{1:n}}|s_{n}=s_{term})\ .$

The marginal probability on the terminal state always transitions to itself,
so that $p(s_{n}=s_{term})=\delta_{s_{term}}(s_{n})$. Thus, a state sequence
solution of length $n$ may end in multiple repetitions of $s_{term}$, such as

$s_{1:n}=\\{s_{1},s_{2},...,s_{n^{\prime}},s_{term},s_{term},...,s_{term}\\}.$

###  4.5. Accuracy Metrics

We applied several state of the art reconstruction algorithms to several
neurons in the brain samples from the MouseLight Project from HHMI Janelia
(Winnubst et al., 2019). In this dataset, projection neurons were sparsely
labeled then imaged with a two-photon microscope at a voxel resolution of
$0.3\times 0.3\times 1\upmu$m. Each axon reconstruction is generated semi-
automatically by two independent annotators. The MouseLight reconstructions
are sampled roughly every $10\upmu$m, so in some cases we retraced the axons
at a higher sampling frequency in order to obtain more precise accuracy
metrics.

We quantified reconstruction accuracy using two metrics, the first of which is
Frechet distance. Frechet distance is commonly described in the setting of dog
walking, where both the dog and owner are following their own predetermined
paths. The Frechet distance between the two paths then is the minimum length
dog leash needed to complete the walk, where both dog and owner are free to
vary their walking speeds but are not allowed to backtrack. In our setting we
compute the Frechet distance between two discrete paths
$P:\\{1,...,L_{p}\\}\rightarrow\mathbb{R}^{3}$,
$Q:\\{1,...,L_{q}\\}\rightarrow\mathbb{R}^{3}$ as defined in Eiter and Mannila
(1994). In this definition, a coupling between $P$ and $Q$ is defined as a
sequence of ordered pairs:

$\displaystyle(P[a_{1}],Q[b_{1}]),(P[a_{2}],Q[b_{2}]),...,(P[a_{K}],Q[b_{K}])$

where the following conditions are put on $\\{a_{k}\\},\\{b_{k}\\}$ to ensure
that they enumerate through the whole sequences $P$ and $Q$:

  * •

$a_{1},b_{1}=1$

  * •

$a_{N}=L_{p}$, $b_{N}=L_{q}$

  * •

$a_{k}=a_{k-1}$ or $a_{k}=a_{k-1}+1$

  * •

$b_{k}=b_{k-1}$ or $b_{k}=b_{k-1}+1$

Then the discrete Frechet distance is defined as:

$\displaystyle\delta_{dF}(P,Q)$ $\displaystyle=\min_{\text{coupling
}\\{a_{k}\\},\\{b_{k}\\}}\max_{k\in\\{1,...,K\\}}\left|P\left[a_{k}\right]-Q\left[b_{k}\right]\right|$

We use the standard Euclidean norm for $|\cdot|$. The discrete Frechet
distance is an upper bound to the continuous Frechet distance between
polygonal curves, and it can be computed more efficiently. Further, if we take
a discrete Frechet distance of zero to be an equivalence relation, then
$\delta_{dF}$ is a metric on this set of equivalence classes and thus is a
natural way to compare non-branching neuronal reconstructions. In this work,
all reconstruction are sampled at at least one point per micron.

Various other performance metrics have been proposed, including an arc-length
based precision and recall (Wang et al., 2011), a critical node matching based
Miss-Extra-Scores (MES) (Xie et al., 2011) and a vertex matching based spatial
distance (SD) (Peng et al., 2010b). We chose to compute SD since it gives a
picture of the average spatial distance between two reconstructions. This
complements the Frechet distance described earlier, which computes the maximum
spatial distance between two reconstructions.

The first step in computing the SD from reconstruction $P$ to reconstruction
$Q$ is, for each point in $P$, finding the distance to the closest point in
$Q$. Directed divergence (DDIV) of $P$ from $Q$ is then defined as the average
of all these distances. Then, SD is computed by averaging the DDIV from $P$ to
$Q$ and the DDIV from $Q$ to $P$.

##  5\. Acknowledgements

We thank the MouseLight team at HHMI Janelia for providing us with access to
this data, and answering our questions about it.

##  6\. Author Contributions

MIM helped to develop the HMM and probabilistic model and DT advised on the
theoretical direction of the manuscript. UM coordinated the data acquisition
for the experiments. TA designed the study, implemented the software, and
managed the manuscript text/figures. All authors contributed to manuscript
revision.

##  7\. Conflict of Interest Statement

MIM owns a significant share of Anatomy Works with the arrangement being
managed by Johns Hopkins University in accordance with its conflict of
interest policies. The remaining authors declare that the research was
conducted in the absence of any commercial or financial relationships that
could be construed as a potential conflict of interest.

The funders had no role in study design, data collection and analysis,
decision to publish, or preparation of the manuscript.

##  8\. Funding

This work is supported by the National Institutes of Health grant RF1MH121539.

## References

  * Acciai et al. (2016) L. Acciai, P. Soda, and G. Iannello.  Automated neuron tracing methods: an updated account.  _Neuroinformatics_ , 14(4):353–367, 2016.
  * Athey et al. (2021) T. L. Athey, J. Teneggi, J. T. Vogelstein, D. J. Tward, U. Mueller, and M. I. Miller.  Fitting splines to axonal arbors quantifies relationship between branch order and geometry.  _Frontiers in Neuroinformatics_ , 2021.
  * Bellman (1958) R. Bellman.  On a routing problem.  _Quarterly of applied mathematics_ , 16(1):87–90, 1958.
  * Berg et al. (2019) S. Berg, D. Kutra, T. Kroeger, C. N. Straehle, B. X. Kausler, C. Haubold, M. Schiegg, J. Ales, T. Beier, M. Rudy, K. Eren, J. I. Cervantes, B. Xu, F. Beuttenmueller, A. Wolny, C. Zhang, U. Koethe, F. A. Hamprecht, and A. Kreshuk.  ilastik: interactive machine learning for (bio)image analysis.  _Nature Methods_ , Sept. 2019.  ISSN 1548-7105.  doi: 10.1038/s41592-019-0582-9.  URL https://doi.org/10.1038/s41592-019-0582-9.
  * Bresenham (1965) J. E. Bresenham.  Algorithm for computer control of a digital plotter.  _IBM Systems journal_ , 4(1):25–30, 1965.
  * Chen et al. (2015) H. Chen, H. Xiao, T. Liu, and H. Peng.  Smarttracing: self-learning-based neuron reconstruction.  _Brain informatics_ , 2(3):135–144, 2015.
  * Choromanska et al. (2012) A. Choromanska, S.-F. Chang, and R. Yuste.  Automatic reconstruction of neural morphologies with multi-scale tracking.  _Frontiers in neural circuits_ , 6:25, 2012.
  * Cohen (1991) L. D. Cohen.  On active contour models and balloons.  _CVGIP: Image Understanding_ , 53(2):211–218, 1991.  ISSN 1049-9660.  doi: https://doi.org/10.1016/1049-9660(91)90028-N.  URL https://www.sciencedirect.com/science/article/pii/104996609190028N.
  * Cover and Thomas (1991) T. M. Cover and J. A. Thomas.  _Elements of information theory_ , volume 2.  Wiley, 1991.
  * Dai et al. (2019) T. Dai, M. Dubois, K. Arulkumaran, J. Campbell, C. Bass, B. Billot, F. Uslu, V. De Paola, C. Clopath, and A. A. Bharath.  Deep reinforcement learning for subpixel neural tracking.  In _International Conference on Medical Imaging with Deep Learning_ , pages 130–150. PMLR, 2019.
  * Dijkstra (1959) E. W. Dijkstra.  A note on two problems in connexion with graphs.  _Numerische mathematik_ , 1(1):269–271, 1959\.
  * Eiter and Mannila (1994) T. Eiter and H. Mannila.  Computing discrete fréchet distance.  Technical report, Citeseer, 1994.
  * Forney (1973) G. D. Forney.  The viterbi algorithm.  _Proceedings of the IEEE_ , 61(3):268–278, 1973\.
  * Friedmann et al. (2020) D. Friedmann, A. Pun, E. L. Adams, J. H. Lui, J. M. Kebschull, S. M. Grutzner, C. Castagnola, M. Tessier-Lavigne, and L. Luo.  Mapping mesoscale axonal projections in the mouse brain using a 3d convolutional network.  _Proceedings of the National Academy of Sciences_ , 117(20):11068–11075, 2020.
  * Kass et al. (1988) M. Kass, A. Witkin, and D. Terzopoulos.  Snakes: Active contour models.  _International journal of computer vision_ , 1(4):321–331, 1988.
  * Khaneja et al. (1998) N. Khaneja, M. I. Miller, and U. Grenander.  Dynamic programming generation of curves on brain surfaces, 1998.
  * Li and Shen (2019) Q. Li and L. Shen.  3d neuron reconstruction in tangled neuronal image with deep networks.  _IEEE transactions on medical imaging_ , 39(2):425–435, 2019.
  * Li et al. (2017) R. Li, T. Zeng, H. Peng, and S. Ji.  Deep learning segmentation of optical microscopy images improves 3-d neuron reconstruction.  _IEEE transactions on medical imaging_ , 36(7):1533–1541, 2017.
  * Li et al. (2019) R. Li, M. Zhu, J. Li, M. Bienkowski, N. N. Foster, H. Xu, T. Ard, I. Bowman, C. Zhou, M. Veldman, X. W. Yang, H. Hintiryan, J. Zhang, and H. Dong.  Precise segmentation of densely interweaving neuron clusters using g-cut.  _Nature Communications_ , 10, 2019.
  * Li et al. (2020) S. Li, T. Quan, H. Zhou, Q. Huang, T. Guan, Y. Chen, C. Xu, H. Kang, A. Li, L. Fu, et al.  Brain-wide shape reconstruction of a traced neuron using the convex image segmentation method.  _Neuroinformatics_ , 18(2):199–218, 2020.
  * Liu et al. (2018) M. Liu, H. Luo, Y. Tan, X. Wang, and W. Chen.  Improved v-net based image segmentation for 3d neuron reconstruction.  In _2018 IEEE International Conference on Bioinformatics and Biomedicine (BIBM)_ , pages 443–448. IEEE, 2018.
  * (22) napari contributors.  napari: a multi-dimensional image viewer for python.  URL doi:10.5281/zenodo.3555620.
  * Peng et al. (2010a) H. Peng, Z. Ruan, D. Atasoy, and S. Sternson.  Automatic reconstruction of 3d neuron structures using a graph-augmented deformable model.  _Bioinformatics_ , 26(12):i38–i46, 2010a.
  * Peng et al. (2010b) H. Peng, Z. Ruan, F. Long, J. H. Simpson, and E. W. Myers.  V3d enables real-time 3d visualization and quantitative analysis of large-scale biological image data sets.  _Nature biotechnology_ , 28(4):348–353, 2010b.
  * Peng et al. (2015) H. Peng, E. Meijering, and G. A. Ascoli.  From diadem to bigneuron, 2015.
  * Peng et al. (2017) H. Peng, Z. Zhou, E. Meijering, T. Zhao, G. A. Ascoli, and M. Hawrylycz.  Automatic tracing of ultra-volumes of neuronal images.  _Nature methods_ , 14(4):332–333, 2017.
  * Quan et al. (2016) T. Quan, H. Zhou, J. Li, S. Li, A. Li, Y. Li, X. Lv, Q. Luo, H. Gong, and S. Zeng.  Neurogps-tree: automatic reconstruction of large-scale neuronal populations with dense neurites.  _Nature methods_ , 13(1):51–54, 2016.
  * Rabiner and Juang (1986) L. Rabiner and B. Juang.  An introduction to hidden markov models.  _ieee assp magazine_ , 3(1):4–16, 1986.
  * Radojevic and Meijering (2019) M. Radojevic and E. Meijering.  Automated neuron reconstruction from 3d fluorescence microscopy images using sequential monte carlo estimation.  _Neuroinformatics_ , 17, 07 2019.  doi: 10.1007/s12021-018-9407-8.
  * Radojević and Meijering (2017) M. Radojević and E. Meijering.  Automated neuron tracing using probability hypothesis density filtering.  _Bioinformatics_ , 33(7):1073–1080, 01 2017.  ISSN 1367-4803.  doi: 10.1093/bioinformatics/btw751.  URL https://doi.org/10.1093/bioinformatics/btw751.
  * Ronneberger et al. (2015) O. Ronneberger, P. Fischer, and T. Brox.  U-net: Convolutional networks for biomedical image segmentation.  In _International Conference on Medical image computing and computer-assisted intervention_ , pages 234–241. Springer, 2015.
  * Scott (2015) D. W. Scott.  _Multivariate density estimation: theory, practice, and visualization_.  John Wiley & Sons, 2015.
  * Snyder and Miller (2012) D. L. Snyder and M. I. Miller.  _Random point processes in time and space_.  Springer Science & Business Media, 2012.
  * Turetken et al. (2013) E. Turetken, F. Benmansour, B. Andres, H. Pfister, and P. Fua.  Reconstructing loopy curvilinear structures using integer programming.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 1822–1829, 2013.
  * Virtanen et al. (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  _Nature Methods_ , 17:261–272, 2020.  doi: 10.1038/s41592-019-0686-2.
  * Wang et al. (2017) C.-W. Wang, Y.-C. Lee, H. Pradana, Z. Zhou, and H. Peng.  Ensemble neuron tracer for 3d neuron reconstruction.  _Neuroinformatics_ , 15(2):185–198, 2017.
  * Wang et al. (2021) H. Wang, C. Zhang, J. Yu, Y. Song, S. Liu, W. Chrzanowski, and W. Cai.  Voxel-wise cross-volume representation learning for 3d neuron reconstruction.  In _International Workshop on Machine Learning in Medical Imaging_ , pages 248–257. Springer, 2021.
  * Wang et al. (2011) Y. Wang, A. Narayanaswamy, C.-L. Tsai, and B. Roysam.  A broadly applicable 3-d neuron tracing method based on open-curve snake.  _Neuroinformatics_ , 9:193–217, 03 2011.  doi: 10.1007/s12021-011-9110-5.
  * Wang et al. (2019) Y. Wang, Q. Li, L. Liu, Z. Zhou, Z. Ruan, L. Kong, Y. Li, Y. Wang, N. Zhong, R. Chai, et al.  Teravr empowers precise reconstruction of complete 3-d neuronal morphology in the whole brain.  _Nature communications_ , 10(1):1–9, 2019.
  * Winnubst et al. (2019) J. Winnubst, E. Bas, T. A. Ferreira, Z. Wu, M. N. Economo, P. Edson, B. J. Arthur, C. Bruns, K. Rokicki, D. Schauder, et al.  Reconstruction of 1,000 projection neurons reveals new cell types and organization of long-range connectivity in the mouse brain.  _Cell_ , 179(1):268–281, 2019.
  * Xiao and Peng (2013) H. Xiao and H. Peng.  App2: automatic tracing of 3d neuron morphology based on hierarchical pruning of a gray-weighted image distance-tree.  _Bioinformatics_ , 29(11), June 2013.
  * Xie et al. (2011) J. Xie, T. Zhao, T. Lee, E. Myers, and H. Peng.  Anisotropic path searching for automatic neuron reconstruction.  _Medical image analysis_ , 15(5):680–689, 2011\.
  * Yang et al. (2019) J. Yang, M. Hao, X. Liu, Z. Wan, N. Zhong, and H. Peng.  Fmst: an automatic neuron tracing method based on fast marching and minimum spanning tree.  _Neuroinformatics_ , 17(2):185–196, 2019.
  * Zhou et al. (2021) H. Zhou, S. Li, A. Li, Q. Huang, F. Xiong, N. Li, J. Han, H. Kang, Y. Chen, Y. Li, et al.  Gtree: an open-source tool for dense reconstruction of brain-wide neuronal population.  _Neuroinformatics_ , 19(2):305–317, 2021.
  * Zhou et al. (2018) Z. Zhou, H.-C. Kuo, H. Peng, and F. Long.  Deepneuron: an open deep learning toolbox for neuron tracing.  _Brain informatics_ , 5(2):1–9, 2018.

## Supplementary Material

##  S1. Computing Endpoints and Tangents of Fragments

As explained in section 4.2, each fragment is a subset of the image $F\subset
D$. Each fragment is assumed to be associated with a segment of an underlying
neuron curve, and we want to estimate the locations of the two endpoints of
the fragment. First, we compute the length of the diagonal of the bounding box
that contains the fragment. We divide this length by two to get $R$. Then,
with each voxel $y$ in the fragment, we associate a set of voxels $N_{y}$
which is the intersection of the voxels in the fragment and the voxels within
distance $R$ from $y$. The voxel with the smallest set $N_{y}$ (by
cardinality) is chosen to be the first endpoint. The second endpoint is the
voxel that has the smallest set $N_{y}$ but is also farther than $R$ away from
the first endpoint.

Currently the tangents at the endpoints is just approximated by the difference
of the endpoints i.e.
$\tau^{0}=\frac{x_{0}-x_{1}}{||x_{0}-x_{1}||},\tau^{1}=-\tau^{1}$. This method
is based on the assumption that fragments are small enough that they are
approximately straight. We also experimented with approximating the endpoint
tangents by computing principal components of voxels near the endpoints, but
found it to be less robust for downstream reconstruction.

##  S2. Curvature Calculation for the Potential in the Markov Chain

The term $\kappa(s_{i-1},s_{i})$ is an approximation of the curvature of the
path that connects $s_{i-1}$ to $s_{i}$. For a curve with the tangent vector
$T(s)$, curvature is defined as
$\kappa(s)=\left|\left|\frac{dT}{ds}\right|\right|$. A finite difference
approximation of curvature is then:

$\displaystyle\kappa(s)^{2}=\left|\left|\frac{dT}{ds}\right|\right|^{2}$
$\displaystyle\approx\left|\left|T(s)-T(s-1)\right|\right|^{2}$
$\displaystyle=\left|\left|T(s)\right|\right|^{2}+\left|\left|T(s-1)\right|\right|^{2}-2T(s)\cdot
T(s-1)$ $\displaystyle=2(1-T(s)\cdot T(s-1))$

Thus,

(8) $\displaystyle\left|\left|\frac{dT}{ds}\right|\right|^{2}$
$\displaystyle\propto 1-T(s+1)\cdot T(s)$

We consider $\tau_{i},\tau_{i-1}$ and the normalized vector between the states
$\tau_{c}\coloneqq\frac{x_{i-1}^{1}-x_{i}^{0}}{|x_{i-1}^{1}-x_{i}^{0}|}$ as
samples of $T(s)$, then use Eq. (8) to estimate the curvature induced by
connecting state $s_{i-1}$ to $s_{i}$:

$\displaystyle(\kappa_{1})^{2}$ $\displaystyle=1-\tau_{i-1}^{1}\cdot\tau_{c}$
$\displaystyle(\kappa_{2})^{2}$ $\displaystyle=1-\tau_{c}\cdot(-\tau_{i}^{0})$
$\displaystyle\kappa(s_{i-1},s_{i})^{2}$
$\displaystyle=\frac{(\kappa_{1})^{2}+(\kappa_{2})^{2}}{2}$ Arithmetic mean

##  S3. Proofs

Recall the likelihood of a complete fragment under the foreground-background
model (Eq. 3b):

$\alpha_{k}(I_{F})\coloneqq\prod_{y\in F}\alpha_{k}(I_{y})\ ,k=0,1\ .$

Lemma 1:

For $n>1$ we have the recursion probability

(9a) $\displaystyle
p(s_{1:n},I_{D})=\left(\frac{\alpha_{1}(I_{F_{n}})}{\alpha_{0}(I_{F_{n}})}\right)^{\delta_{D\setminus
F_{1:n-1}}(F_{n})}p(s_{n}|s_{n-1})p(s_{1:n-1},I_{D})$

implying the factored probability:

(9b) $\displaystyle
p(s_{1:n},I_{D})=\prod_{i=2}^{n}\left(\frac{\alpha_{1}(I_{F_{i}})}{\alpha_{0}(I_{F_{i}})}\right)^{\delta_{D\setminus
F_{1:i-1}}(F_{i})}p(s_{i}|s_{i-1})p(s_{1},I_{D})$

###### Proof.

Factor the event $I_{D}=(I_{F_{1:n}},I_{D\setminus F_{1:n}})$, then

(10) $\displaystyle p(s_{1:n},I_{F_{1:n}},I_{D\setminus F_{1:n}})$
$\displaystyle=$ $\displaystyle
p(I_{F_{n}}|s_{1:n},I_{F_{1:n-1}},I_{D\setminus
F_{1:n}})p(s_{n}|s_{1:n-1},I_{D\setminus
F_{1:n}})p(s_{1:n-1},I_{F_{1:n-1}},I_{D\setminus F_{1:n}})$ $\displaystyle=$
$\displaystyle p(I_{F_{n}}|s_{n})^{\delta_{D\setminus
F_{1:n-1}}(F_{n})}p(s_{n}|s_{n-1})p(s_{1:n-1},I_{F_{1:n-1}},I_{D\setminus
F_{1:n}})$ $\displaystyle=$
$\displaystyle\alpha_{1}(I_{F_{n}})^{\delta_{D\setminus
F_{1:n-1}}(F_{n})}p(s_{n}|s_{n-1})p(s_{1:n-1},I_{F_{1:n-1}},I_{D\setminus
F_{1:n}})$

We rewrite the last term using the splitting property:

$\displaystyle p(s_{1:n-1},I_{F_{1:n-1}},I_{D\setminus F_{1:n}})$
$\displaystyle=$ $\displaystyle p(I_{F_{1:n-1}}|s_{1:n-1},I_{D\setminus
F_{1:n}})p(s_{1:n-1},I_{D\setminus F_{1:n}})$ $\displaystyle=$ $\displaystyle
p(I_{F_{1:n-1}}|s_{1:n-1})p(I_{D\setminus F_{1:n}}|s_{1:n-1})p(s_{1:n-1})$
$\displaystyle=$ $\displaystyle
p(I_{F_{1:n-1}}|s_{1:n-1})\frac{p(I_{D\setminus
F_{1:n-1}}|s_{1:n-1})}{p(I_{F_{n}}|s_{1:n-1})^{\delta_{D\setminus
F_{1:n-1}}(F_{n})}}p(s_{1:n-1})$ $\displaystyle=$
$\displaystyle\frac{1}{\alpha_{0}(I_{F_{n}})^{\delta_{D\setminus
F_{1:n-1}}(F_{n})}}p(s_{1:n-1},I_{F_{1:n-1}},I_{D\setminus F_{1:n-1}})$

with the last substitution following from the background model. Substituting
into 10 yields the probability written as a recursion 9a:

$p(s_{1:n},I_{D})=\left(\frac{\alpha_{1}(I_{F_{n}})}{\alpha_{0}(I_{F_{n}})}\right)^{\delta_{D\setminus
F_{1:n-1}}(F_{n})}p(s_{n}|s_{n-1})p(s_{1:n-1},I_{D})$

$\sqcap$$\sqcup$

Statement 1: Define the most probable solution $s_{1:n}\in\mathcal{S}^{n}$ by
the joint probability

${\mathrm{argmax}}_{s_{1:n}\in\mathcal{S}^{n}}p(s_{1:n},I_{F_{1:n}})$. Then we
have

(11)
$\max_{s_{1:n}\in\mathcal{S}^{n}}p(s_{1:n},I_{F_{1:n}})=\max_{s_{1:n}\in\mathcal{S}^{n}}\prod\nolimits_{i=2}^{n}\left(\alpha_{1}(I_{F_{i}})\right)^{\delta_{D\setminus
F_{1:i-1}}(F_{i})}p(s_{i}|s_{i-1})\,p(s_{1},I_{F_{1}})\ .$

Further, if $\alpha_{1}(I_{y})\leq 1$ for all $y$, then the globally optimal
solution to the fixed start and end point problem is a nonrepeating state
sequence and can be obtained by computing the shortest path in a directed
graph where the vertices are the states, and the edge weight from state
$s_{i-1}$ to $s_{i}$ is given by:

(12) $\displaystyle e(s_{i-1},s_{i})=-\log\alpha_{1}(I_{F_{i}})-\log
p(s_{i}|s_{i-1})$

###### Proof.

First we will demonstrate the factorization in Eq. 11:

$\displaystyle p(s_{1:n},I_{F_{1:n}})$
$\displaystyle=p(s_{1:n},I_{F_{1:n-1}},I_{F_{n}})$
$\displaystyle=p(I_{F_{n}}|s_{1:n},I_{F_{1:n-1}})p(s_{1:n},I_{F_{1:n-1}})$
$\displaystyle=p(I_{F_{n}}|s_{1:n},I_{F_{1:n-1}})p(I_{F_{1:n-1}}|s_{1:n})p(s_{1:n})$
$\displaystyle=p(I_{F_{n}}|s_{1:n},I_{F_{1:n-1}})p(I_{F_{1:n-1}}|s_{1:n-1})p(s_{n}|s_{1:n-1})p(s_{1:n-1})$
$\displaystyle=(\alpha_{1}(I_{F_{n}}))^{\delta_{D\setminus
F_{1:n-1}}(F_{n})}p(s_{n}|s_{n-1})p(I_{F_{1:n-1}}|s_{1:n-1})p(s_{1:n-1})\ .$

Applying this to factor the conditional probability , $I_{F_{1:i}}$,
$i=n-1,\dots,1$ gives the joint product. Next, we will show how this can be
solved using shortest path algorithms if $\alpha_{1}(I_{y})\leq 1$ for all
$y$:

First, we want to show that the globally optimal sequence is nonrepeating.
This is clear because if $\alpha_{1}(I_{y})\leq 1$, then every term in the
product in equation (6) is bounded by 1. Thus, for any state sequence that
repeats states, if we remove all elements between the two instances of the
repeated states, then this new sequence will have at least as much probability
$p(s_{1:n},I_{F_{1:n}})$.

For nonrepeating state sequences, the probability $p(s_{1:n},I_{F_{1:n}})$ of
(6) can be simplified:

$\displaystyle
p(s_{1:n},I_{F_{1:n}})=\prod\nolimits_{i=2}^{n}\left(\alpha_{1}(I_{F_{i}})\right)p(s_{i}|s_{i-1})p(s_{1},I_{F_{1}})\
.$

Taking the logarithm yields a sequentially additive cost function that can be
solved with shortest path algorithm on a graph with edge weights given by
Equation (12).

$\sqcap$$\sqcup$

##  S4. Viterbi Algorithm Counter-example

Here we present a simple counter example demonstrating that the indicator
function in 5 cannot be ignored, implying the globally optimal solution cannot
be obtained efficiently with the Viterbi algorithm. The details of the state
likelihoods and transition probabilities are given in Figure S1a.

Figure S1. Example showing that the Viterbi algorithm is invalid in our
setting. a) Transition probabilities, and image likelihood ratios of a two
state Markov model where, in our probabilistic model, the globally optimal
state sequence from state $1$ to state $2$ with length $4$ cannot be computed
with the classic Viterbi algorithm. Panels b) and c) depict paths through the
state space. The nodes represent the states and the line segments represent
state transitions, as in Figure 8 of Forney [1973]. Panel b) shows true
shortest path from state $1$ to state $2$ with length $4$. Panel c)
illustrates three iterations of the Viterbi algorithm, with the lines
representing the paths that are stored at each iteration. In this example, the
blue path at the bottom of panel c) is identified by naive application of the
Viterbi algorithm, but it differs from the true shortest path in panel b).
Each transition has a number above it indicating its contribution to the joint
probability -
$\left(\frac{\alpha_{1}(I_{F_{i}})}{\alpha_{0}(I_{F_{i}})}\right)^{\delta_{D\setminus
F_{1:i-1}}(F_{i})}p(s_{i}|s_{i-1})$ from Eqn. (5). The number at the end of
the path is the product of all such terms - the joint probability of the given
path.

Here, we remove the indicator function in 5 and apply the Viterbi algorithm in
the usual way. After three iterations, the Viterbi algorithm stores
$(1,1,1,2)$ as the highest probability path of length 4 that ends in state 2
(Figure S1c). The joint probability of this path is $1/800$. However, there is
an alternative path that the algorithm missed, $(1,2,1,2)$, which has higher
joint probability $1/400$ (Figure S1b). The algorithm did not identify the
globally optimal sequence because it did not “see” that the cost of
transitioning from state 1 to 2 would drop from $1/200$ to $1/2$ after
visiting state 2 the first time.

##  S5. Other Figures

Figure S2. The most common failure mode is an inadequate fragment space due to
signal dropout. a) shows maximally probable reconstructions (blue lines) that
deviate from the manual reconstructions (red lines). b) shows an example with
an inadequate fragment space in more detail. The full image subvolume is shown
in panel i) with the blue line depicting the most probable reconstruction and
the red line depicting the manual reconstruction. The close up views in panels
ii,iii) show that the image signal under the manual reconstruction is weak
(red dotted box). Panel iv) shows the fragments in different colors. Indeed,
there are only a few fragments generated in the region delineated by the red
dotted box, where the manual reconstruction is. Panel v) shows the manual
reconstruction in red and the algorithm’s reconstruction in blue. The
algorithm starts on the right path, then veers away from the true path at the
red arrow.  Figure S3. GTree produced sensible reconstructions on our dataset,
but often fell short of fully reconstructing the axon segment in question. We
observed two common failure modes with this algorithm. First, the
reconstruction would sometimes fail to extend far enough down the axon. An
example of this is shown in a) where the neuron reconstruction (brown overlay)
does not extend up the axon identified with the red arrow. The other failure
mode is when an axon would be split into several components such as in b).
Each colored overlay in b), including green, blue, and purple represents a
distinct reconstruction object produced by GTree. In other words, GTree
severed the neuron (emphasized by the red circle) into several components and
failed to fully reconstruct the axon. Figure S4. The Snake tracing algorithm
produced incoherent reconstructions on our dataset. a) shows one of the image
subvolumes, and the brown overlay in b) shows a reconstruction by Snake which
is a dense, convoluted set of tangled paths. This algorithm was applied to 10
different subvolumes and the results were similar each time. There were no
options to change any parameters in the Vaa3D implementation of this
algorithm, and the algorithm took over an hour to process any subvolume larger
than 50 megabytes, so we terminated its use after only 10 subvolumes. Figure
S5. The Advantra algorithm, as implemented in Vaa3D, did not produce coherent
reconstructions on the test dataset. Even under different hyperparameter
settings, the algorithm would often produce a reconstruction composed of only
a single point, like the one shown here. Figure S6. We applied the ViterBrain
pipeline to an image subvolume with dimensions $3332\times 3332\times 1000$
voxels, which encompasses a cubic millimeter of tissue. Pictured is a set of
maximum intensity projections of a downsampled version of this subvolume. a)
shows the image, b) shows the image with a color overlay depicting the neuron
fragments and c) shows a trace in red that follows an axon. The scale bar
represents $125$ microns.