 Research article
 Open Access
 Published:
Neuron anatomy structure reconstruction based on a sliding filter
BMC Bioinformatics volume 16, Article number: 342 (2015)
Abstract
Background
Reconstruction of neuron anatomy structure is a challenging and important task in neuroscience. However, few algorithms can automatically reconstruct the full structure well without manual assistance, making it essential to develop new methods for this task.
Methods
This paper introduces a new pipeline for reconstructing neuron anatomy structure from 3D microscopy image stacks. This pipeline is initialized with a set of seeds that were detected by our proposed Sliding Volume Filter (SVF), given a noncircular crosssection of a neuron cell. Then, an improved open curve snake model combined with a SVF external force is applied to trace the full skeleton of the neuron cell. A radius estimation method based on a 2D sliding band filter is developed to fit the real edge of the crosssection of the neuron cell. Finally, a surface reconstruction method based on nonparallel curve networks is used to generate the neuron cell surface to finish this pipeline.
Results
The proposed pipeline has been evaluated using publicly available datasets. The results show that the proposed method achieves promising results in some datasets from the DIgital reconstruction of Axonal and DEndritic Morphology (DIADEM) challenge and new BigNeuron project.
Conclusion
The new pipeline works well in neuron tracing and reconstruction. It can achieve higher efficiency, stability and robustness in neuron skeleton tracing. Furthermore, the proposed radius estimation method and applied surface reconstruction method can obtain more accurate neuron anatomy structures.
Background
Neuron morphology and structure information is critical for neuroscience research. Hence, reconstructing the entire anatomy structure of a neuron is an essential task in the field of neuron informatics [1, 2]. However, reconstructing the anatomy structure of a neuron artificially is labor intensive. Efficient, advanced methods for anatomy structure reconstruction of neurons are greatly demanded in this field. Specifically, with the rapid development of microscopic imaging technology, a wide range of scales of bioimages can be obtained, which is helpful for us to develop new methods and algorithms to meet the needs in neuroscience research [3, 4]. The reconstructed digital neuron structure, including axons and dendrites as well as thickness information, can be used in conjunction with electrophysiological simulations to determine the complex mechanisms of the nervous system [5, 6].
The computeraided manual neuron reconstruction method was first proposed in 1965 and was achieved by a biologist using a microscope [7]. Following this milestone, numerous algorithms and open softwares were introduced to reduce manual labor consisting of the boring task of tracing and analysis [8–11], but most of them were still limited to semiautomation and required manual validation by experts to achieve accurate reconstruction of whole neurons. Hence, the lack of powerful and effective computational tools for automatically reconstructing neuron cells has emerged as a major technical bottleneck in neuroscience research. This problem motivated the DIgital reconstruction of Axonal and DEndritic Morphology (DIADEM) challenge [12] and BigNeuron project [13, 14], which began in 2010 and 2015 respectively. They provided an opensource platform for researchers from all over the world and aimed to promote the development of computer algorithms for reconstructing the full anatomy structure of neurons. The data sets from DIADEM are most widely used in the domain of neuron reconstruction to date. However, the BigNeuron proposed some new challenges for the further research in the field of neuron reconstruction.
Generally speaking, before the DIADEM project, the neuron tracing methods were categorized into several types: shortest path methods [15, 16], minimum spanning tree methods [17, 18], sequential tracing methods [19, 20], skeletonization methods [21, 22], neuromuscular projection fiber tracing methods [23–25] and active contourbased tracing methods [26–28]. Based on these methods, some new improved methods were proposed [29]. The DIADEM final listed five wellperformed algorithms: the modelbased method [30], geometrybased method [31], probabilistic approachbased method [32], open snakebased method [33] and cost minimization treesbased approach [34]. In the modelbased method, Myers’s team employed the idea of shortest paths to refine local tracing, which is based on the model of AlKofahi [19] and a formal tube model. This pipeline can reconstruct the neuron from raw or preprocessed images [30]. In the geometrybased method, Erdogmus’s team introduced a principal curve to represent the skeleton of axons, and they then extracted the topology information using a recursive principal curve tracing method [31]. In the probabilistic approachbased method, Gonzalez’s team built a set of candidate trees to choose the best one by a global objective function, which combined geometric priors from image evidence [32]. In the open snakebased method, Roysam’s team proposed a three dimensional open curve snake model that was initiated automatically by a set of skeletons from binary images generated by the 2D graph cut presegmentation method, and the snake curve could be stretched bidirectionally along the centerline to trace the neuron cell structure [33]. Stepanyants’s team proposed treesbased method, which can merge individual branches into trees based on a cost minimization strategy [34]. After the DIADEM final, Liu’s group proposed a 3D neuronal morphology reconstruction method based on the augmented ray burst sampling method [35]. This method consisted of a single step to achieve the tracing and reconstruction, in which the centerline extraction or the extra radius estimation was unnecessary but the first seed must be set artificially. Peng’s team proposed series of efficient methods for neuron reconstruction, such as an anisotropic path searching method [36], an allpath pruning method [37], a hierarchicalpath pruning method based on a grayweighted image distancetree [38], an automatic distancefield neuron tracing method based on global threshold foreground extraction [39], a smart tracing method based on machine learning [40] and a method based on reverse mapping and assembling of 2D projections [41]. These methods can work well with the neuron center lines tracing under the complex and noisy background. Kakadiaris’s team proposed a learning 3D tubular modelsbased method, which can use a morphologyguided deformable model to extract the dendritic centerline and use minimum shapecost tree to represent the neuron morphology [42]. In addition, to achieve more accurate neuron tracing results, some open source softwares have been developed, such as flNeuronTool [35], FarSight [33], V3D [10], and Vaa3D [43, 44]. Along with all the existing algorithms, these open source softwares also promote the development of neuron reconstruction.
Despite the large number of proposed neuron tracing algorithms mentioned above, few methods can automatically reconstruct the complete and detailed neuron morphology, including complex dendritic and axonal arbors and variable thickness information. Moreover, because of the limited computer power, the automatic and accurate reconstruction of neuron anatomy structure is still a significant challenge.
In this paper, we propose a new 3D seed detection method based on Sliding Volume Filter (SVF) to initialize our framework, and we designed an open curve snake model combined with a SVF external force for centerline extraction and tracing. This open curve snake model has higher efficiency in the convergence of endpoints and detection of branch collision. In addition, radius estimation is another critical problem in neuron reconstruction, and accurate radius estimation can benefit simulation and functional research. Hence, this paper also proposes a new radius estimation method based on a 2D sliding band to estimate the radius of a neuron. The proposed radius estimation method can fit the real edges of neuron noncircular crosssections better than previous methods. Finally, a surface reconstruction method based on contour lines is adopted to reconstruct detailed neuron morphology.
Methods
As shown in Fig. 1, some critical steps, such as seeding, tracing, radius estimating and surface reconstruction, are included in the pipeline of our protocol. The details of every critical step will be explained.
Seed detection
Seed detection is a critical procedure in the open snakebased tracing protocol, and an ideal seed list can ensure tracing accuracy. The proposed seeding method includes the following two stages:

