CGAL 5.3 - Surface Mesh Topology
User Manual

Author
Guillaume Damiand and Francis Lazarus

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.

Introduction

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.

Shortest Non-Contractible Curve

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:

  • Given a surface mesh \(\cal{M}\) and a dart handle dh, compute a shortest non-contractible combinatorial curve passing through the source vertex of dh,
  • Given a surface mesh \(\cal{M}\), compute a shortest non-contractible combinatorial curve without the previous vertex requirement. When all the edges have the same unit length, the length of a shortest non-contractible curve is known as the edge width of the surface,
  • Given a surface mesh \(\cal{M}\), compute a shortest non-contractible topological curve. It can be assumed that this curve does not cross the edges of \(\cal{M}\) and only passes through the vertices. It follows that the curve can be described by a circular sequence of traversed faces alternating with the vertices it passes through. The length of this curve (i.e., the number of vertices it passes through) is known as the face width of the surface.

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.

Homotopy Test

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:

  • Given a combinatorial surface \(\cal{M}\) and a closed combinatorial curve specified as a sequence of edges of \(\cal{M}\), decide if the curve is contractible on \(\cal{M}\),
  • Given a combinatorial surface \(\cal{M}\) and two closed combinatorial curves on \(\cal{M}\), decide if the two curves are related by a continuous transformation,
  • Given a combinatorial surface \(\cal{M}\) and two, non-necessarily closed, combinatorial curves on \(\cal{M}\), decide if the two curves are related by a continuous transformation that fixes the curve extremities. The curves should have common endpoints, otherwise the answer to the query is trivially negative.

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.

free-vs-fixed-endpoints.svg
Figure 83.1 On the upper left surface the green curve is contractible. The red and blue curves share the same (green) endpoint. (Being closed, their two endpoints coincide.) Although these last two curves are not homotopic with fixed endpoints they are freely homotopic as shown by the suggested continuous transformation of the blue curve.

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].

Simplicity Test

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:

  • Given a combinatorial surface \(\cal{M}\) and a closed combinatorial curve specified as a sequence of edges of \(\cal{M}\), decide if the curve is homotopic to a simple one on \(\cal{M}\).

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.

API Description

Specifying the Input Surface and Curves

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.

  • Simply push the next halfedge using the push_back() member function. One can also specify if this halfedge should have its direction flipped so as to satisfy the condition of a Path_on_surface (see the description of `can_be_pushed()` below) This can be done even when the path is empty,
  • The user may push the index of the next halfedge instead of the halfedge itself with the member function push_back_by_index(). This may however be at the cost of an overhead computation mapping the index to the actual dart,
  • The path may be extended with the member function extend_positive_turn() by specifying the next halfedge thanks to a number of positive turns with respect to the previous dart/halfedge in the path. Calling this previous halfedge 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,
  • Finally, when the input surface is a model of 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.

Polygonal Schema Helper

Polygonal Schema

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.

incremental-builder.svg
Figure 83.2 Left, a surface described by a single facet with eight edges pairwise identified. The resulting (topological) surface is shown in Figure 83.9. Middle, a surface described by three labeled quadrilaterals. Right, a single labeled facet. The corresponding surface is topologically equivalent to the middle example.

The code for creating the above left and middle examples appear in the polygonal schema examples below. The class provides the following functionalities.
  • add_facet(s) adds a polygon to the current collection of polygons. If the polygon has "n" sides, "s" is a sequence of "n" edge labels possibly preceded by a minus and separated by blanks.
  • alternatively, the user can add a facet by adding edge labels one at a time using the member functions init_facet(), add_edges_to_facet() and finish_facet()

Polygonal Schema with Boundary

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.

  • perforate_facet (h) perforates the facet identified by the halfedge h. If s is the label of the oriented edge corresponding to h, one may equivalently call perforate_facet (s).
  • Similarly, fill_facet (h or s) turns a facet into a plain one.
  • The member functions get_dart_labeled(s) and get_label(h) allow to easily pass from a halfedge to its label and vice versa.

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}\).

Curves on Polygonal Schema

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.

Compute Shortest Non-contractible Cycle

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:

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.

Compute Face Width

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.

