A Variational Loop Shrinking Analogy for Handle and Tunnel Detection and Reeb Graph Construction on Surfaces

05/27/2021 ∙ by Alexander Weinrauch, et al. ∙ Max Planck Society TU Graz 0

The humble loop shrinking property played a central role in the inception of modern topology but it has been eclipsed by more abstract algebraic formalism. This is particularly true in the context of detecting relevant non-contractible loops on surfaces where elaborate homological and/or graph theoretical constructs are favored in algorithmic solutions. In this work, we devise a variational analogy to the loop shrinking property and show that it yields a simple, intuitive, yet powerful solution allowing a streamlined treatment of the problem of handle and tunnel loop detection. Our formalization tracks the evolution of a diffusion front randomly initiated on a single location on the surface. Capitalizing on a diffuse interface representation combined with a set of rules for concurrent front interactions, we develop a dynamic data structure for tracking the evolution on the surface encoded as a sparse matrix which serves for performing both diffusion numerics and loop detection and acts as the workhorse of our fully parallel implementation. The substantiated results suggest our approach outperforms state of the art and robustly copes with highly detailed geometric models. As a byproduct, our approach can be used to construct Reeb graphs by diffusion thus avoiding commonly encountered issues when using Morse functions.



There are no comments yet.


page 1

page 5

page 6

page 7

page 8

page 9

page 10

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

Probably one of the most inspiring achievements of topology is the classification theorem for surfaces, which establishes equivalence classes based on the Euler characteristic and orientability. This development triggered a formidable effort focalized on the so called Poincaré conjuncture aiming at extending classification to higher dimensions and in particular to the manifold setting and captured the imagination of generations of mathematicians and the general public alike. Loosely speaking the classification theorem for surfaces (2-manifolds) tells us that any orientable surface is equivalent to a sphere with a certain number of "handles" sewn onto it. In this respect, the surface can be constructed from a sphere by topological surgery, which can be understood as a set of cutting, stitching and deforming operations.

Although a sound theoretical foundation of the subject matter has been laid out in algebraic topology, see e.g. (Munkres, 1984), localization of geometrically meaningful handles on surfaces remains a highly challenging topic. Over the last decades, several algorithms for detecting non trivial loops on general graphs have been proposed in computational geometry, see e.g. (Erickson, 2012) and the references therein, but most of these results remain theoretical and have not turned into efficient implementations. The topic of extracting such topological features is not only relevant as a theoretical pursuit but has fundamental applications in geometry processing, covering tasks such as mesh parametrization, mesh repair, and feature recognition, and spills beyond to fields such as biotechnology and bioinformatics, see e.g. (Brezovsky et al., 2013; Voss and Gerstein, 2010), therefore there is need for efficient practical solutions. The pioneering efforts of Dey and coworkers (2007; 2013) formalized the notion of geometrically meaningful topological features such as handle and tunnel loops and presented practical solutions using intermediate structures relying on persistence homology and Reeb graphs. Nonetheless, the cost of all intermediate constructs and their limitations impedes performance and intensifies memory usage.

Our goal is to provide a simple, intuitive, yet efficient strategy for detection of handle and tunnel loops. To this end, we re-examine the problem in the light of the humble loop shrinking property which lies behind Poincaré’s intuition and we seek to develop a variational analog to it which would allow us to evolve freely on the surface. In doing so, we avoid the difficulties encountered by existing methods in their efforts to marry homotopy classes, elaborate graph constructs, and practical geometric requirements.

Consider a person walking on a given surface while holding a sufficiently long thread from both ends, as illustrated in Figure 2. Had the person been on a sphere, she would be able to spool it back. On the other hand, if the person is on a torus, it won’t be able to re-spool the thread because it passes through the 2-dimensional hole of the torus.

Figure 2. Loop shrinking property illustrated on simple geometric objects.

To capture the essence of the loop shrinking property, we consider a diffusion process starting from an arbitrary location on the surface. If we are on the surface of a sphere, the advancing front will grow steadily but eventually it will start shrinking and vanish to one point similar to the thread losing hold. On the other hand, if we are on the surface of a torus, the growing patch will initially have a single boundary (Figure 3-left-top), and eventually, it will meet itself as it folds like a cannoli yielding a tunnel like region with two separate boundaries (Figure 3

-left-middle). At this point, we can only confirm that we are evolving on a tubular structure but we cannot infer the existence of a handle. Only at the moment when the two advancing fronts meet (Figure 

3-left-bottom), then a handle loop is detected. A tunnel loop can be obtained as a streamline by tracing back from the handle loop along the diffusion gradient, see Figure 4. Clearly here, the number of possible tunnel loops is large and all of them will pass through the initial point , so there is is no reason to expect them to exhibit some geometric optimality. This can be remedied by performing a diffusion initiated from the ring of the handle loop (Figure 3-left-bottom). Once the diffusion converges we can trace multiple streamlines sampled along the handle loop, thus obtaining a set of tunnel loops and we can then select the shortest.

This idea extends naturally to higher genus settings. For instance, let us consider the case of the double torus in Figure 3-right-top, the diffusion from the point marked on the surface will behave initially in a similar manner to the torus case, but as the two fronts merge (Figure 3-right-middle), the resulting single front will eventually split into two fronts (Figure 3-right-bottom) which would meet each other again yielding the second handle loop.

Figure 3. Handle loop discovery through splitting and collision of diffusion fronts on a simple torus (left) and a double torus (right).