(1)
We used the proposed SVF based method to select coarse seeding points in the interior of neuron cells.

(2)
The ridge criterion was used to achieve the further filter to obtain better seeding points, which are always near the center of the neuron cell.
A. SVFBased seeding
In the field of computer vision and image processing, the convex region is defined as follows:

a)
A rounded convex region is a region with higher intensity in the center than the edge, and the gradient vectors of this region point to its center.

b)
A tubelike convex region is a region with higher intensity along its centerline than the edge, and the gradient vectors point to the centerline from the edge.
Quelhas’s group proposed a 2D Sliding Band Filter (SBF) for cell nucleus detection based on the characteristic of a rounded convex region [45]. In the data sets of microscopic imaging, a 3D neuron cell has not only a tubelike convex region but also a noncircular crosssection. Given these two characteristics, we extended the SBF into 3D space and designed a Sliding Volume Filter (SVF) to enhance the tubelike convex region for seed detection of neuron volume data.
To explain the calculation of SVF, we first explained the Voxel Convergence Index (VCI). As shown in Fig. 2a, O is an interested voxel in 3D volume datasets with its coordinate located at (x, y, z). A sphere support region R is located around the center O, and P are the voxels in the support region R except at O, whose coordinate is (i, j, k). ϕ(i, j, k) is the angle between PO and the gradient vector direction. The VCI of P is defined as follows:
Figure 2b and c show the calculation scheme of the sliding volume filter in a support region R whose radius is rad. To finish the discretization computation efficiently, the polar coordinate is introduced into this scheme, and the SVF is calculated as
with
where M is the number of support region lines radiating from the center pixel O(x, y, z), ρ denotes the radial coordinate, a and b stand for the angular coordinates, d is the thickness of sliding volume, r is the center position of the sliding volume in the support region line ranging from R _{ min } to R _{ max }, Q is the points between [r−d/2,r + d/2], and φ(qx _{ ρ }, qy _{ ρ }, qz _{ ρ }) is the angle between the gradient vector at Q and the direction of QO. Additionally, the angles a ∈ [0, 2π] and b ∈ [0, π] are divided into 2 L parts and L parts, respectively. Thus, M = 2 L ^{2}. Specially, the number of parts of L determines the accuracy and efficiency of computation.
After the SVF was applied to the neuron volume data for seed detection voxel by voxel, we selected the voxels as the raw seeds whose SVF response values are higher than the threshold T. Notably, there are more gradient vectors that point to the center of a tubelike structure in the marginal regions than in the other regions [45]. Hence, the sliding volumes of support regions of interior points are more likely to converge in the marginal regions. As shown in Fig. 3a, the voxel A in the interior of the nerve is more likely to be selected as a raw seed than the external voxel B. Because A has a higher SVF response value than voxel B, the orientations of gradient vectors in the sliding volumes of support region of A are more likely to point to A. However, the orientations of gradient vectors in the sliding volumes of the support region of B are not consistent and sometimes point away from B. Moreover, a nerve cell is not a uniform tubelike structure but instead has variable thickness. Therefore, SVF is the proper filter for raw seed selection.
B. Ridge criterion
The Aylward’s ridge criterion method was applied to the raw seeds for the final seed choice [20]. As shown in Fig. 3a, I is the volume data set, ∇I(p) is the gradient vector at voxel p with its coordinate (x, y, z) in I, and ev _{1}, ev _{2} and ev _{3} are the eigenvectors computed from Hessian matrix of I. ev _{1}(p) is the principle direction along the center lines of the tubelike structure, and ev _{2}(p) and ev _{3}(p) are the other two orthogonal eigenvectors. The seeds near the center of the tubelike structcure meet the condition of Eq. 4.
The raw seeding points were further chosen according to the ridge criterion Eq. 4. As shown in Fig. 3b , after the steps of SVF and ridge criterion, the proper seeds near the center line of the nerve are chosen and included in the seed list, in which the seed points are sorted by the response values. Simultaneously, response values from SVF are used to enhance the intensity of voxels in the original data, which benefits the deformation of the open curve model in the following step. The SVF volume enhancement method is denoted as
where I _{ SVF }(p) is the intensity of point p after SVF enhancement, I(p) is the intensity of point p before SVF enhancement, and SVF(p) is the SVF value of point p.
SEFOCS Neuron tracing
Tracing the full neuron skeleton is still a challenging task in neuron science, although many methods have been proposed. In this section, a new tracing model is proposed called an SVF external force open curve snake (SEFOCS, SEFOpen Curve Snake). The open curve snake model was initially applied to automated actin filament segmentation and tracking [33, 46]. Extended the application of the open curve snake model to neuron tracing. However, the computation was tedious in the tracing framework of [33], especially in the step of branch detection. The proposed SEFOCS includes three parts: open curve deformation, curve extension, and collision detection.
A. Open curve deformation
This model is a parametric open curve model, and the total snake energy can be defined as
E _{ Total } is the total image energy combined with internal energy and external energy. This model is a traditional deformable model, which resembles previous work in [16]. The open snake model is a parametric curve, c(s) = (x(s), y(s), z(s)), s∈[0, 1], and the snake internal and external energy are defined as follows:
with
In Eq. 7, α and β are the “elasticity coefficient” and “stiffness coefficient”, respectively, in internal energy, and they can control the regularity of the curve in the process of evolution. In Eq. 8, the external energy term is used to make the snake deform near the center line of the neuron and stretch the endpoints to the tail of the neuron. ∇E _{ im } is the negative normalized Gradient Vector Flow (GVF) of the volume data enhanced by SVF, p signifies point (x(s), y(s), z(s)) on the open curve, and I _{ SVF } is the volume after SVF enhancement in this paper. Instead of the original 3D image GVF, we calculated the GVF of I _{ SVF }. The SVF can enhance the tubelike convex region to smooth the GVF. As shown in Fig. 4a, the blue arrows show examples of gradient vectors from the volume enhanced by SVF. The vectors point toward the centers of neurons, which can make the seed points (the yellow points in Fig. 4a) move to the center position (the position of the red points in Fig. 4a). Specifically, the stretching force ∇E _{ str }(c(s)) is only implemented to the final endpoints c (0) and c (1). The c _{s}(s)/c _{s}(s) denotes the direction of the stretching force. The value \( \left\Vert \sqrt{\lambda_2{\lambda}_3}{\scriptscriptstyle \raisebox{1ex}{$$}\!\left/ \!\raisebox{1ex}{$$}\right.}\left({\lambda}_1+0.01\right)1\right\Vert \) is used to measure the tubelike level around the end point. When a curve reaches the end of a neuron, the end points will lose the tubelike characteristic. Hence, ∇E _{ str }(c(s)) approaches zero, and the open active curve converges. According to a large number of experiments, this strategy is not only efficient and reliable but also can avoid the leakage of the neuron boundary. To minimize the energy function E _{ Total }, the points on the snake curves are updated as:
where the parameters t and γ control the iteration numbers and size of the step at each iteration, respectively. The iterations are stopped when t reaches the threshold of the max iteration number.
B. Curve extension
The initial open snake curve is formed by three points (fewer than three points will not be traced as a branch of neuron). The first point p has the best response value in the seed list, and the other two points are generated by extending along the first principal direction to ev _{1}(p) and −ev _{1}(p). As shown in Fig. 4b, along with the open snake curve moving to the center of neuron, it also extends toward the two inverse tangential directions, c _{s}(p _{0}) and −c _{s}(p _{1}), in which the p _{0} and p _{1} are the two temporary endpoints. During the procedure of extension, the seed points belonging to one curve were labelled with new values (the default value is zero) in accord with the ID of the curve. For example, in Fig. 4b the yellow points and green points belong to different curves.
C. Collision detection
Neurons have many branches, especially in the dendrite region. Hence, detecting branching points and handling collision are essential. In the proposed scheme, two types of collision exist in the collision region and are shown in Fig. 4b. The first collision is branching point collision, which occurs when the open snakes reach a seed point whose value is not zero, and this point is recorded as the branching point (pink point in Fig. 4b). This branching point detection strategy is based on labelling seeds and is highly efficient. It also can handle the second type of collision, contour lines collision. The contour lines coming from the following step of radius estimation are the foundation of neuroanatomy reconstruction. However, due to the ambiguity of radius estimation in the collision region, the contour lines from two curves easily intersect. In Fig. 4b, this situation is illustrated in the imaginary pink circle and the embedded image, which is an experimental result in the branching region. This collision will influence the accuracy of the following reconstruction algorithm. Hence, a backoff strategy is proposed to avoid contour line collision. First, radius estimation in the branching points will not be executed. Second, if an extending curve reaches the branching point, it will be cut back the length of D, which is usually set as double the average estimated radius of the current curve.
In other words, the imaginary pink circle is not necessary in radius estimation because the following reconstruction algorithm would interpolate the information using triangular meshes automatically. Finally, the tracing algorithm ends when all the seed points are traversed.
The entire tracing algorithm procedure is shown as follows:
In summary, compared to the open snake method in [33], we improved this model in the following three aspects. First, the volume after SVF enhancement has more straightforward gradient vectors, which point to the center line of the neuron and can be used in driving the initial lines to the center of the neuron. Second, the proposed method can cut down the computation of the stretching force of end nodes. Third, in the step of collision detection, compared to the method based on labelling neighbor voxels, the method based on labelling seeds has higher detection efficiency and benefits the following reconstruction procedure.
Radius estimation
Radius estimation is another critical task in neuron anatomy reconstruction, and it can provide more quantitative information for neuroscience research. Peng, Aylward, and Wang had proposed some radius estimation methods [16, 20, 33], but most of them are based on the assumption that the neuron have a uniform tubelike structure, whose crosssections are regular circles. However, the real crosssections are not regular circles, as shown in the embedded image of Fig. 5. To reconstruct the neuron morphological structure more accurately, fitting the real edge of the neuron cell is achieved by a new proposed radius estimation method based on a 2D Sliding Band Filter (SBF) [45]. The SBF can converge on the real edge of a neuron crosssection that has the rounded convex region in [45].
Figure 5 shows the scheme of the radius estimation method based on a 2D sliding band filter. We could obtain the crosssection according to the normal vector v _{ 1 }, which points to the tangential direction of the open curve. Additionally, the v _{ 2 } and v _{ 3 } are the orthogonals in the crosssection. To obtain an accurate estimation of neuron crosssection, n radiuses radiating from the center point S on the snake curve are estimated as different lengths. The radius lengths are equal to r in the Eq. 10 with the maximum SBF response value.
with,
where B is the boundary points on the crosssection, which are at the centers of sliding bands and will be used to fit the real edge. (x _{ r } ^{n} , y _{ r } ^{n} , z _{ r } ^{n} ) are the spatial coordinates of S. The computation method of SCI in point P, which is in the range of [r−d/2, r + d/2], has been introduced in Eq. 4. d is the width of the sliding band, r is the distance between B and the center point S, and it can slide in the range of [R _{ min } , R _{ max }] to obtain the optimal position of B with the maximum SBF response value. The boundary points B can be connected clockwise to fit the edge of the neuron cell.
In the proposed method, the parameter n is related to the accuracy of radius estimation. As shown in the embedded image of Fig. 5, the larger n is, the more accurate the crosssection fitting will be. However, considering the efficiency and accuracy in the actual application, the parameter n should be adjusted flexibly.
Neuron surface reconstruction
Most of the traditional neuron reconstruction methods were based on the fast marching method and some supplemental processes for connecting different fragments [33, 47]. However, in this paper, Liu’s nonparallel contour lines surface reconstruction method is employed for surface reconstruction [48], considering the nonparallel characteristic of circles generated from previous steps. On the premise of an accurate description of the entire neuron anatomy structure, this method is efficient. Although this method had been widely used in other biological models, it has rarely been used in neuron model reconstruction. The generated mesh model of the neuron can benefit the future finite element mesh subdivision and simulation.
Figure 6 shows the scheme of Liu’s method. First, it constructs medial axes (MA) between adjacent contour lines (Fig. 6a). Second, the points and lines from different contours are projected on the MA (Fig. 6b). Third, triangular meshes are used to connect the curve networks to their projection points on the MA (Fig. 6c). Finally, the surface meshes, which are connected with different contour lines, are formed as the boundaries between neighboring compartments [48]. To obtain a smoother neuron surface, we use a surface diffusion smoothing algorithm to minimize the curvature of the model surface to obtain a smooth 3D model [49]. As shown in Fig. 6, the initial neuron surface (Fig. 6e) is formed by contour lines (Fig. 6d) which are obtained by radius estimation, and the final smooth neuron surface is shown in Fig. 6f. In addition, the most outstanding advantage of Liu’s method is that it can automatically handle branch reconstruction, especially of circles without intersections in the branching region (the intersection problem was resolved through removing the collisions of circles in the SEFOCS neuron tracing step). As shown in Fig. 7, in the branching region, two branches could be automatically reconstructed with different label colors.
Results and discussion
Parameters
We validated and evaluated the steps of seeding, tracing, radius estimation and neuron reconstruction in the proposed method using synthetic data and real data from the DIADEM challenge [12] and parts of datasets from BigNeuron project [13, 14]. All of the experiments were performed on an ordinary computer (Intel Core i5 3.2 CPU, NVIDIA GeForce GTX 960, 8 GB RAM, Windows 7). The proposed algorithm was developed using C++ language. In addition, to compare the other methods equally, we did not adopt any manual interactive operations shown in the Fig. 1, such as preprocessing, picking and expending seeds, checking and validating data, tracing editing, branch refining, and rooting setting, although these operation can improve the result of neuron tracing.
According to the image data sets (DIADEM challenge and BigNeuron project) chosen in this paper, Table 1 shows the parameters selected for the following experiments. Some experimental parameters such as d, T, L, γ, ɑ, and β remain constant for the following experiments. Other parameters such as rad, R _{ max }, R _{ min } , n and t could be adjusted for optimal results.
Seeding
Generally speaking, an ideal seed point is located in the neuron body of the foreground, and its position is near the center of the neuron cell’s crosssection. To quantify the performance of SVF seed detection, we evaluated the seeding method using an artificial helix dataset and the real dataset of the DIADEM challenge [12], according to following point deviation measurement principle:
where P_{g} denotes point sets in the gold standard, P_{s} denotes point sets generated by tested methods and N is the number of points in P_{s}. d _{ min } is the distance between a seeding point and its nearest point in the gold standard.
We compared the proposed seeding methods with the two most widely used seeding point detection methods, the global threshold method [18] and the LoG threshold method [50]. We set the parameters rad = 20, R _{ max= }16, and R _{ min= }5 in this experiments.
In the proposed SVF filter method, compared with two other methods, most of the seed points are in the interior of the neuron body. Table 2 shows the seed deviation results of the artificial helix body and OP_1, in which the SVF seeding method can achieve the lowest deviation.
As shown in Fig. 8, some seeds fall outside of the artificial helix body, which are detected by traditional threshold methods and highlighted with arrows. Figure 9 shows the seeding results in drosophila olfactory axonal datasets (OP_1 of DIADEM challenge), which are generated by the three mentioned methods. These results also suggest that our method is better than the two other seed detection methods.
We also compared the enhancement results from the three methods, and we chose the same crosssection of the neuron volume image to demonstrate that the SVF seeding method can enhance the region around the center line and simultaneously save the contrast information of the tubelike volume. The results are shown in Fig. 10, in which the red part has a higher response value than the blue part. This result suggests that the SVF method can extract the center region better than the other two methods and can enhance the original volume data properly. Furthermore, the better results in both seed detection and SVF volume enhancement can benefit neuron tracing.
Skeleton tracing
A. Tracing accuracy
We adopted the drosophila olfactory axonal datasets (OP_1 to OP_9 from DIADEM challenge) and neocortical layer 1 axons subset 1 datasets (NC_1 to NC_6 from DIADEM challenge) to evaluate the performance of the proposed neuron skeleton tracing method in the term of accuracy. Meanwhile, we compared the SEFOCS tracing method with some startoftheart algorithms, such as the Open Curve Snake tracing method (OCS) [33], Neural Circuit Tracer method (NCT) [34], allpath pruning method (APP) [37], improved allpath pruning method (APP2) [38], distancefield based method (DF) [39], and 3D tubular models based method (TM) [42].
To compare with these methods fairly, we conducted the experiments without any manual interactions and corrections, using the widely used accuracy principle to measure the test results. Similarly, we set the parameters rad = 20, R _{ max= }16, R _{ min= }5, and t = 10 in the proposed method, and we chose the better parameters for other six methods according to the features of different datasets. The measured principle is defined as:
where Precision is measured as the proportion of the length of correct traces to the total length of the traces generated by the tested methods, and Recall is the proportion of the length of correct traces to the length of the gold standard of adopted datasets.
Figure 11a and b show the reconstruction results of OP_1, and Fig. 11c and d show the results of OP_4. All of these results were generated by our proposed SEFOpen curve snake method. To illustrate the higher performance of our proposed method, we choose different colors to indicate the differences; the blue lines are the gold standard, and the green lines are the skeleton reconstructed by our method. Additionally, more tracing results of the other datasets are shown in the Additional file 1.
Table 3 summarizes the comparisons between the OCS and other six methods in terms of precision and recall. We can see that SEFOCS is far better than other six methods in most datasets in terms of accuracy and recall. In addition, the SEFOCS outperforms other six methods in average accuracy and recall. We also conducted the DIADEM score test [51] to evaluate the proposed method in the precision of reconstructed neuron topology and compare with the other methods. To the best of our knowledge, due to the various features of different datasets, no methods can get higher DIADEM score in all the datasets automatically. Hence a box plot is adopted to show the DIADEM score distribution of different methods tested in the different datasets. As we can see from Fig. 12, our method can achieve higher DIADEM score in most of the tested datasets and outperforms other six methods in average value (0.87 ± 0.001), median (0.86) and minimum (0.81). This results also proved that the proposed method has higher stability. In order to evaluate the automaticity of the proposed method, we used fixed parameters to test our method in this paper. Actually, The changed parameters can also be tried to get more meaningful tracing results. For instance, when the neuron data includes a big cell body, the bigger rad parameter is needed. Additionally, some other methods also can be tried to trace neuron according to the features of different neuron cells. For example, the APP, APP2 and DF methods can achieve better results sometimes.
B. Tracing robustness
To verify the robustness of our method, we designed three kinds of experiments. Firstly, the datasets with different levels of signal attenuation were tested. Secondly, the datasets with deferent levels of Gaussian noise were tested. Thirdly, the datasets (checked6_frog_scrippts, checked6_human_culturedcell_Cambridge_in_vitro_confocal_GFP,checked6_human_allen_confocal and checked6_fruitfly_larvae_gmu) from BigNeuron project were also tested using the proposed method.
Firstly, we compared the length of the traced skeleton with OCS in handling image signal reduction. Compared with the traditional robustness test method, which always added Gaussian noise to the original volume, this paper’s test method has a special meaning. Unbalanced light will lead to different levels of signal attenuation in the process of capturing images from a microscope. The OP_1 data set is chosen as an example, and we reduced the image signal from 10 % to 40 %. The content of the neuron images with different degrees of reduction is shown in Fig. 13. In Fig. 14, we list the lengths of the neuron skeleton traced by the two methods for comparison. With the image information reducing from 10 % to 40 %, our method can trace more information than the open curve snake in skeleton length. This result conveys that our method has higher robustness upon image signal reduction. All these results suggest that the proposed method performs better than the OCS method.
Secondly, we tested the robustness of our method using the datasets with different levels of Gaussian noise (The mean is 0 and the variances are 0.01, 0.02, 0.03 and 0.04 respectively.). As we can see from Fig. 15, the blue lines represent the tracing results of the proposed method. The tracing results are tolerable even when the variance is 0.04 and the major branches of neurons are not missed.
Thirdly, the datasets from the BigNeuron project were also tested. To my best knowledge, the BigNeuron will be a new trend in this field and most of datasets are challenging. The Fig. 16 shows some tracing results of the datasets of the BigNeuron project. Similarly, the blue lines represent the tracing results of our method. The results are tolerable even these datasets are complex and sometimes include a big cell body (In the Fig. 16, the big cell body is highlighted using the red rectangle). However, rad parameter must be set bigger (we set the parameters rad = 35, Rmax = 31, Rmin = 5, and t = 10 in the proposed method) to get better results when the datasets contain a big cell body. Additionally, the APP2 from Vaa3D [43, 44] can also achieve good results automatically when a big cell body exists in the datasets.
Radius estimation and surface reconstruction
Adaptive radius estimation and surface reconstruction methods can improve the neuron model, which has branches of varying widths. In radius estimation, the proposed method can fit the edge of the crosssection of the neuron cell. Furthermore, the credibility of radius estimation can be adjusted by the parameter n. As we can see from Fig. 5, the higher parameter n, the greater the credibility of radius estimation and the higher the computation intensity. Though a higher credibility of radius estimation must be achieved by higher computation intensity, the proposed method has solved the noncircular radius estimation problem mentioned in [33]. Generally, the parameter n = 16 is sufficient for most applications. Hence, we set n = 16 in the efficiency comparison experiments. Figure 17b shows the radius estimation results of the entire neuron using the proposed method based on the skeleton of the original volume, which is shown in Fig. 17a.
The adopted reconstruction method can interpolate meshes based on contour lines from radius estimation and can also handle branching reconstruction problems efficiently. To illustrate the performance of the proposed framework, the reconstruction result of the most complex data (OP_1) in the OP series data sets is shown in Fig. 17. The morphology of a reconstructed neuron cell of OP_1 was obtained, and the different branches were labelled with different colors during reconstruction.
Computational efficiency
In the term of automatic computation efficiency, we tested every step of the proposed pipeline and compared with other methods to evaluate the 3D neuron reconstruction efficiency. In addition, we set the parameters rad = 20, Rmax = 16, Rmin = 5, and t = 10 in the proposed method.
As we can see from Table 4, the proposed method is more efficient than the OCS framework, especially in the seeding step because seeding in the OCS framework is based on a complex procedure including graphcut segmentation and skeleton and seeding point selection. In contrast, our seeding method is more concise and efficient. The higher tracing efficiency also demonstrated the improvements in stretching force in the open snake model and the collision point detection strategy. Additionally, Table 4 also shows that the radius estimation and reconstruction are more efficient in the SEFOCS framework than in the OCS framework. These results prove that the proposed radius estimation and surface reconstruction methods outperform the corresponding methods in the OCS framework.
To test the ability of parallel computation, we also developed our CUDA implementation for the four main steps. The computation time is shown in the Table 4. We can see that the proposed method is faster than OCS method in the same parallel environment. In addition, the average speedup ratio of SEFOCS can achieve 63.94, which is higher than OCS’s average speedup ratio 48.8. This also demonstrates that the proposed method has higher parallel computation ability.
Table 5 shows the efficiency comparison results with some other methods from Table 3. Actually, most of neuron reconstruction methods don’t contain the surface reconstruction procedure. Hence, we conducted the comparison experiments in another way. We cut down the computational cost of the surface reconstruction step in order to carry out the comparison among different methods fairly. What’s more, the experiments are conducted three times for every method corresponding to every data set to avoid the errors from the operation system environment. The average results of three times are shown in Table 5. The comparison results show that the proposed method achieve the lowest average computational cost. We also can see that the computational cost of our method mainly depend on the size of the data sets and will realize higher efficiency with the development of computation parallel capacity.
Conclusions
Neuron cell anatomy structure reconstruction plays a very important role in the field of neurology. In this paper, we have developed a new neuron tracing framework, which is based on a sliding filter. We improved every step of the traditional framework compared to the OCS framework. First, given a noncircular crosssection of a neuron, the sliding filter method was introduced to the proposed seeding method (SVF) and radius estimation method (SBF), which is critical for accurately tracing skeletons and reconstructing real morphology. Second, on the basis of better seeding results, the traditional open curve snake model was improved by introducing a new external force to aid the curve evolution for neuron skeleton tracing and a new strategy for collision detection. Finally, a surface reconstruction method based on contour lines was used to generate whole neuron morphology.
A series of experiments have proved that the proposed framework has higher efficiency, stability and robustness in tracing accuracy. In addition, the proposed estimation method and adopted neuron reconstruction method can obtain more accurate neuron morphology, which is meaningful for future works such as simulation and analysis of neuron function in the field of neuroscience research.
Availability of supporting data
The source code can be available in the website [52]. The OP and NC datasets come from DIADEM challenge project, it can be downloaded from [53]. The checked6_frog_scrippts datasets, the checked6_human_culturedcell_Cambridge_in_vitro_confocal_GFP datasets, the checked6_human_allen_confocal datasets and the checked6_fruitfly_larvae_gmu datasets are available in the BigNeuron project whose website is [54].
Abbreviations
 VCI:

