CGAL 5.6 - Triangulated Surface Mesh Simplification
User Manual

Authors
Fernando Cacciola, Mael Rouxel-Labbé, and Baskın Şenbaşlar
Illustration-Simplification-ALL.jpg

Introduction

Surface mesh simplification is the process of reducing the number of faces used in a surface mesh while keeping the overall shape, volume and boundaries preserved as much as possible. It is the opposite of subdivision.

The algorithm presented here can simplify any oriented 2-manifold surface, with any number of connected components, with or without boundaries (border or holes) and handles (arbitrary genus), using a method known as edge collapse. Roughly speaking, the method consists of iteratively replacing an edge with a single vertex, removing 2 triangles per collapse.

Edges are collapsed according to a priority given by a user-supplied cost function, and the coordinates of the replacing vertex are determined by another user-supplied placement function. The algorithm terminates when a user-supplied stop predicate is met, such as reaching the desired number of edges.

The algorithm implemented here is generic in the sense that it does not require the surface mesh to be of a particular type but to be a model of the MutableFaceGraph and HalfedgeListGraph concepts. We give examples for Surface_mesh, Polyhedron_3, and OpenMesh.

The design is policy-based (https://en.wikipedia.org/wiki/Policy-based_design), meaning that you can customize some aspects of the process by passing a set of policy objects. Each policy object specifies a particular aspect of the algorithm, such as how edges are selected and where the replacement vertex is placed. All policies have a sensible default. Furthermore, the API uses the so-called named-parameters technique which allows you to pass only the relevant parameters, in any order, omitting those parameters whose default is appropriate.

Overview of the Simplification Process

The free function that implements the simplification algorithm takes not only the surface mesh and the desired stop predicate but a number of additional parameters which control and monitor the simplification process. This section briefly describes the process in order to set the background for the discussion of the parameters to the algorithm.

There are two slightly different "edge" collapse operations. One is known as edge-collapse while the other is known as halfedge-collapse. Given an edge e joining vertices w and v, the edge-collapse operation replaces e,w and v for a new vertex r, while the halfedge-collapse operation pulls v into w, eliminating e and leaving w in place. In both cases the operation removes the edge e along with the 2 triangles adjacent to it.

This package uses the halfedge-collapse operation, which is implemented by removing, additionally, 1 vertex (v) and 2 edges, one per adjacent triangle. It optionally moves the remaining vertex (w) into a new position, called placement, in which case the net effect is the same as in the edge-collapse operation.

Naturally, the surface mesh that results from an edge collapse deviates from the initial surface mesh by some amount, and since the goal of simplification is to reduce the number of triangles while retaining the overall look of the surface mesh as much as possible, it is necessary to measure such a deviation. Some methods attempt to measure the total deviation from the initial surface mesh to the completely simplified surface mesh, for example, by tracking an accumulated error while keeping a history of the simplification changes. Other methods, like the ones implemented in this package, attempt to measure only the cost of each individual edge collapse (the local deviation introduced by a single simplification step) and plan the entire process as a sequence of steps of increasing cost.

Global error tracking methods produce highly accurate simplifications but take up a lot of additional space. Cost-driven methods, like the ones in this package, produce slightly less accurate simplifications but take up much less additional space, even none in some cases.

This package provides two cost-driven methods. The first cost-driven method implemented in this package, namely Lindstrom&Turk, is mainly based on [4], [5], with contributions from [3], [2] and [1]. The second cost-driven method implemented in this package, namely Garland&Heckbert, is mainly based on [2], with enhancements from Trettner and Kobbelt [7].

The algorithm proceeds in two stages. In the first stage, called collection stage, an initial collapse cost is assigned to each and every edge in the surface mesh. Then in the second stage, called collapsing stage, edges are processed in order of increasing cost. Some processed edges are collapsed while some are just discarded. Collapsed edges are replaced by a vertex and the collapse cost of all the edges now incident on the replacement vertex is recalculated, affecting the order of the remaining unprocessed edges.

Not all edges selected for processing are collapsed. A processed edge can be discarded right away, without being collapsed, if it does not satisfy certain topological and geometric conditions.

The algorithm presented in [2] contracts (collapses) arbitrary vertex pairs and not only edges by considering certain vertex pairs as forming a pseudo-edge and proceeding to collapse both edges and pseudo-edges in the same way as in [4], [5]. However, contracting an arbitrary vertex-pair may result in a non-manifold surface mesh, but the current state of this package can only deal with manifold surface meshes, thus, it can only collapse edges. That is, this package cannot be used as a framework for vertex contraction. Therefore, our implementation of [2] only collapses edges.

Cost Strategy

The specific way in which the collapse cost and vertex placement is calculated is called the cost strategy. The user can choose different strategies in the form of policies and related parameters, passed to the algorithm.

The current version of the package provides a set of policies implementing three strategies: the Lindstrom-Turk strategy, which is the default, the Garland-Heckbert family of strategies, and a strategy consisting of an edge-length cost with an optional midpoint placement (much faster but less accurate).

Lindstrom-Turk Cost and Placement Strategy

The main characteristic of the strategy presented in [4], [5] is that the simplified surface mesh is not compared at each step with the original surface mesh (or the surface mesh at a previous step) so there is no need to keep extra information, such as the original surface mesh or a history of the local changes. Hence the name memoryless simplification.

At each step, all remaining edges are potential candidates for collapsing and the one with the lowest cost is selected. The cost of collapsing an edge is given by the position chosen for the vertex that replaces it.

The replacement vertex position is computed as the solution to a system of 3 linearly-independent linear equality constraints. Each constraint is obtained by minimizing a quadratic objective function subject to the previously computed constraints. There are several possible candidate constraints and each is considered in order of importance. A candidate constraint might be incompatible with the previously accepted constraints, in which case it is rejected and the next constraint is considered. Once 3 constraints have been accepted, the system is solved for the vertex position. The first constraints considered preserves the shape of the surface mesh boundaries (in case the edge profile has boundary edges). The next constraints preserve the total volume of the surface mesh. The next constraints, if needed, optimize the local changes in volume and boundary shape. Lastly, if a constraint is still needed (because the ones previously computed were incompatible), a third (and last) constraint is added to favor equilateral triangles over elongated triangles.

The cost is then a weighted sum of the shape, volume and boundary optimization terms, where the user specifies the unit weighting unit factor for each term.

The local changes are computed independently for each edge using only the triangles currently adjacent to it at the time when the edge is about to be collapsed, that is, after all previous collapses. Thus, the transitive path of minimal local changes yields at the end a global change reasonably close to the absolute minimum.

Garland-Heckbert Cost and Placement Strategy

As in the case of the Lindstrom-Turk strategy, the Garland-Heckbert strategy introduced in [2] does not compare the resulting mesh with the original mesh and does not depend on an history of local changes. Instead, it encodes approximate distance to the original mesh by using quadric matrices that it assigns to each vertex.

In its classic version, a quadric matrix \( Q \) is assigned to each vertex \( v \) and encodes the total squared distance that any point \( p \) has to \( v \)'s neighboring faces, which is given by the matrix product \( p'Qp \). At each step, the edge which minimizes the cost of collapsing is selected for the collapse operation. The cost of collapsing an edge is calculated by minimizing the error function \( p'Qp \) where \( Q \) is the combined quadric matrices of the edge's endpoints and \( p \) is the point that minimizes the cost (i.e., decision variables). The point \( p \) that minimizes the objective function is picked as the new placement point. Since the error function is quadratic, finding its minimum can be done by simply calculating the gradient and equating it to zero. If this fails due to a singularity, the optimal point and the cost are found on the edge. After placing the new vertex, a new quadric matrix is assigned to it by simply summing the quadric matrices of the two extremities of the edge that has been collapsed. Additional pseudo-faces that are perpendicular to neighboring faces are added to quadric matrices of border vertices in order to preserve sharp borders of the mesh as much as possible.

An extension of the Garland-Heckbert was proposed by Trettner and Kobbelt [7], who proposed the concept of probabilistic quadrics. In this new approach, the energy minimization is no longer performed using distances to input planes or polygons as in the classic version; instead, input geometry is made uncertain by introducing Gaussian noise in the input (vertex positions and face normals). This variance naturally deteriorates the tightness of the result, but on the other hand it enables creating more uniform triangulations and the approach is more tolerant to noise, while still maintaining feature sensitivity.

Figure 71.1 Sappho's Head model (leftmost, 34882 vertices). From left to right, simplified output (1745 vertices) and symmetric Hausdorff distance for the four Garland-Heckbert variations: plane (0.217912), probabilistic plane (0.256801), triangle (0.268872), and probabilistic triangle (0.490846).

Cost Strategy Policies

The cost strategy used by the algorithm is selected by means of three policies: GetPlacement, GetCost, and Filter.

The GetPlacement policy is called to compute the new position for the remaining vertex after the halfedge-collapse. It returns an optional value, which can be absent if the edge should not be collapsed.

The GetCost policy is called to compute the cost of collapsing an edge. This policy uses the placement to compute the cost (which is an error measure) and determines the ordering of the edges.

The algorithm maintains an internal data structure (a mutable priority queue) which allows each edge to be processed in increasing cost order. Such a data structure requires some per-edge additional information, such as the edge's cost. If the record of per-edge additional information occupies N bytes of storage, simplifying a surface mesh of 1 million edges (a normal size) requires 1 million times N bytes of additional storage. Thus, to minimize the amount of additional memory required to simplify a surface mesh only the cost is attached to each edge and nothing else.

But this is a trade-off: the cost of a collapse is a function of the placement (the new position chosen for the remaining vertex) so before GetCost is called for each and every edge, GetPlacement must also be called to obtain the placement parameter to the cost function. But that placement, which is a 3D point, is not attached to each and every edge since that would easily triple the additional storage requirement. On the one hand, this dramatically saves on memory but on the other hand is a processing waste because when an edge is effectively collapsed, GetPlacement must be called again to know were to move the remaining vertex.

Earlier prototypes showed that attaching the placement to the edge, thus avoiding one redundant call to the placement function after the edge collapsed, has little impact on the total running time. This is because the cost of an each edge is not just computed once but changes several times during the process so the placement function must be called several times just as well. Caching the placement can only avoid the very last call, when the edge is collapsed, but not all the previous calls which are needed because the placement (and cost) changes.

Finally, we explain the PlacementFilter policy. While the cost is a scalar that is used in the priority queue, there may be additional criteria coming in to decide if the edge collapse shall be performed or not. While such a criterion could be easily integrated into the cost function, that is setting the cost to infinity in order not to be considered a candidate for a collapse, we test a criterion only for an edge when it is the next edge to collapse. This makes the mesh simplification faster in case the computation of the criterion is expensive, for example when we check if the simplified mesh is in a tolerance envelope of the input mesh.

API

API Overview

Since the algorithm is free from robustness issues there is no need for exact predicates nor constructions and Simple_cartesian<double> can be used safely. In the current version, 3.3, the LindstromTurk policies are not implemented for homogeneous coordinates, so a Cartesian kernel must be used.

The simplification algorithm is implemented as the free template function Surface_mesh_simplification::edge_collapse(). The function has two mandatory and several optional parameters.

Mandatory Parameters

There are two main parameters to the algorithm: the surface mesh to be simplified (in-place) and the stop predicate.

The surface mesh to simplify must be a model of the MutableFaceGraph and HalfedgeListGraph concepts.

The stop predicate is called after each edge is selected for processing, before it is classified as collapsible or not (thus before it is collapsed). If the stop predicate returns true the algorithm terminates.

Optional Named Parameters

The notion of named parameters was also introduced in the BGL. You can read about it in [6] or the following site: https://www.boost.org/libs/graph/doc/bgl_named_params.html. Named parameters allow the user to specify only those parameters which are really needed, by name, making the parameter ordering unimportant.

Say there is a function f() that takes 3 parameters called name, age and gender, and you have variables n,a and g to pass as parameters to that function. Without named parameters, you would call it like this: f(n,a,g), but with named parameters, you call it like this: f(name(n).age(a).gender(g)).

That is, you give each parameter a name by wrapping it into a function whose name matches that of the parameter. The entire list of named parameters is really a composition of function calls separated by a dot ( \( .\)). Thus, if the function takes a mix of mandatory and named parameters, you use a comma to separate the last non-named parameter from the first named parameters, like this:

f(non_named_par0, non_named_pa1, name(n).age(a).gender(g))

When you use named parameters, the ordering is irrelevant, so this: f(name(n).age(a).gender(g)) is equivalent to this: f(age(a).gender(g).name(n)), and you can just omit any named parameter that has a default value.

Sample Call

/*
surface_mesh : the surface_mesh to simplify
stop_predicate : policy indicating when the simplification must finish
vertex_index_map(vimap) : property-map giving each vertex a unique integer index
edge_index_map(eimap) : property-map giving each edge a unique integer index
edge_is_constrained_map(ebmap): property-map specifying whether an edge is a constrained edge or not
get_cost(cf) : function object computing the cost of a collapse
get_placement(pf) : function object computing the placement for the remaining vertex
filter(filter) : function object to reject a candidate chosen for the next edge collapse
visitor(vis) : function object tracking the simplification process
*/
int r = edge_collapse(surface_mesh, stop_predicate,
CGAL::parameters::vertex_index_map(vimap)
.edge_index_map(eimap)
.edge_is_border_map(ebmap)
.get_cost(cf)
.get_placement(pf)
.filter(filter)
.visitor(vis));

Examples

Example Using a Surface_mesh

The following example illustrates the simplification of a Surface_mesh. The unspecified cost strategy defaults to Lindstrom-Turk.


File Surface_mesh_simplification/edge_collapse_surface_mesh.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Surface_mesh_simplification/edge_collapse.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_count_ratio_stop_predicate.h>
#include <chrono>
#include <fstream>
#include <iostream>
typedef Kernel::Point_3 Point_3;
typedef CGAL::Surface_mesh<Point_3> Surface_mesh;
int main(int argc, char** argv)
{
Surface_mesh surface_mesh;
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/cube-meshed.off");
std::ifstream is(filename);
if(!is || !(is >> surface_mesh))
{
std::cerr << "Failed to read input mesh: " << filename << std::endl;
return EXIT_FAILURE;
}
if(!CGAL::is_triangle_mesh(surface_mesh))
{
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
std::chrono::steady_clock::time_point start_time = std::chrono::steady_clock::now();
// In this example, the simplification stops when the number of undirected edges
// drops below 10% of the initial count
double stop_ratio = (argc > 2) ? std::stod(argv[2]) : 0.1;
SMS::Edge_count_ratio_stop_predicate<Surface_mesh> stop(stop_ratio);
int r = SMS::edge_collapse(surface_mesh, stop);
std::chrono::steady_clock::time_point end_time = std::chrono::steady_clock::now();
std::cout << "\nFinished!\n" << r << " edges removed.\n" << surface_mesh.number_of_edges() << " final edges.\n";
std::cout << "Time elapsed: " << std::chrono::duration_cast<std::chrono::milliseconds>(end_time - start_time).count() << "ms" << std::endl;
CGAL::IO::write_polygon_mesh((argc > 3) ? argv[3] : "out.off", surface_mesh, CGAL::parameters::stream_precision(17));
return EXIT_SUCCESS;
}

Example Using a Default Polyhedron

The following example illustrates the simplification of a Polyhedron_3 with default vertices, halfedges, and faces. The unspecified cost strategy defaults to Lindstrom-Turk.


File Surface_mesh_simplification/edge_collapse_polyhedron.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_3.h>
// Simplification function
#include <CGAL/Surface_mesh_simplification/edge_collapse.h>
// Stop-condition policy
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_count_stop_predicate.h>
#include <iostream>
#include <fstream>
typedef CGAL::Polyhedron_3<Kernel> Surface_mesh;
int main(int argc, char** argv)
{
Surface_mesh surface_mesh;
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/small_cube.off");
std::ifstream is(filename);
if(!is || !(is >> surface_mesh))
{
std::cerr << "Failed to read input mesh: " << filename << std::endl;
return EXIT_FAILURE;
}
if(!CGAL::is_triangle_mesh(surface_mesh))
{
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
// This is a stop predicate (defines when the algorithm terminates).
// In this example, the simplification stops when the number of undirected edges
// left in the surface mesh drops below the specified number (1000)
const std::size_t edge_count_treshold = (argc > 2) ? std::stoi(argv[2]) : 1000;
SMS::Edge_count_stop_predicate<Surface_mesh> stop(edge_count_treshold);
// This the actual call to the simplification algorithm.
// The surface mesh and stop conditions are mandatory arguments.
// The index maps are needed because the vertices and edges
// of this surface mesh lack an "id()" field.
std::cout << "Collapsing edges of Polyhedron: " << filename << ", aiming for " << edge_count_treshold << " final edges..." << std::endl;
int r = SMS::edge_collapse(surface_mesh, stop,
CGAL::parameters::vertex_index_map(get(CGAL::vertex_external_index, surface_mesh))
.halfedge_index_map(get(CGAL::halfedge_external_index, surface_mesh)));
std::cout << "\nFinished!\n" << r << " edges removed.\n"
<< (surface_mesh.size_of_halfedges()/2) << " final edges.\n";
std::ofstream os(argc > 3 ? argv[3] : "out.off");
os.precision(17);
os << surface_mesh;
return EXIT_SUCCESS;
}

Example Using an Enriched Polyhedron

The following example is equivalent to the previous example but using an enriched polyhedron whose halfedges support an id field to store the edge index needed by the algorithm.


File Surface_mesh_simplification/edge_collapse_enriched_polyhedron.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_3.h>
// Extended polyhedron items which include an id() field
#include <CGAL/Polyhedron_items_with_id_3.h>
#include <CGAL/Surface_mesh_simplification/edge_collapse.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_count_ratio_stop_predicate.h>
#include <iostream>
#include <fstream>
typedef Kernel::Point_3 Point;
// Setup an enriched polyhedron type which stores an id() field in the items
typedef boost::graph_traits<Surface_mesh>::vertex_descriptor vertex_descriptor;
typedef boost::graph_traits<Surface_mesh>::halfedge_descriptor halfedge_descriptor;
int main(int argc, char** argv)
{
Surface_mesh surface_mesh;
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/small_cube.off");
std::ifstream is(filename);
if(!is || !(is >> surface_mesh))
{
std::cerr << "Failed to read input mesh: " << filename << std::endl;
return EXIT_FAILURE;
}
if(!CGAL::is_triangle_mesh(surface_mesh))
{
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
// The items in this polyhedron have an "id()" field
// which the default index maps used in the algorithm
// need to get the index of a vertex/edge.
// However, the Polyhedron_3 class doesn't assign any value to
// this id(), so we must do it here:
int index = 0;
for(halfedge_descriptor hd : halfedges(surface_mesh))
hd->id() = index++;
index = 0;
for(vertex_descriptor vd : vertices(surface_mesh))
vd->id() = index++;
// In this example, the simplification stops when the number of undirected edges
// drops below xx% of the initial count
const double ratio = (argc > 2) ? std::stod(argv[2]) : 0.1;
SMS::Edge_count_ratio_stop_predicate<Surface_mesh> stop(ratio);
// The index maps are not explicitelty passed as in the previous
// example because the surface mesh items have a proper id() field.
// On the other hand, we pass here explicit cost and placement
// function which differ from the default policies, omitted in
// the previous example.
std::cout << "Collapsing edges of mesh: " << filename << ", aiming for " << 100 * ratio << "% of the input edges..." << std::endl;
int r = SMS::edge_collapse(surface_mesh, stop);
std::cout << "\nFinished!\n" << r << " edges removed.\n"
<< (surface_mesh.size_of_halfedges()/2) << " final edges.\n";
std::ofstream os((argc > 3) ? argv[3] : "out.off");
os.precision(17);
os << surface_mesh;
return EXIT_SUCCESS;
}

Example for Simplification of OpenMesh

The following example shows how the mesh simplification package can be applied on a mesh data structure which is not part of CGAL, but a model of FaceGraph.

What is particular in this example is the property map that allows to associate 3D CGAL points to the vertices.


File Surface_mesh_simplification/edge_collapse_OpenMesh.cpp

#include <OpenMesh/Core/IO/MeshIO.hh>
#include <OpenMesh/Core/Mesh/PolyMesh_ArrayKernelT.hh>
#include <CGAL/boost/graph/graph_traits_PolyMesh_ArrayKernelT.h>
// Simplification function
#include <CGAL/Surface_mesh_simplification/edge_collapse.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_count_stop_predicate.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_length_cost.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Midpoint_placement.h>
#include <iostream>
#include <fstream>
typedef OpenMesh::PolyMesh_ArrayKernelT</* default traits*/> Surface_mesh;
typedef boost::graph_traits<Surface_mesh>::edge_descriptor edge_descriptor;
class Constrained_edge_map
{
public:
typedef boost::read_write_property_map_tag category;
typedef bool value_type;
typedef bool reference;
typedef edge_descriptor key_type;
Constrained_edge_map(Surface_mesh& sm)
: sm_(sm)
{
sm_.add_property(constraint);
}
inline friend value_type get(const Constrained_edge_map& em, key_type e)
{
bool b = em.sm_.property(em.constraint,em.sm_.edge_handle(e.idx()));
return b;
}
inline friend void put(const Constrained_edge_map& em, key_type e, value_type b)
{
em.sm_.property(em.constraint,em.sm_.edge_handle(e.idx())) = b;
}
private:
Surface_mesh& sm_;
OpenMesh::EPropHandleT<bool> constraint;
};
int main(int argc, char** argv)
{
Surface_mesh surface_mesh;
Constrained_edge_map constraints_map(surface_mesh);
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/cube-meshed.off");
OpenMesh::IO::read_mesh(surface_mesh, filename);
if(!CGAL::is_triangle_mesh(surface_mesh)){
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
// For the purpose of the example we mark 100 edges as constrained edges
int count = 0;
for(edge_descriptor e : edges(surface_mesh))
put(constraints_map, e, (count++ < 100));
// This is a stop predicate (defines when the algorithm terminates).
// In this example, the simplification stops when the number of undirected edges
// left in the surface mesh drops below the specified number (1000)
const std::size_t stop_n = (argc > 2) ? std::stoi(argv[2]) : 1000;
SMS::Edge_count_stop_predicate<Surface_mesh> stop(stop_n);
// This the actual call to the simplification algorithm.
// The surface mesh and stop conditions are mandatory arguments.
std::cout << "Collapsing edges of mesh: " << filename << ", aiming for " << stop_n << " final edges..." << std::endl;
int r = SMS::edge_collapse(surface_mesh, stop,
CGAL::parameters::halfedge_index_map(get(CGAL::halfedge_index,surface_mesh))
.vertex_point_map(get(boost::vertex_point, surface_mesh))
.edge_is_constrained_map(constraints_map));
surface_mesh.garbage_collection();
std::cout << "\nFinished!\n" << r << " edges removed.\n"
<< num_edges(surface_mesh) << " final edges.\n";
OpenMesh::IO::write_mesh(surface_mesh, "out.off");
return EXIT_SUCCESS;
}

Example with Edges Marked as Non-Removable

The following example shows how to use the optional named parameter edge_is_constrained_map to prevent edges from being removed. Edges marked as constrained are guaranteed to be in the final surface mesh. However, the vertices of the constrained edges may change and the placement may change the points. The wrapper CGAL::Surface_mesh_simplification::Constrained_placement guarantees that these points are not changed.


File Surface_mesh_simplification/edge_collapse_constrained_border_surface_mesh.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
// Simplification function
#include <CGAL/Surface_mesh_simplification/edge_collapse.h>
// Midpoint placement policy
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Midpoint_placement.h>
//Placement wrapper
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Constrained_placement.h>
// Stop-condition policy
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_count_stop_predicate.h>
#include <iostream>
#include <fstream>
typedef Kernel::Point_3 Point_3;
typedef CGAL::Surface_mesh<Point_3> Surface_mesh;
typedef boost::graph_traits<Surface_mesh>::halfedge_descriptor halfedge_descriptor;
typedef boost::graph_traits<Surface_mesh>::edge_descriptor edge_descriptor;
// BGL property map which indicates whether an edge is marked as non-removable
struct Border_is_constrained_edge_map
{
const Surface_mesh* sm_ptr;
typedef edge_descriptor key_type;
typedef bool value_type;
typedef value_type reference;
typedef boost::readable_property_map_tag category;
Border_is_constrained_edge_map(const Surface_mesh& sm) : sm_ptr(&sm) {}
friend value_type get(const Border_is_constrained_edge_map& m, const key_type& edge) {
return CGAL::is_border(edge, *m.sm_ptr);
}
};
// Placement class
typedef SMS::Constrained_placement<SMS::Midpoint_placement<Surface_mesh>,
Border_is_constrained_edge_map > Placement;
int main(int argc, char** argv)
{
Surface_mesh surface_mesh;
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/mesh_with_border.off");
std::ifstream is(filename);
if(!is || !(is >> surface_mesh))
{
std::cerr << "Failed to read input mesh: " << filename << std::endl;
return EXIT_FAILURE;
}
if(!CGAL::is_triangle_mesh(surface_mesh))
{
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
Surface_mesh::Property_map<halfedge_descriptor, std::pair<Point_3, Point_3> > constrained_halfedges;
constrained_halfedges = surface_mesh.add_property_map<halfedge_descriptor,std::pair<Point_3, Point_3> >("h:vertices").first;
std::size_t nb_border_edges=0;
for(halfedge_descriptor hd : halfedges(surface_mesh))
{
if(CGAL::is_border(hd, surface_mesh))
{
constrained_halfedges[hd] = std::make_pair(surface_mesh.point(source(hd, surface_mesh)),
surface_mesh.point(target(hd, surface_mesh)));
++nb_border_edges;
}
}
// Contract the surface mesh as much as possible
SMS::Edge_count_stop_predicate<Surface_mesh> stop(0);
Border_is_constrained_edge_map bem(surface_mesh);
// This the actual call to the simplification algorithm.
// The surface mesh and stop conditions are mandatory arguments.
std::cout << "Collapsing as many edges of mesh: " << filename << " as possible..." << std::endl;
int r = SMS::edge_collapse(surface_mesh, stop,
CGAL::parameters::edge_is_constrained_map(bem)
.get_placement(Placement(bem)));
std::cout << "\nFinished!\n" << r << " edges removed.\n"
<< surface_mesh.number_of_edges() << " final edges.\n";
CGAL::IO::write_polygon_mesh((argc > 2) ? argv[2] : "out.off", surface_mesh, CGAL::parameters::stream_precision(17));
// now check!
for(halfedge_descriptor hd : halfedges(surface_mesh))
{
if(CGAL::is_border(hd,surface_mesh))
{
--nb_border_edges;
if(constrained_halfedges[hd] != std::make_pair(surface_mesh.point(source(hd, surface_mesh)),
surface_mesh.point(target(hd, surface_mesh))))
{
std::cerr << "oops. send us a bug report\n";
}
}
}
assert(nb_border_edges==0);
return EXIT_SUCCESS;
}

Example with Bounded Change of Face Normals

The surface mesh simplification does not guarantee that the resulting surface has no self intersections. Even the rather trivial mesh shown in Figure 71.2 results in a self intersection when one edge is collapsed using the Lindstrom-Turk method.

SMS-selfintersection.png
Figure 71.2 Simple mesh before and after the collapse of edge v-w into vertex w. While the normals of f1 and f2 are almost equal, they are opposed after the edge collapse.

The class Surface_mesh_simplification::Bounded_normal_change_filter checks if a placement would invert the normal of a face around the stars of the two vertices of an edge that is candidate for an edge collapse. It then rejects this placement by returning boost::none.

Note
This filter class replaces the usage of the class Surface_mesh_simplification::Bounded_normal_change_placement. Using the filter is faster as it is only performed on the edge to be collapsed next, and not during the update of all edges incident to the vertex that is the result of the edge collapse.


File Surface_mesh_simplification/edge_collapse_bounded_normal_change.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Timer.h>
#include <CGAL/Surface_mesh_simplification/edge_collapse.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_count_stop_predicate.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_cost.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_placement.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounded_normal_change_filter.h>
#include <iostream>
#include <fstream>
struct Dummy_placement {
template <typename Profile>
boost::optional<typename Profile::Point> operator()(const Profile&) const
{
return boost::none;
}
template <typename Profile>
boost::optional<typename Profile::Point> operator()(const Profile&, const boost::optional<typename Profile::Point>& op) const
{
return op;
}
};
int main(int argc, char** argv)
{
Surface_mesh surface_mesh;
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/fold.off");
std::ifstream is(filename);
if(!is || !(is >> surface_mesh))
{
std::cerr << "Failed to read input mesh: " << filename << std::endl;
return EXIT_FAILURE;
}
if(!CGAL::is_triangle_mesh(surface_mesh))
{
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
// This is a stop predicate (defines when the algorithm terminates).
// In this example, the simplification stops when the number of undirected edges
// left in the surface mesh drops below the specified number
const std::size_t stop_n = (argc > 2) ? std::stoi(argv[2]) : num_edges(surface_mesh)/2 - 1;
SMS::Edge_count_stop_predicate<Surface_mesh> stop(stop_n);
typedef SMS::LindstromTurk_placement<Surface_mesh> Placement;
CGAL::Timer t;
t.start();
// This the actual call to the simplification algorithm.
// The surface mesh and stop conditions are mandatory arguments.
// The index maps are needed because the vertices and edges
// of this surface mesh lack an "id()" field.
std::cout << "Collapsing edges of mesh: " << filename << ", aiming for " << stop_n << " final edges..." << std::endl;
SMS::Bounded_normal_change_filter<> filter;
SMS::edge_collapse(surface_mesh, stop,
CGAL::parameters::get_cost(SMS::LindstromTurk_cost<Surface_mesh>())
.filter(filter)
.get_placement(Placement()));
std::cout << t.time() << " sec" << std::endl;
CGAL::IO::write_polygon_mesh((argc > 3) ? argv[3] : "out.off", surface_mesh, CGAL::parameters::stream_precision(17));
return EXIT_SUCCESS;
}

Example with Polyhedral Envelope

The surface mesh simplification can be done in a way that the simplified mesh stays inside an envelope of the input mesh. This makes use of the class Polyhedral_envelope which enables to check whether a query point, segment, or triangle lies within a polyhedral envelope, which consists of the union of inflated triangles. While the user gives a tolerance \(\epsilon\), the check is conservative, that is there may be triangles which would be inside a surface obtained as the Minkowski sum envelope with a sphere of radius \(\epsilon\), but which are outside the polyhedral envelope.


File Surface_mesh_simplification/edge_collapse_envelope.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Surface_mesh_simplification/edge_collapse.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_count_stop_predicate.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_cost.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_placement.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounded_normal_change_filter.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Polyhedral_envelope_filter.h>
//bbox
#include <CGAL/Polygon_mesh_processing/bbox.h>
#include <iostream>
#include <fstream>
typedef Kernel::Point_3 Point_3;
typedef SMS::LindstromTurk_cost<Surface> Cost;
typedef SMS::LindstromTurk_placement<Surface> Placement;
typedef SMS::Polyhedral_envelope_filter<Kernel,SMS::Bounded_normal_change_filter<> > Filter;
int main(int argc, char** argv)
{
Surface mesh;
std::ifstream is(argc > 1 ? argv[1] : CGAL::data_file_path("meshes/helmet.off"));
is >> mesh;
SMS::Edge_count_stop_predicate<Surface> stop(0); // go as far as you can while in the envelope
Point_3 cmin = (bbox.min)();
Point_3 cmax = (bbox.max)();
const double diag = CGAL::approximate_sqrt(CGAL::squared_distance(cmin, cmax));
std::cout << "eps = " << 0.01*diag << std::endl;
Placement placement;
Filter filter(0.01*diag);
SMS::edge_collapse(mesh, stop, CGAL::parameters::get_cost(Cost()).filter(filter).get_placement(placement));
std::ofstream out("out.off");
out << mesh << std::endl;
out.close();
return EXIT_SUCCESS;
}

Example with Visitor

The last example shows how to use a visitor with callbacks that are called at the different steps of the simplification algorithm.


File Surface_mesh_simplification/edge_collapse_visitor_surface_mesh.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
// Simplification function
#include <CGAL/Surface_mesh_simplification/edge_collapse.h>
// Visitor base
#include <CGAL/Surface_mesh_simplification/Edge_collapse_visitor_base.h>
// Stop-condition policy
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_count_ratio_stop_predicate.h>
#include <iostream>
#include <fstream>
typedef Kernel::Point_3 Point_3;
typedef CGAL::Surface_mesh<Point_3> Surface_mesh;
typedef boost::graph_traits<Surface_mesh>::halfedge_descriptor halfedge_descriptor;
typedef boost::graph_traits<Surface_mesh>::vertex_descriptor vertex_descriptor;
typedef SMS::Edge_profile<Surface_mesh> Profile;
// The following is a Visitor that keeps track of the simplification process.
// In this example the progress is printed real-time and a few statistics are
// recorded (and printed in the end).
//
struct Stats
{
std::size_t collected = 0;
std::size_t processed = 0;
std::size_t collapsed = 0;
std::size_t non_collapsable = 0;
std::size_t cost_uncomputable = 0;
std::size_t placement_uncomputable = 0;
};
struct My_visitor : SMS::Edge_collapse_visitor_base<Surface_mesh>
{
My_visitor(Stats* s) : stats(s) {}
// Called during the collecting phase for each edge collected.
void OnCollected(const Profile&, const boost::optional<double>&)
{
++(stats->collected);
std::cerr << "\rEdges collected: " << stats->collected << std::flush;
}
// Called during the processing phase for each edge selected.
// If cost is absent the edge won't be collapsed.
void OnSelected(const Profile&,
boost::optional<double> cost,
std::size_t initial,
std::size_t current)
{
++(stats->processed);
if(!cost)
++(stats->cost_uncomputable);
if(current == initial)
std::cerr << "\n" << std::flush;
std::cerr << "\r" << current << std::flush;
}
// Called during the processing phase for each edge being collapsed.
// If placement is absent the edge is left uncollapsed.
void OnCollapsing(const Profile&,
boost::optional<Point> placement)
{
if(!placement)
++(stats->placement_uncomputable);
}
// Called for each edge which failed the so called link-condition,
// that is, which cannot be collapsed because doing so would
// turn the surface mesh into a non-manifold.
void OnNonCollapsable(const Profile&)
{
++(stats->non_collapsable);
}
// Called after each edge has been collapsed
void OnCollapsed(const Profile&, vertex_descriptor)
{
++(stats->collapsed);
}
Stats* stats;
};
int main(int argc, char** argv)
{
Surface_mesh surface_mesh;
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/small_cube.off");
std::ifstream is(filename);
if(!is || !(is >> surface_mesh))
{
std::cerr << "Failed to read input mesh: " << filename << std::endl;
return EXIT_FAILURE;
}
if(!CGAL::is_triangle_mesh(surface_mesh))
{
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
// In this example, the simplification stops when the number of undirected edges
// drops below xx% of the initial count
const double ratio = (argc > 2) ? std::stod(argv[2]) : 0.1;
SMS::Edge_count_ratio_stop_predicate<Surface_mesh> stop(ratio);
Stats stats;
My_visitor vis(&stats);
// The index maps are not explicitelty passed as in the previous
// example because the surface mesh items have a proper id() field.
// On the other hand, we pass here explicit cost and placement
// function which differ from the default policies, omitted in
// the previous example.
int r = SMS::edge_collapse(surface_mesh, stop, CGAL::parameters::visitor(vis));
std::cout << "\nEdges collected: " << stats.collected
<< "\nEdges processed: " << stats.processed
<< "\nEdges collapsed: " << stats.collapsed
<< std::endl
<< "\nEdges not collapsed due to topological constraints: " << stats.non_collapsable
<< "\nEdge not collapsed due to cost computation constraints: " << stats.cost_uncomputable
<< "\nEdge not collapsed due to placement computation constraints: " << stats.placement_uncomputable
<< std::endl;
std::cout << "\nFinished!\n" << r << " edges removed.\n"
<< surface_mesh.number_of_edges() << " final edges.\n";
CGAL::IO::write_polygon_mesh((argc > 3) ? argv[3] : "out.off", surface_mesh, CGAL::parameters::stream_precision(17));
return EXIT_SUCCESS;
}

Example Using Garland-Heckbert Policies

Each Garland-Heckbert simplification strategy is implemented with a single class that regroups both the cost and the placement policies, which must be used together as they share vertex quadric data. The classic strategy of using plane-based quadric error metrics is implemented with the class Surface_mesh_simplification::GarlandHeckbert_plane_policies. Although both policies must be used together, it is still possible to wrap either policy using behavior modifiers such as Surface_mesh_simplification::Bounded_normal_change_placement.


File Surface_mesh_simplification/edge_collapse_garland_heckbert.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_count_ratio_stop_predicate.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounded_normal_change_placement.h>
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/GarlandHeckbert_policies.h>
#include <CGAL/Surface_mesh_simplification/edge_collapse.h>
#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
#include <chrono>
#include <fstream>
#include <iostream>
#include <vector>
typedef Kernel::FT FT;
typedef Kernel::Point_3 Point_3;
typedef CGAL::Surface_mesh<Point_3> Surface_mesh;
typedef SMS::GarlandHeckbert_plane_policies<Surface_mesh, Kernel> Classic_plane;
typedef SMS::GarlandHeckbert_probabilistic_plane_policies<Surface_mesh, Kernel> Prob_plane;
typedef SMS::GarlandHeckbert_triangle_policies<Surface_mesh, Kernel> Classic_tri;
typedef SMS::GarlandHeckbert_probabilistic_triangle_policies<Surface_mesh, Kernel> Prob_tri;
template <typename GHPolicies>
void collapse_gh(Surface_mesh& mesh,
const double ratio)
{
std::chrono::steady_clock::time_point start_time = std::chrono::steady_clock::now();
SMS::Edge_count_ratio_stop_predicate<Surface_mesh> stop(ratio);
// Garland&Heckbert simplification policies
typedef typename GHPolicies::Get_cost GH_cost;
typedef typename GHPolicies::Get_placement GH_placement;
typedef SMS::Bounded_normal_change_placement<GH_placement> Bounded_GH_placement;
GHPolicies gh_policies(mesh);
const GH_cost& gh_cost = gh_policies.get_cost();
const GH_placement& gh_placement = gh_policies.get_placement();
Bounded_GH_placement placement(gh_placement);
int r = SMS::edge_collapse(mesh, stop,
CGAL::parameters::get_cost(gh_cost)
.get_placement(placement));
std::chrono::steady_clock::time_point end_time = std::chrono::steady_clock::now();
std::cout << "Time elapsed: "
<< std::chrono::duration_cast<std::chrono::milliseconds>(end_time - start_time).count()
<< "ms" << std::endl;
std::cout << "\nFinished!\n" << r << " edges removed.\n" << edges(mesh).size() << " final edges.\n";
}
// Usage:
// ./command [input] [ratio] [policy] [output]
// policy can be "cp" (classic plane), "ct" (classic triangle), "pp" (probabilistic plane), "pt" (probabilistic triangle)
int main(int argc, char** argv)
{
Surface_mesh mesh;
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/cube-meshed.off");
if(!CGAL::IO::read_polygon_mesh(filename, mesh))
{
std::cerr << "Failed to read input mesh: " << filename << std::endl;
return EXIT_FAILURE;
}
{
std::cerr << "Input geometry is not triangulated." << std::endl;
return EXIT_FAILURE;
}
std::cout << "Input mesh has " << num_vertices(mesh) << " nv "
<< num_edges(mesh) << " ne "
<< num_faces(mesh) << " nf" << std::endl;
const double ratio = (argc > 2) ? std::stod(argv[2]) : 0.2;
std::cout << "Collapsing edges of mesh: " << filename << ", aiming for " << 100 * ratio << "% of the input edges..." << std::endl;
const std::string policy = (argc > 3) ? argv[3] : "cp";
if(policy == "cp")
collapse_gh<Classic_plane>(mesh, ratio);
else if(policy == "ct")
collapse_gh<Classic_tri>(mesh, ratio);
else if(policy == "pp")
collapse_gh<Prob_plane>(mesh, ratio);
else
collapse_gh<Prob_tri>(mesh, ratio);
CGAL::IO::write_polygon_mesh((argc > 4) ? argv[4] : "out.off", mesh, CGAL::parameters::stream_precision(17));
return EXIT_SUCCESS;
}

Note that these policies depend on the third party Eigen library.

Design and Implementation History

The core of the package, as well as most of the simplification strategies, are the work of Fernando Cacciola, between 2006 and 2009.

Andreas Fabri added the Surface_mesh_simplification::Bounded_normal_change_placement functionality in CGAL 4.11.

The implementation of the Garland-Heckbert simplification policies is the result of the work of Baskın Şenbaşlar (Google Summer of Code 2019), and Julian Komaromy (Google Summer of Code 2021) They both were mentored by Mael Rouxel-Labbé, who also contributed to the code and to the documentation.