healthmedicinet

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

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 45]. 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 region R except at O, whose coordinate is (i, j, k). ?(i, j, k) is the angle between PO and the gradient vector direction. The VCI of P is defined as follows:

(1)

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 R whose radius is rad. To finish the discretization computation efficiently, the polar coordinate is introduced
into this scheme, and the SVF is calculated as

(2)

with

(3)

where M is the number of support region lines radiating from the center pixel O(x, y, z), ? denotes the radial coordinate, a and b stand for the angular coordinates, d is the thickness of sliding volume, r is the center position of the sliding volume in the support region line ranging from
R _min
to R _max
, Q is the points between [r?d/2,r?+?d/2], and ?(qx _?
,?qy _?
,?qz _?
) is the angle between the gradient vector at Q and the direction of QO. Additionally, the angles a???[0,?2?] and b???[0,??] are divided into 2?L parts and L parts, respectively. Thus, M?=?2 L ²
. 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 45]. Hence, the sliding volumes of support regions of interior points are more likely
to converge in the marginal regions. As shown in Fig. 3a, the voxel A in the interior of the nerve is more likely to be selected as a raw seed than the
external voxel B. Because A has a higher SVF response value than voxel B, the orientations of gradient vectors
in the sliding volumes of support region of A are more likely to point to A. However,
the orientations of gradient vectors in the sliding volumes of the support region
of B are not consistent and sometimes point away from B. Moreover, a nerve cell is
not a uniform 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 20]. As shown in Fig. 3a, I is the volume data set, ?I(p) is the gradient vector at voxel p with its coordinate (x, y, z) in I, and ev₁
, ev₂
and ev₃
are the eigenvectors computed from Hessian matrix of I. ev₁
(p) is the principle direction along the center lines of the tube-like structure, and
ev₂
(p) and ev₃
(p) are the other two orthogonal eigenvectors. The seeds near the center of the tube-like
structcure meet the condition of Eq. 4.

(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

(5)

where I _SVF
(p) is the intensity of point p after SVF enhancement, I(p) is the intensity of point p before SVF enhancement, and SVF(p) is the SVF value of point p.

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 33], 46]. Extended the application of the open curve snake model to neuron tracing. However,
the computation was tedious in the tracing framework of 33], especially in the step of branch detection. The proposed 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

(6)

E _Total
is the total image energy combined with internal energy and external energy. This
model is a traditional deformable model, which resembles previous work in 16]. The open snake model is a parametric curve, c(s)?=?(x(s), y(s), z(s)), s?[0, 1], and the snake internal and external energy are defined as follows:

(7)

(8)

with

In Eq. 7, ? and ? are the “elasticity coefficient” and “stiffness coefficient”, respectively, in internal
energy, and they can control the regularity of the curve in the process of evolution.
In Eq. 8, the external energy term is used to make the snake deform near the center line of
the neuron and stretch the endpoints to the tail of the neuron. ?E _im
is the negative normalized Gradient Vector Flow (GVF) of the volume data enhanced
by SVF, p signifies point (x(s), y(s), z(s)) on the open curve, and I _SVF
is the volume after SVF enhancement in this paper. Instead of the original 3D image
GVF, we calculated the GVF of I _SVF
. The SVF can enhance the 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 endpoints c?(0) and c?(1). The c_s
(s)/||c_s
(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,
?E _str
(c(s)) approaches zero, and the open active curve converges. According to a large number
of experiments, this strategy is not only efficient and reliable but also can avoid
the leakage of the neuron boundary. To minimize the energy function E _Total
, the points on the snake curves are updated as:

(9)

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 ev₁
(p) and ?ev₁
(p). As shown in Fig. 4b, along with the open snake curve moving to the center of neuron, it also extends
toward the two inverse tangential directions, c_s
(p₀
) and ?c_s
(p₁
), in which the p₀
and p₁
are the two temporary endpoints. During the procedure of extension, the seed points
belonging to one curve were labelled with new values (the default value is zero) in
accord with the ID of the curve. For example, in Fig. 4b the yellow points and green points belong to different curves.

C. Collision detection

Neurons have many branches, especially in the dendrite region. Hence, detecting branching
points and handling collision are essential. In the proposed scheme, two types of
collision exist in the collision region and are shown in Fig. 4b. The first collision is branching point collision, which occurs when the open snakes
reach a seed point whose value is not zero, and this point is recorded as the branching
point (pink point in Fig. 4b). This branching point detection strategy is based on labelling seeds and is highly
efficient. It also can handle the second type of collision, contour lines collision.
The contour lines coming from the following step of radius estimation are the foundation
of neuroanatomy reconstruction. However, due to the ambiguity of radius estimation
in the collision region, the contour lines from two curves easily intersect. In Fig. 4b, this situation is illustrated in the imaginary pink circle and the embedded image,
which is an experimental result in the branching region. This collision will influence
the accuracy of the following reconstruction algorithm. Hence, a backoff strategy
is proposed to avoid contour line collision. First, radius estimation in the branching
points will not be executed. Second, if an extending curve reaches the branching point,
it will be cut back the length of D, which is usually set as double the average estimated radius of the current curve.

In other words, the imaginary pink circle is not necessary in radius estimation because
the following reconstruction algorithm would interpolate the information using triangular
meshes automatically. Finally, the tracing algorithm ends when all the seed points
are traversed.

The entire tracing algorithm procedure is shown as follows:

In summary, compared to the open snake method in 33], we improved this model in the following three aspects. First, the volume after SVF
enhancement has more straightforward gradient vectors, which point to the center line
of the neuron and can be used in driving the initial lines to the center of the neuron.
Second, the proposed method can cut down the computation of the stretching force of
end nodes. Third, in the step of collision detection, compared to the method based
on labelling neighbor voxels, the method based on labelling seeds has higher detection
efficiency and benefits the following reconstruction procedure.