Voxel Convergence Index
 SBF:

Sliding Band Filter
 SVF:

Sliding Volume Filter
 SEFOCS:

SVF external force open curve snake
 OCS:

Open Curve Snake
 MA:

Medial axes
 NCT:

Neural Circuit Tracer method
 APP:

Allpath pruning method
 APP2:

Improved allpath pruning method
 DF:

Distancefield based method
 TM:

3D tubular models based method
 S:

Seeding
 T:

Tracing
 RE:

Radius estimation
 SR:

Surface reconstruction
References
 1.
Roysam B, Shain W, Ascoli GA. The central role of neuroinformatics in the national academy of engineering’s grandest challenge: reverse engineer the brain. Neuroinformatics. 2009;7(1):1–5.
 2.
van Pelt J, van Ooyen A, Uylings HBM. The need for integrating neuronal morphology databases and computational environments in exploring neuronal structure and function. Anat Embryol. 2001;204(4):255–65.
 3.
Peng HC. Bioimage informatics: a new area of engineering biology. Bioinformatics. 2008;24(17):1827–36.
 4.
Peng HC, Bateman A, Valencia A, Wren JD. Bioimage informatics: a new category in bioinformatics. Bioinformatics. 2012;28(8):1057–57.
 5.
Arbib MA, Bonaiuto JJ, BornkesselSchlesewsky I, Kemmerer D, MacWhinney B, Nielsen FA, et al. Action and language mechanisms in the brain: data, models and neuroinformatics. Neuroinformatics. 2014;12(1):209–25.
 6.
