Coverage Axis: Inner Point Selection for 3D Shape Skeletonization

Given a closed, oriented and bounded two-manifold surface $\mathcal{S}$ in $\mathbb{R}^{3}$, medial axis is defined as the locus of the center of maximally inscribed spheres that are tangent to $\mathcal{S}$ at two or more points. The medial axis $\mathcal{M}$ of $\mathcal{S}$, together with its radius function $\mathcal{R}$, forms the medial axis transform (MAT), denoted by a pair $\big{(}\mathcal{M},\mathcal{R}\big{)}$.

Voronoi diagram is a partition of the domain $\Omega\subset\mathbb{R}^{d}$ into regions, close to a set of points called generators $\{\mathbf{x}_{i}\in\Omega\}_{i=1}^{n}$. Each region, also called cell, is defined as

A well-known MAT initialization technique is to use the Voronoi diagram inside a model generated by sampled points on the surface.

Power diagrams [Aur87] can be viewed as an extension of Voronoi diagrams, where each generator $\mathbf{x}_{i}$ is equipped with a weight $w_{i}$ to control its influence. By defining the power distance $d^{\text{pow}}(x,x_{i})$ between $x$ and the weighted generator $x_{i}$ to be $\|\mathbf{x}-\mathbf{x}_{i}\|^{2}-w_{i}$, the cell associated with $\mathbf{x}_{i}$ is defined by

A generator with a larger weight is more dominant, and when all the weights are equal, the Power diagram is reduced to a Voronoi diagram. In [ACK01], the approximate MAT structure is built based on regular triangulation, the dual of Power diagram, with the weights as squared radius of inside and outside poles (a subset of Voronoi vertices). Examples of a Voronoi diagram and a Power diagram are shown in Figure 1 respectively.

Our method starts with the generation of candidate points located inside a given 3D shape. The input to our system can be either surface meshes or point clouds with normals.

For a mesh input, following a common initialization technique, we compute the Voronoi diagram of a set of points sampled on the mesh and extract the inner points, which become the candidate points of the skeleton. The inside Voronoi diagram becomes the initial coverage axis structure.

For a point cloud input (with normals), similar to a mesh input, we first compute the Voronoi diagram w.r.t. surface samples. Then, we extract the inner points as follows. We first compute its Delaunay triangulation which is the dual of the Voronoi diagram w.r.t. the input point cloud. Consider a Voronoi vertex $p$ and the vertices of its dual tetrahedron $p^{\prime}_{0},p^{\prime}_{1},p^{\prime}_{2}$ and $p^{\prime}_{3}$. The candidate $p$ is considered as an inner point only if we have $\overrightarrow{pp^{\prime}_{i}}\cdot\overrightarrow{n}(p^{\prime}_{i})>0$, $\forall i=0,1,2,3$. Here $\overrightarrow{n}(p^{\prime}_{i})$ is the input normal of $p^{\prime}_{i}$ and $\cdot$ is dot product. See more details in Appendix A. Note other inside-outside query methods such as [BDS*18] are also applicable.

The Voronoi-based point generation is based on principled geometric transform which can generate candidates with good quality. Nevertheless, note our method does not rely on strictly defined MAT for initialization. We can also randomly sample a set of points inside the shape and use them as skeleton candidates (see Sec. 6.1 for a detailed discussion). This property significantly enhances the flexibility of our method for handling various inputs in different applications.

Given a set of elements (called universe) and its subsets whose union equals the universe, the set cover problem (SCP) is to identify the minimum sub-collection of these subsets whose union equals the universe. Inspired by this formulation, we now introduce our Coverage Axis. In the last step, we obtain a set of inner point candidates $P=\{p_{i}\}$. We estimate the radii of these points $R=\{r_{i}\}$ as their closest distances to the boundary surface, resulting in candidate balls. Now we slightly dilate all the balls by adding a small value to their radii, i.e., $r_{i}^{\prime}=r_{i}+\delta_{r}$, leading to a set of dilated balls $B=\{(p_{i},r_{i}^{\prime})\}$. The dilation makes our algorithm robust to insignificant details and noise, for which we give detailed discussions later in Sec. 6.2.