Testing Homotopy

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:

  • is_contractible( \(p_1\)) returns true if the closed curve \(p_1\) is contractible,
  • are_freely_homotopic( \(p_1\), \(p_2\)) returns true if the closed curves \(p_1\) and \(p_2\) are freely homotopic,
  • are_homotopic_with_fixed_endpoints( \(p_1\), \(p_2\)) returns 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,

Note
The user must not modify the input surface as long as homotopy tests are performed with this 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.

Testing Simplicity

Given a Surface_mesh_topology::Path_on_surface \(p\), the class Surface_mesh_topology::Curves_on_surface_topology provides the following function:

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.

Note
The user must not modify the input surface as long as simplicity tests are performed with this 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.

Examples

Compute Shortest Non-contractible Cycle

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

#include <CGAL/Linear_cell_complex_for_combinatorial_map.h>
#include <CGAL/Linear_cell_complex_constructors.h>
#include <CGAL/Curves_on_surface_topology.h>
#include <CGAL/Path_on_surface.h>
#include <CGAL/draw_face_graph_with_paths.h>
#include <CGAL/Random.h>
#include <iostream>
#include <cstdlib>
double cycle_length(const LCC_3& lcc, const Path_on_surface& cycle)
{ // Compute the length of the given cycle.
double res=0;
for (std::size_t i=0; i<cycle.length(); ++i)
{ res+=std::sqrt
(CGAL::squared_distance(lcc.point(cycle[i]),
lcc.point(lcc.other_extremity(cycle[i])))); }
return res;
}
void display_cycle_info(const LCC_3& lcc, const Path_on_surface& cycle)
{ // Display information about the given cycle.
if (cycle.is_empty()) { std::cout<<"Empty."<<std::endl; return; }
std::cout<<"Root: "<<lcc.point(cycle[0])<<"; "
<<"Number of edges: "<<cycle.length()<<"; "
<<"Length: "<<cycle_length(lcc, cycle)<<std::endl;
}
int main(int argc, char* argv[])
{
std::string filename(argc==1?"data/3torus.off":argv[1]);
bool draw=(argc<3?false:std::string(argv[2])=="-draw");
LCC_3 lcc;
if (!CGAL::load_off(lcc, filename.c_str())) // Load the off file.
{
std::cout<<"Cannot read file '"<<filename<<"'. Exiting program"<<std::endl;
return EXIT_FAILURE;
}
std::cout<<"File '"<<filename<<"' loaded. Finding shortest non contractible cycle..."<<std::endl;
LCC_3::Dart_const_handle root=lcc.dart_handle
(CGAL::get_default_random().get_int(0, static_cast<int>(lcc.number_of_darts()))); // One dart of the mesh
Path_on_surface cycle1=
cst.compute_shortest_non_contractible_cycle_with_base_point(root);
Path_on_surface cycle2=
cst.compute_shortest_non_contractible_cycle_with_base_point(root, wf);
std::cout<<"Cycle 1 (pink): "; display_cycle_info(lcc, cycle1);
std::cout<<"Cycle 2 (green): "; display_cycle_info(lcc, cycle2);
if (draw) { CGAL::draw(lcc, {cycle1, cycle2}); }
return EXIT_SUCCESS;
}

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

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Curves_on_surface_topology.h>
#include <CGAL/Path_on_surface.h>
#include <CGAL/squared_distance_3.h>
#include <CGAL/draw_face_graph_with_paths.h>
#include <fstream>
double cycle_length(const Mesh& mesh, const Path_on_surface& cycle)
{ // Compute the length of the given cycle.
double res=0;
for (std::size_t i=0; i<cycle.length(); ++i)
{ res+=std::sqrt
(CGAL::squared_distance(mesh.point(mesh.vertex(mesh.edge(cycle[i]), 0)),
mesh.point(mesh.vertex(mesh.edge(cycle[i]), 1)))); }
return res;
}
void display_cycle_info(const Mesh& mesh, const Path_on_surface& cycle)
{ // Display information about the given cycle.
if (cycle.is_empty()) { std::cout<<"Empty."<<std::endl; return; }
std::cout<<"Root: "<<mesh.point(mesh.vertex(mesh.edge(cycle[0]), 0))<<"; "
<<"Number of edges: "<<cycle.length()<<"; "
<<"Length: "<<cycle_length(mesh, cycle)<<std::endl;
}
int main(int argc, char* argv[])
{
std::string filename(argc==1?"data/3torus.off":argv[1]);
bool draw=(argc<3?false:(std::string(argv[2])=="-draw"));
Mesh sm;
if(!CGAL::IO::read_polygon_mesh(filename, sm))
{
std::cout<<"Cannot read file '"<<filename<<"'. Exiting program"<<std::endl;
return EXIT_FAILURE;
}
std::cout<<"File '"<<filename<<"' loaded. Finding edge-width of the mesh..."<<std::endl;
Path_on_surface cycle1=cst.compute_edge_width(true);
Path_on_surface cycle2=cst.compute_shortest_non_contractible_cycle(wf, true);
std::cout<<"Cycle 1 (pink): "; display_cycle_info(sm, cycle1);
std::cout<<"Cycle 2 (green): "; display_cycle_info(sm, cycle2);
if (draw) { CGAL::draw(sm, {cycle1, cycle2}); }
return EXIT_SUCCESS;
}