Despite the inherent simplicity of the process described above, special attention needs to be paid to front tracking, bookkeeping, memory usage and performance in order to build a practical solution suitable for large data sets and complex surface layouts. A possible way forward is to use the level set method for carrying out the front tracking of the diffusion on the surface, however the sharp interface would require substantial bookkeeping and costly computations and checks on top of commonly known issues of the method pertaining to the distance function evaluation and re-initialization. Instead, we rely on the recent layered fields approach proposed as an alternative model to the ubiquitous Voronoi diagrams for modeling natural tessellation modeling (Zayer et al., 2018). We use a narrow band to represent the advancing front to avoid discontinuities across the front. We regard the propagation as field evolving on a layer on the surface. The splitting of the interface corresponds to the splitting of the field into two different fields evolving on different layers. For this purpose, we develop a dynamic approach to layer creation and layer collision which allows us to keep track of all the handles and tunnels on non trivial higher genus surfaces.

We represent the different layers as rows of a sparse matrix which both allows for carrying numerical diffusion computations and the tracking of topological features. As our geometry processing operations are channeled through linear algebra we take full advantage of parallel linear algebra primitives. In particular, our implementation operates fully on graphics hardware. A windfall of our approach is the ability to generate a Reeb graph simply by performing diffusion initiated at the detected handle loops. In this way, the commonly encountered shortcomings of working with height based morse functions are avoided in the first place.

In summary, this work makes the following contributions:

  • Novel variational abstraction for the loop shrinking property

  • Succinct diffuse interface Model for identifying handle and tunnel loops and capturing front propagation and collision

  • Dynamic sparse matrix representation for bookkeeping and processing

  • Fully parallel algorithmic realization on graphics hardware

  • Simple and intuitive Reeb graph extraction

Previous work

The body of work on extracting topological primitives on surfaces spans efforts across computational geometry and mesh processing. Theoretical algorithms for extracting non-trivial loops on surfaces have been proposed, e.g., (Erickson and Whittlesey, 2005; Kutz, 2006; de Verdière and Erickson, 2006)

. However there are no existing realizations of these theoretical efforts and no guarantees that the resulting loops are geometrically meaningful handles or tunnels. In the context of mesh simplification, several heuristics for identifying small handles have been proposed,

e.g.(El-Sana and Varshney, 1997; Guskov and Wood, 2001). The use of intermediate graphs for non-trivial loop identification has been explored in terms of Reeb graphs in (Shattuck and Leahy, 2001; Wood et al., 2004) and in terms of medial axis in (Zhou et al., 2007).

More closely related to our work are the efforts directly targeting the localization of handle and tunnel loops. Dey and coworkers addressed the problem by relying on persistence homology in (Dey et al., 2008), and then later on by making use of Reeb graphs (Dey et al., 2013). In the former, the intermediate tetrahedral tessellation required for carrying out homology computations poses several challenges both to feasibility and performance and bloats the problem size unnecessarily. Dropping persistence homology, in the latter, in favor of Reeb graphs improved both performance and scope, but computational cost and memory requirements remain considerable.

A data structure based on tree-cotree decomposition which allows for constructing generators of fundamental groups of surfaces was proposed by Eppstein (2003). Adjusting edge weights of the tree-cotree decomposition based on principal curvature directions, (Diaz-Gutierrez et al., 2009)

attempt to inject geometric meaning into those graph cycles by to steering them to align with those directions. These methods are prone to produce redundant non-trivial cycles which require post processing. An iterative approach alternating between principal directions to find “good” cycles was proposed in 

(Chen et al., 2018). It should be noted that curvature directions on general surfaces do not necessarily match the targeted cycle directions as the surface can be heavily decorated, see for instance, the dragons models in Figure 12 and this limits the scope to overly smooth meshes.

Without loss of generality, the above discussed methods do not operate directly on the D manifold but require one or multiple intermediate structures. Subsequently, performance and results quality depend on the quality of these underlaying structures. For instance, in (Dey et al., 2013), the Reeb graph quality depends among other thing on the choice of the initial direction as usually the associated morse function is some height function. Likewise, (Diaz-Gutierrez et al., 2009; Chen et al., 2018)

depend on principal directions estimations which are generally local, sensitive to surface fluctuations and noise. More importantly, the resulting algorithmic pipeline is not streamlined and therefore pose further challenges to code vectorization in view of the ubiquitous parallelism in modern hardware infrastructure.

2. Mathematical representation and data structure

General setting.

A majority of existing methods operate in the combinatorial manifold setting where the computation of shortest homology generators is relatively simpler than the more geometrically meaningful piecewise linear manifold. While this simplification offers certain algorithmic guarantees and computational advantages, the results depends largely on the quality of the underlying mesh edges. The challenges of working directly on the piecewise linear manifold have long been recognized, for instance, (Erickson and Whittlesey, 2005) discuss using level-sets for speeding up shortest path computations but acknowledge that “even the simpler problem of finding the shortest non-separating cycle in a piecewise-linear manifold appears to be open”.

Similar to (Dey et al., 2013), our work seeks to find geometrically meaningful generating pairs, meaningful is to be understood in the sense of shortest loops. In our illustrations, e.g., Figure 3, as in most common objects, the handle loop is way shorter than the tunnel loop (think a mug handle for instance), so for the sake of presentation flow, we presume in the ensuing discussion that a handle loops are found before their tunnel counterparts. Of course, it is possible to imagine scenarios where the opposite happens, e.g., deep screw holes in the drill in Figure 13) but this does not affect the pair extraction per se and the distinction between handle and tunnel loops will be discussed aside.

Diffuse interface.