Our goal is to find the minimum number of dilated balls that cover all the sampled points $S=\{s_{j}\}$ on the surface. For this purpose, we introduce a coverage matrix $\mathbf{D}\in\{0,1\}^{m\times n}$, where $m$ and $n$ are total numbers of sampled surface points and candidate skeletal points, respectively. Each element $d_{ji}\in\{0,1\}$ of $\mathbf{D}$ indicates if a surface point $s_{j}$ is covered by the dilated ball $(p_{i},r_{i}^{\prime})$:

Let $\mathbf{v}\in\{0,1\}^{n\times 1}$ be a decision vector, where the $i$-th element $v_{i}\in\{0,1\}$ indicates if a candidate skeletal point $p_{i}$ is selected. Now we can derive the formulation of our 0-1 integral optimization problem:

where $\mathds{1}$ is a vector of ones and $\geq$ is applied element-wise to the vector entries. Here we minimize the norm of the decision vector $\mathbf{v}$ under a certain constraint. The constraint enforces each sampled point to be covered by the selected skeletal balls, which ensures consistency between the selected skeleton and the input shape. With this simple formulation, our method allows a global consideration of the overall shape structure and effectively selects the most predominant inner candidates, giving an expressive and compact skeletal representation to faithfully capture the fundamental geometry and topology. We solve Eq. 2 using Mixed-integer linear programming (MILP) solver in Matlab (MathWorks 2021). The pipeline of inner point selection is demonstrated in Figure 2. We further discuss a series of properties of our method, e.g., centrality, robustness to noise, etc., in Sec. 6.

It is notable that the constraint in Eq. 2 does not always lead to a solution. For example, if there are points not covered by any dilated balls (assume the dilation is tiny), no solution can be found. Nevertheless, this barely happens and the fixed parameters always lead to a solution for the extensive experiments as we will show later. This is owing to the proper initialization for generating inside candidates: Voronoi Diagram or densely sampled points. Voronoi Diagram produces balls maximally inscribed to the surface, thus already giving tight-fitting candidates. The randomly generated inner points are considerably dense, of which a subset approximates the vertices of the Voronoi Diagram, and the excessive balls provide abundant candidates for solving SCP.

Once the skeletal points are determined, in this step, we analyze their connectivity to form a structured mesh. Generally, if the input is a surface mesh, the topological connections of the skeletal points are easier to be derived from the surface. However, this is difficult when the input is a point cloud that only consists of unorganized points.

Considering our method flexibly allows different input types (i.e., surface mesh and point cloud) with different candidate generation strategies (i.e., Voronoi-based and random sampling), we propose to use different ways to build the connections of skeletal points. More details can be found in Appendix B.

We conduct comparisons with representative MAT simplification methods including Q-MAT [LWS*15] and SAT [MGP10]. The errors are scaled by the length of the diagonal of the input model bounding box.

For point cloud input, we compare with two closely relevant methods: Deep Point Consolidation (DPC) [WHG*15] and Point2Skeleton [LLL*21], given that these two methods can also predict MAT-based skeletons directly from point clouds. Note Point2Skeleton is a deep learning-based method that requires training from a large amount of data.

A good skeleton requires each skeletal point to be located at the center of its corresponding local geometry. Our method is able to effectively select the local centers for generating the skeletal representation. This property also enables our algorithm to handle randomly generated inner points instead of relying on computing the strictly defined Voronoi diagram as an input. In this section, we elaborate on this centrality property of our point selection strategy.