Radius estimation

Radius estimation is another critical task in neuron anatomy reconstruction, and it
can provide more quantitative information for neuroscience research. Peng, Aylward,
and Wang had proposed some radius estimation methods 16], 20], 33], but most of them are based on the assumption that the neuron have a uniform 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) 45]. The SBF can converge on the real edge of a neuron cross-section that has the rounded
convex region in 45].

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 ₁
is the tangential direction in S _i
, v ₂
and v ₃
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 ₁
, which points to the tangential direction of the open curve. Additionally, the v ₂
and v ₃
are the orthogonals in the cross-section. To obtain an accurate estimation of neuron
cross-section, 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.

(10)

with,

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 _r
,?y _r
,?z _r
) are the spatial coordinates of S. The computation method of SCI in point P, which is in the range of [r?d/2, r?+?d/2], has been introduced in Eq. 4. d is the width of the sliding band, r is the distance between B and the center point S, and it can slide in the range of [R _min , R _max
] to obtain the optimal position of B with the maximum SBF response value. The boundary points B can be connected clockwise to fit the edge of the neuron cell.

In the proposed method, the parameter is related to the accuracy of radius estimation. As shown in the embedded image of
Fig. 5, the larger is, the more accurate the cross-section fitting will be. However, considering the
efficiency and accuracy in the actual application, the parameter should be adjusted flexibly.

Neuron surface reconstruction

Most of the traditional neuron reconstruction methods were based on the fast marching
method and some supplemental processes for connecting different fragments 33], 47]. However, in this paper, Liu’s non-parallel contour lines surface reconstruction
method is employed for surface reconstruction 48], 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 48]. To obtain a smoother neuron surface, we use a surface diffusion smoothing algorithm
to minimize the curvature of the model surface to obtain a smooth 3D model 49]. As shown in Fig. 6, the initial neuron surface (Fig. 6e) is formed by contour lines (Fig. 6d) which are obtained by radius estimation, and the final smooth neuron surface is
shown in Fig. 6f. In addition, the most outstanding advantage of Liu’s method is that it can automatically
handle branch reconstruction, especially of circles without intersections in the branching
region (the intersection problem was resolved through removing the collisions of circles
in the 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 48]. b Projection of points and lines on MA 48]. c Triangulation of adjacent contour lines 48]. d The contour lines of neuron cell. e The initial surface from triangulation of adjacent contour lines. f The final surface model after smoothing

Fig. 7. Process of neuron reconstruction in branching region. a The input branching data from 48]. b The reconstruction result of input data of (a). c The input data of neuron contour lines. d The reconstruction result of (c), in which different branches are labelled by different colors