While our approach is driven by diffusion, a pure diffusion alone is not sufficient for tracking different advancing fronts and resolving front collisions. In fact, it would lead simply to a smeared interface which would eventually end up as a flat curve (merged fronts). This issue gets further amplified on higher genus surfaces as multiple branchings yield multiple concurrent fronts. In order to resolve these issues and to avoid the numerical problems commonly encountered with sharp interface such as derivative discontinuities, we adopt a diffuse interface approach where the advancing front is not a sharp line as in the level set method but a narrow band.

In this way the growing patch can be defined by a function define dover the whole mesh and which takes value inside the patch and outside the cell and is smoothly varying on the narrow a narrow band where . In this way, it is easy to localize the narrow band by simple inspection of field values. Furthermore, at any time, it is possible to extract a sharp boundary as a level set of the narrow band.

Mathematical characterization.

In our context, front splitting due to branching, see e.g. Figure 3, leads to the creation of new separate regions and front collision requires colliding fronts to freeze thus blocking each other. These behavior patterns are similar to the natural tessellation model proposed in (Zayer et al., 2018) which describes the evolution of the interface of multiple cells growing simultaneously on a surface and yields Voronoi-like surface tessellations.

Within this setting, we can regard the diffusion of individual regions or cells as fields evolving on different layers and their interactions are governed by a set of requirements and weighted energy terms which drive the overall cell interactions:

  1. Partition of unity (the different fields sum up to at any given location on the surface).

  2. Integrity and locality of the individual cells in terms of the weighted product . The interaction does not cause a cell to disintegrate or break.

  3. Distinguishability of the individual narrow bands in terms of the weighted gradients . Narrow bands corresponding to each cell do not merge and flatten out.

  4. Behavior of narrow bands under contact which we model by a certain function which will be discussed shortly.

We account for the partition of unity constraint by introducing the Lagrange multiplier . We define the global energy of the model on the surface as a combination of three remaining terms discussed above. The ’s can be treated as independent variables, and the arising Lagrangian then reads


where the integral spans the surface and the summation is over the number of cells or regions .

Following (Zayer et al., 2018), the minimization of the Lagrangian associated with the energy arising from the above considerations yields the following time dependent equation which describes the evolution of the individual fields


As observed in the above equation, We do not need to define explicitly, we only need to define the difference term . The function encodes the desired behavior when two boundary bands meet each other. In the simple case, where only two bands meet, we have by virtue of partition of unity, . In practice , we seek a function that produces a symmetric behavior in the interval and vanishes at and and does not change sign. One such function is , in this way, we have

In the above equation, we can recognize the pure diffusion term characterized by the Laplacian which guarantees the smoothness of the growing interface. The first term in the equation, ensures a stable and stationary interface and therefore balances the effect of the diffusion term. The last term characterizes the driving force of the interface between the th and the th cell.

The growth rate is controlled by the mobility term . The solution can then be carried by a simple explicit Euler stepping scheme . In all our experiments we used a constant time step. For a discussion of the effect of the different time step the reader is referred to  (Zayer et al., 2018) or standard numerical methods textbooks.

Implementation details.

In our discretization on surface meshes, we uses the standard linear finite element approximation. The different fields can be thought of as individual layers on the surface and can be numerically encoded as a sparse matrix where each column represents a vertex of the original input mesh and each row represents a field , i.e. a region. This representation is convenient in the sense that we can precompute the Laplacian matrix of the mesh and at each time step the evaluation of the Laplacian of the fields amount to a sparse matrix matrix multiplication with . The Laplacian of is just the th row in the resulting product. This formulation offers a clear advantage in view of efficient parallel implementation. Furthermore, checking if two cells share a boundary amounts their corresponding rows in and looking for common nonzeros along the columns (vertices). To steer the field evolution, an additional base layer is initialized to at all vertices. When a region is initialized, we set the corresponding vertices inside to and to zero outside. To help enforce partition of unity, especially across borders, the field is normalized by dividing each column by the sum of its components. This normalization is further enforced after each time step. Furthermore, instead of a costly explicit tracking of cross-cell interactions at each iteration, we use the base layer for carrying this task implicitly. When two boundary bands move towards each other (or multiple bands at a junction), the base layer between them erodes (the corresponding vertices in are no longer ). With this in mind, both and the mobility term are set to zero except when either or index the base layer. We used and resp. in our current implementation. The values of the gradient energy coefficients are set to , and the penalty term to for , and to otherwise.

Figure 4. In the initial pass (top), the field interface (rim of the orange region) splits into two different components (top-left), both of which propagate till they reach each other (top-middle) yielding an initial handle estimate (top-right). In the second pass (middle), a diffusion initiated leftwards of the handle (middle-left) progresses till it reaches the opposite side of the handle (middle-middle). A streamline tracing against the diffusion gradient field yields the tunnel loop (middle-right). Similarly, in the third pass (bottom), a diffusion from the tunnel loop upwards (bottom-left) reaches the opposite side of the tunnel (bottom-middle) yielding a refined handle loop (bottom-right).

3. Variational Loop Shrinking

From the outline given in the introduction, our method can be regarded as sequence of three distinct passes. The initial pass registers the number of handle and tunnel pairs and gives an initial estimate for each handle loop. In the second pass, tunnel loops corresponding to each initial handle are created. In the third pass, the tunnel loop is used to generate a refined handle loop. Figure 4 depicts these three stages on the simple case of a torus. Throughout these passes, we mitigate the effect of the random seeding of the diffusion process in the initial pass and we improve the quality of the reported loops. If the focus is only on producing topologically correct loops without any geometric considerations, the final pass may be skipped and the estimates from the initial pass can be used as the handle loops.