In these two examples, the mesh and the cycles can be visualized if CGAL_Qt5 is enabled.

Compute Face Width

The following example computes the face width, and visualizes it if CGAL_Qt5 is enabled.


File Surface_mesh_topology/facewidth.cpp

#include <CGAL/Linear_cell_complex_for_combinatorial_map.h>
#include <CGAL/Linear_cell_complex_constructors.h>
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <CGAL/Curves_on_surface_topology.h>
#include "draw_facewidth.h"
using Dart_const_handle=LCC_3::Dart_const_handle;
int main(int argc, char* argv[])
{
std::cout<<"Program facewidth_on_unweighted_map started."<<std::endl;
std::string filename(argc==1?"data/double-torus.off":argv[1]);
bool draw=(argc<3?false:std::string(argv[2])=="-draw");
std::ifstream inp(filename);
if (inp.fail())
{
std::cout<<"Cannot read file '"<<filename<<"'. Exiting program"<<std::endl;
return EXIT_FAILURE;
}
LCC_3 lcc;
CGAL::load_off(lcc, inp);
std::cout<<"File '"<<filename<<"' loaded. Finding the facewidth..."<<std::endl;
CST cst(lcc, true);
std::vector<Dart_const_handle> cycle=cst.compute_face_width(true);
if (cycle.size()==0)
{ std::cout<<" Cannot find such cycle."<<std::endl; }
else
{
std::cout<<" Number of faces: "<<cycle.size()<<std::endl;
if (draw) { draw_facewidth(lcc, cycle); }
}
return EXIT_SUCCESS;
}

Basic Homotopy Test

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

#include <CGAL/Linear_cell_complex_for_combinatorial_map.h>
#include <CGAL/Linear_cell_complex_constructors.h>
#include <CGAL/Curves_on_surface_topology.h>
#include <CGAL/Path_on_surface.h>
#include <CGAL/draw_face_graph_with_paths.h>
void create_path_1(Path_on_surface<LCC_3_cmap>& p)
{
p.push_back_by_index(14); // Its starting dart
for (int i=0; i<7; ++i)
{ p.extend_positive_turn(2); } // Extend the path
}
void create_path_2(Path_on_surface<LCC_3_cmap>& p)
{ p.push_back_by_index({202, 206, 335, 317, 322, 69, 62, 414}); }
void create_path_3(Path_on_surface<LCC_3_cmap>& p)
{
p.push_back_by_index(470); // Its starting dart
for (int i=0; i<13; ++i)
{ p.extend_positive_turn(2); } // Extend the path
}
int main(int argc, char** argv)
{
bool draw=(argc>1?std::string(argv[1])=="-draw":false);
LCC_3_cmap lcc;
if (!CGAL::load_off(lcc, "data/double-torus.off"))
{
std::cout<<"ERROR reading file data/double-torus.off"<<std::endl;
exit(EXIT_FAILURE);
}
Path_on_surface<LCC_3_cmap> p1(lcc), p2(lcc), p3(lcc);
create_path_1(p1);
create_path_2(p2);
create_path_3(p3);
bool res1=cst.is_contractible(p1);
std::cout<<"Path p1 (pink) "<<(res1?"IS":"IS NOT")
<<" contractible."<<std::endl;
bool res2=cst.are_freely_homotopic(p1, p2);
std::cout<<"Path p1 (pink) "<<(res2?"IS":"IS NOT")
<<" homotopic with path p2 (green)."<<std::endl;
bool res3=cst.are_freely_homotopic(p1, p3);
std::cout<<"Path p1 (pink) "<<(res3?"IS":"IS NOT")
<<" homotopic with path p3 (orange)."<<std::endl;
if (draw)
{ CGAL::draw(lcc, {p1, p2, p3}); }
return EXIT_SUCCESS;
}

