Triangle mesh skeletonization using non-deterministic voxel thinning and graph spectrum segmentation

. In the context of shape processing, the estimation of the medial axis is relevant for the simplification and re-parameterization of 3D bodies. The currently used methods are based on (1) General fields, (2) Geometric methods and (3) voxel-based thinning. They present shortcomings such as (1) overrepresentation and non-smoothness of the medial axis due to high frequency nodes and (2) biased-skeletons due to skewed thinning. To partially overcome these limitations, this article presents a non-deterministic algorithm for the estimation of the 1D skeleton of triangular B-Reps or voxel-based body representations. Our method articulates (1) a novel randomized thinning algorithm that avoids possible skewings in the final skeletonization, (2) spectral-based segmentation that eliminates short dead-end branches, and (3) a maximal excursion method for reduction of high frequencies. The test results show that the randomized order in the removal of the instantaneous skin of the solid region eliminates bias of the skeleton, thus respecting features of the initial solid. An Alpha Shape-based inversion of the skeleton encoding results in triangular boundary Representations of the original body, which present reasonable quality for fast non-minute scenes. Future work is needed to (a) tune the spectral filtering of high frequencies off the basic skeleton and (b) extend the algorithm to solid regions whose skeletons mix 1D and 2D entities.


Introduction
The medial axis () of a compact 3D region Ω ⊂  3 , is defined as the set of all points  ∈ Ω such that the closest point in the boundary Ω is not unique.Such a medial axis is unique, defined as the finite union of compact 1-and 2-manifold sets.An approximation of the (Ω) is called an skeleton (Ω) .Some applications for skeletonization algorithms include shape decomposition, animation, medical images, virtual navigation and computational mechanics [18,20].This manuscript presents a skeletonization algorithm based on non-deterministic voxel thinning and graph spectral decomposition which articulates: (1) a novel non-deterministic thinning algorithm for computing non-skewed 1-manifold (Ω), (2) spectral graph segmentation for clearing small death branches, and (3) application of a maximal excursion algorithm to smooth the resulting skeleton and reduce high frequencies.Finally, our algorithm applies Alpha-Shape inversion on the computed skeleton (Ω), producing a retrieved surface   that is an approximation of the input body surface  = (, ) .The remainder of this article is organized as follows: Section 2 discusses the skeletonization methods already present in the literature.Section 3 discusses our skeletonization method.Section 4 presents skeletonization and alpha-shape inversion results for different triangle meshes.Finally, Section 5 presents the conclusions and introduces what remains for future work.

Literature review
Skeletonization methods are algorithms that produce an estimation of the Medial Axis of a 3D solid body that provides a lower dimensional (1D and 2D) representation of the original model.Depending on the approach for this dimensional reduction, skeletonization algorithms can be classified as: (1) general-field methods, (2) geometric methods and (3) voxel-thinning methods.The state of the art for the afore mentioned methodologies is briefly discussed below.For a more detailed survey, we refer the reader to [22].
General Field methods estimate the medial axis by computing a map field on Ω.The skeleton of Ω is extracted by connecting singularities on the map field [11,22].Refs [2,23] compute a distance field map that have proven to be efficient and becoming suitable for large datasets [4].Refs [10,22] compute potential fields that improves the robustness of the methods in the presence of input noisy shapes.Other methods use projection planes and inscribed balls to find the 2D curve skeleton.Afterwards, a convex hull operation is applied to generate the 1D skeleton [14].General field methods can be applied to mesh and voxel data inputs [14,22].In addition, they can produce also 2-manifold skeletons [2,22].However, they (1) do not ensure the connectivity of the resulting skeleton [10,23] and, (2) they produce non-smooth skeletons with small dead-branches [19,20,22].
Geometric methods find the symmetries and geometric properties in the boundary of a shape region to generate the medial axis [19,20].Voronoi diagrams [7], mesh contraction [25], octree-graphs [3], and union-of-balls [16] are popular geometric methods to compute the skeleton.These methods are applied to mesh or point-cloud representations [20], preserving the connectivity of the resulting skeleton [25] and generating non-skewed skeletons [25].However, these methods still produce skeletons with small dead-branches and high frequencies [19,20,22].
Voxel-based thinning methods iteratively remove the boundary set of a voxel representation of the original region [6,11,18].These algorithms usually preserve the connectivity of the skeleton by defining explicit rules such as: simple point criterion [15,17] or removing and recovering templates [12].In addition, some thinning methods are able to produce 2D (surface) skeletons [11].Current thinning-based methods present shortcomings such as: (1) computational expensive algorithms, (2) biased skeletons [13,27], and (3) nonsmooth skeletons with small dead-branches [19,20,22].Some authors have reduced the computational cost of these algorithms by applying parallel programming [6,12,27].However, the resulting skeleton does not preserve the connectivity.

