CGAL 5.1 - 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 82.1 below illustrates the three types of queries.

free-vs-fixed-endpoints.svg
Figure 82.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 [2], providing a linear time algorithm for the above homotopy tests. This is a simplified version of the linear time algorithm by Lazarus and Rivaud [3].

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 82.2 Left, a surface described by a single facet with eight edges pairwise identified. The resulting (topological) surface is shown in Figure 82.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.

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 <fstream>
#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>
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"));
std::ifstream inp(filename);
if (inp.fail())
{
std::cout<<"Cannot read file '"<<filename<<"'. Exiting program"<<std::endl;
return EXIT_FAILURE;
}
Mesh sm;
inp>>sm;
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;
}

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 82.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 82.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 82.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 82.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 82.4, according to the number of darts of the two input paths.

computation-time-path-homotopy.svg
Figure 82.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 82.5. The free homotopy test takes 17 seconds for paths of length 10,000,000.

computation-time-polygonal-schema.svg
Figure 82.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 82.6.
    spanning_tree.svg
    Figure 82.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 82.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 82.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 82.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 82.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 82.11.
    bracket.svg
    Figure 82.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 82.12.
    spur.svg
    Figure 82.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 82.13
    push_right.svg
    Figure 82.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.

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.