Peng HC, Roysam B, Ascoli GA. Automated image computing reshapes computational neuroscience. BMC Bioinformatics. 2013;14:293. doi:10.1186/1471210514293.
 7.
Glaser EM, Vanderloos H. A semiautomatic computermicroscope for the analysis of neuronal morphology. IEEE Trans Biomed Eng. 1965;12:22–31.
 8.
Fordholevinski TS, Dahlberg TA, Agranoff BW. A microcomputerbased image analyzer for quantitating neurite outgrowth. Brain Res. 1986;368(2):339–46.
 9.
Ascoli GA. Neuroinformatics grand challenges. Neuroinformatics. 2008;6(1):1–3.
 10.
Peng HC, Ruan ZC, Long FH, Simpson JH, Myers EW. V3D enables realtime 3D visualization and quantitative analysis of largescale biological image data sets. Nat Biotechnol. 2010;28(4):348–53.
 11.
Peng HC, Long FH, Myers EW. VANO: a volumeobject image annotation system. Bioinformatics. 2009;25(5):695–7.
 12.
Brown KM, Barrionuevo G, Canty AJ, De Paola V, Hirsch JA, Jefferis GSXE, et al. The DIADEM data Sets: representative light microscopy images of neuronal morphology to advance automation of digital reconstructions. Neuroinformatics. 2011;9(23):143–57.
 13.