Basic Simplicity Test

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

#include <CGAL/Linear_cell_complex_for_combinatorial_map.h>
#include <CGAL/Linear_cell_complex_constructors.h>
#include <CGAL/Curves_on_surface_topology.h>
#include <CGAL/Path_on_surface.h>
#include <CGAL/draw_face_graph_with_paths.h>
void create_path(Path_on_surface<LCC_3_cmap>& p)
{
p.push_back_by_index(682); // Its starting dart
for (int i=0; i<11; ++i)
{ p.extend_positive_turn(2); } // Extend the path
for (int i=0; i<5; ++i)
for (int i=0; i<2; ++i)
for (int i=0; i<5; ++i)
for (int i=0; i<2; ++i)
for (int i=0; i<2; ++i)
for (int i=0; i<8; ++i)
for (int i=0; i<4; ++i)
for (int i=0; i<5; ++i)
for (int i=0; i<5; ++i)
for (int i=0; i<3; ++i)
for (int i=0; i<11; ++i)
}
int main(int argc, char** argv)
{
bool draw=(argc>1?std::string(argv[1])=="-draw":false);
LCC_3_cmap lcc;
if (!CGAL::load_off(lcc, "data/double-torus.off"))
{
std::cout<<"ERROR reading file data/double-torus.off"<<std::endl;
exit(EXIT_FAILURE);
}
create_path(p);
bool res=cst.is_homotopic_to_simple_cycle(p);
std::cout<<"Path p (pink) "<<(res?"IS":"IS NOT")
<<" simple."<<std::endl;
if (draw)
{ CGAL::draw(lcc, {p}); }
return EXIT_SUCCESS;
}

Polygonal Schema

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

#include <CGAL/Polygonal_schema.h>
#include <CGAL/Path_on_surface.h>
#include <CGAL/Curves_on_surface_topology.h>
#include <iostream>
#include <cstdlib>
int main()
{
PS ps;
ps.add_facet("a b -a -b c d -c -d");
Path_on_surface<PS> p1(ps); p1.push_back_by_label("a");
Path_on_surface<PS> p2(ps); p2.push_back_by_label("b c a -c -b");
bool res1=cst.are_freely_homotopic(p1, p2);
std::cout<<"Paths p1 and p2 "<<(res1?"ARE":"ARE NOT")
<<" freely homotopic."<<std::endl;
bool res2=cst.are_base_point_homotopic(p1, p2);
std::cout<<"Paths p1 and p2 "<<(res2?"ARE":"ARE NOT")
<<" base point homotopic."<<std::endl;
return EXIT_SUCCESS;
}

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

#include <CGAL/Polygonal_schema.h>
#include <CGAL/Path_on_surface.h>
#include <CGAL/Curves_on_surface_topology.h>
#include <iostream>
#include <cstdlib>
int main()
{
PS ps;
ps.add_facet("a b -a c"); // First facet, giving directly its sequence of edges
ps.add_facet("d -c e -b"); // Second facet
ps.init_facet(); // Third facet
ps.add_edges_to_facet("f"); // Here, each edge is added one at a time
ps.add_edges_to_facet("-d");
ps.add_edges_to_facet("-f");
ps.add_edges_to_facet("-e");
ps.finish_facet();
ps.perforate_facet("f");
std::cout<<"Number of cells of the combinatorial maps: ";
ps.display_characteristics(std::cout)<<std::endl;
p.push_back_by_label("a b -a e -b d");
bool res=cst.is_contractible(p);
std::cout<<"Path "<<(res?"IS":"IS NOT")<<" contractible."<<std::endl;
return EXIT_SUCCESS;
}

