1Research Center of Perception and Computing, School of Computer Science and Technology, Harbin Institute of Technology, Harbin, China
2Department of Computer Science and Engineering, Incheon National University, Incheon, Korea
BMC Bioinformatics 2015, 16:342 doi:10.1186/s12859-015-0780-0
The electronic version of this article is the complete one and can be found online at:http://www.biomedcentral.com/1471-2105/16/342
|Received:||14 May 2015|
|Accepted:||16 October 2015|
|Published:||24 October 2015|
© 2015 Luo et al.
Open AccessThis 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.
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.
This paper introduces a new pipeline for reconstructing neuron anatomy structure from 3-D microscopy image stacks. This pipeline is initialized with a set of seeds that were detected by our proposed Sliding Volume Filter (SVF), given a non-circular cross-section 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 cross-section of the neuron cell. Finally, a surface reconstruction method based on non-parallel curve networks is used to generate the neuron cell surface to finish this pipeline.
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.
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.
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 , . 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 bio-images can be obtained, which is helpful for us to develop new methods and algorithms to meet the needs in neuroscience research , . 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 , .
The computer-aided manual neuron reconstruction method was first proposed in 1965 and was achieved by a biologist using a microscope . Following this milestone, numerous algorithms and open softwares were introduced to reduce manual labor consisting of the boring task of tracing and analysis –, but most of them were still limited to semi-automation 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  and BigNeuron project , , which began in 2010 and 2015 respectively. They provided an open-source 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 , , minimum spanning tree methods , , sequential tracing methods , , skeletonization methods , , neuromuscular projection fiber tracing methods – and active contour-based tracing methods –. Based on these methods, some new improved methods were proposed . The DIADEM final listed five well-performed algorithms: the model-based method , geometry-based method, probabilistic approach-based method , open snake-based method  and cost minimization trees-based approach . In the model-based method, Myers’s team employed the idea of shortest paths to refine local tracing, which is based on the model of Al-Kofahi  and a formal tube model. This pipeline can reconstruct the neuron from raw or preprocessed images. In the geometry-based 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 . In the probabilistic approach-based 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 . In the open snake-based 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 2-D graph cut pre-segmentation method, and the snake curve could be stretched bi-directionally along the centerline to trace the neuron cell structure . Stepanyants’s team proposed trees-based method, which can merge individual branches into trees based on a cost minimization strategy . After the DIADEM final, Liu’s group proposed a 3D neuronal morphology reconstruction method based on the augmented ray burst sampling method . 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 , an all-path pruning method , a hierarchical-path pruning method based on a gray-weighted image distance-tree, an automatic distance-field neuron tracing method based on global threshold foreground extraction , a smart tracing method based on machine learning  and a method based on reverse mapping and assembling of 2D projections . 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 models-based method, which can use a morphology-guided deformable model to extract the dendritic centerline and use minimum shape-cost tree to represent the neuron morphology . In addition, to achieve more accurate neuron tracing results, some open source softwares have been developed, such as flNeuronTool , FarSight , V3D , and Vaa3D, . 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 non-circular cross-sections better than previous methods. Finally, a surface reconstruction method based on contour lines is adopted to reconstruct detailed neuron morphology.
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.
Fig. 1. Pipeline of neuron anatomy structure reconstruction
Seed detection is a critical procedure in the open snake-based 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. SVF-Based 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 tube-like 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 . In the data sets of microscopic imaging, a 3D neuron cell has not only a tube-like convex region but also a non-circular cross-section. Given these two characteristics, we extended the SBF into 3D space and designed a Sliding Volume Filter (SVF) to enhance the tube-like 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 regionR 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:
Fig. 2. Scheme of Spatial Convergence Index. a The model of 3D spatial convergence index. b The model of 3D sliding volume filter in y-z plate section. c the discretization calculation of SVF using the polar coordinates
Figure 2b and c show the calculation scheme of the sliding volume filter in a support region Rwhose radius is rad. To finish the discretization computation efficiently, the polar coordinate is introduced into this scheme, and the SVF is calculated as
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 minto 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 tube-like structure in the marginal regions than in the other regions . 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 voxelA 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 tube-like structure but instead has variable thickness. Therefore, SVF is the proper filter for raw seed selection.
Fig. 3. Scheme of computation of Sliding Volume Filter, as well as the selection of seed points. a The procedure of seed points filtrate, after the SVF and ridge criterion, proper seed points are chosen. b The seeding points selection after step1 and step2
B. Ridge criterion
The Aylward’s ridge criterion method was applied to the raw seeds for the final seed choice . 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 ev1 , ev2 and ev3 are the eigenvectors computed from Hessian matrix of I. ev1 (p) is the principle direction along the center lines of the tube-like structure, and ev2 (p) and ev3 (p) are the other two orthogonal eigenvectors. The seeds near the center of the tube-like 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 pbefore SVF enhancement, and SVF(p) is the SVF value of point p.
SEF-OCS 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 (SEF-OCS, SEF-Open Curve Snake). The open curve snake model was initially applied to automated actin filament segmentation and tracking , . Extended the application of the open curve snake model to neuron tracing. However, the computation was tedious in the tracing framework of , especially in the step of branch detection. The proposed SEF-OCS 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 . 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:
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 tube-like 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 endpointsc (0) and c (1). The cs (s)/||cs (s)|| denotes the direction of the stretching force. The value is used to measure the tube-like level around the end point. When a curve reaches the end of a neuron, the end points will lose the tube-like characteristic. Hence, ∇Estr (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:
Fig. 4. Scheme of SEF open curve snake model. a The open curve is driven to the center of neuron by external force in the volume after SVF.b The procedure of open snake curve extension and collision detection in the branching region
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 ev1 (p) and −ev1 (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, cs (p0 ) and −cs (p1 ), in which the p0 and p1 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 , 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 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 , , , but most of them are based on the assumption that the neuron have a uniform tube-like structure, whose cross-sections are regular circles. However, the real cross-sections 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) . The SBF can converge on the real edge of a neuron cross-section that has the rounded convex region in .
Fig. 5. Illustration of radius estimation of the neuron cross-section. The left embedded image shows the real cross-section of neuron, and the estimation result with different parameters. And in the right image v 1 is the tangential direction in S i , v 2 and v 3 are the orthogonal vectors which define the cross-section
Figure 5 shows the scheme of the radius estimation method based on a 2D sliding band filter. We could obtain the cross-section 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 cross-section. To obtain an accurate estimation of neuron cross-section, 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.
where B is the boundary points on the cross-section, which are at the centers of sliding bands and will be used to fit the real edge. (x rn , y rn , z rn ) 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. dis 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 cross-section 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 , . However, in this paper, Liu’s non-parallel contour lines surface reconstruction method is employed for surface reconstruction , considering the non-parallel 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 . 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 . 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 SEF-OCS neuron tracing step). As shown in Fig. 7, in the branching region, two branches could be automatically reconstructed with different label colors.
Fig. 6. Process of contour reconstruction. a Construction of projective plate MA . b Projection of points and lines on MA . c Triangulation of adjacent contour lines . d The contour lines of neuron cell. e The initial surface from triangulation of adjacent contour lines. f The final surface model after smoothing
Results and discussion
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  and parts of datasets from BigNeuron project , . 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.
Table 1. Parameter selection
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 cross-section. 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 , 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  and the LoG threshold method . 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.
Table 2. Comparisons of different seeding methods on test datasets
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.
Fig. 8. Comparison of seeding method on test dataset; a The seeds detection result of global threshold method. b The seeds detection result of LoG threshold method. c The seeds detection result of SVF method
Fig. 9. Comparison of seeding method on drosophila olfactory axonal data sets (OP_1). a The seeds detection result of global threshold method. b The seeds detection result of LoG threshold method. c The seeds detection result of SVF method and the enlargement image in intensive region of dendrites
We also compared the enhancement results from the three methods, and we chose the same cross-section 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 tube-like 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.
Fig. 10. Comparisons among seeding methods for image enhancement. The results of global threshold method lose the contrast information between center and edge; The results of LoG threshold method extend the center region exorbitantly; The results of SVF seeding method can enhance the original volume properly
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 SEF-OCS tracing method with some start-of-the-art algorithms, such as the Open Curve Snake tracing method (OCS) , Neural Circuit Tracer method (NCT) , all-path pruning method (APP) , improved all-path pruning method (APP2) , distance-field based method (DF) , and 3D tubular models based method (TM) .
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 SEF-Open 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.
Fig. 11. The tracing results of OP_1 and OP_4. a The tracing result of OP_1 from multi-view. b The magnified result of branch parts of OP_1. cThe tracing result of OP_4 from multi-view; d The magnified result of branch parts of OP_4
Table 3 summarizes the comparisons between the OCS and other six methods in terms of precision and recall. We can see that SEF-OCS is far better than other six methods in most datasets in terms of accuracy and recall. In addition, the SEF-OCS outperforms other six methods in average accuracy and recall. We also conducted the DIADEM score test  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 radparameter 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.
Table 3. Comparisons among different methods with different image datasets in tracing accuracy
Fig. 12. The box plot of DIADEM scores of the different methods for different datasets (OP_1 to OP_9 and NC_1 to NC_6 datasets). The median is the middle pink bar. The box indicates the lower quartile (splits 25 % of lowest data) and the upper quartile (splits 75 % of highest data). The red bar and blue bar are the maximum and minimum values. The blue diamond denotes the mean value of the scores
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.
Fig. 14. Changes of skeleton length with signal reduction (Unit: Volex)
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.
Fig. 15. The tracing results of NC_2 dataset with different levels of Gaussian noise. a The dataset with a variance of 0.01. b The dataset with a variance of 0.02. c The dataset with a variance of 0.03. d The dataset with a variance of 0.04. To prove the robustness of our method clearly, we follow the similar design of experiment in 
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 ,  can also achieve good results automatically when a big cell body exists in the datasets.
Fig. 16. The tracing results of the datasets from the BigNeuron project. aThe tracing result of the done_err_Recon112012no2-2 data of checked6_frog_scrippts. b The tracing result of image 7 data of checked6_human_culturedcell_Cambridge_in_vitro_confocal_GFP. c The tracing result of in_house1 data of checked6_human_allen_confocal. dThe tracing result of done_1_CL-I_X_OREGON_R_ddaD_membrane-GFP data of checked6_fruitfly_larvae_gmu
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 cross-section 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 non-circular radius estimation problem mentioned in . 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.
Fig. 17. Radius estimation and anatomical reconstruction of OP_1. a The original result by volume rendering method. b The result of radius estimation by 2D sliding band method, in which the blue contour lines are used to fit the edges of neuron cross-sections. c The anatomical reconstruction result based on the contour lines, in which the different branches are labelled with different colors
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.
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 graph-cut 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 SEF-OCS 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.
Table 4. Comparisons with OCS method in the efficiency of the proposed pipeline
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 SEF-OCS 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.
Table 5. Comparisons in the efficiency of neuron tracing
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 non-circular cross-section 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 . The OP and NC datasets come from DIADEM challenge project, it can be downloaded from . 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 .
VCI: Voxel Convergence Index
SBF: Sliding Band Filter
SVF: Sliding Volume Filter
SEF-OCS: SVF external force open curve snake
OCS: Open Curve Snake
MA: Medial axes
NCT: Neural Circuit Tracer method
APP: All-path pruning method
APP2: Improved all-path pruning method
DF: Distance-field based method
TM: 3D tubular models based method
RE: Radius estimation
SR: Surface reconstruction
The authors declare that they have no competing interests.
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 1:. This document includes additional figures not included in the paper. Some other tracing results are shown in the supplementary material (12 figures). S1-S7 are some tracing results of BigNeuron datasets. S8-S9 are some other tracing results of NC datasets. S10-S12 are some other tracing results of OP datasets. (PDF 1439 kb)
Format: PDF Size: 1.4MB Download file
This file can be viewed with: Adobe Acrobat Reader
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.
- 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. PubMed Abstract | Publisher Full Text
- 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. PubMed Abstract | Publisher Full Text
- Peng HC. Bioimage informatics: a new area of engineering biology. Bioinformatics. 2008; 24(17):1827-36. PubMed Abstract | Publisher Full Text
- Peng HC, Bateman A, Valencia A, Wren JD. Bioimage informatics: a new category in bioinformatics. Bioinformatics. 2012; 28(8):1057-57. PubMed Abstract |Publisher Full Text
- Arbib MA, Bonaiuto JJ, Bornkessel-Schlesewsky 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. PubMed Abstract |Publisher Full Text
- Peng HC, Roysam B, Ascoli GA. Automated image computing reshapes computational neuroscience. BMC Bioinformatics. 2013;14:293. doi:10.1186/1471-2105-14-293.
- Glaser EM, Vanderloos H. A semi-automatic computer-microscope for the analysis of neuronal morphology. IEEE Trans Biomed Eng. 1965; 12:22-31. PubMed Abstract |Publisher Full Text
- Fordholevinski TS, Dahlberg TA, Agranoff BW. A microcomputer-based image analyzer for quantitating neurite outgrowth. Brain Res. 1986; 368(2):339-46. Publisher Full Text
- Ascoli GA. Neuroinformatics grand challenges. Neuroinformatics. 2008; 6(1):1-3.PubMed Abstract | Publisher Full Text
- Peng HC, Ruan ZC, Long FH, Simpson JH, Myers EW. V3D enables real-time 3D visualization and quantitative analysis of large-scale biological image data sets. Nat Biotechnol. 2010; 28(4):348-53. PubMed Abstract | Publisher Full Text
- Peng HC, Long FH, Myers EW. VANO: a volume-object image annotation system.Bioinformatics. 2009; 25(5):695-7. PubMed Abstract | Publisher Full Text
- 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(2-3):143-57. PubMed Abstract | Publisher Full Text
- Peng HC, Meijering E, Ascoli GA. From DIADEM to BigNeuron. Neuroinformatics. 2015;13(3):259-60. PubMed Abstract | Publisher Full Text
- Peng HC, Hawrylycz M, Roskams J, Hill S, Spruston N, Meijering E et al.. BigNeuron: large-scale 3D neuron reconstruction from optical microscopy images. Neuron. 2015;87(2):252-6. PubMed Abstract | Publisher Full Text
- 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. Publisher Full Text
- Peng HC, Ruan ZC, Atasoy D, Sternson S. Automatic reconstruction of 3D neuron structures using a graph-augmented deformable model. Bioinformatics. 2010;26(12):i38-46. PubMed Abstract | Publisher Full Text
- Yuan XS, Trachtenberg JT, Potter SM, Roysam B. MDL constrained 3-D grayscale skeletonization algorithm for automated extraction of dendrites and spines from fluorescence confocal images. Neuroinformatics. 2009; 7(4):213-32. PubMed Abstract |Publisher Full Text
- Gonzalez G, Turetken E, Fleuret F, Fua P. Delineating trees in noisy 2D images and 3D image-stacks. Proc Cvpr IEEE. 2010:2799-806. doi:10.1109/CVPR.2010.5540010.
- Al-Kofahi KA, Lasek S, Szarowski DH, Pace CJ, Nagy G, Turner JN et al.. Rapid automated three-dimensional tracing of neurons from confocal image stacks. IEEE T Inf Technol B. 2002; 6(2):171-87. Publisher Full Text
- 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. PubMed Abstract | Publisher Full Text
- Cohen AR, Roysam B, Turner JN. Automated tracing and volume measurements of neurons from 3-D confocal fluorescence microscopy data. J Microsc-Oxford. 1994;173:103-14. Publisher Full Text
- He W, Hamilton TA, Cohen AR, Holmes TJ, Pace C, Szarowski DH et al.. Automated three-dimensional tracing of neurons in confocal and brightfield images. Microsc Microanal. 2003; 9(4):296-310. PubMed Abstract | Publisher Full Text
- Cai HM, Xu XY, Lu J, Lichtman J, Yung SP, Wong STC. Using nonlinear diffusion and mean shift to detect and connect cross-sections of axons in 3D optical microscopy images. Med Image Anal. 2008; 12(6):666-75. PubMed Abstract | Publisher Full Text
- Lu J, Fiala JC, Lichtman JW. Semi-automated reconstruction of neural processes from large numbers of fluorescence images. Plos One. 2009; 4(5):e5655. PubMed Abstract |Publisher Full Text
- Srinivasan R, Li Q, Zhou XB, Lu J, Lichtman J, Wong STC. Reconstruction of the neuromuscular junction connectome. Bioinformatics. 2010; 26(12):i64-70.PubMed Abstract | Publisher Full Text
- Schmitt S, Evers JF, Duch C, Scholz M, Obermayer K. New methods for the computer-assisted 3-D reconstruction of neurons from confocal image stacks. Neuroimage. 2004; 23(4):1283-98. PubMed Abstract | Publisher Full Text
- 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. PubMed Abstract | Publisher Full Text
- Vasilkoski Z, Stepanyants A. Detection of the optimal neuron traces in confocal microscopy images. J Neurosci Meth. 2009;178(1):197–204.
- Halavi M, Hamilton KA, Parekh R, Ascoli GA. Digital reconstructions of neuronal morphology: three decades of research trends. Front Neurosci. 2012; 6:49.PubMed Abstract | Publisher Full Text
- 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(2-3):247-61. PubMed Abstract | Publisher Full Text
- Bas E, Erdogmus D. Principal curves as skeletons of tubular objects. Neuroinformatics. 2011; 9(2-3):181-91. PubMed Abstract | Publisher Full Text
- Turetken E, Gonzalez G, Blum C, Fua P. Automated reconstruction of dendritic and axonal trees by global optimization with geometric priors. Neuroinformatics. 2011;9(2-3):279-302. PubMed Abstract | Publisher Full Text
- Wang Y, Narayanaswamy A, Tsai CL, Roysam B. A broadly applicable 3-D neuron tracing method based on open-curve snake. Neuroinformatics. 2011; 9(2-3):193-217.PubMed Abstract | Publisher Full Text
- Chothani P, Mehta V, Stepanyants A. Automated tracing of neurites from light microscopy stacks of images. Neuroinformatics. 2011; 9(2-3):263-78. PubMed Abstract |Publisher Full Text
- 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. PubMed Abstract | Publisher Full Text
- 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.
- Peng HC, Long FH, Myers G. Automatic 3D neuron tracing using all-path pruning.Bioinformatics. 2011; 27(13):I239-47. PubMed Abstract | Publisher Full Text
- Xiao H, Peng HC. APP2: automatic tracing of 3D neuron morphology based on hierarchical pruning of a gray-weighted image distance-tree. Bioinformatics. 2013;29(11):1448-54. PubMed Abstract | Publisher Full Text
- Yang JZ, Gonzalez-Bellido PT, Peng HC. A distance-field based automatic neuron tracing method. BMC Bioinformatics. 2013; 14(1):93. PubMed Abstract |BioMed Central Full Text
- Chen HB, Xiao H, Liu TM, Peng HC. SmartTracing: self-learning-based neuron reconstruction. Brain Inform. 2015; 1(1):1.
- 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/s12021-015-9278-1.
- Santamaria-Pang A, Hernandez-Herrera 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.PubMed Abstract | Publisher Full Text
- 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.PubMed Abstract | Publisher Full Text
- Peng HC, Tang JY, Xiao H, Bria A, Zhou JL, Butler V et al. Virtual finger boosts three-dimensional imaging and microsurgery as well as terabyte volume image visualization and analysis. Nat Commun. 2014;5:4342. doi:10.1038/ncomms5342.
- 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. PubMed Abstract | Publisher Full Text
- 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.
- Sethian JA. A review of level set and fast marching methods for image processing.Nato Sci Ser Ii Math. 2002; 75:365-96.
- Liu L, Bajaj C, Deasy JO, Low DA, Ju T. Surface reconstruction from non-parallel curve networks. Comput Graph Forum. 2008; 27(2):155-63. Publisher Full Text
- Xu GL, Pan Q, Bajaj CL. Discrete surface modelling using partial differential equations. Comput Aided Geom D. 2006; 23(2):125-45. Publisher Full Text
- 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(105-07):949-60. PubMed Abstract | Publisher Full Text
- Gillette TA, Brown KM, Ascoli GA. The DIADEM metric: comparing multiple reconstructions of the Same Neuron. Neuroinformatics. 2011; 9(2-3):233-45.PubMed Abstract | Publisher Full Text
- Neuron reconstruction project. http://biometrics.hit.edu.cn/projects/neuron-reconstruction. Accessed 10 Oct 2015.
- DIADEM challenge datasets. http://diademchallenge.org/data_sets.html. Accessed 1 Feb 2012.
- BigNeuron project datasets. http://alleninstitute.org/bigneuron/data/. Accessed 1 Jan 2015.