Peng HC, Meijering E, Ascoli GA. From DIADEM to BigNeuron. Neuroinformatics. 2015;13(3):259–60.
 14.
Peng HC, Hawrylycz M, Roskams J, Hill S, Spruston N, Meijering E, et al. BigNeuron: largescale 3D neuron reconstruction from optical microscopy images. Neuron. 2015;87(2):252–6.
 15.
Meijering E, Jacob M, Sarria JCF, Steiner P, Hirling H, Unser M. Design and validation of a tool for neurite tracing and analysis in fluorescence microscopy images. Cytom Part A. 2004;58A(2):167–76.
 16.
Peng HC, Ruan ZC, Atasoy D, Sternson S. Automatic reconstruction of 3D neuron structures using a graphaugmented deformable model. Bioinformatics. 2010;26(12):i38–46.
 17.
Yuan XS, Trachtenberg JT, Potter SM, Roysam B. MDL constrained 3D grayscale skeletonization algorithm for automated extraction of dendrites and spines from fluorescence confocal images. Neuroinformatics. 2009;7(4):213–32.
 18.
Gonzalez G, Turetken E, Fleuret F, Fua P. Delineating trees in noisy 2D images and 3D imagestacks. Proc Cvpr IEEE. 2010:2799806. doi:10.1109/CVPR.2010.5540010.
 19.