We first generate the point candidates by randomly sampling inside the shape and then applying our point selection algorithm introduced in Sec 4.2. The number of randomly generated inner points is $10000$, and the offset is set to $0.02$. A selection of results from these randomly generated points are shown in Figure 8. Thanks to the idea of optimal set coverage, it can be observed that the selected points are located at the centers of local geometry, even if the initial points are derived by random sampling. An evaluation of the influence of different randomly generated candidate inner point numbers $|P|$ for centrality is given in Figure 9. We find different numbers of candidate inner points do not show a large influence on the centrality.

In the following, we give an analysis of the rationale of the centrality property of our algorithm. Recall in Eq. 2, the constraint of the optimization is that the union of dilated balls needs to cover the sampled surface points. As explained in Sec. 6.4, the dilated radius should not deviate too much from the original radius, which means $r^{\prime}\approx r$. Consider a local shape shown in Figure 10. The coverage goal of one candidate point $P_{a}$ is defined by its farthest distance to the local geometry, e.g., $d_{a}$, given that all the surface samples are covered only if the farthest sample is covered. According to the definition of inner balls, the coverage ability is determined by its radius $r_{a}$ (note that we have $r_{a}^{\prime}\approx r_{a}$), which is actually the distance to the nearest surface samples of the local geometry.

The goal of finding the minimum number of dilated balls to cover all surface samples can be viewed as selecting those points with the most powerful coverage ability (cover all the target samples). That to say, its radius $r_{a}$ should be as equal as possible to the distance to farthest point $d_{a}$ so that we have $r_{a}^{\prime}\geq d_{a}$ which satisfies the coverage constraint. Those balls whose radii are close to the distance to the farthest point that needs to be covered will be preferred during optimization. Note that $r_{a}\approx d_{a}$ is consistent with the definition of the medial axis, which is intrinsically a set of centers of maximally inscribed spheres. In this case, the inner point $P_{b}$ is preferred for the local geometry coverage since $d_{b}-r_{b}<d_{a}-r_{a}$. Note that we always have $d_{b}\geq r_{b}$. Finally, as shown in Figure 10, $P_{b}$ is a better approximation than $P_{a}$ of skeleton point and will be selected by our strategy.

We evaluate the robustness of our method to surface noise. We compare with Q-MAT [LWS*15], which has the state-of-the-art performance in terms of approximation accuracy and goodness of structure (i.e., compactness and shape awareness). The number of vertices of the skeletal points are set the same for both methods. As shown in Figure 11, given noisy input, our method can still preserve the main features and achieve relatively high-precision shape approximation with a compact representation. The robustness mainly owes to the ball dilation strategy as well as the formulation based on global coverage. The dilated balls cover the insignificant details cause by local perturbations, while the overall coverage formulation gives a high-level optimization rather than focusing on local details. In contrast, Q-MAT [LWS*15], which uses step-by-step edge collapse based on local errors, cannot capture the fine details when the input surface has great noise.

As aforementioned, the relationship between inner points and surface is built by the coverage, which is achieved by dilating the medial radius. There are two typical ways for inner ball dilation. One is by adding a specified value $\delta_{r}$ to the radius of each ball. The other is to scale the ball by multiplying the radius by a factor $\sigma_{r}$. We now discuss the effect of these two dilation strategies. To facilitate the discussion, we refer to the first way as offset and the second as scaling.

We take the model Octopus and Ant-3 as examples. Let $r$ and $r^{\prime}$ denote the original medial radius and dilated medial radius respectively. For the offset manner, we set $r^{\prime}=r+\delta_{r}$ where $\delta_{r}=0.02$, while for the scaling manner, we set $r^{\prime}=r\times\sigma_{r}$ with $\sigma_{r}=1.5$.

As shown in Figure 12, the offset manner tends to generate selected points that are distributed more evenly, while the scaling manner can distinguish parts that have different radii, which leads to relatively large gap around the joints of different parts (as shown in the box region). The reason behind this is, adding a constant value is independent of the original radius, but the scaling makes the larger ball dilate more and control more area, which results in part differentiation.