Conclusions of the literature review
Based on the literature review, the algorithms for skeletonization can be classified in: (1) General field methods are efficient for large datasets, robust to noisy shapes and some methods produce curve and surface skeletons but can present problems like discontinuities, (2) Geometric method conserve connectivity and generate well centered skeletons, and (3) Voxel-based thinning can be robust to noisy shapes, preserve connectivity and some methods produce curve and surface skeletons but are computationally expensive and the skeletons are skewed.Also, the surveys show that the most part of the algorithms are not smoothed and/or produce small death branches.To partially overcome this problem, this article presents (1) a random labelling that neutralizes the bias, (2) deletion of the small death branches using a spectral segmentation of the graph skeleton to improve the centricity and smoothness of the skeleton and (3) the use of maximal excursion algorithm to improve the smoothness.We also perform an alpha-shape-based inversion that retrieve an approximation of the original body.

Methodology
This section discusses the method used to compute the skeleton that approximates the Medial Axis (Ω) of a compact and connected 3D solid Ω ⊂  3 .Our method requires a triangle mesh  = (, ) as input.The mesh M discretizes the boundary ∂Ω of the solid body Ω. M is a closed 2-manifold defined by the set of points  and the set of triangles .The proposed method proceeds as follows: (1) Our algorithm converts  into a voxel representation (Ω) 0 , where a point  ∈ Ω is called a black voxel and a point  ∈  3 \Ω is called white voxel.(2) A voxel-based thinning algorithm is applied to compute (Ω).In order to neutralize the skeleton bias of Voxel-based algorithms, a random labelling is applied to the voxelization in this step.The Graph Laplacian  of (Ω) and its corresponding spectral decomposition are computed.The spectral decomposition enables a segmentation of the skeleton in the frequency domain, which characterizes the large branches of (Ω). (4) Aided by the spectral segmentation, a shortest path algorithm is applied to (Ω), resulting in the suppression of the small death-branches.(5) A maximal excursion algorithm is applied to reduce high frequencies and improve the smoothness of (Ω). (6) Finally, an alpha-shape-based inversion is used to retrieve the original body Ω. Figure 1 shows the flowchart of the proposed algorithm.The remainder of this section describes the details of the skeletonization algorithm.

Voxel-based thinning
In order to apply a thinning algorithm to the solid body Ω, the triangle mesh  = (, ) is converted into a voxel representation (Ω) 0 .For this purpose, the voxelization algorithm presented in [1] is implemented.Thinning algorithms iteratively detect the boundary voxels (Ω) of (Ω) 0 , and delete the voxels in (Ω) that comply with an elimination criterion.To understand the implemented thinning algorithm, it is necessary to describe: (1) The neighbours of a voxel, (2) the simple point criterion, and (3) the end point criterion.
The neighbours of a voxel  ∈ Ω are the 26 voxels that surround such a voxel.This neighborhood is sub-classified as 6-18-and 26-neighbors, as follows: (a) 6-neighbor voxels share one face with the central voxel p.(b) 18-neighbor voxels share at least one edge with the central voxel p. (c) 26-neighbor voxels share at least one vertex with the central voxel [17].In order to preserve the connectivity of the skeleton, the simple point criterion is implemented [15,17].The simple point criterion, requires a voxel p to satisfy the following conditions: 1.
The central voxel p has a white voxel in its 6-neighbor.

2.
The central voxel p has a black voxel in its 26-neighbor.

3.
The set of black 26-neighbors of the central voxel p is 26-connected.4.
The set of white 6-neighbors of the central voxel p is 6-connected in the set of white 18-neighbors.
On the other hand, the end point criterion is implemented to preserve the shape of the input body.A voxel satisfies the end point criterion if it has a single voxel neighbor in its 26-neighbor set.
The implemented thinning algorithm follows the next steps: (1) To voxelize the input mesh  = (, ) generating the voxel set (Ω) 0 .(2) The boundary voxels (Ω)  are detected by evaluating the simple point and the end point criteria.(3) The voxel set is compared with the boundary voxel set.If the voxel set is equal to the boundary ( (Ω)  = (Ω)  ?), the thinning procedure finishes.(4) the random labelling ((Ω)  ) is applied.(5) The voxels in the skin Ω that comply with the simple point criterion are removed.Finally, (6) steps 2-5 are repeated until step 4 is satisfied, generating the resulting skeleton (Ω).