AlKofahi KA, Lasek S, Szarowski DH, Pace CJ, Nagy G, Turner JN, et al. Rapid automated threedimensional tracing of neurons from confocal image stacks. IEEE T Inf Technol B. 2002;6(2):171–87.
 20.
Aylward SR, Bullitt E. Initialization, noise, singularities, and scale in height ridge traversal for tubular object centerline extraction. IEEE Trans Med Imaging. 2002;21(2):61–75.
 21.
Cohen AR, Roysam B, Turner JN. Automated tracing and volume measurements of neurons from 3D confocal fluorescence microscopy data. J MicroscOxford. 1994;173:103–14.
 22.
He W, Hamilton TA, Cohen AR, Holmes TJ, Pace C, Szarowski DH, et al. Automated threedimensional tracing of neurons in confocal and brightfield images. Microsc Microanal. 2003;9(4):296–310.
 23.
Cai HM, Xu XY, Lu J, Lichtman J, Yung SP, Wong STC. Using nonlinear diffusion and mean shift to detect and connect crosssections of axons in 3D optical microscopy images. Med Image Anal. 2008;12(6):666–75.
 24.
Lu J, Fiala JC, Lichtman JW. Semiautomated reconstruction of neural processes from large numbers of fluorescence images. Plos One. 2009;4(5):e5655. doi:10.1371/journal.pone.0005655.
 25.
