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 (2manifolds) 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 reexamine 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 respool the thread because it passes through the 2dimensional hole of the torus.
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 3lefttop), and eventually, it will meet itself as it folds like a cannoli yielding a tunnel like region with two separate boundaries (Figure 3
leftmiddle). 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
3leftbottom), 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 3leftbottom). 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 3righttop, 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 3rightmiddle), the resulting single front will eventually split into two fronts (Figure 3rightbottom) which would meet each other again yielding the second handle loop.
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 reinitialization. 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 nontrivial 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., (ElSana and Varshney, 1997; Guskov and Wood, 2001). The use of intermediate graphs for nontrivial 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 treecotree decomposition which allows for constructing generators of fundamental groups of surfaces was proposed by Eppstein (2003). Adjusting edge weights of the treecotree decomposition based on principal curvature directions, (DiazGutierrez 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 nontrivial 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, (DiazGutierrez 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 levelsets for speeding up shortest path computations but acknowledge that “even the simpler problem of finding the shortest nonseparating cycle in a piecewiselinear 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 Voronoilike 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:

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

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

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

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
(1) 
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
(2) 
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 crosscell 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.
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.
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 4topleft and 5left. 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 vertexonly 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 4topmiddle and 5secondleft. 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 4topright, 5right, and 6. Since the number of handles equals the number of involved advancing fronts minus one, one boundary loop has to be ignored. Figure 4topright 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.
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 4middle. 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.
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
7top. 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 4bottom. 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 7bottom. Please note that although our approach is purely heuristic, it is well justified and avoids any costly numerical optimization3.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 crosssection thickness alongside the arms. Therefore, there are multiple candidate locations which produce meaningful handle loop.
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 multithreaded 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.
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.
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 
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 
6. Results
All our experiment were carried on a Intel(R) Core(TM) i77700 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 nontrivial 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.
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.
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 longestablished 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 wellknown 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.
7. Conclusion
We have reabstracted 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.
References
 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 9783540444381 Cited by: §5.
 Software tools for identification, visualization and analysis of protein tunnels and channels. Biotechnology advances 31 (1), pp. 38–49. Cited by: §1.
 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.
 Tightening nonsimple paths and cycles on surfaces. In Proceedings of the Seventeenth Annual ACMSIAM Symposium on Discrete Algorithm, SODA ’06, USA, pp. 192–201. External Links: ISBN 0898716055 Cited by: §1.
 An efficient computation of handle and tunnel loops via reeb graphs. ACM Trans. Graph. 32 (4). External Links: ISSN 07300301, Link, Document Cited by: §1, §1, §1, §2, Table 1, Table 2, §5, §6, §6, §6.
 Computing geometryaware handle and tunnel loops in 3d models. ACM Trans. Graph. 27 (3), pp. 1–9. External Links: ISSN 07300301, Link, Document Cited by: §1.
 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.
 Curvature Aware Fundamental Cycles. Computer Graphics Forum. External Links: ISSN 14678659, Document Cited by: §1, §1.
 Efficient algorithms for computing reeb graphs. Comput. Geom. Theory Appl. 42 (6–7), pp. 606–616. External Links: ISSN 09257721, Link, Document Cited by: §5.
 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.
 Dynamic generators of topologically embedded graphs. In Proceedings of the Fourteenth Annual ACMSIAM Symposium on Discrete Algorithms, SODA ’03, USA, pp. 599–608. External Links: ISBN 0898715385 Cited by: §1.
 Greedy optimal homotopy and homology generators. In Proceedings of the Sixteenth Annual ACMSIAM Symposium on Discrete Algorithms, SODA ’05, USA, pp. 1038–1046. External Links: ISBN 0898715857 Cited by: §1, §2.
 Combinatorial optimization of cycles and bases. Advances in Applied and Computational Topology 70, pp. 195–228. Cited by: §1.
 Topological noise removal. In Proceedings of Graphics Interface 2001, GI ’01, CAN, pp. 19–26. External Links: ISBN 0968880800 Cited by: §1.
 Topology matching for fully automatic similarity estimation of 3d shapes. pp. 203–212. External Links: Document Cited by: §6.
 Computing shortest nontrivial cycles on orientable surfaces of bounded genus in almost linear time. In Proceedings of the TwentySecond Annual Symposium on Computational Geometry, SCG ’06, New York, NY, USA, pp. 430–438. External Links: ISBN 1595933409, Link, Document Cited by: §1.
 Elements of algebraic topology. Advanced book classics, Perseus Publishing. External Links: ISBN 9780201054873, LCCN lc84006250 Cited by: §1.
 Robust online 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.
 Automated graphbased analysis and correction of cortical volume topology. IEEE Transactions on Medical Imaging 20 (11), pp. 1167–1177. External Links: Document Cited by: §1.
 Surface coding based on morse theory. IEEE Computer Graphics and Applications 11 (5), pp. 66–78. External Links: Document, ISSN 15581756 Cited by: §5.
 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.
 The Topology ToolKit. IEEE Transactions on Visualization and Computer Graphics (Proc. of IEEE VIS). Cited by: §5.
 3V: cavity, channel and cleft volume calculator and extractor. Nucleic acids research 38 (suppl_2), pp. W555–W562. Cited by: §1.
 Removing excess topology from isosurfaces. ACM Trans. Graph. 23 (2), pp. 190–208. External Links: ISSN 07300301, Link, Document Cited by: §1.
 Layered fields for natural tessellations on surfaces. ACM Trans. Graph. 37 (6). External Links: ISSN 07300301, Link, Document Cited by: §1, §2, §2, §2.
 Topology repair of solid models using skeletons. IEEE Transactions on Visualization and Computer Graphics 13 (4), pp. 675–685. External Links: ISSN 10772626, Link, Document Cited by: §1.
Comments
There are no comments yet.