We further discuss the time efficiency of our point selection strategy. The statistics is summarized in Table 3. We use Mixed-integer linear programming (MILP) solver in MATLAB (MathWorks 2021) on the platform explained in Sec 5.

As aforementioned, the experiments are conducted with mesh input using $1500$ surface point samples, candidate inner points generated by $4000$ surface samples. The dilation offset is set to $0.02$. It can be viewed from Table 3 that the point selection only takes a few seconds for most shapes. For some cases, the complexity of the searching space of the Set Coverage problem may increase, which leads to one to two minutes of running time. The average running time for all the samples is $16.4$ s.

Table 3 also reveals the difference in running time for various shapes, where the difference is because the complexity of the solution space in SCP varies for different geometries. Generally, solving for the coverage of plate-like geometries is more complex than tube-like ones, since the plate-like models have more eligible candidates that satisfy the coverage constraints, which leads to a larger search space during optimization. Since many shapes consist of both tube-like and plate-like structures, there is a variation in running times. We further report running time comparison in Table 4.

SCP is a known NP-Hard problem [Har82] and its running time complexity is $O(|P|^{2}2^{|S|})$ where $|P|$ and $|S|$ are the number of candidate inner points and surface points, respectively. To investigate the influence of $|S|$ and $|P|$ on running time efficiency and the resulting number of selected skeletal points $|V|$, we conduct additional experiments of which results are summarized in Table 1 and Table 2 respectively. We adopt the Octopus-1 model for evaluation. We fix $|C|=4000$ in Table 1 and $|S|=1500$ in Table 2. From the results, both $|S|$ and $|P|$ do not show a significant influence on the selected vertex number. However, $|S|$ has a larger impact on the running time, while $|P|$ shows less effect, which is consistent with the aforementioned analysis of time complexity.

To evaluate the time efficiency for different algorithm settings, we conduct an additional experiment where the initial candidate points are generated by random sampling. Here, we randomly sample $10000$ points inside a shape and test the time efficiency using two different dilation strategies and different dilation factors. The results are given in Table  7. It can be observed the running time is stable in different configurations.

Since the set coverage problem is NP-hard, the action space grows exponentially with respect to the number of surface points. Fortunately, an interesting feature of the proposed method is the flexibility in handling large complex models. The overall insights lie in the combination of divide and conquer strategy and our set coverage formulation. In this case, we take a mesh model with Voronoi initialization for explanation.

The first step is decomposing a mesh into smaller and meaningful sub-meshes with some existing methods, e.g., [YL21, HHF*19]. Then we have surface samples divided into several subsets $S_{1},S_{2},...,S_{n}$ where $n$ is the number of surface components. According to the correspondence between the vertices on the inside Voronoi diagram and the surface samples, we partition the inside Voronoi diagram into $P_{1},P_{2},...,P_{n}$ accordingly. For each component $i$, with surface samples $S_{i}$ and inner point candidates $P_{i}$, we solve the SCP and obtain the selection result $P^{+}_{i}$ embedded on the sub-Voronoi diagram. This process reduces the size of the problem and speeds up our algorithm by divide and conquer.

Finally, we merge each Voronoi diagram as well as the selected points. Then, we have all the selected points embedded on the original inside Voronoi diagram so that the model’s simplification can be achieved by removing redundant components as we described in Sec. 4.3. We demonstrate the results in Figure 14 (a-b). Note that the complex model is with $30000$ surface samples which are further used for generating candidate inner points.

Our method allows users to specify the ignorable surface parts that can be neglected. To indicate the ignorable surface points, the constraint in Eq 2 are reformulated as $\mathbf{Dv}\geq\mathbf{B}$, where the corresponding values for the ignorable surface points in $\mathbf{B}$ are set to $0$ and the others to $1$. In this way, the optimization tends to remove the skeletal points in the specified parts to favor parsimony. Users can also directly indicate if a skeletal point should be preserved or removed, for which we just simply set the corresponding values in $\mathbf{v}$ to 0 or 1 after the point selection.