Srinivasan R, Li Q, Zhou XB, Lu J, Lichtman J, Wong STC. Reconstruction of the neuromuscular junction connectome. Bioinformatics. 2010;26(12):i64–70.
 26.
Schmitt S, Evers JF, Duch C, Scholz M, Obermayer K. New methods for the computerassisted 3D reconstruction of neurons from confocal image stacks. Neuroimage. 2004;23(4):1283–98.
 27.
Cai HM, Xu XY, Lu J, Lichtman JW, Yung SP, Wong STC. Repulsive force based snake model to segment and track neuronal axons in 3D microscopy image stacks. Neuroimage. 2006;32(4):1608–20.
 28.
Vasilkoski Z, Stepanyants A. Detection of the optimal neuron traces in confocal microscopy images. J Neurosci Meth. 2009;178(1):197–204.
 29.
Halavi M, Hamilton KA, Parekh R, Ascoli GA. Digital reconstructions of neuronal morphology: three decades of research trends. Front Neurosci. 2012;6:49. doi:10.3389/fnins.2012.00049.
 30.
Zhao T, Xie J, Amat F, Clack N, Ahammad P, Peng HC, et al. Automated reconstruction of neuronal morphology based on local geometrical and global structural models. Neuroinformatics. 2011;9(23):247–61.
 31.