Volume estimate
To retrieve the original body Ω from the skeleton (Ω), it is necessary to estimate the distance of each voxel in (Ω) to the boundary Ω .For this purpose, the next two methods can be used: 1. Arcelli et.al. [2] use the method of reverse distance transform to retrieve the original region.This method find the radius inscribing the sphere with the minimum radius that is tangent to ∂Ω in one or more points.
2. Munoz et.al. [18] use a measure of mass (Computed using the diffusive advection proposed by Jalba et.al. [11]) in each voxel of SK(Ω) to approximate a cylinder in each segment and considering a constant density.

Graph skeleton representation
To generate the spectral representation, it is necessary to convert the voxelized skeleton (Ω) to its graph representation  = (, ), where  is a connected and non-oriented graph defined by a set of vertices  and a set of edges .In graph representation, each point   ∈  is represented by the center point of each voxel in (Ω), and each edge   ∈  is computed by connecting the adjacent black-voxels in each 26-neighborhood.A function   :  (  )−>  is assigned to each edge   ∈ .  stores the radius of every point contained in the line segment  (  ), defined by the edge   .In this process, any 3-vertex closed cycle (triangle) is deleted and replaced with its respective barycenter.

Spectral segmentation
The spectral decomposition of the graph  = (, ) is a representation of the skeleton (Ω) in the frequency domain.As such, this spectral decomposition preserves the topology of .However, the low frequencies of this skeleton can be extracted from the first 3 eigenvectors (i.e. with lowest non-zero eigenvalue) of the Laplacian operator on graphs [24,28].The graph Laplacian L and its corresponding spectrum is computed as follows: 1.
Find the    weighted Laplacian matrix  of , defined as follows: 2.

Small death branches removal
With the identification of the large features  = ( 1 ,  2 , … ), the next step is use the shortest path algorithm of Dijkstra [8] to remove the small death branches in  = (, ).This process compute a new graph   = (  ,   ) without the small death branches where   ⊂ , the vertices   ⊂  and the edges   ⊂ .

High frequencies removal
To reduce the oscillations in the skeleton, a maximal excursion algorithm (edge linearization) is used [18].This reduce the high frequencies and improve the smoothness of   = (  ,   ), generating a smoothed graph   = (  ,   ).In each iteration of this algorithm, is evaluated if the maximum distance of each new segment is greater than a predefined threshold.The algorithm stops when the maximum distance of all new segments is less than the threshold.The mass of a new segment is the sum of the mass of all the vertices in the new segment created by the old skeleton.

Retrieval of the original region
To retrieve the original body Ω, an alpha-shape based inversion is used [9], generating a 2manifold triangular mesh   = (  ,   ) that is an estimation to Ω, defined by a set of points   ⊂ Ω and a set of triangles   .For this propose, a set of discs with different radius are created in each segment of   = (  ,   ).The points of these discs are the input of the alpha-shape.

Conclusions and future work
This manuscript presents a novel method to improve the centeredness centricity and the smoothness of the skeletons generated by thinning algorithms.The presented method articulates: (a) a random labelling applied to the voxel-based algorithm that neutralize the skewed produced by the order of analysis in the skeletons, (b) the deletion of the small death branches using a spectral segmentation, and (c) the smoothing of the skeleton reducing the high frequencies with a maximal excursion algorithm (edge linearization).The results show that using the spectral segmentation in graph of skeleton is an effective method to produce a smoothed skeleton and delete undesirable noise.Also, the retrieved models show that our method can produce a good approximation of the original body.Future work includes (1) the tune the spectral filtering of high frequencies off the basic skeleton and (b) extend the algorithm to solid regions whose skeletons mix 1D and 2D entities.
), 3(b) and 3(c) shows the voxel skeletons (Ω) with their respective mass distribution.Figures3(d),3(e) and 3(f) shows graph skeletons with the set of computed radius using the reverse distance transform method.In this step SK(Ω) presents high frequencies (e.g.figure3(b)) and small death branches (figure3(c)).

Fig. 3 .
Fig. 3. Resulting skeletons from our unskewed thinning algorithm and the respective node masses.The skeletons present high frequencies and short dead branches.(a) Truss (b) Snake and (c) Horse skeleton with the distribution of mass.(d) Truss (e) Snake and (f) Horse graph skeleton with the computed radius.