Open Path

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

#include <CGAL/Surface_mesh.h>
#include <CGAL/Curves_on_surface_topology.h>
#include <CGAL/Path_on_surface.h>
#include <CGAL/draw_face_graph_with_paths.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point_3;
void create_path_1(Path_on_surface<SM>& p)
{
p.push_back_by_index(88); // Its starting dart
for (int i=0; i<3; ++i)
{ p.extend_positive_turn(2); } // Extend the path
}
void create_path_2(Path_on_surface<SM>& p)
{
p.push_back_by_index(300); // Its starting dart
for (int i=0; i<3; ++i)
{ p.extend_negative_turn(2); } // Extend the path
}
void create_path_3(Path_on_surface<SM>& p)
{
p.push_back_by_index(87); // Its starting dart
p.extend_positive_turn(1); // Extend the path
for (int i=0; i<3; ++i)
}
int main(int argc, char* argv[])
{
bool draw=(argc>1?std::string(argv[1])=="-draw":false);
std::ifstream in("data/double-torus.off");
if (!in.is_open())
{
std::cout<<"ERROR reading file data/double-torus.off"<<std::endl;
exit(EXIT_FAILURE);
}
SM sm;
in>>sm;
Path_on_surface<SM> p1(sm), p2(sm), p3(sm);
create_path_1(p1);
create_path_2(p2);
create_path_3(p3);
bool res1=cst.are_base_point_homotopic(p1, p2);
std::cout<<"Path p1 (pink) "<<(res1?"IS":"IS NOT")
<<" base point homotopic with path p2 (green)."<<std::endl;
bool res2=cst.are_base_point_homotopic(p2, p3);
std::cout<<"Path p2 (green) "<<(res2?"IS":"IS NOT")
<<" base point homotopic with path p3 (orange)."<<std::endl;
if (draw)
{ CGAL::draw(sm, {p1, p2, p3}); }
return EXIT_SUCCESS;
}

Benchmarks

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.

Combinatorial Surface Topology Computation Time

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.

computation-time-reduce-surface.svg
Figure 83.3 Computation time of Surface_mesh_topology::Curves_on_surface_topology constructions, according to the number of darts or the input surface.

Path Homotopy Tests

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.

computation-time-path-homotopy.svg
Figure 83.4 Computation time of 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.

computation-time-polygonal-schema.svg
Figure 83.5 Computation time of 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.

Implementation Details

Compute Shortest Non-Contractible Cycle

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.

Compute Face Width

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.

Reducing to a Quadrangulation

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.

  1. A spanning tree of the graph composed of the vertices and edges of \(\cal{M}\) is computed. See Figure 83.6.
    spanning_tree.svg
    Figure 83.6 Left, a combinatorial map with three faces (red, yellow, green). Right, a spanning tree of its graph.

  2. The edges of the spanning are contracted. The resulting surface has the same topology as \(\cal{M}\) and has a single vertex.
    contract_tree.svg
    Figure 83.7 The contraction of a spanning tree merges all the vertices into a single vertex.

  3. The faces are merged into a single face by iteratively erasing edges incident to distinct faces. Those edges corresponds to a spanning tree of the dual combinatorial map.
    merge_faces.svg
    Figure 83.8 The green, red and yellow faces are merged by removing edges. The resulting reduced surface has a single vertex and a single face.

    Cutting through the graph of the reduced surface, we obtain a face that can be flattened into the plane.
    cut-open.svg
    Figure 83.9 If \(\cal{M}\) is obtained by gluing \(g\) tori, i.e. \(\cal{M}\) has genus \(g\), the cut-open reduced surface has \(4g\) sides.

  4. A vertex is introduced in the middle of this unique face and joined by new edges to the corners of this face to form a triangulation with \(4g\) triangles. Gluing back along the (old) edges of the reduced surface and deleting them amounts to merge the triangles by pairs. We thus obtain a quadrangulated surface \(\cal{Q}\) with \(2g\) quadrilaterals, \(4g\) edges and 2 vertices.
    quad_mesh.svg
    Figure 83.10 Triangles of the same color are merged into quadrilaterals. All the light blue vertices correspond to a same vertex on the glued surface.

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}\).

