Chapter 49
Triangulated Surface Mesh Simplification

Fernando Cacciola

Table of Contents

49.1   Introduction

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 to be of a particular type. Instead, it defines the concept of a EdgeCollapsableMesh, which presents the surface as being a halfedge data structure, and any surface that is a model of that concept can be simplified. The concept is defined not in terms of a monolithic class, but in terms of a set of functions and traits, making it easy to adapt any concrete surface type, even if it is not a halfedge data structure at all. In particular, the concept definition follows the design of the Boost Graph Library (Bgl) [SLL02].

The design is policy-based (http://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.

49.2   Overview of the Simplification Process

The free function that implements the simplification algorithm takes not only the surface 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', dissapearing '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 that results from an edge collapse deviates from the initial surface by some amount, and since the goal of simplification is to reduce the number of triangles while retaining the overall look of the surface as much as possible, it is necessary to measure such a deviation. Some methods attempt to measure the total deviation from the initial surface to the completely simplified surface, for example, by tracking an accumulated error while keeping a history of the simplification changes. Other methods, like the one 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 one in this package, produce slightly less accurate simplifications but take up much less additional space, even none in some cases.

The cost-driven method implemented in this package is mainly based on [LT98, LT99], with contributions from [HDD⁺93], [GH97] and [DEGN99].

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. 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 doesn't satisfy certain topological and geometric conditions.

The algorithm presented in [GH97] 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 [LT98, LT99] ( which is the algorithm implemented here). However, contracting an arbitrary vertex-pair may result in a non-manifold surface, but the current state of this package can only deal with manifold surfaces, thus, it can only collapse edges. That is, this package cannot be used as a framework for vertex contraction.

49.3   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 two strategies: the Lindstrom-Turk strategy, which is the default, and a strategy consisting of an edge-length cost with an optional midpoint placement (much faster but less accurate).

49.3.1   Lindstrom-Turk Cost and Placement Strategy

The main characteristic of the strategy presented in [LT98, LT99] is that the simplified surface is not compared at each step with the original surface (or the surface at a previous step) so there is no need to keep extra information, such as the original surface 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 boundaries (in case the edge profile has boundary edges). The next constraints preserve the total volume of the surface. 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.

49.3.2   Cost Strategy Policies

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

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, in which case the remaining vertex is kept in its place.

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 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 only the cost is attached to each edge and nothing else.

But this is a tradeoff: 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 shown 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.

49.4   API

49.4.1   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. ¹

The simplification algorithm is implemented as the free template function CGAL::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 to be simplified (in-place) and the stop predicate.

The surface to simplify must be a model of the EdgeCollapsableMesh concept. Many concrete surface types, such as CGAL::Polyhedron_3 with only triangular faces, become models of that concept via a technique known as external adaptation, which is described in [SLL02] and this Bgl web page: http://www.boost.org/libs/graph/doc/leda_conversion.html

External adaptation is a way to add an interface to an object without coercing the type of the object (which happens when you adapt it by means of a wrapper). That is, the formal parameter to the edge_collapse function that implements the simplification is the concrete surface object itself, not an adaptor which delegates the functionality to the concrete type.

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 [SLL02] or the following site: http://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                  : the surface 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_border_map(ebmap): property-map specifying whether an edge is a border 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 
visitor(vis)             : function object tracking the simplification process 
*/
int r = edge_collapse(surface
                     ,stop_predicate
                     ,vertex_index_map(vimap)
                       .edge_index_map(eimap)
                       .edge_is_border_map(ebmap)
                       .get_cost(cf)
                       .get_placement(pf)
                       .visitor(vis)
                     );

49.4.2   Examples

Example Using a Default Polyhedron

The following example illustrates the simplest of the cases. It uses an ordinary polyhedron and only one of the optional parameters. The unspecified cost strategy defaults to Lindstrom-Turk.

File: examples/Surface_mesh_simplification/edge_collapse_polyhedron.cpp

#include <iostream>
#include <fstream>

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/IO/Polyhedron_iostream.h>

// Adaptor for Polyhedron_3
#include <CGAL/Surface_mesh_simplification/HalfedgeGraph_Polyhedron_3.h>

// Simplification function
#include <CGAL/Surface_mesh_simplification/edge_collapse.h>

// Stop-condition policy
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Count_stop_predicate.h>

typedef CGAL::Simple_cartesian<double> Kernel;
typedef CGAL::Polyhedron_3<Kernel> Surface; 

namespace SMS = CGAL::Surface_mesh_simplification ;

int main( int argc, char** argv ) 
{
  Surface surface; 
  
  std::ifstream is(argv[1]) ; is >> surface ;

  // 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 drops below the specified number (1000)
  SMS::Count_stop_predicate<Surface> stop(1000);
     
  // This the actual call to the simplification algorithm.
  // The surface and stop conditions are mandatory arguments.
  // The index maps are needed because the vertices and edges
  // of this surface lack an "id()" field.
  int r = SMS::edge_collapse
            (surface
            ,stop
            ,CGAL::vertex_index_map(boost::get(CGAL::vertex_external_index,surface)) 
                  .edge_index_map  (boost::get(CGAL::edge_external_index  ,surface)) 
            );
  
  std::cout << "\nFinished...\n" << r << " edges removed.\n" 
            << (surface.size_of_halfedges()/2) << " final edges.\n" ;
        
  std::ofstream os( argc > 2 ? argv[2] : "out.off" ) ; os << surface ;
  
  return 0 ;      
}

// EOF //

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. It also shows how to explicitly specify a cost strategy (other than the default) and how to use a visitor object to track the simplification process.

File: examples/Surface_mesh_simplification/edge_collapse_enriched_polyhedron.cpp

#include <iostream>
#include <fstream>

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/IO/Polyhedron_iostream.h>

// Adaptor for Polyhedron_3
#include <CGAL/Surface_mesh_simplification/HalfedgeGraph_Polyhedron_3.h>

// Simplification function
#include <CGAL/Surface_mesh_simplification/edge_collapse.h>

// Visitor base
#include <CGAL/Surface_mesh_simplification/Edge_collapse_visitor_base.h>

// Extended polyhedron items which include an id() field
#include <CGAL/Polyhedron_items_with_id_3.h>

// Stop-condition policy
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Count_ratio_stop_predicate.h>

// Non-default cost and placement policies
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Midpoint_and_length.h> 

typedef CGAL::Simple_cartesian<double> Kernel;

typedef Kernel::Point_3 Point ;

//
// Setup an enriched polyhedron type which stores an id() field in the items
//
typedef CGAL::Polyhedron_3<Kernel,CGAL::Polyhedron_items_with_id_3> Surface; 

typedef Surface::Halfedge_handle Halfedge_handle ;
typedef Surface::Vertex_handle   Vertex_handle ;

namespace SMS = CGAL::Surface_mesh_simplification ;

typedef SMS::Edge_profile<Surface> 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 My_visitor : SMS::Edge_collapse_visitor_base<Surface>
{
  My_visitor() 
    : collected(0)
    , processed(0)
    , collapsed(0)
    , non_collapsable(0)
    , cost_uncomputable(0) 
    , placement_uncomputable(0) 
  {} 

  // Called during the collecting phase for each edge collected.
  virtual void OnCollected( Profile const&, boost::optional<double> const& ) const
  {
    ++ collected ;
    std::cerr << "\rEdges collected: " << collected << std::flush ;
  }                
  
  // Called during the processing phase for each edge selected.
  // If cost is absent the edge won't be collapsed.
  virtual void OnSelected(Profile const&          
                         ,boost::optional<double> cost
                         ,std::size_t             initial
                         ,std::size_t             current
                         ) const
  {
    ++ processed ;
    if ( !cost )
      ++ 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.
  virtual void OnCollapsing(Profile const&          
                           ,boost::optional<Point>  placement
                           ) const
  {
    if ( !placement )
      ++ 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 into a non-manifold.
  virtual void OnNonCollapsable( Profile const& ) const
  {
    ++ non_collapsable;
  }                
  
  // Called AFTER each edge has been collapsed
  virtual void OnCollapsed( Profile const&, Vertex_handle hv ) const
  {
    ++ collapsed;
  }                
  
  mutable std::size_t collected ;
  mutable std::size_t processed ;
  mutable std::size_t collapsed ;
  mutable std::size_t non_collapsable ;
  mutable std::size_t cost_uncomputable  ;
  mutable std::size_t placement_uncomputable ; 
} ;


int main( int argc, char** argv ) 
{
  Surface surface; 
  
  std::ifstream is(argv[1]) ; is >> surface ;

  // 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( Surface::Halfedge_iterator eb = surface.halfedges_begin()
     , ee = surface.halfedges_end()
     ; eb != ee
     ; ++ eb
     ) 
    eb->id() = index++;

  index = 0 ;
  for( Surface::Vertex_iterator vb = surface.vertices_begin()
     , ve = surface.vertices_end()
     ; vb != ve
     ; ++ vb
     ) 
    vb->id() = index++;
    
  // In this example, the simplification stops when the number of undirected edges
  // drops below 10% of the initial count
  SMS::Count_ratio_stop_predicate<Surface> stop(0.1);

  My_visitor vis ;
    
  // The index maps are not explicitelty passed as in the previous
  // example because the surface items have a proper id() field.
  // On the other hand, we pass here explicit cost and placement
  // function which differ from the default policies, ommited in
  // the previous example.
  int r = SMS::edge_collapse
           (surface
           ,stop
           ,CGAL::get_cost     (SMS::Edge_length_cost  <Surface>())
                 .get_placement(SMS::Midpoint_placement<Surface>())
                 .visitor      (vis)
           );
  
  std::cout << "\nEdges collected: " << vis.collected
            << "\nEdges proccessed: " << vis.processed
            << "\nEdges collapsed: " << vis.collapsed
            << std::endl
            << "\nEdges not collapsed due to topological constrians: " 
            << vis.non_collapsable
            << "\nEdge not collapsed due to cost computation constrians: " 
            << vis.cost_uncomputable 
            << "\nEdge not collapsed due to placement computation constrians: " 
            << vis.placement_uncomputable 
            << std::endl ; 
            
  std::cout << "\nFinished...\n" << r << " edges removed.\n" 
            << (surface.size_of_halfedges()/2) << " final edges.\n" ;
        
  std::ofstream os( argc > 2 ? argv[2] : "out.off" ) ; os << surface ;
  
  return 0 ;      
}

// EOF //

Example with edges marked as non-removable

The following example shows how to use the optional named parameter edge_is_border_map to prevent edges from being removed even if they are not really borders.

File: examples/Surface_mesh_simplification/edge_collapse_constrained_polyhedron.cpp

#include <iostream>
#include <fstream>

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/IO/Polyhedron_iostream.h>

// Adaptor for Polyhedron_3
#include <CGAL/Surface_mesh_simplification/HalfedgeGraph_Polyhedron_3.h>

// Simplification function
#include <CGAL/Surface_mesh_simplification/edge_collapse.h>

// Stop-condition policy
#include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Count_stop_predicate.h>

// Map used to mark edges as fixed
#include <CGAL/Unique_hash_map.h>

typedef CGAL::Simple_cartesian<double> Kernel;
typedef CGAL::Polyhedron_3<Kernel> Surface; 

namespace SMS = CGAL::Surface_mesh_simplification ;

//
// BGL property map which indicates whether an edge is border OR is marked as non-removable
//
class Constrains_map : public boost::put_get_helper<bool,Constrains_map>
{
public:

  typedef boost::readable_property_map_tag                    category;
  typedef bool                                                value_type;
  typedef bool                                                reference;
  typedef boost::graph_traits<Surface const>::edge_descriptor key_type;

  Constrains_map() : mConstrains(false) {}

  reference operator[](key_type const& e) const { return e->is_border() || is_constrained(e) ; }
  
  void set_is_constrained ( key_type const& e, bool is )  { mConstrains[e]=is; }
  
  bool is_constrained( key_type const& e ) const { return mConstrains.is_defined(e) ? mConstrains[e] : false ; }
  
private:
  
  CGAL::Unique_hash_map<key_type,bool> mConstrains ;
  
};

int main( int argc, char** argv ) 
{
  Surface surface; 
  
  std::ifstream is(argv[1]) ; is >> surface ;

  // 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 drops below the specified number (1000)
  SMS::Count_stop_predicate<Surface> stop(10);
     
  Constrains_map constrains_map ;
     
     
  // This example marks ALL edges as non-removable, but a real world application would mark only selected ones.   
  for( Surface::Halfedge_iterator eb = surface.halfedges_begin(), ee = surface.halfedges_end() ; eb != ee ; ++ eb ) 
    constrains_map.set_is_constrained(eb,true);
     
  // This the actual call to the simplification algorithm.
  // The surface and stop conditions are mandatory arguments.
  // The index maps are needed because the vertices and edges
  // of this surface lack an "id()" field.
  int r = SMS::edge_collapse
            (surface
            ,stop
            ,CGAL::vertex_index_map(boost::get(CGAL::vertex_external_index,surface)) 
                  .edge_index_map  (boost::get(CGAL::edge_external_index  ,surface)) 
                  .edge_is_border_map(constrains_map)
            );
  
  std::cout << "\nFinished...\n" << r << " edges removed.\n" 
            << (surface.size_of_halfedges()/2) << " final edges.\n" ;
        
  std::ofstream os( argc > 2 ? argv[2] : "out.off" ) ; os << surface ;
  
  return 0 ;      
}

// EOF //

Footnotes