Figure 4 (
Figure 4(a) shows the skewed Horse graph skeleton produced using the algorithm without the proposed randomization.Figure 4(b) shows the unskewed Horse graph skeleton produced using the algorithm with the randomization.The comparison between the red squares in figures 4(a) and 4(b) shows that the small death branches are skewed in figure 4(a) due the direction of analysis, while the branches are unskewed in the randomized version (our algorithm).Figures 5(a) and 5(b) shows the skewed and unskewed retrieved model Snake (beige color) computed using the reverse distance transform method and compared with the input mesh  = (, ) (blue color).Compared with the figure 5(b) (unskewed), figure 5(a) (skewed) have significant less volume that the original body, indicating a deviation grater from the medial axis.

Figure 4 (
Figure 4(a) shows the skewed Horse graph skeleton produced using the algorithm without the proposed randomization.Figure 4(b) shows the unskewed Horse graph skeleton produced using the algorithm with the randomization.The comparison between the red squares in figures 4(a) and 4(b) shows that the small death branches are skewed in figure 4(a) due the direction of analysis, while the branches are unskewed in the randomized version (our algorithm).Figures 5(a) and 5(b) shows the skewed and unskewed retrieved model Snake (beige color) computed using the reverse distance transform method and compared with the input mesh  = (, ) (blue color).Compared with the figure 5(b) (unskewed), figure 5(a) (skewed) have significant less volume that the original body, indicating a deviation grater from the medial axis.

Figure 6 (
Figure 6(a) shows the graph skeleton  = (, ) of the Horse dataset with the small death branches and high frequencies.Figure 6(b) shows the eigenvalues with the lowest frequencies of  were the small death branches are squashed and allows the manual graph segmentation  = (_1, _2, . . .).Figure 6(c) shows  = (, ) after the application of the shortest path algorithm using the   = ( 1 ,  2 , . .), that deletes the small death branches.Figure 6(d) shows  = (, ) after the application of the maximal excursion algorithm that reduce the high frequencies.The qualitative comparison shows that the final model (figure 6(c)) is smoothed in comparison with the initial graph skeleton (figure 6(a)).Also, the number of points that represent the model decrease significantly, from 735 points (figure 6(a)) to 25 points (figure 6(d)).

Figure 6 (
Figure 6(a) shows the graph skeleton  = (, ) of the Horse dataset with the small death branches and high frequencies.Figure 6(b) shows the eigenvalues with the lowest frequencies of  were the small death branches are squashed and allows the manual graph segmentation  = (_1, _2, . . .).Figure 6(c) shows  = (, ) after the application of the shortest path algorithm using the   = ( 1 ,  2 , . .), that deletes the small death branches.Figure 6(d) shows  = (, ) after the application of the maximal excursion algorithm that reduce the high frequencies.The qualitative comparison shows that the final model (figure 6(c)) is smoothed in comparison with the initial graph skeleton (figure 6(a)).Also, the number of points that represent the model decrease significantly, from 735 points (figure 6(a)) to 25 points (figure 6(d)).

Figure 6 (
Figure 6(a) shows the graph skeleton  = (, ) of the Horse dataset with the small death branches and high frequencies.Figure 6(b) shows the eigenvalues with the lowest frequencies of  were the small death branches are squashed and allows the manual graph segmentation  = (_1, _2, . . .).Figure 6(c) shows  = (, ) after the application of the shortest path algorithm using the   = ( 1 ,  2 , . .), that deletes the small death branches.Figure 6(d) shows  = (, ) after the application of the maximal excursion algorithm that reduce the high frequencies.The qualitative comparison shows that the final model (figure 6(c)) is smoothed in comparison with the initial graph skeleton (figure 6(a)).Also, the number of points that represent the model decrease significantly, from 735 points (figure 6(a)) to 25 points (figure 6(d)).

Figure 6 (
Figure 6(a) shows the graph skeleton  = (, ) of the Horse dataset with the small death branches and high frequencies.Figure 6(b) shows the eigenvalues with the lowest frequencies of  were the small death branches are squashed and allows the manual graph segmentation  = (_1, _2, . . .).Figure 6(c) shows  = (, ) after the application of the shortest path algorithm using the   = ( 1 ,  2 , . .), that deletes the small death branches.Figure 6(d) shows  = (, ) after the application of the maximal excursion algorithm that reduce the high frequencies.The qualitative comparison shows that the final model (figure 6(c)) is smoothed in comparison with the initial graph skeleton (figure 6(a)).Also, the number of points that represent the model decrease significantly, from 735 points (figure 6(a)) to 25 points (figure 6(d)).

Fig. 6 .
Fig. 6.Geometric simplification of the graph skeleton.(a) Horse graph skeleton with small death branches and high frequencies.(b) Lowest frequency eigenfunction of Horse graph.(c) Horse graph skeleton after deleting small death branches.(d) Horse graph skeleton after reducing high frequencies.

Figure 7
Figure7shows the retrieval model using an alpha shape-based inversion.Figures7(a), 7(b) and 7(c) shows the retrieval after the deletion of the small death branches that is a good approximation to Ω.

Fig. 7 .
Fig. 7. Retrieval of the original model using alpha-shape based inversion.Retrieval of the model (a) Truss (b) Snake and (c) Horse.