|
CGAL 5.3 - Surface Mesh Topology
|
|
CGAL 5.3 - Surface Mesh Topology
|
This package provides a toolbox for manipulating curves on a combinatorial surface from the topological point of view. Two main functionalities are proposed. One is the computation of shortest curves that cannot be continuously deformed to a point. This includes the computation of the so-called edge width and face width of the vertex-edge graph of a combinatorial surface. The other functionality is the homotopy test for deciding if two given curves on a combinatorial surface can be continuously deformed one into the other.
All the computations in this package, either for the shortest non-contractible curve or for the homotopy tests, are performed on a input surface represented as a model of CombinatorialMap or any model of FaceGraph. Note that combinatorial maps are based on darts and FaceGraphs are based on halfedges. To avoid repetitions we use the terms darts and halfedges interchangeably in the sequel. The input surface is supposed to be connected and orientable.
Given a combinatorial surface, one may consider combinatorial curves, which are described as sequences of edges, or topological curves, which are continuous curves on the topological surface underlying the combinatorial surface. The length of a combinatorial curve is the sum of the lengths of its edges. Here, we measure the length of a topological curve as the number of crossings of this curve with the vertex-edge graph of the combinatorial surface. A closed curve, either topological or combinatorial, that cannot be continuously deformed to a point on the topological surface is said non-contractible. This package offers the following functions:
dh, compute a shortest non-contractible combinatorial curve passing through the source vertex of dh,It is important to clarify how we compare the lengths of two combinatorial curves in order to compute the shortest one. "Shortest" can be understood as "having the least amount of edges" or "having the smallest total length of its edges". In the former case, we consider that the mesh is unweighted; in the latter case, we consider that the mesh is weighted, and one must specify how the weight, or length, of each edge is calculated (see concept WeightFunctor). When the vertices of the mesh have Euclidean coordinates, the Euclidean distance between two connected vertices defines a natural weight for the corresponding edge. A weight functor Euclidean_length_weight_functor is provided for this purpose.
The algorithm to find a shortest non-contractible curve through a specified vertex is based on the paper by Cabello et al. [1]. The time complexity is linear, though in the weighted case it is raised by a logarithmic factor, assuming that the weight computation takes constant time per edge. Computing the edge width takes quadratic time by running the first function on each vertex, and its complexity is also raised by a logarithmic factor when considering a weighted map. Computing the face width consists of constructing the radial graph of the original mesh and computing the edge width of the radial graph. It thus takes quadratic time. Computing face width on weighted maps is currently not supported.
Given a curve drawn on a surface one can ask if the curve can be continuously deformed to a point (i.e. a zero length curve). In other words, does there exist a continuous sequence of curves on the surface that starts with the input curve and ends to a point? Curves that deform to a point are said contractible. Any curve on a sphere is contractible but this is not true for all curves on a torus or on a surface with more complicated topology. The algorithms in this section are purely topological and do not assume any geometry on the input surface. In particular, the surface is not necessarily embedded in a Euclidean space.
The algorithm implemented in this package builds a data structure to efficiently answer queries of the following forms:
The second query asks if the curves are freely homotopic while the third one asks if the curves are homotopic with fixed endpoints. The three queries are globally referred to as homotopy tests. Figure 83.1 below illustrates the three types of queries.
The algorithms used are based on a paper by Erickson and Whittlesey [3], providing a linear time algorithm for the above homotopy tests. This is a simplified version of the linear time algorithm by Lazarus and Rivaud [4].
Given a cycle drawn on a surface one can ask if the cycle can be continously deformed to a cycle that does not intersect with itself. Any contractible cycle deforms to a simple cycle but this is not true for more complicated cycles. The algorithm in this section is purely topological and do not assume any geometry on the input surface.
The algorithm implemented in this package builds a data structure to efficiently answer queries of the following forms:
The algorithm used is based on a paper by Despré and Lazarus [2], providing a \(O(n + l\log{l})\)-time algorithm where \(n\) is the complexity of \(\cal{M}\) and \(l\) is the length of the path.
The main class for this package is Surface_mesh_topology::Curves_on_surface_topology. Its constructor takes the input surface. An internal representation of the surface (described below) is computed the first time an homotopy test is called.
Each combinatorial curve on this surface is contained in an instance of the class Surface_mesh_topology::Path_on_surface. An object in this class behaves as a list. This list is initially empty and the halfedges corresponding to the sequence of consecutive oriented edges of an input curve should be pushed back in this list. The class provides four ways for extending a nonempty path.
Path_on_surface (see the description of `can_be_pushed()` below) This can be done even when the path is empty,h, extending by a positive one turn is thus equivalent to extend the path with next(h). An analogous member function extend_negative_turn() is provided for convenience,PolygonalSchema, which is a model of GenericMap with labeled edges as explained in section Polygonal Schema Helper, the user may push the label of the next halfedge instead of the halfedge itself with the member function push_back_by_label().In the first two cases, let A be the source vertex of the added dart or the target vertex if the added dart is flipped, let B be the target vertex of the last dart in the path or the source vertex if the last dart is flipped: A and B should coincide. The user is responsible for ensuring this condition. The member functions can_be_pushed(), can_be_pushed_by_index() and can_be_pushed_by_label() return true if and only if the condition is satisfied.
Specifying a path on a combinatorial surface might be a tedious task. Indeed, knowing in advance the pointer, index or turn of each consecutive halfedge in a path is not always easy. In order to facilitate this task, we provide an intuitive model of CombinatorialMap called Surface_mesh_topology::Polygonal_schema_with_combinatorial_map, a model of the PolygonalSchema concept. In this model, a surface is viewed as a collection of clockwise oriented polygonal facets with labeled boundary (oriented) edges. Boundary edges with the same label are glued together to form a surface. Each label should appear at most twice in the collection and a label that appears only once corresponds to a boundary edge. The label of the opposite of an oriented edge is preceded by a minus. For example, the opposite of 'a1' is '-a1'. Since we are dealing with orientable surfaces only, each label that appears twice must appear once with a minus. The user can add facets to the surface one at a time. Each facet is specified by the sequence of its oriented edge labels given as a string where the labels are words (any sequence of characters, except space) separated by blank spaces. In the next figure we see three examples of combinatorial maps described by a collection of facets with labeled edges.
As noted in the previous section Polygonal Schema, every label that appears only once in a polygonal schema corresponds to a boundary edge. The Polygonal Schema helper offers another mechanism to create surfaces with boundary. Each facet already added to a polygonal schema may be perforated to create a hole in the surface. The edges of a perforated facet thus becomes boundary edges. The reverse operation consists in filling the supposedly perforated facet. The class provides the following interface.
h. If s is the label of the oriented edge corresponding to h, one may equivalently call perforate_facet (s).As an example, one may perforate all the facets of a polygonal schema \(\cal{M}\) to obtain a "skeleton" surface equivalent to a thickening of the graph composed of the edges of \(\cal{M}\).
A Surface_mesh_topology::Curves_on_surface_topology can be constructed with a Surface_mesh_topology::Polygonal_schema_with_combinatorial_map as input surface. In this case, every halfedge has a label (possibly preceded by a minus) and a path can be specified by the sequence of labels corresponding to its halfedge sequence. A repeated call to the function push_back_by_label() allows the user to specify the path in this way.
Since the data structures to represent a surface are edge-centralized, in order to specify a vertex where the curve is computed, the user can use any dart belonging to this vertex. We use the term cycle as a synonymous of closed curve.
The class Curves_on_surface_topology provides the following three functions:
compute_shortest_non_contractible_cycle_with_base_point(dh, weight_functor) : Compute a shortest non-contractible cycle going through the source vertex of dh,compute_shortest_non_contractible_cycle(weight_functor) : Very similar to the previous function, except that one does not specify a vertex. It computes a shortest non-contractible cycle through every vertex and returns the shortest cycle among them,compute_edge_width() : Compute the edge width of the mesh, equivalent to compute_shortest_non_contractible_cycle( Unit_weight_functor()).The above functions return an instance of Path_on_surface . The optional argument weight_functor is used to calculate the length of the edges. If not given, all the edge lengths are set to 1, i.e. the mesh is unweighted.
Face width returns a shortest non-contractible topological curve described as a circular sequence of traversed faces alternating with the vertices it passes through. Each face or vertex in this sequence is identified by one of its dart handles. In practice, we choose the same dart handle for a face as for the next vertex it passes through. This way, we only need to return one dart handle for a face and its following vertex in the sequence.
The function compute_face_width() computes the sequence of dart handles as described above and returns an std::vector of dart handles, where each dart represents a traversed face followed by an incident vertex.
Given two Surface_mesh_topology::Path_on_surface \(p_1\) and \(p_2\), the class Surface_mesh_topology::Curves_on_surface_topology provides the following three functions:
true if the closed curve \(p_1\) is contractible,true if the closed curves \(p_1\) and \(p_2\) are freely homotopic,true if the paths \(p_1\) and \(p_2\) are homotopic with fixed endpoints. This call is equivalent to is_contractible( \(p_1\cdot \overline{p_2}\)), where \(p_1\cdot \overline{p_2}\) is the concatenation of \(p_1\) and the reverse of \(p_2\).A common first step in the homotopy test algorithms is to simplify the input combinatorial surface. This preprocessing step is done once and for all for a given mesh, the first time an homotopy test is called. The simplified surface is a quadrangulation, every face of which is a quadrilateral, stored in a Surface_mesh_topology::Curves_on_surface_topology. It has 2 vertices and \(2g\) quadrilaterals where \(g\) is the genus of the input surface. This is otherwise independent of the size of input surface,
Surface_mesh_topology::Curves_on_surface_topology.Each time a Surface_mesh_topology::Path_on_surface is provided for a homotopy test, it is first transformed to an equivalent path in the quadrangulation stored by the Surface_mesh_topology::Curves_on_surface_topology. This transformation is transparent to the user who has never access to the quadrangulation.
Given a Surface_mesh_topology::Path_on_surface \(p\), the class Surface_mesh_topology::Curves_on_surface_topology provides the following function:
is_homotopic_to_simple_cycle(p) returns true if the closed curve \(p\) is homotopic to some simple cycle.Like homotopy tests, the first step is to simplify the input combinatorial surface. The algorithm will share the surface with homotopy tests and invoke the simplification if the preprocessing has not been done yet.
Surface_mesh_topology::Curves_on_surface_topology.Each time a Surface_mesh_topology::Path_on_surface is provided for a simplicity test, it is first transformed to an equivalent path in the quadrangulation stored by the Surface_mesh_topology::Curves_on_surface_topology. This transformation is transparent to the user who has never access to the quadrangulation.
In the next two examples, we present various ways to compute shortest non-contractible cycles.
One can store the original mesh in a Combinatorial_map instance and run the algorithm without regarding the geometric distances, i.e. the unweighted case (first call to compute_shortest_non_contractible_cycle_with_base_point). Alternatively, one can take the geometric distances into consideration by providing a weight functor to calculate the weight of the edge containing the given dart (second call to compute_shortest_non_contractible_cycle_with_base_point). Note that the time complexity is raised by a logarithmic factor.
File Surface_mesh_topology/shortest_noncontractible_cycle.cpp
In order to find the edge width of the surface, one can make use of the routine compute_edge_width as in the following example. The weighted shortest non contractible cycle is also computed (calling compute_shortest_non_contractible_cycle). In this example, a Surface_mesh is used to store the mesh.
File Surface_mesh_topology/edgewidth_surface_mesh.cpp
In these two examples, the mesh and the cycles can be visualized if CGAL_Qt5 is enabled.
The following example computes the face width, and visualizes it if CGAL_Qt5 is enabled.
File Surface_mesh_topology/facewidth.cpp
The following example shows how to load an off file and how to create three closed paths on this surface. Contractibility and free homotopy tests are then performed. The example also shows how to use the CGAL viewer if CGAL_Qt5 is enabled.
File Surface_mesh_topology/path_homotopy_double_torus.cpp
The following example shows how to test the simplicity of a closed path on a double torus. The original path is visualized if CGAL_Qt5 is enabled.
File Surface_mesh_topology/path_simplicity_double_torus_2.cpp
Here, we show with two examples how to create a surface from a list of faces specified by edge label sequences. In this first example, we build a genus two torus surface from a single face, also called a polygonal schema. See left Figure 83.2 for an illustration. Two closed paths are then created. The paths are freely homotopic but not homotopic with fixed endpoint.
File Surface_mesh_topology/path_homotopy_with_symbols.cpp
In this second example, we build a genus two torus surface from a set of three squares. See middle Figure 83.2 for an illustration. The first two faces are added each with a single call to the member function add_facet(). The third face is build incrementally by adding its edge labels one at a time. We then create a contractible closed path.
File Surface_mesh_topology/path_homotopy_with_symbols_2.cpp
In this third example, we create non closed paths on the same mesh as in the first example and perform homotopy tests with fixed endpoints. Here, a Surface_mesh is used as an alternative to a CombinatorialMap.
File Surface_mesh_topology/open_path_homotopy.cpp
The machine used is a PC running Ubuntu 18.04 with an Intel CPU Core i7-4790 CPU clocked at 3.60GHz with 32GB of RAM.
The first time an homotopy test is called, we build a special quadrangulation of the surface as internal representation (as explained in Section Implementation Details). The complexity of this operation is linear in the number of darts of the input surface, as we can see in Figure 83.3.
For this benchmark, we computed 38 Surface_mesh_topology::Curves_on_surface_topology objects for different input surfaces with different number of darts (between 8,000 and 30,000,000). We show in the figure the computation time of the quadrangulation according to the number of darts of the input surface. We remind that this computation needs be done only once if you want to perform several path homotopy tests on the same surface.
Surface_mesh_topology::Curves_on_surface_topology constructions, according to the number of darts or the input surface. In this second benchmark, we use two surfaces as input (with respectively 543,652 vertices, 1,631,574 edges and 1,087,716 faces, and 55,498 vertices, 167,106 edges and 111,404 faces). We generate 100 random pairs of closed paths on each surface. The first path is generated randomly, with a lower bound for its length given by a random number between 100 and 10,000: passed the lower bound, the path is randomly extended until it returns to its origin vertex. The second path is generated from a sequence of elementary deformations of the first path, so that the two paths are homotopic. The number of elementary deformations is a random number between 100 and 10,000.
The computation time of the 200 are_freely_homotopic() tests are given in Figure 83.4, according to the number of darts of the two input paths.
are_freely_homotopic() tests, according to the number of darts or the input paths. The third benchmark is similar to the previous one, except that we use a genus 5 surface with one vertex, 10 edges and 1 face. The length of the paths to be tested is varied between 10 and 30,000,000. The computation time of the are_freely_homotopic() tests are given in Figure 83.5. The free homotopy test takes 17 seconds for paths of length 10,000,000.
are_freely_homotopic() tests, according to the number of darts or the input paths for random paths on the canonical reduced surface of genus five. Note that the linear time implementation of the algorithms is remarkably observed in each benchmark.
The algorithm to find shortest non-contractible cycle through a vertex in [1] can be summarized as follows.
A mesh \(\cal{M}\) consists of components such as vertices, edges, and faces, thus it can be seen as a graph \(G\) embedded in the surface \(\Sigma\).
Let \(T\) be a spanning tree of \(G\). Let \(C^*\) be the subgraph of the dual graph \(G^*\) of \(G\) with the same vertex set as \(G^*\) and the edge set be \(E(G^*)\backslash E(T)^*\). Repeatedly remove from \(C^*\) the edges with an incident vertex of degree one, the remaining set of edges is denoted as \(E_{nc}(T)^*\). It has been proven that for any edge \(ab\in E_{nc}(T)\), the concatenation of the path from a vertex \(v\) following \(T\) to \(a\), the edge \(ab\), and the path from \(b\) following \(T\) back to \(v\) is a closed path (denoted as \(\tau(T, v, ab)\)) and is non-contractible. Furthermore, if \(T\) is a BFS tree (or Dijkstra tree in the weighted case) rooted at \(v\), the shortest cycle found among \(\tau(T, v, e)\) for any \(e\in E_{nc}(T)\) is the shortest non-contractible cycle through \(v\).
Although it is said in Specifying the Input Surface and Curves that the given mesh should be closed (no dart is 2-free), the algorithm to find shortest non-contractible cycles also works if the surface contains boundaries.
The reader is recommended to read the section Compute Shortest Non-Contractible Cycle before reading this section.
The face width is the minimum number of intersection points between \(G\) and any non-contractible cycle of \(\Sigma\) (these cycles do not necessarily follow the edges of \(G\)). As a result, the face width of \(G\) is half the edge width of \(R(G)\), where \(R(G)\) denotes the radial graph of \(G\). The radial graph of \(G\) is a bipartite graph, constructed as follows. Start with the original graph. For every face of the original graph, add a vertex in the interior of the face and connect this vertex to all the vertices of the face. After doing this for all faces, remove all edges from the original graph.
A quadrangulation is a combinatorial map whose faces are quadrilaterals, i.e. have four sides. For efficiency of the homotopy test, the input combinatorial surface \(\cal{M}\) is first turned into a quadrangulation with only two vertices. The transformation is performed as follows.
This quadrangulation \(\cal{Q}\) is stored in a Surface_mesh_topology::Curves_on_surface_topology. In order to perform a homotopy test, each input curve \(C\) is transformed into a (closed) path in \(\cal{Q}\) as follows. If an edge of \(C\) is part of the contracted spanning tree, we simply ignore that edge. Otherwise the edge can be replaced by two consecutive edges of \(\cal{Q}\) to obtain a new path \(C'\) in the vertex-edge graph of \(\cal{Q}\) so that \(C'\) is a continuous deformation of \(C\). Hence, deciding if \(C\) is contractible in \(\cal{M}\) is equivalent to test if \(C'\) is contractible in \(\cal{Q}\).
In order to test if two input curves \(C\) and \(D\) in \(\cal{M}\) are homotopic they are first replaced by curves \(C'\) and \(D'\) in \(\cal{Q}\) as above. Those curves are further transformed into canonical forms that only depend on their homotopy classes. The transformation to canonical form relies on three basic operations that we now describe.
The canonical form of a curve is obtained by flattening its brackets, removing its spurs and pushing its right L-shapes to the right until the curve has no more brackets, spurs or L-shapes. This can be done in time proportional to the number of edges of the curve. Note that the above three operations preserve the homotopy class of the curve.
It can be proven that the canonical form is uniquely defined and only depends on the homotopy class of the curve. Hence, the curves \(C'\) and \(D'\) in \(\cal{Q}\) are homotopic if and only if their canonical forms are equal. Since each curve is defined as a sequence of (oriented) edges up to a cyclic permutation, we resort to the Knuth-Morris-Pratt algorithm to decide in linear time if the canonical forms are the same up to a cyclic permutation.
The simplicity test relies on the fact that a closed curve is homotopic to a simple one if and only if its canonical form can be made intersection free via infinitesimal perturbations together with some homotopy preserving operations. One can imagine each edge in the quadrangulation has a width and each vertex in the quadrangulation has an area so that paths visiting the same vertex/edge multiple times can avoid intersection after perturbation within a vertex or an edge. See Figure 83.14 for an example.
Such a perturbation can be encoded as a transverse ordering of the edges from the canonical form which traverse the same edge in the quadrangulation. The idea of the algorithm is to traverse the canonical form, one edge at a time, to inductively build such orderings and try to avoid intersection as best as we can.
There is an easy case where we know for sure that a closed curve cannot be deformed to simple one: If the canonical form can be expressed as concatenation of two or more copies of the same path. So the first step of the algorithm is to detect repetition.
Let \(P\) be a path and let \(+\) be the operator of concatenation. It can be shown that \(P\) contains no repetition if and only if there are only 2 matchings of \(P\) in \(P+P\) (matching the first and the second copy). The algorithm resorts to the Knuth-Morris-Pratt algorithm to decide in linear time if the canonical form contains repetitions.
Apart from applying perturbation, the algorithm also tries to avoid crossings using a homotopy-preserving operation called as a switch. A switch is triggered whenever an intersection could be avoided by turning a left L-shaped subpath into a right L-shaped subpath. See Figure 83.15 for an example.
As the algorithm inductively builds orderings, it has to determine a relative order between the edge being processed and the edges that have been ordered. There are three cases to consider.
The transverse orderings are stored in red-black trees, one for each edge of the quadrangulation. So each insertion or search takes \(O(\log{l})\) time, where \(l\) is the length of the closed curve.
After computing a tentative ordering within the edges of the path, we have to verify that such an ordering could result in an intersection free arrangement. Since there is no intersection within an edge, we only need to verify this for each vertex in the quadrangulation. Each vertex is naturally associated with a circular ordering of the incident path edges by concatenating clockwise the orderings computed for every incident edge in the quadrangulation. We consider the two consecutive edges composing a turn (one going in the vertex, one going out of the vertex) at the vertex being verified as a pair. The ordering at the vertex is intersection free if and only if there is no two pairs crossing each other according to the clockwise ordering around the vertex. In other words, for any two pairs \((a, a')\) and \((b, b')\), none of the subsequences \(a, b, a', b'\) or \(a, b', a', b\) should appear in the clockwise ordering. This is very similar to verifying balanced parentheses in a string. We traverse clockwise at each vertex and use a stack-based algorithm to verify in linear time that the ordering produces a cycle without self-intersection.
The code was developed in 2018 by Guillaume Damiand and Francis Lazarus. Felix Castillon contributed to the extension of the homotopy test to the case of surfaces with boundaries. Thien Hoang added methods to compute shortest non-contractible cycles, edge width and face width as part of the program Google Summer of Code 2019. Shuhao Tan added methods to test simplicity of a closed curve as part of the program Google Summer of Code 2020.
1.8.13