Figure 5. Handle detection on the double torus during the initial pass. The creation of new colors, reflect dynamic changes to multi-layered field representation such as layer creation or or merging as the advancing fronts split or merge.

3.1. Initial Pass

The diffusion process is initialized from a random point on the surface and carried out by a simple explicit Euler stepping scheme. After each update, we check the state of the boundary of each advancing front (by inspecting nonzeros along the columns of which is much cheaper than computing isolines at each iteration. Isolines are computed only when a collision or branching are detected). We encounter a branching on the surface when there is more than one boundary loop (narrow band), see Figures 4-top-left and 5-left. In this case, each boundary loop is transferred to a new layer (row) in . As is modeled as a compressed sparse row (CSR) matrix, the transfer means we only have to adjust the column index for each vertex part of a specific connected component. The parent layer which gave birth to these new fronts can then be ignored in future iterations, because it has no space left to diffuse further, see e.g. Figure 3 for an illustration. Please note that for visualizations purposes, see Figures, 4 and 5

, the transfer is visualized in a vertex-only coloring and not interpolated, therefore it looks jagged; the transfer itself is fully integrated in the smooth diffusion process.

A handle is detected when multiple advancing fronts meet each other, Figures 4-top-middle and 5-second-left. This can be inferred from by inspecting which advancing fronts (rows) contribute field values to the same vertices (columns). Therefore, if a vertex has contributions from two or more layers, those layers should be merged. The number of handles detected is the number of layers involved minus one, e.g., two layers yield one handle. Inversely, to the scattering of multiple advancing fronts to multiple layers, now we gather all involved layers into a single new layer. This time the energy values per vertex of all involved layers are accumulated into the energy value of the last involved layer, clearing all the others. Afterwards the row index of the accumulation result is shifted to represent the new layer. This process avoids any expensive sparsity pattern changes and is fast to apply. The handle loops are formed by the involved advancing fronts before they are merged into a single layer, see Figures 4-top-right,  5-right, and 6. Since the number of handles equals the number of involved advancing fronts minus one, one boundary loop has to be ignored. Figure 4-top-right visualizes the formed handle loop as a green line. In all of our experiments, the longest boundary loop is ignored as this choice seems to give the best visual results. In general, this is not required and our final results are independent of this choice, as long as one layer is ignored.

Figure 6. Visualization of the progression of the initial pass on the ball mesh. Note that some handle loops are part of multiple tunnel loops. The red rectangle highlights and example where a naive tunnel loop pass will not report the correct tunnel. The book keeping challenges arising in scenarios like this one highlight the advantage of using a diffuse interface.

3.2. Tunnel loop pass

It is possible to generate tunnel loops by backtracing streamlines against the diffusion gradient in the initial pass. However this turns out to be a poor algorithmic choice as some tunnel loops may have to pass by the initial random diffusion seed. In order to produce well behaved tunnel loops we take a more judicious approach. We setup a diffusion process to start on one side of the handle loop whereas the other acts as a border preventing propagation in the opposite side. In this way, the diffusion process can get from one side of the handle loop to other, see Figure 4-middle. If the other side is reached, we can generate Streamlines back to the start points. The first point reached on the other side is used as the start point for the streamline trace. Due to the assumption that the diffusion happens at a constant speed on the surface, this will give us an approximation of the shortest path from one side to the other. This naive approach, however, does not work for complex geometries, because there may be multiple tunnel candidates for a handle estimate. The ball shown in figure 6 visualizes that problem. Whenever a handle estimate lies on a tubular connection between two holes, an example is highlighted by the red rectangle in figure 6, there will be two possibilities, one on each side, for a tunnel. Without further restricting the diffusion process we cannot guarantee which tunnel will be found by the Streamlines. Consequently, the same tunnel may be reported by different handle estimates, if they lie along the same tunnel. To fix this problem, we introduce an implicit ordering of the reported tunnels. The second diffusion process for each handle estimate is restricted to the surface area which was already covered at the time the handle estimate was detected. This ensures that a handle estimate finds the shortest tunnel fully covered by the initial diffusion process and no tunnel is reported more than once.

Figure 7. Examples for improved handle loops compared to the estimates of the initial pass. Estimates are rendered as green tubes, tunnel loops are in blue and improved handle loops as orange.

3.3. Handle refinement pass

There is a priori no reason for the handle loop generated in the initial pass to be optimally placed or shaped, for instance, perfectly round at the thinnest section of handle. Often times, it will be skewed or twisted, see examples in Figure 

7-top. By applying the same steps to generate the tunnel loop we can generate a refined handle loop. We start from one side of the tunnel loop and after we reach the other side, we trace backwards starting from the first point reached, see Figure 4-bottom. This time the diffusion process is restricted by the corresponding tunnel loop only. The streamline tracing will yield an improved, ideally placed handle loop, see Figure 7-bottom. Please note that although our approach is purely heuristic, it is well justified and avoids any costly numerical optimization

Figure 8. Detected handles and tunnel loops for different initializations. The purple point shows the start point of the diffusion process.

3.4. Robustness to initial seed placement

As discussed above, the multiple pass strategy makes our approach robust to the initial random seed placement. This is confirmed empirically in our test results. A typical scenario is shown in Figure 8. Starting from different locations on the surface of the mother and child model, our method generates nearly identical results. Please note that the handle loops, in orange, slightly differ. This is expected because there are multiple location having similar cross-section thickness alongside the arms. Therefore, there are multiple candidate locations which produce meaningful handle loop.

Figure 9. Distinguishing handle and tunnel loops (left) by offsetting them along their corresponding surface normals (right). As illustrated, tunnel loops tend to shrink whereas handle loops expand.

