M-AMST: an automatic 3D neuron tracing method based on mean shift and adapted minimum spanning tree
© The Author(s). 2017
Received: 7 June 2016
Accepted: 11 March 2017
Published: 29 March 2017
Understanding the working mechanism of the brain is one of the grandest challenges for modern science. Toward this end, the BigNeuron project was launched to gather a worldwide community to establish a big data resource and a set of the state-of-the-art of single neuron reconstruction algorithms. Many groups contributed their own algorithms for the project, including our mean shift and minimum spanning tree (M-MST). Although M-MST is intuitive and easy to implement, the MST just considers spatial information of single neuron and ignores the shape information, which might lead to less precise connections between some neuron segments. In this paper, we propose an improved algorithm, namely M-AMST, in which a rotating sphere model based on coordinate transformation is used to improve the weight calculation method in M-MST.
Two experiments are designed to illustrate the effect of adapted minimum spanning tree algorithm and the adoptability of M-AMST in reconstructing variety of neuron image datasets respectively. In the experiment 1, taking the reconstruction of APP2 as reference, we produce the four difference scores (entire structure average (ESA), different structure average (DSA), percentage of different structure (PDS) and max distance of neurons’ nodes (MDNN)) by comparing the neuron reconstruction of the APP2 and the other 5 competing algorithm. The result shows that M-AMST gets lower difference scores than M-MST in ESA, PDS and MDNN. Meanwhile, M-AMST is better than N-MST in ESA and MDNN. It indicates that utilizing the adapted minimum spanning tree algorithm which took the shape information of neuron into account can achieve better neuron reconstructions. In the experiment 2, 7 neuron image datasets are reconstructed and the four difference scores are calculated by comparing the gold standard reconstruction and the reconstructions produced by 6 competing algorithms. Comparing the four difference scores of M-AMST and the other 5 algorithm, we can conclude that M-AMST is able to achieve the best difference score in 3 datasets and get the second-best difference score in the other 2 datasets.
We develop a pathway extraction method using a rotating sphere model based on coordinate transformation to improve the weight calculation approach in MST. The experimental results show that M-AMST utilizes the adapted minimum spanning tree algorithm which takes the shape information of neuron into account can achieve better neuron reconstructions. Moreover, M-AMST is able to get good neuron reconstruction in variety of image datasets.
KeywordsM-AMST Neuron reconstruction Mean shift Sphere model Coordinate transformation
Understanding how the brain works from the aspect of cognition and structure is one of the greatest challenges for modern science . On one hand, it is meaningful to systematically investigate human information processing mechanisms from both macro and micro points of views by cooperatively using experimental, computational, cognitive neuroscience. On the other hand, acquiring knowledge of the neuron’ morphological structure is also of particular importance to simulate the electrophysiological behavior which intricately links with cognitive function and promotes our understanding of brain. Based on the above views, several research journals (Brain Informatics , Bioimage Informatics  and Neuroinformatics ), worldwide neuron reconstruction contest (DIADEM ) and bench testing project (BigNeuron ) have been launched. One of their basic tasks is how to extract the neuronal morphology from the molecular and cellular microscopic images, namely neuron reconstruction or neuron tracing.
In order to get the neuron tracing algorithms with high performance as many as possible, the BigNeuron project aims at gathering a worldwide community to define and advance the state-of-the-art of single neuron reconstruction. The primary method to achieve that goal is to bench test as many varieties of automated neuron reconstruction methods as possible against as many neuron datasets as possible following standardized data protocols . So far, varieties of neuron reconstruction methods based on image segmentation theories such as fuzzy set , level set [8, 9], active contour model [10–12], graph theory , and clustering [14, 15] have been contributed to the project. For example, APP2 algorithm based on level set theory can generate reliable tree morphology of neuron with the fastest tracing speed . A neuron tracing algorithm named Micro-Optical Sectioning Tomography ray-shooting (MOST for short) achieves a good result in terabytes 3D datasets of the whole mouse brain . Additionally, a neuron tracing algorithm named SIMPLE is a DT-based method and can produce better reconstruction in dragonfly thoracic ganglia neuron images than other methods . A neuron tracing method based on graph theory, namely neuron tracing minimum spanning tree (N-MST for short), also gets reasonable reconstructions for several neuron image datasets. Due to the spatial nature of image, the methods mentioned above are all take the spatial information into account. However, in some segmentation scenarios, the objects of interest may be reasonably characterized by an intensity distribution. In the 3D image, the more voxels distributed in an image region, the region has a higher voxel density. For such a situation, it is important to integrate intensity information into a spatial algorithm. The neuron tracing method based on clustering is the algorithm which adopts spatial information and intensity distribution of neuron simultaneously. Moreover, because the clustering algorithms are intuitive and, some of them, easy to implement, they are very popular and widely used in image segmentation. For instance, mean shift is a nonparametric density gradient estimation using a generalized kernel approach and is one of the most powerful clustering techniques. Cai et al. proposed a cross-sections of axons detection and connection method using nonlinear diffusion and mean shift . The automatic method can shift the centroids of cross-section on slice A iteratively until the sample mean convergence on slice B. They concluded the centroid on slice A and the centroid on slice B correspond to the same axon. Comaniciu et al. proposed a robust approach for the analysis of a complex multimodal feature space and to delineate arbitrarily shaped clusters in it . The basic computational module of the technique is the mean shift. They proved for discrete data the convergence of a recursive mean shift procedure to the nearest stationary point of the underlying density function and, thus, its utility in detecting the modes of the density. They also claimed that the mean shift algorithm is a density estimation-based non-parametric clustering approach that the data space can be regarded as the empirical probability density function (p.d.f) of the represented parameter. As we know, dense region in the data space corresponds to local maxima of the p.d.f, that is, to the modes of the unknown density. Once the location of a mode is determined, the cluster associated with it is delineated based on the local structure of the data space. As it happens, the neuron image generated by fluorescent probes has the characteristics of spatial distribution, intensity discretization and the portions around the neuron skeleton have a higher voxel density. Based on this, we developed a neuron tracing algorithm based on mean shift and minimum spanning tree as a contribution to the BigNeuron project in the beginning. Specifically, the algorithm can move each voxel to the local mean until automatically get the convergence region which has the local maxima of the p.d.f. Meanwhile, the voxels in the convergence regions can also be considered as the classification voxels which indicate the modes of unknown density, other voxels which shift toward and finally locate at the regions after several iterations can be marked as the same classification as the corresponding classification voxel. They also can be regarded as the voxels subordinated to the classification voxels. Based on the subordinate voxels belong to the different classification voxel, the local structure of the neuron can be delineated. In this method, not only the information of voxel density distribution of neuron image can be captured correctly, but also got the sufficient voxels belong to several modes to delineate the whole neuron topological structure. In the basis of the classification voxels and their own subordinate voxels, we can use the minimum spanning tree (MST) algorithm to reconstruct the neuron. It is worth noting that with MST, the information of spatial distribution of neuron is also adopted to get the neuron reconstruction.
However, the experimental result of M-MST shows that although the M-MST algorithm can reconstruct 120 images successfully, it generates less precise connection between some neuron segments since the MST just takes the spatial distance as the weight of edge to build the paternity between nodes. The node is defined as the voxel with the property of spatial coordinates, radius, node type and parent compartment. The situation about less precise connection between some neuron segments caused by MST can be illustrated as the following:
Node A and node B belong to the different neuron segments and have the nearest distance or the smallest edge weight in the neuron image. According to the topological structure of original neuron image, there is a gap between them. In this case, it is not suitable to use the minimum spanning tree algorithm to build paternity of the two nodes directly. If we can detect the gap between the two nodes and set their weight of edge in MST as a high value which is exceed a lot than the real spatial distance, the MST will choose other pair of nodes which have no gap between them to form the neuron segment. The pair of nodes which have no gap between them are more likely subordinate to the same neuron segment. Once the gap is detected, the shape information can be captured. Therefore, the weight calculation method considers the shape information of neuron will help for achieving more precise neuron reconstruction.
Topological structure segmentation
A. Voxel clustering using mean shift algorithm
On the whole, we follow the sequence of neuron reconstruction operations: binarization, skeletonization, rectification and graph representation . In the binarization operation, we firstly define the voxel whose intensity is less than a threshold as the dark spot and otherwise the bright spot. For each voxel, the number of the dark spots and the bright spots among the 26 surrounding voxels are calculated respectively. And then, we calculate the ratio of the number of the dark spots and the bright spots, the ratio is compared with a threshold to determine whether the voxel is a foreground or not. In the skeletonization operation, we use mean shift algorithm to extract the neuron skeleton.
Mean shift involves shifting a kernel iteratively to a region with higher density until convergence. We shift the 3D coordinate of each voxel using a Gaussian kernel described as the following:
Assume a sphere centered on voxel P and with radius r. Using X-axis as example, the x in the formula (1) can be calculated by
The new coordinate value of the sphere center in X-axis is calculated by
B. Covered node pruning
- (1)For a pair of marks, prune the covered mark according to the distance between them and their own radiuses. Figure 2 gives three covering situations of mark. The red and purple dots represent two different marks and their own radiuses are r1 and r2 respectively. We define two marks should all be kept (Fig. 2a) if the difference between their Euclidean distance D and the sum of their radiuses is greater than a threshold. Conversely, prune one of marks (Fig. 2b and c) without defining a particular pruning priority. The threshold should be set greater than 3 according to the prior knowledge which claimed that the human eye can tell the detail variation above 3 voxels.
Remark the subordinate nodes of the removing mark to the corresponding keeping mark. Specifically, due to some covered or overlapped marks are pruned, the nodes which are subordinate to the removing marks should be re-subordinate to the keeping marks. We deal with this using a two-fold QMap data structure, the keys of first and second QMap mean the order number of keeping mark and the removing mark in the original mark set respectively. The values of second QMap mean the subordinate node set belongs to the corresponding removing mark. After remarking, we get the keeping marks and their own subordinate node set.
For each node set, prune the subordinate nodes overlapped or covered by others but always keep the marks. The pruning method is consistent with the first step.
Pathway extraction using a sphere model
Based on the marks and their own subordinate node set, we use MST algorithm to reconstruct the neuron and finish the graph representation step. The main principle of MST is to select a pathway with minimum total weights to connect all vertices in a connected and undirected graph. Taking the spatial distance between each pair of nodes as the weight is the most obvious weight calculation method. In M-MST, we build an undirected graph by connecting all extracted nodes firstly and then calculate the spatial distance between each pair of nodes as the weight of edge. However, building the paternity between each pair of nodes according to their spatial distance could easily lead to the less precise connection between some neuron segments. The weight calculation method combines the shape and spatial information of neuron will help for achieving more precise neuron reconstruction. Therefore, in order to get an accurate neuronal topological reconstruction, a rotating sphere model based on coordinate transformation is proposed to improve the weight calculation method of edge in M-MST. The core idea of the rotating sphere model is to move the sphere centered on each node in the node set to progressively approach other nodes. For each pair of nodes, we define the two nodes as starting node and terminal node respectively. A pathway could be extracted between the starting node and the terminal node if there is no gap between them or one node is not far away from the other. The main principle of the rotating sphere model based on coordinate transformation is described as the following.
Neuron reconstruction adopting MST algorithm
For each rotating step, we take the sphere center chose in the previous rotating step as the parent node of the current sphere center. A node list can be got if the starting node approaches the terminal node successfully. The node list is also a list which is composed of voxels with paternity. We defined the node list as the pathway between the starting node and the terminal node. Comparing with using the Euclidean distance as the weight of edge of the MST, it is more reasonable to adopt the accumulating distance of the node list as the weight of edge. Specifically, the accumulating distance can be calculated by summing up the distance between each pair of the parent node and the child node in the list. And then, take the accumulating Euclidean distance as the weight of edge of the MST. Notable, the weight will be set to a numerical value which is far greater than the real spatial distance if the sphere rotating model fails to extract the pathway. The reason for the sphere rotating model failed can be concluded as two aspects, the one is a gap exists between the two nodes, the other is the rotating times beyond the pre-configured threshold since one node is far away from the other. The calculation method of the weight in M-AMST descripted as the following:
if pathway exist,
then W E = ∑ i = 1 n − 1 Dis(p i , p (i + 1));
else W E = M ∗ Dis(S, T),
where WE is the accumulating Euclidean distance of the pathway between the corresponding pair of nodes, pi and pi+1 mean the child node and parent node in the node list respectively, Dis indicates the Euclidean distance between the two nodes, M means a positive integer value which is greater than 10.
After that, we use MST algorithm to build a graph representation in SWC format. Specifically, the weights between all pairs of nodes can form a weight diagonal matrix, which can be acted as the input of MST algorithm. In order to reconstruct the neuron image more elaborate, based on the neuron reconstruction by MST, we fulfill the node list into the pathway between the corresponding pair of nodes.
Parameters in the implementation
Binarization. We define the voxel whose intensity which is less than a threshold as the dark spot and otherwise the bright spot. For each voxel, the number of the dark spots and the bright spots among the 26 surrounding voxels are calculated respectively. And then, we calculate the ratio of the number of the dark spots and the bright spots, the ratio is compared with a threshold to determine whether the voxel is a foreground or not. The intensity threshold and ratio threshold are set to be 30 and 0.3 respectively.
Skeletonization. For each foreground voxel, the mean shift algorithm defines a sphere with certain spatial range and takes the voxel as the center of the sphere. Then it shifts the voxel to the local mean iteratively until getting the convergence region. In order to avoid the endless loop, the number of shifting times of the foreground voxel is set to 100. The radius of the sphere is set to 5.
Rectification. We prune the nodes yielded from skeletonization step using a node pruning method. For a pair of nodes, we define two nodes should all be kept if the difference between their Euclidean distance D and the sum of their radiuses is greater than a certain voxel distance. The voxel distance is set to 3.
Graph representation. In the implementation of the rotating sphere model, the radius of sphere is set to 2 and the sphere was split into 4 quadrants. We extract the pathway between each pair of nodes and calculate the length of it by accumulating the Euclidean distance between the parent node and child node in the node list. A weight diagonal matrix can be generated to act as the input of MST algorithm. Based on the neuron reconstruction by MST, we fulfill the node list into the pathway between the corresponding pair of nodes. After that, in order to remove the unnecessary or redundant spurs in the tree structure, we adopt a hierarchical pruning method utilized in APP2. The code can be downloaded from www.github.com/Vaa3D/vaa3d_tools/tree/master/bigneuron_ported/APP2_ported.
A. Neuron reconstruction efficiency comparison between M-AMST and M-MST
B. Running time comparison between M-AMST and M-MST
Comparison of running time (seconds) of M-AMST and M-MST on few images
Running time (s)
C. Comparison with other reconstruction algorithms
We compared the 120 neuron reconstructions generated by M-AMST with another four tracing algorithms which listed as follows: M-MST, MOST, SIMPLE, N-MST. Notably, the four tracing algorithms are also developed by the member groups of BigNeuron and the corresponding code can be downloaded from www.github.com/Vaa3D/vaa3d_tools/tree/master/bigneuron_ported. Due to the APP2 tracing algorithm is the fastest tracing algorithm and is reliable in generating tree morphology of neuron among the existing methods, we select the reconstructions generated by APP2 as the reference to compare the effect of the five algorithms. Moreover, we calculate four difference scores (ESA, DSA, PDS and MDNN) of the reconstructions produced by APP2 and the five tracing algorithms. Correspondingly, the four difference scores measure the overall average spatial divergence between two reconstructions, the spatial distance between different structures in two reconstruction, and the percentage of the neuron structure that noticeably varies in independent reconstruction, as well as the maximum distance to the nearest reconstruction elements between two reconstructions. The smaller value the four difference scores get, the neuron reconstruction effect of the competing algorithm is closer to the reference algorithm. To make a fair comparison, the reported results of the competing algorithms correspond to the default parameters called by respective plugins.
D. Tracing different neuron images and comparing with the gold standard reconstruction
Several neuron segments with low intensity are not reconstructed successfully. The reason can be concluded into two points. On one hand, the foreground identification method used in the binarization step might ignore the nodes with low intensity. On the other hand, the voxels with low intensity located around the neuron skeleton will be moved in several iterations and cannot be identified as the classification marks.
After analyzing the reconstructions generated by M-AMST thoroughly, we find that M-AMST did not thoroughly solve the problem of connecting the different segments by crossing the gap. The main reason for this is the pathway extraction method failed when there is a gap between the pair of nodes. In this case, although a numerical value which is far greater than the real spatial distance is assigned to the weight of edge between them, MST might still select this pathway if there is no smaller weight can be chose.
The node pruning method probably degrades the reconstruction accuracy of M-AMST. For each pair of nodes, we define one of the pair of nodes should be deleted if the difference between their Euclidean distance D and the sum of their radiuses is less than a threshold. A higher threshold could cause the method cannot keep the sufficient nodes to reconstruct the whole topology. Conversely, a lower threshold could cause the reconstruction with redundant neuron segments due to the method keeps the excessive nodes.
In response to limitations mentioned above, the implementation details of M-AMST algorithm will be further refined. In the near future, we will keep working on developing the neuron image de-noising and tracing algorithms based on the machine learning method so as to continue making contributions to the BigNeuron project.
We develop a pathway extraction method using a rotating sphere model based on coordinate transformation to improve the weight calculation approach in MST. The corresponding experiment shows that utilizing the adapted minimum spanning tree algorithm which took the shape information of neuron into account can achieve better neuron reconstructions. Moreover, the adoptability of M-AMST in reconstructing variety of neuron images is also testified. The result indicates that comparing with the gold standard reconstruction, the neuron reconstruction generated by the M-AMST is able to achieve good difference scores in variety of neuron datasets.
Improved all-path pruning method
Mean shift and adapted minimum spanning tree
Mean shift and minimum spanning tree
Micro-Optical Sectioning Tomography ray-shooting
Neuron tracing minimum spanning tree
probability density function
The authors thank the BigNeuron project and Dr. Hanchuan Peng for providing the testing image data used in this article and many discussions.
This work is partially supported by the National Basic Research Program of China (No. 2014CB744600), National Natural Science Foundation of China (No. 61420106005), Beijing Natural Science Foundation (No. 4164080), and Beijing Outstanding Talent Training Foundation (No. 2014000020124G039). The four projects are all funded the first workshop of BigNeuron hold in Beijing. The workshop aims at gathering Chinese neuron scientists to study and share the released neuron dataset, analyze and reconstruct the neuron image and give the interpretation for neuron reconstruction. The funding for writing and publishing the manuscript is provided by the National Natural Science Foundation of China.
Availability of data and materials
The 120 confocal images of single neurons in the Drosophila brain used in this article are provided by the BigNeuron project and also can be downloaded from the flycircuit.tw database. Besides the neuron datasets, the code of APP2 and other several neuron reconstruction algorithms can be also downloaded from the Internet. The code of APP2 can be downloaded from www.github.com/Vaa3D/vaa3d_tools/tree/master/bigneuron_ported/APP2_ported. Other algorithms include M-MST, MOST, SIMPLE and N-MST are developed by the member groups of BigNeuron and the corresponding code can be downloaded from www.github.com/Vaa3D/vaa3d_tools/tree/master/bigneuron_ported.
WZ carried out the design and implementation of M-AMST algorithm and finished the work of writing the draft. HY and HM helped to generate the neuron reconstructions of 120 confocal images using 5 algorithms, they also compared the neuron reconstructions of 5 algorithms respectively and get the statistical results of four scores. WZ generated the neuron reconstruction of 7 datasets using the M-AMST, the reconstructions of other 5 algorithms were results of the bench-test finished by Dr. Hanchuan Peng’s group. YJ provided the suggestion of manuscript revision. ZN agreed to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. All authors read and approved the final manuscript.
The authors declare that they have no competing interests.
Consent for publication
Ethics approval and consent to participate
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
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.
- Meijering E. Neuron tracing in perspective. Cytometry A. 2010;77 A(7):693–704.View ArticleGoogle Scholar
- Zhong N, Yau SS, Ma J, Shimojo S, Just M, Hu B, Wang G, Oiwa K, Anzai Y. Brain Informatics-Based Big Data and the Wisdom Web of Things. IEEE Intell Syst. 2015;30(5):2–7.View ArticleGoogle Scholar
- Peng H. Bioimage informatics: a new area of engineering biology. Bioinformatics. 2008;24(17):1827–36.View ArticlePubMedPubMed CentralGoogle Scholar
- Ascoli GA. Neuroinformatics Grand Challenges. Neuroinformatics. 2008;6(1):1–3.View ArticlePubMedGoogle Scholar
- Peng H, Meijering E, Ascoli GA. From DIADEM to BigNeuron. Neuroinformatics. 2015;13(3):259–60.View ArticlePubMedPubMed CentralGoogle Scholar
- Peng H, Hawrylycz M, Roskams J, Hill S, Spruston N, Meijering E, Ascoli Giorgio A. BigNeuron: large-scale 3d neuron reconstruction from optical microscopy images. Neuron. 2015;87(2):252–6.View ArticlePubMedPubMed CentralGoogle Scholar
- Pal SK. Fuzzy sets in image processing and recognition. In: Fuzzy Systems, 1992, IEEE International Conference on: 8-12 Mar 1992 1992. 119-126.
- Malladi R, Sethian JA. Level set and fast marching methods in image processing and computer vision. In: Image Processing, 1996 Proceedings, International Conference on: 16-19 Sep 1996 1996. 489-492 vol.481.
- Xiao H, Peng H. APP2: automatic tracing of 3D neuron morphology based on hierarchical pruning of a gray-weighted image distance-tree. Bioinformatics. 2013;29(11):1448–54.View ArticlePubMedPubMed CentralGoogle Scholar
- Peng H, Long F, Myers G. Automatic 3D neuron tracing using all-path pruning. Bioinformatics. 2011;27(13):i239–47.View ArticlePubMedPubMed CentralGoogle Scholar
- Wang Y, Narayanaswamy A, Tsai C-L, Roysam B. A broadly applicable 3-D neuron tracing method based on open-curve snake. Neuroinformatics. 2011;9(2):193–217.View ArticlePubMedGoogle Scholar
- Cai H, Xu X, Lu J, Lichtman JW, Yung SP, Wong ST. Repulsive force based snake model to segment and track neuronal axons in 3D microscopy image stacks. Neuroimage. 2006;32(4):1608–20.View ArticlePubMedGoogle Scholar
- Peng H, Ruan Z, Atasoy D, Sternson S. Automatic reconstruction of 3D neuron structures using a graph-augmented deformable model. Bioinformatics. 2010;26(12):i38–46.View ArticlePubMedPubMed CentralGoogle Scholar
- Cai H, Xu X, Lu J, Lichtman J, Yung SP, Wong ST. 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.View ArticlePubMedGoogle Scholar
- Oliver A, Munoz X, Batlle J, Pacheco L, Freixenet J. Improving Clustering Algorithms for Image Segmentation using Contour and Region Information. In: 2006 IEEE International Conference on Automation, Quality and Testing, Robotics: 25-28 May 2006 2006. 315-320.
- Wu J, He Y, Yang Z, Guo C, Luo Q, Zhou W, Chen S, Li A, Xiong B, Jiang T, et al. 3D BrainCV: Simultaneous visualization and analysis of cells and capillaries in a whole mouse brain with one-micron voxel resolution. NeuroImage. 2014;87:199–208.View ArticlePubMedGoogle Scholar
- Yang J, Gonzalez-Bellido PT, Peng H. A distance-field based automatic neuron tracing method. BMC Bioinformatics. 2013;14(1):1–11.View ArticleGoogle Scholar
- Comaniciu D, Meer P. Mean shift: a robust approach toward feature space analysis. IEEE Trans Pattern Anal Mach Intell. 2002;24(5):603–19.View ArticleGoogle Scholar
- Peng H, Ruan Z, Long F, Simpson JH, Myers EW. V3D enables real-time 3D visualization and quantitative analysis of large-scale biological image data sets. Nat Biotech. 2010;28(4):348–53.View ArticleGoogle Scholar