Canonical Form

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.

  1. A bracket in a curve is a subsequence of edges along a row of quadrilaterals, surrounded by two edges along the end sides of the row. A bracket can be flattened by replacing the corresponding subpath with a shorter subpath going along the other long side of the row. See Figure 83.11.
    bracket.svg
    Figure 83.11 Left, a blue curve in a quadrangulation (for clarity the quadrangulation has more than two vertices). Middle, a bracket of the blue curve. Right, the bracket has been flattened.

  2. A spur in a curve is a subsequence of two opposite edges. A spur can be deleted to shorten the curve. See Figure 83.12.
    spur.svg
    Figure 83.12 Removing a spur.

  3. A right L-shape in a curve is a subsequence of edges going along the exterior side of a sequence of quadrilaterals forming an L, with the interior of the L to its right. This notion takes into account the traversal direction of the curve. A right L-shape subpath can be pushed to the right by replacing it with the other side of the L-shaped sequence of quadrilaterals. See Figure 83.13
    push_right.svg
    Figure 83.13 Pushing an L-shaped subpath to its right.

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.

Homotopy Test

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.

Simplicity Test

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.

perturbation_sample.svg
Figure 83.14 Applying a perturbation to remove 2 intersections between the red and the blue subpaths.

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.

Detect Repetition

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.

Avoid Crossing by Switching

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.

switch_sample.svg
Figure 83.15 Top-left, an intersection between the red subpath and the green subpath at the start of the left L-shaped red subpath. Top-right, switch the left L-shape to right L-shape to avoid the intersection. Bottom-left or right, no switch is performed because no intersection can be avoided.

Decide Relative Order

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.

  1. The current edge and its predecessor, say \(e\), form a subpath of length two so that a parallel copy of this subpath has already been processed, and \(e\) is adjacent to its parallel copy in the previously constructed ordering. In this case, the current edge must be adjacent to its copy in the transverse ordering. See Figure 83.16.
    relative_order_corner.svg
    Figure 83.16 The red edge is being processed. The blue edge is adjacent to the green edge (the predecessor \(e\) of the red edge) in the previously constructed ordering and forms the same turn (blue-pink turn) as the green-red turn. The red edge should be right next to the pink edge in the ordering so as to avoid crossings.

  2. When the previous situation does not occur, the current edge has to be compared against every parallel copy already processed. In this case, the predecessors of the copy and of the current edge form a Y shape with the current edge. The circular order of the edges in this Y can be used to determine the relative order between the current edge and its copy. See Figure 83.17.
    relative_order_normal.svg
    Figure 83.17 The red edge is being processed and is compared against the pink edge. Since the green edge (the predecessor of the red edge) is to the right of the blue edge around the vertex, the red edge must be to the right of the pink edge in the transverse ordering.

  3. There is one special case where the comparison of the current edge with a parallel copy cannot be deduced form previous computations: When this copy happens to be the very first edge processed in the traversal. Indeed, the predecessor of the first edge (aka the last edge of the path) has not been processed yet. If the last edge runs parallel to the predecessor of the current edge, we cannot determine their relative order. So the idea is to keep following the predecessors of the current edge and of the first edge until they diverge, at which point the relative order around the vertex can be used to determine the relative order. See Figure 83.18. This can be precomputed by finding all the longest common suffixes of the path against its circular shifts. A modified Knuth-Morris-Pratt algorithm is applied to preprocess the path in linear time.
    relative_order_first.svg
    Figure 83.18 The red edge is being processed and is compared against the pink edge which is the first edge of the path. The blue and green edges are the first diverging pair when tracing backward. The dashed line means that edges have not been processed yet. Since the green edge lies to the right of the blue edge around the vertex, the red edge must be to the right of the pink edge in the ordering.

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.

Verify Ordering

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.

History

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.