3.5. Distinguishing handle and tunnel loops

So far, we have assumed that loop handles are found first. This is not always the case, especially for deep holes as for instance in the drill in Figure 13. Our algorithmic flow does not depend on this distinction therefore, we can do identification at the end. Keeping to the simplicity of our overall approach, basic checks such as slightly offsetting the final loop along its corresponding surface normals and checking the changes in its chord length, as illustrated in Figure 9, was sufficient for making the distinction in our extensive test cases. Tunnel loops favor hyperbolic parts and will tend to shrink after offsetting whereas handle loops tend to expand. It might be possible to engineer scenarios where this basic distinction fails however we believe this a reasonable price for a tradeoff between fast practical solution and costly bullet proof robustness.

4. Performance Considerations

While it would have been easier to implement our current algorithmic pipeline in a multi-threaded fashion, our aim it to harness the tight memory requirement on the gpu and fine grained parallelism pertaining to modern graphics hardware. The two major challenges we faced are finding ways to parallelize the whole algorithm as well as keeping the communication and synchronization between CPU and GPU at a minimum.

As discussed earlier, the time stepping is dominated by Laplacian evaluation which is encoded as sparse matrix matrix multiplication, i.e. SpGEMM, and can readily take advantage of existing numerical kernels.

In all times, our method needs to keep an eye on the evolution of the advancing fronts. The extraction of the advancing fronts on a layer requires a list of all triangles which are part of the advancing fronts. A triangle is part of the advancing front, if one or two vertices of that triangle have an energy level above a given threshold and below . By applying this condition while iterating over all triangles in parallel, we atomically build the list of all triangles part of the advancing front.

Identifying the number of separate advancing fronts, while an easy serial exercise is not straightforward in parallel. It can then be modeled as a connected component labeling problem on an unstructured grid. Each marked triangle models a node in the grid and the triangle connectivity describes the edges in the grid. The connectivity can be extracted from the sparsity pattern of the Laplacian matrix used by the diffusion process. Connected component labeling is a challenging task to implement efficiently on the GPU. The nature of this problems requires scattered memory accesses as well as multiple passes, which both are problematic in terms of performance on the GPU. Our method uses the algorithm described in the work by Soman et al. (Soman et al., 2010), which minimizes the communication between CPU and GPU and reduces the scattered write operations to improve performance.

Figure 10. Geometrically embedded Reeb graphs on a selection of relevant test cases.
Figure 11. Reeb graph creation on the genus- deer model

The tunnel loop pass (subsection 3.2) and handle refinement pass (subsection 3.3) only differ in the initialization of the diffusion process and are discussed as one. Compared to the diffusion process during the initial pass, in these passes, the diffusion process has exactly one single layer. No scattering or merging is performed. This allows us to replace the sparse resizable system matrix by a dense vector. The dense vector contains the energy value for each vertex on the single layer. The base layer, required by the diffusion process, is implicitly modeled as one minus the energy value of the first layer. Furthermore, the SpGEMM during the diffusion update is replaced by a simple sparse matrix vector multiplication (SpMV). This performance optimization of the diffusion process is particularly important because each tunnel and handle loop requires a separate diffusion process. To fully utilize the GPU, multiple gradient and handle refinement passes run in parallel. Due to the constant memory consumption by the dense vector, we can calculate the exact number of possible parallel passes based on the available GPU memory.

Figure 12. Typical results of our approach on multiple data sets: From top left, The Pegasus, dancing children, chair, Grayloc, Thai statue, Dragon Phoenix, Dragon Ball, Metal Key, Casting, and Botijo Jar.

5. Reeb Graph Extraction

While our objective has been the detection of handle tunnel and loops, looking at Figure 5 it is clear that the diffusion process somewhat capture all the branching in the surface. While the topological information encoded by a function defined on given surface by tracking the connected components of it level sets is generally captured in the notion of Reeb graphs. Interest into Reeb graphs in computer graphics goes back to the work of (Shinagawa et al., 1991) which opened the door for a steady research effort on the topic. This includes performance oriented solutions, e.g. (Pascucci et al., 2007), mainstream implementations such as the on available from within the Topology Toolkit (Tierny et al., 2017) , improvements to the algorithmic complexity, eg. (Doraiswamy and Natarajan, 2009), and extensions and generalizations e.g. (Biasotti et al., 2000).

In this digression, our goal is not produce a state of the art Reeb graph algorithm but to show the versatility of our approach, in the sense, that by the same token we get the relevant handle and tunnel loops and meaningful geometric embedding of the Reeb graph associated with our diffusion on the surface.

Just using the diffusion initialized at random point as described in the subsection 3.1 would naturally yield topologically correct graph, however as we seek to associate a geometric embedding to the surface, we use the second pass from subsection 3.2 to improve the location of the junctions. Typical results of our approach are shown in Figure 10. Clearly for surfaces without any handle or tunnel loops it is sufficient to run the initial pass to track all the relevant branching, as illustrated in the case of the deer model in Figure 11 . The location of the connecting points across the geometric embedding are obtained by simply taking the mean of the narrow band.

An estimate of the performance can be inferred by comparing the cost of (initial and secondary pass) from table 1 to one of the top performing methods (Pascucci et al., 2007) since it is the method used in the Reeb graph step in (Dey et al., 2013).

