The 3D reconstruction of neuronal morphology is a powerful technique for investigating nervous systems. Due to the noises in optical microscopic images, the automated reconstruction of neuronal morphology has been a challenging problem. We propose a novel automatic neuron reconstruction algorithm, Rivulet, to target the challenges raised by the poor quality of the optical microscopic images. After the neuron images being de-noised with an anisotropic filter, the Rivulet algorithm combines multi-stencils fast-marching and iterative back-tracking from the geodesic farthest point on the segmented foreground. The neuron segments are dumped or merged according to a set of criteria at the end of each iteration. The proposed Rivulet tracing algorithm is evaluated with data provided from the BigNeuron Project. The experimental results demonstrate that Rivulet outperforms the compared state-of-the-art tracing methods when the images are of poor quality

Cai, Weidong

Feng, Dagan

Liu, Sidong

Liu, Siqi

Peng, Hanchuan

Zhang, Donghao

English

Sydney eScholarship

Reconstruction of 3D Neuron Morphology Using Rivulet Back-TrackingDonghao Zhang?, Siqi Liu?, Sidong Liu?, Dagan Feng?, Hanchuan Peng†, Weidong Cai?? School of Information Technologies, University of Sydney, Australia† Allen Institute for Brain Science, Seattle, WA, USAABSTRACTThe 3D reconstruction of neuronal morphology is a pow-erful technique for investigating nervous systems. Due tothe noises in optical microscopic images, the automated re-construction of neuronal morphology has been a challengingproblem. We propose a novel automatic neuron reconstruc-tion algorithm, Rivulet, to target the challenges raised by thepoor quality of the optical microscopic images. After theneuron images being de-noised with an anisotropic filter, theRivulet algorithm combines multi-stencils fast-marching anditerative back-tracking from the geodesic farthest point onthe segmented foreground. The neuron segments are dumpedor merged according to a set of criteria at the end of eachiteration. The proposed Rivulet tracing algorithm is evaluatedwith data provided from the BigNeuron Project. The ex-perimental results demonstrate that Rivulet outperforms thecompared state-of-the-art tracing methods when the imagesare of poor quality.Index Terms— Neuron morphology, curvilinear structuretracing, fast marching, vesselness filtering1. INTRODUCTIONThe digital reconstruction of neurons is an important tool inacquiring the morphological information of single neurons,which is critical to understand their functions within neuralcircuits [1]. Single neurons can be captured from 3D opti-cal microscopic images with various labelling visualizationtechniques, such as bulk dye loading, intracellular injections,immunolabeling and genetic labelling [2]. The single neu-ron morphology can be reconstructed from the microscopicimages and represented as 3D tubular models. The major-ity of the neuronal morphology datasets were collected bycomputer-aided manual tracing [2]. However, manual neu-ron reconstruction is a labour intensive task that is difficult toscale to large datasets.Automatic/semi-automatic neuron reconstruction algo-rithms have been developed in recent years to enable thelarge-scale acquisition of neuron morphological data. InThis work was supported in part by ARC grants. Siqi Liu(siqi.liu@sydney.edu.au) and Weidong Cai (tom.cai@sydney.edu.au) are thecorresponding authors of this paper.Neuron Studio [3], the tracing is initialised from a manuallychosen source point and shoots rays from a spherical model(Rayburst) at new traced locations. Neuron Tube [4] andMOST algorithm [5] improve the robustness of the Rayburstalgorithm by considering more shape specifics. Snake [6]traces from various seeding locations following a prepro-cessed gradient vector flow (GVF) map. All Path Pruningalgorithms (APP1 and APP2) [7, 8] prune the redundantbranches from a neuronal tree that results from an initialover-reconstruction based on fast-marching.Although these powerful algorithms have been proposed,the automatic reconstruction of single neurons remains abottleneck in the processing pipeline of microscopic neuronimages due to the varying image qualities. The broken andfuzzy neurite segments cause major difficulties in neurontracing. Most existing algorithms are incapable of bridgingthe gaps between the broken segments without introducingmiswiring. The noises from labelling techniques and mi-croscopies cannot be trivially removed by the median filteror Gaussian filter that are widely used in processing generalimages. These noises cause over-tracing in most seed de-pendent algorithms [3, 5, 6], and also increase the risk ofbranch miswiring in post-processing. In addition, the imper-fect shapes captured in real microscopic images frequentlyviolate the shape-based assumption relied on by many tracingalgorithms.In this study, we propose a novel tracing algorithm to tar-get the above mentioned challenges. The curvilinear struc-tures are segmented using an anisotropic filter. The multi-stencils fast-marching is applied to a distance field to producea time-crossing map. Inspired by the all path pruning algo-rithms [7, 8], the longest remaining neuron segments are dis-covered iteratively. To achieve this, the Rivulet back-trackingis performed based on the gradient of the time-crossing mapto trace from the geodesic farthest point on the segmentedforeground to the soma location. After each iteration, the dis-covered path is erased from the time-crossing map to avoidduplicated tracing. The whole Rivulet back-tracking processonly stops when a pre-defined proportion of the foregroundarea is covered; thus, it is highly robust to under-tracing. TheRivulet back-tracking is also robust to small structural gapsand holes caused by a low signal to noise ratio (SNR). Weuse the percentile of back-tracking locations traced withinforeground voxels to measure the confidence of each tracedsegment. Since the fast-marching is performed on a distancetransform map which has brighter intensities near the centre-line of the neuron segments, Rivulet is insensitive to tubularstructures with non-smooth boundaries.2. METHODS2.1. Hessian Based Neuron SegmentationThe image is firstly cropped to remove the empty backgroundregions. Since noise points and irrelevant non-curvilinearstructures would cause errors in neuron tracing, we applyanisotropic filter [9] to remove the image noise and preservethe curvilinear neuronal structures. The Hessian matrix H(x)is obtained by calculating the second-order gradient of theGaussian filtered image. The function g(x) is defined asg(x) =i=3∑i=1αikiλi (λ1 ≈ 0, λ1  λ2, λ1  λ3)0 otherwise(1)where λ1, λ2 and λ3 are the eigenvalues of H(x); αi are pre-defined as α1 = 0.5, α2 = 0.5, α3 = 25; ki is the normal-ising term and ki = exp(−λ2i /∑i=3i=1 λi). The foregroundbinary image is segmented by applying a threshold of 0 onthe filtered image f(x) = exp(−|∇x|2)g(x). Though theanisotropic filter is capable of eliminating a number of noisepoints without corrupting the neuron segments of interest, forhighly noisy images noise points may remain in the binaryimage Ib(x) ∈ {0, 1} which are considered in the Rivulettracing.2.2. Multi-Stencils Fast-MarchingThe first step of Rivulet tracing is to generate time-crossingmap T based on multi-stencils fast-marching (MSFM) bysolving the Eikonal equation Eq. (2):F =dxdT, |∇T |F = 1, T (Γ0) = 0 (2)where F is the speed evolution function of fast marching, Γ0is the starting condition [10]. MSFM solves Eq. (2) alongmultiple stencils, thus avoiding the numerical errors of con-ventional fast-marching methods along diagonal directions.The speed image F is computed as a boundary distance trans-form with sub-voxel precision that has bright intensities nearthe centre-line of the neuron segments. To compute F , theboundary of the segmented foreground map B(Ib) is firstlyobtained by subtracting the Ib by the dilated foreground Ib ⊕K, where K is a 3D dilation kernel of size 3 × 3 × 3. ThenMSFM is applied at all the points on the boundary [11]. Thevoxels within the foreground region of F represent the closestdistance to any voxels in B(Ib). All voxels in the backgroundregion of F are set to zero. Thus, the time crossing map T (x)is defined as:T (x) = minP(s,x)∫ L0C(P (t))dt (3)C(P (t)) =(F (P (t))maxF)4(4)where P(s,x) is a set of possible paths L from the source points to x; C(.) is the cost function for back-tracking procedureto avoid cutting corners. T is then used for the Rivulet back-tracking procedure to determine the centre-lines of neuronsegments.2.3. Rivulet Back-TrackingThe source point s of Rivulet back-tracking is chosen at thelocation with maxF . s is usually the soma position, there-fore this step is equivalent to soma detection. Then the back-tracking procedure starts from the point P (t0) in the seg-mented foreground with the maximum value in T (x). In otherwords, P (t0) is the point with longest accumulative Ham-ming distance along the centreline from s. P (t) is tracked onT (x) by the gradient descent with fourth-order Runge-Kuttamethod (RK4).We define a gap G > 0 when back-tracking steps on thebackground voxels. G increments by 1 for each continuousstep on the background and reset to 0 when a foregroundvoxel is stepped. When G is greater than a predefined thresh-old, the back-tracking of P (t) stops. G is designed to stopthe back-tracking procedure to avoid miswiring. The back-tracking also stops when the Euclidean distance betweenP (ti) and P (ti−15) is less than one voxel size or P (ti) iswithin the voxel of s.We introduce a simple criterion to evaluate the confidenceof each traced branch P (t), as:P (t) =1|P (t)|∫ 10Ib(P (t))dt (5)where |P (t)| is the length of P (t). |P (t)| can be illustrated asthe proportion of gradient descent steps made on foregroundvoxels. The obtained branches with low confidences arediscarded. P (t) makes the Rivulet tracing procedure morerobust to the noise points. The first branch is directly addedto the neuron tree and considered as the initial trunk and thelast point of the branch is set as root. When a new branchis tracked, we apply time sequential wiring to determinewhether the branch is connected to the tree. The branch is notconnected if the Euclidean distance of termini p and the treeTr, D(p, Tr) < R×(rp+3) orD(p, Tr) < R×(rmin+3),where R is the connecting threshold; rp is the radius of pand rmin is the radius of the node with the shortest dis-tance to p. A spherical region ωi produced by the each newGround Truth Neuron Tube Snake APP2 RivuletFig. 1. The first row shows the reconstructions of a fly neuron with non-uniform distributed noises; the second row shows thereconstructions of a fruit-fly neuron with complex arborisation; the third row displays the reconstructions of a noise corruptedfrog neuron with serious discontinuity. The 3D neuron reconstructions shown here were all produced automatically withoutmanual correction.branch point is erased from T (x). The regions on T (x)covered by the tracked branches cannot be revisited in fol-lowing back-tracking iterations. The foreground point withhighest value in the erased T (x) is used as the starting pointof next back-tracking iteration. After each iteration, a cov-erage percentage Ψ is increased by∑x⊂ΩIb(x)/∑Ib(x),where Ω = {ω1 ∪ ω2 ∪ · · · ∪ ωn}. A higher Ψ revealsmore details of small branches. The back-tracking processonly terminates when Ψ reaches a predefined threshold (eg.98%). The time complexity of Rivulet tracing is O(n log n),where n is the number of the foreground voxels. Since thetime-crossing map erasing avoids duplicated fast-marchingand gradient calculation, the O(n log n) term comes with thefast-marching method alone.3. EXPERIMENTS AND RESULTSThe data we used to evaluate the proposed Rivulet algorithmwere provided by the BigNeuron Project [12]1. All the trac-ing results were visualised with Vaa3D [13]. Rivulet has been1http://alleninstitute.org/bigneuron/about/implemented in an open-sourced GUI Matlab toolbox2. Allhyper-parameters of the evaluated tracing algorithms weretuned by exhaust searches and validated with visual check.NeuTube Snake APP2 RivuletSD 6.63 8.35 4.71 3.59SSD 12.07 14.25 10.19 5.59SSD% 0.37 0.45 0.41 0.57Table 1. The mean SD, SSD and SSD% scores [7] of thecompared methods in Fig. 1.We compared Rivulet with the state-of-the-art automaticneuron tracing algorithms. The results are shown in Fig. 1.In the first row, Rivulet was able to trace the majority of theneuron arbors without being affected by the noises. Rivuletwas also shown to trace the image in the second row with anysmall gaps and complex arborisation. The image in the thirdrow has a low contrast which resulted in blurry boundaries be-tween the neuron and the background. Since the neuron seg-mentation is not used as the boundary wall in Rivulet tracing,2https://github.com/lsqshr/Rivulet-Neuron-Tracing-Toolbox(a) A neuron image with dis-continuous structures(b) A complete reconstruc-tion of the corrupted imageFig. 2. The discontinuous and fuzzy neurons images wereconstructed by Rivulet.(a) A noisy confocal micro-scopic neuron image(b) The reconstruction resultby RivuletFig. 3. The Rivulet reconstruction result of an extremelynoisy image.Rivulet was shown to be more robust to dark branches anddiscontinuous branches. In the quantitative analysis shown inTable 1, Rivulet was also shown to outperform the other com-pared methods regarding the spatial distance (SD), substantialspatial distance scores (SSD) [7]. The slightly higher meanpercentage of substantial spatial distance (SSD%) of Rivuletis due to the larger number of nodes of Rivulet.The reconstruction results of a fly neuron with discon-tinuous and fuzzy structures by Rivulet are shown in Fig. 2.The discontinuous structures were caused by imperfect stain-ing and excitation power during image acquisition. The gapthreshold and time sequential wiring enable Rivulet to tracediscontinuous neuron arbours and resistant to over-tracing atthe same time. The unevenly distributed fluorescent mark-ers of fly neuron lead to extremely noisy images. For thenoisy image shown in Fig. 3, Rivulet was the only algorithmamong the evaluated ones that reconstructed meaningful re-sults, mainly because of the gap threshold and the confidencedumping strategy.4. CONCLUSIONSIn this study, we proposed a novel automatic tracing algo-rithm, Rivulet, based on iterative back-tracking. Evaluatedwith data provided by the BigNeuron project, Rivulet showsgreat robustness to discontinuous and noisy neuronal struc-tures in confocal microscopic images. Rivulet also outper-formed other state-of-the-art neuron tracing algorithms on thechallenging images with low SNR, discontinuous arbours andcomplex topologies.References[1] H Peng and M Hawrylycz et al., “BigNeuron: Large-scale 3D neuron reconstruction from optical microscopyimages,” Neuron, vol. 87, no. 2, pp. 252 – 256, 2015.[2] R Parekh, G Ascoli, “Neuronal morphology goes dig-ital: a research hub for cellular and system neuro-science,” Neuron, vol. 77, no. 6, pp. 1017–1038, 2013.[3] S.L. Wearne, A Rodriguez, et al., “New techniques forimaging, digitization and analysis of three-dimensionalneural morphology on multiple scales,” Neuroscience,vol. 136, no. 3, pp. 661–680, 2005.[4] L Feng, T Zhao, et al., “Neutube 1.0: a new designfor efficient neuron reconstruction software based on theswc format,” eneuro, vol. 2, no. 1, pp. 49, 2015.[5] X Ming and A Li et al., “Rapid reconstruction of 3d neu-ronal morphology from light microscopy images withaugmented rayburst sampling,” PLoS ONE, vol. 8, no.12, pp. e84557, 2013.[6] Y Wang and A Narayanaswamy et al., “A broadly ap-plicable 3D neuron tracing method based on open-curvesnake,” Neuroinformatics, vol. 9, no. 2-3, pp. 193–217,2011.[7] H Peng, F Long, and G Myers, “Automatic 3D neurontracing using all-path pruning,” Bioinformatics, vol. 27,no. 13, pp. i239–i247, 2011.[8] H Xiao, H Peng, “APP2: automatic tracing of 3D neu-ron morphology based on hierarchical pruning of a gray-weighted image distance-tree,” Bioinformatics, vol. 29,no. 11, pp. 1448–1454, 2013.[9] J Yang, P.T. Gonzalez-Bellido, et al., “A distance-fieldbased automatic neuron tracing method,” BMC Bioin-formatics, vol. 14, no. 1, pp. 93, 2013.[10] D Adalsteinsson and J.A. Sethian, “A fast level setmethod for propagating interfaces,” Journal of compu-tational physics, vol. 118, no. 2, pp. 269–277, 1995.[11] S.M. Hassouna, A Farag, et al., “Multistencils fastmarching methods: A highly accurate solution to theeikonal equation on cartesian domains,” PAMI, IEEETransactions on, vol. 29, no. 9, pp. 1563–1574, 2007.[12] H Peng, E Meijering, and G Ascoli, “From DIADEMto BigNeuron,” Neuroinformatics, vol. 13, no. 3, pp.259–260, 2015.[13] H Peng, A Bria, et al., “Extensible visualization andanalysis for multidimensional images using Vaa3D,”Nature Protocols, vol. 9, no. 1, pp. 193–208, 2014.

Reconstruction of 3D neuron morphology using Rivulet back-tracking

https://ses.library.usyd.edu.au/bitstream/2123/20514/2/Cai07_ISBI-2016_Donghao-Zhang.pdf

Reconstruction of 3D neuron morphology using Rivulet back-tracking

Abstract

Similar works

Full text

Available Versions

Sydney eScholarship