Bas E, Erdogmus D. Principal curves as skeletons of tubular objects. Neuroinformatics. 2011;9(23):181–91.
 32.
Turetken E, Gonzalez G, Blum C, Fua P. Automated reconstruction of dendritic and axonal trees by global optimization with geometric priors. Neuroinformatics. 2011;9(23):279–302.
 33.
Wang Y, Narayanaswamy A, Tsai CL, Roysam B. A broadly applicable 3D neuron tracing method based on opencurve snake. Neuroinformatics. 2011;9(23):193–217.
 34.
Chothani P, Mehta V, Stepanyants A. Automated tracing of neurites from light microscopy stacks of images. Neuroinformatics. 2011;9(23):263–78.
 35.
Ming X, Li AA, Wu JP, Yan C, Ding WX, Gong H et al. Rapid reconstruction of 3D neuronal morphology from light microscopy images with augmented rayburst sampling. Plos One. 2013;8(12):e84557. doi:10.1371/journal.pone.0084557.
 36.
Xie J, Zhao T, Lee T, Myers E, Peng HC. Anisotropic path searching for automatic neuron reconstruction. Med Image Anal. 2011;15(5):680–9.
 37.
Peng HC, Long FH, Myers G. Automatic 3D neuron tracing using allpath pruning. Bioinformatics. 2011;27(13):I239–47.
 38.
Xiao H, Peng HC. APP2: automatic tracing of 3D neuron morphology based on hierarchical pruning of a grayweighted image distancetree. Bioinformatics. 2013;29(11):1448–54.
 39.
Yang JZ, GonzalezBellido PT, Peng HC. A distancefield based automatic neuron tracing method. BMC Bioinformatics. 2013;14(1):93.
 40.
Chen HB, Xiao H, Liu TM, Peng HC. SmartTracing: selflearningbased neuron reconstruction. Brain Inform. 2015;1(1):1.
 41.
Zhou Z, Liu XX, Long B, Peng HC. TReMAP: Automatic 3D neuron reconstruction based on tracing, reverse mapping and assembling of 2D projections. Neuroinformatics. 2015. doi:10.1007/s1202101592781.
 42.
SantamariaPang A, HernandezHerrera P, Papadakis M, Saggau P, Kakadiaris I. Automatic morphological reconstruction of neurons from multiphoton and confocal microscopy images using 3D tubular models. Neuroinformatics. 2015;13(3):297–320.
 43.
Peng HC, Bria A, Zhou Z, Iannello G, Long FH. Extensible visualization and analysis for multidimensional images using Vaa3D. Nat Protoc. 2014;9(1):193–208.
 44.
Peng HC, Tang JY, Xiao H, Bria A, Zhou JL, Butler V et al. Virtual finger boosts threedimensional imaging and microsurgery as well as terabyte volume image visualization and analysis. Nat Commun. 2014;5:4342. doi:10.1038/ncomms5342.
 45.
Quelhas P, Marcuzzo M, Mendonca AM, Campilho A. Cell nuclei and cytoplasm joint segmentation using the sliding band filter. IEEE Trans Med Imaging. 2010;29(8):1463–73.
 46.
Li HS, Shen T, Smith MB, Fujiwara I, Vavylonis D, Huang XL. Automated actin filament segmentation, tracking and tip elongation measurements based on open active contour models. 2009 IEEE Int Symp Biomed Imaging. 2009;1 And 2:1302–5.
 47.
Sethian JA. A review of level set and fast marching methods for image processing. Nato Sci Ser Ii Math. 2002;75:365–96.
 48.
Liu L, Bajaj C, Deasy JO, Low DA, Ju T. Surface reconstruction from nonparallel curve networks. Comput Graph Forum. 2008;27(2):155–63.
 49.
Xu GL, Pan Q, Bajaj CL. Discrete surface modelling using partial differential equations. Comput Aided Geom D. 2006;23(2):125–45.
 50.
Byun JY, Verardo MR, Sumengen B, Lewis GP, Manjunath BS, Fisher SK. Automated tool for the detection of cell nuclei in digital microscopic images: application to retinal images. Mol Vis. 2006;12(10507):949–60.
 51.
Gillette TA, Brown KM, Ascoli GA. The DIADEM metric: comparing multiple reconstructions of the Same Neuron. Neuroinformatics. 2011;9(23):233–45.
 52.
Neuron reconstruction project. http://biometrics.hit.edu.cn/projects/neuronreconstruction. Accessed 10 Oct 2015.
 53.
DIADEM challenge datasets. http://diademchallenge.org/data_sets.html. Accessed 1 Feb 2012.
 54.
BigNeuron project datasets. http://alleninstitute.org/bigneuron/data/. Accessed 1 Jan 2015.
Acknowledgements
This work was supported by the Incheon National University International Cooperative Research Grant in 2013. Our deepest gratitude goes to the anonymous reviewers for their careful work and constructive suggestions that have helped us to improve this paper substantially. We also appreciate the support of softwares and datasets from the team of the BigNeuron project.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
GL and DS carried out the program of the neuron tracing and reconstruction and finish the work of writing the draft, KW designed the workflow of neuron cell tracing and reconstruction pipeline. JC provided the suggestion of manuscript revision and pipeline optimization. All authors read and approved the final manuscript.
Additional file
Additional file 1:
This document includes additional figures not included in the paper. Some other tracing results are shown in the supplementary material (12 figures). S1S7 are some tracing results of BigNeuron datasets. S8S9 are some other tracing results of NC datasets. S10S12 are some other tracing results of OP datasets. (PDF 1439 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Neuron anatomy structure reconstruction
 Radius estimation
 Sliding filter
 Open curve snake model