Ours (Dey et al., 2013) Length ratio
Mesh(#faces, Genus) Total Initial Pass Tunnel Pass Handle Pass Total Reeb graph Map/Link Edge annot. Shorten. ours/(Dey et al., 2013)
Torus (3.1k, 1) 0.16 0.12 0.00 0.00 0.02 0.01 0 0 0.01 1.005
Double Torus (25k, 2) 0.48 0.42 0.00 0.00 0.18 0.05 0.01 0.02 0.1 1.010
MotherWithChild(28K,4) 0.42 0.29 0.08 0.04 0.52 0.05 0.05 0.02 0.40 0.988
Botijo(82K,5) 1.78 0.65 0.72 0.41 1.81 0.21 0.17 0.11 1.32 0.960
Casting(93K,9) 1.31 0.70 0.30 0.30 3.10 0.40 0.10 0.16 2.44 0.999
Happy Buddha(98K,8) 1.13 0.61 0.28 0.24 2.04 0.24 0.10 0.15 1.55 0.988
Ball (184K, 120) 17.65 13.01 1.97 1.16 401.77 0.28 34.76 3.16 363.57 0.918
Metal Key(390k,10) 10.48 6.87 2.11 1.50 36.52 10.35 9.52 0.52 16.13 0.984
Wooden Chair(400K,7) 7.49 4.45 2.28 0.77 31.21 2.77 18.55 0.43 9.46 0.982
Dragon Phoenix (750K, 11) 12.89 8.75 2.15 1.99 100.83 9.77 50.66 1.73 38.67 0.942
Grayloc (921k, 5) 22.25 14.23 4.11 3.92 DNF
V745 Sco nova (923K, 184) 273.99 104.21 86.32 82.34 DNF
Metal Table(950K, 198) 151.33 92.31 36.02 21.32 3,732.92 4.49 35.74 16.68 3,676.01 0.954
Dancing children(1.3M,8) 20.57 13.29 4.01 3.13 102.54 38.26 7.60 3.16 53.52 0.989
Drill (1.5m, 13) 18.77 13.92 2.67 2.18 90.32 13.52 1.14 2.74 72.92 0.846
Dragon Tamer(2M, 335) 429.57 254.93 93.90 78.63 >20,000.00 16.47 885.50 249.81 Crash
Dragon Ball (2.4M, 5, 40) 37.29 21.88 9.40 6.01 91.70 38.14 1.28 2.84 49.44 0.975
Napoleon (6.5M, 25) 705.34 376.50 170.86 157.99 >2,970.16 2730.67 239.49 Crash
Statue (10M, 3) 167.71 152.18 7.41 8.12 491.69 253.45 19.12 7.91 211.21 0.976
Table 1. Runtime comparison for various test models. Timings are measured in seconds. The length ratio compares the total handle and tunnel loops length of our method against (Dey et al., 2013).
Ours (Dey et al., 2013)
Mesh(#faces, Genus, Mem. Size) Inital Pass T/H. Pass T/H. Pass (#threads) Total
Torus (3.1k, 1, 0.05) 86 2 2 (1) 10
Double Torus (25k, 2, 0.4) 92 4 6 (2) 42
MotherWithChild (28K, 4, 0.48) 92 6 18 (4) 76
Botijo (82K, 5, 1.5) 108 64 108 (5) 193
Casting(93K,9, 1.6) 114 18 134 (9) 210
Happy Buddha(98K,8,1.8) 97 16 112 (8) 189
Ball (184K, 120, 3.15) 152 152 458 (10) 2135
Metal Key (390k,10, 6.6) 210 82 552 (10) 810
Wooden Chair (400K,7, 6.8) 214 82 472 (7) 921
Dragon Phoenix (750K, 11, 12) 340 106 1062 (10) 1540
Grayloc (921k, 5, 16) 366 354 1155 (5) 0
V745 Sco nova (923K, 184, 16) 502 187 1402 (10) DNF
Metal Table (950K,198, 16) 336 166 882 (10) 16891
Dancing children (1.3M,8, 21) 462 240 612 (8) 4992
Drill (1.5m, 13, 25) 496 338 1844 (13) 4142
Dragon Tamer (2M, 335, 34) 636 290 1410 (10) >32562
Dragon Ball (2.4M, 5, 40) 728 292 1452 (5) 4557
Napoleon (6.5M, 25, 111) 1844 1566 7744 (10) >15000
Statue (10M, 3, 171) 496 338 1844 (13) 4142
Table 2. Memory consumption for various test models. For our method peak GPU memory usage is reported for the initial pass and the tunnel and handle (H./T.) pass combined. For the tunnel and handle passes, we report the memory consumption with and without running multiple tunnel/handle pairs in parallel. The initial pass and the gradients pass do run one after another, therefore, their maximum is the maximum of our whole method. For (Dey et al., 2013) peak CPU memory consumption is reported.

6. Results

All our experiment were carried on a Intel(R) Core(TM) i7-7700 with 32Gb system memory and a Nvidia RTX 2080ti 11gb as GPU. As initialization point for the initial diffusion process we used the first vertex of the mesh.

Figure 12 shows a collection of non-trivial test case meshes featuring: large loops as in the pegasus model and the dancing children, object with thin features such as the key and the chairs and the mechanical parts. In all these cases, our approaches detect the handles and tunnels correctly. The handles and tunnels in the dragon models are heavily decorated with ornaments. The reported tunnel loops correctly avoid following the ornaments which face inwards compared to the tunnel loop. Following those ornaments would increase the length of the loop. The handle loops are located at the smallest parts despite the ornaments.

Figure 13 shows the tunnel and handles loops computed for the napoleon mesh. This mesh features very large tunnels, e.g. the tunnel covering the stone plate, two legs and the body of the horse, and also very tiny ones, e.g. at the top part of the foot rest or at the reins.

The drill shown in Figure 13 features very small tunnels combined with large handles. All tunnels and handles are reported correctly, because our method is not based on any presumptions on the length of handle or tunnel loops. Note that not all screw holes feature a tunnel in the mesh, since it is a 3D scan.

Figure 13. Evaluation of our approach on the Napoleon model (with zoomed views) and the Drill model (double sided view).

Table 1 compares the runtime of our method with the state of the art approach (Dey et al., 2013). The results show that our method is slower for relatively small meshes but scales significantly better with larger meshes. This can be explained by the constant overhead of GPU initialization as well as memory transfers to and from the GPU to obtain the results. Our experiments included five meshes with a high genus: ball, metal table, V745 Sco nova and dragon tamer. The results of these meshes are visualized in figure 14 and figure 1. In these cases our method outperformed the other method by over an order of magnitude. V745 Sco nova and dragon tamer, which has the highest genus of all tested meshes, did not complete with their method. Also the metal table required multiple runs until the method completed without an error, we report only the runtime of the successful run.

The number of pairs detected by both methods is the same. As it is cumbersome to match the pairs of both methods automatically for comparison, we compare the total length of all loops. The length ratio in Table 1 suggest that our method yields shorter loops in general. Please not that the torus and the double torus are regularly meshed models where the optimal tunnel and handles loops are well aligned with the mesh edges. Since the method of  (Dey et al., 2013) performs the shortening in combinatorial fashion it lands directly on those edges. Still the difference with our method is not significant.

Table 2 contains the memory consumptions obtained during our experiments. Again, except for very small meshes, our method uses less memory even when multiple handle and tunnel refinement passes run in parallel. We report the minimum memory requirement for the handle and tunnel pass as well as the maximum when running multiple passes in parallel. The memory consumption does not directly scale by the number of parallel runs because some read only information, such as the mesh, can be shared between these passes. Our initial tests showed that running more than 10 passes in parallel did not improve the performance further. This is probably caused by driver overhead scheduling the parallel kernel runs on the GPU.

Figure 14. Very high genus test cases, a complicated ball, a metal table, and a model of the V745 Sco nova (courtesy of NASA).

Limitations and discussion.

The performance of (Dey et al., 2013) depends largely on the fast Reeb graph computations of  (Pascucci et al., 2007) which in turn uses “free” height functions (coordinates) as Morse functions. Certainly, other functions can be used, the “streaming meshes” format has been used in the latter but then the cost of the underlying Fiedler vector computation should be factored in. Similarly, when geodesic or harmonic functions, e.g., (Hilaga et al., 2001) are used, their cost will impede performance. Therefore the comparisons we are providing are representative.

Our approach operates on the piecewise linear manifold setting, and cannot offer the guarantees of well and long-established combinatorial approaches but empirical evidence shows that our method produces the same number of pairs as [Dey2013]. This is expected since the method relies on first principals in topology embodied in the loop shrinking property. Our method builds on the ability to conduct standard numerical simulations on surface meshes and is therefore bound by fairly well-known requirements on methods such as the finite element method. We do not see this as a limitation but rather a motivation for exploring common grounds between combinatorial and variational approaches.

There can be scenarios where tiny and coarsely meshed handles can affect the quality of our results or even escape the width of our advancing front. The only test case where we encountered such a scenario is depicted in Figure 15. Although the handle is only spanned by three triangles, our method still captures it but the loop is clearly suboptimal. Although this does not respect the smoothness assumptions of the linear manifold required in our variational setting, it can be addressed by adaptive refinement along the front.

Figure 15. The Happy Buddha model (left) and a zoom-in on the area marked by the arrow showing the effect of a tiny coarse handle on the quality of our result.

7. Conclusion

We have re-abstracted the shrinking loop property as a continuous growth process steered by diffusion. The diffusion field was modeled using diffuse interface and the dynamics of the advancing fronts were captured by assuming a multilayered field representation where the creation and and merging of layers are respectively steered by front splitting and collision. Our approach is simple and versatile allowing the detection of handle and tunnel loops as well as the creation of Reeb graphs. Furthermore the time stepping nature of our approach allows for a close control of the whole process by simple adjustment of the field and regions of diffusion. We foresee that this control will help extend our approach to other application such topological surgery, mesh segmentation, and rigging. We hope that adopting our methodology will help break the deadlock in performance and bring new machinery to help understanding and harnessing geometric problems.


  • S. Biasotti, B. Falcidieno, and M. Spagnuolo (2000) Extended reeb graphs for surface understanding and description. In Discrete Geometry for Computer Imagery, G. Borgefors, I. Nyström, and G. S. di Baja (Eds.), Berlin, Heidelberg, pp. 185–197. External Links: ISBN 978-3-540-44438-1 Cited by: §5.
  • J. Brezovsky, E. Chovancova, A. Gora, A. Pavelka, L. Biedermannova, and J. Damborsky (2013) Software tools for identification, visualization and analysis of protein tunnels and channels. Biotechnology advances 31 (1), pp. 38–49. Cited by: §1.
  • J. Chen, J. Jester, and M. Gopi (2018) Fast computation of tunnels in corneal collagen structure. In Proceedings of Computer Graphics International 2018, CGI 2018, New York, NY, USA, pp. 57–65. External Links: ISBN 9781450364010, Link, Document Cited by: §1, §1.
  • É. C. de Verdière and J. Erickson (2006) Tightening non-simple paths and cycles on surfaces. In Proceedings of the Seventeenth Annual ACM-SIAM Symposium on Discrete Algorithm, SODA ’06, USA, pp. 192–201. External Links: ISBN 0898716055 Cited by: §1.
  • T. K. Dey, F. Fan, and Y. Wang (2013) An efficient computation of handle and tunnel loops via reeb graphs. ACM Trans. Graph. 32 (4). External Links: ISSN 0730-0301, Link, Document Cited by: §1, §1, §1, §2, Table 1, Table 2, §5, §6, §6, §6.
  • T. K. Dey, K. Li, J. Sun, and D. Cohen-Steiner (2008) Computing geometry-aware handle and tunnel loops in 3d models. ACM Trans. Graph. 27 (3), pp. 1–9. External Links: ISSN 0730-0301, Link, Document Cited by: §1.
  • T. K. Dey, K. Li, and J. Sun (2007) On computing handle and tunnel loops. In Proceedings of the 2007 International Conference on Cyberworlds, CW ’07, USA, pp. 357–366. External Links: ISBN 0769530052 Cited by: §1.
  • P. Diaz-Gutierrez, D. Eppstein, and M. Gopi (2009) Curvature Aware Fundamental Cycles. Computer Graphics Forum. External Links: ISSN 1467-8659, Document Cited by: §1, §1.
  • H. Doraiswamy and V. Natarajan (2009) Efficient algorithms for computing reeb graphs. Comput. Geom. Theory Appl. 42 (6–7), pp. 606–616. External Links: ISSN 0925-7721, Link, Document Cited by: §5.
  • J. El-Sana and A. Varshney (1997) Controlled simplification of genus for polygonal models. In Proceedings of the 8th Conference on Visualization ’97, VIS ’97, Washington, DC, USA, pp. 403–ff.. External Links: ISBN 1581130112 Cited by: §1.
  • D. Eppstein (2003) Dynamic generators of topologically embedded graphs. In Proceedings of the Fourteenth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’03, USA, pp. 599–608. External Links: ISBN 0898715385 Cited by: §1.
  • J. Erickson and K. Whittlesey (2005) Greedy optimal homotopy and homology generators. In Proceedings of the Sixteenth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’05, USA, pp. 1038–1046. External Links: ISBN 0898715857 Cited by: §1, §2.
  • J. Erickson (2012) Combinatorial optimization of cycles and bases. Advances in Applied and Computational Topology 70, pp. 195–228. Cited by: §1.
  • I. Guskov and Z. J. Wood (2001) Topological noise removal. In Proceedings of Graphics Interface 2001, GI ’01, CAN, pp. 19–26. External Links: ISBN 0968880800 Cited by: §1.
  • M. Hilaga, Y. Shinagawa, T. Komura, and T. Kunii (2001) Topology matching for fully automatic similarity estimation of 3d shapes. pp. 203–212. External Links: Document Cited by: §6.
  • M. Kutz (2006) Computing shortest non-trivial cycles on orientable surfaces of bounded genus in almost linear time. In Proceedings of the Twenty-Second Annual Symposium on Computational Geometry, SCG ’06, New York, NY, USA, pp. 430–438. External Links: ISBN 1595933409, Link, Document Cited by: §1.
  • J.R. Munkres (1984) Elements of algebraic topology. Advanced book classics, Perseus Publishing. External Links: ISBN 9780201054873, LCCN lc84006250 Cited by: §1.
  • V. Pascucci, G. Scorzelli, P. Bremer, and A. Mascarenhas (2007) Robust on-line computation of reeb graphs: simplicity and speed. In ACM SIGGRAPH 2007 Papers, SIGGRAPH ’07, New York, NY, USA. External Links: Link, Document Cited by: §5, §5, §6.
  • D. W. Shattuck and R. M. Leahy (2001) Automated graph-based analysis and correction of cortical volume topology. IEEE Transactions on Medical Imaging 20 (11), pp. 1167–1177. External Links: Document Cited by: §1.
  • Y. Shinagawa, T. L. Kunii, and Y. L. Kergosien (1991) Surface coding based on morse theory. IEEE Computer Graphics and Applications 11 (5), pp. 66–78. External Links: Document, ISSN 1558-1756 Cited by: §5.
  • J. Soman, K. Kothapalli, and P. J. Narayanan (2010) Some gpu algorithms for graph connected components and spanning tree. Parallel Processing Letters 20 (04), pp. 325–339. External Links: Document, Link, https://doi.org/10.1142/S0129626410000272 Cited by: §4.
  • J. Tierny, G. Favelier, J. A. Levine, C. Gueunet, and M. Michaux (2017) The Topology ToolKit. IEEE Transactions on Visualization and Computer Graphics (Proc. of IEEE VIS). Cited by: §5.
  • N. R. Voss and M. Gerstein (2010) 3V: cavity, channel and cleft volume calculator and extractor. Nucleic acids research 38 (suppl_2), pp. W555–W562. Cited by: §1.
  • Z. Wood, H. Hoppe, M. Desbrun, and P. Schröder (2004) Removing excess topology from isosurfaces. ACM Trans. Graph. 23 (2), pp. 190–208. External Links: ISSN 0730-0301, Link, Document Cited by: §1.
  • R. Zayer, D. Mlakar, M. Steinberger, and H. Seidel (2018) Layered fields for natural tessellations on surfaces. ACM Trans. Graph. 37 (6). External Links: ISSN 0730-0301, Link, Document Cited by: §1, §2, §2, §2.
  • Q. Zhou, T. Ju, and S. Hu (2007) Topology repair of solid models using skeletons. IEEE Transactions on Visualization and Computer Graphics 13 (4), pp. 675–685. External Links: ISSN 1077-2626, Link, Document Cited by: §1.