CGAL 4.12.1 - 3D Surface Subdivision Methods
User Manual

# Introduction

Subdivision methods are simple yet powerful ways to generate smooth surfaces from arbitrary polygonal meshes. Unlike spline-based surfaces (e.g NURBS) or other numeric-based modeling techniques, users of subdivision methods do not need the mathematical knowledge of the subdivision methods. The natural intuition of the geometry suffices to control the subdivision methods.

Subdivision_method_3 works for the classes Polyhedron_3 and Surface_mesh, as they are models of the concept MutableFaceGraph, and it aims to be easy to use and to extend. Subdivision_method_3 is not a class, but a namespace which contains four popular subdivision methods and their refinement functions. These include Catmull-Clark, Loop, Doo-Sabin and $$\sqrt{3}$$ subdivisions. Variations of these methods can be easily extended by substituting the geometry computation of the refinement host.

# Subdivision Method

In this chapter, we explain some fundamentals of subdivision methods. We focus only on the topics that help you understand the design of the package. Details on subdivision methods can be found in [6]. Some terminology introduced in this section will be used again in later sections. If you are only interested in using a specific subdivision method, Section A Quick Example: Catmull-Clark Subdivision gives a quick tutorial on using a Catmull-Clark subdivision.

A subdivision method recursively refines a coarse mesh and generates an ever closer approximation to a smooth surface. The coarse mesh can have arbitrary shape, but it has to be a 2-manifold. In a 2-manifold, every interior point has a neighborhood homeomorphic to a 2D disk. Subdivision methods on non-manifolds have been developed, but are not considered in Subdivision_method_3. The chapter teaser shows the steps of Catmull-Clark subdivision on a CAD model. The coarse mesh is repeatedly refined by a quadrisection pattern, and new points are generated to approximate a smooth surface.

Many refinement patterns are used in practice. Subdivision_method_3 supports the four most popular patterns, and each of them is used by Catmull-Clark[1], Loop[4], Doo-Sabin[2] and $$\sqrt{3}$$ subdivision[3] (left to right in the figure below). We name these patterns by their topological characteristics instead of the associated subdivision methods. PQQ indicates the Primal Quadtrateral Quadrisection. PTQ indicates the Primal Triangle Quadrisection. DQQ indicates the Dual Quadtrateral Quadrisection. $$\sqrt{3}$$ indicates the converging speed of the triangulation toward the subdivision surface.

The figure demonstrates these four refinement patterns on the 1-disk of a valence-5 vertex/face. Refined meshes are shown below the source meshes. Points on the refined mesh are generated by averaging neighbor points on the source mesh. A graph, called stencil, determines the source neighborhood whose points contribute to the position of a refined point. A refinement pattern usually defines more than one stencil. For example, the PQQ refinement has a vertex-node stencil, which defines the 1-ring of an input vertex; an edge-node stencil, which defines the 1-ring of an input edge; and a face-node stencil, which defines an input face. The stencils of the PQQ refinement are shown in the following figure. The blue neighborhoods in the top row indicate the corresponding stencils of the refined nodes in red.

Stencils with weights are called geometry masks. A subdivision method defines a geometry mask for each stencil, and generates new points by averaging source points weighted by the mask. Geometry masks are carefully chosen to meet requirements of certain surface smoothness and shape quality. The geometry masks of Catmull-Clark subdivision are shown below.

The weights shown here are unnormalized, and $$n$$ is the valence of the vertex. The generated point, in red, is computed by a summation of the weighted points. For example, a Catmull-Clark face-node is computed by the summation of $$1/4$$ of each point on its stencil.

A stencil can have an unlimited number of geometry masks. For example, a face-node of PQQ refinement may be computed by the summation of $$1/5$$ of each stencil node instead of $$1/4$$. Although it is legal in Subdivision_method_3 to have any kind of geometry mask, the result surfaces may be odd, not smooth, or not even exist. [6] explains the details on designing masks for a quality subdivision surface.

# A Quick Example: Catmull-Clark Subdivision

Assuming that you are familiar with Surface_mesh, you can integrate Subdivision_method_3 into your program without much effort.

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/boost/graph/graph_traits_Surface_mesh.h>
#include <CGAL/Timer.h>
#include <boost/lexical_cast.hpp>
#include <iostream>
#include <fstream>
using namespace std;
using namespace CGAL;
namespace params = CGAL::parameters;
int main(int argc, char** argv) {
if (argc > 4) {
cerr << "Usage: CatmullClark_subdivision [d] [filename_in] [filename_out] \n";
cerr << " d -- the depth of the subdivision (default: 1) \n";
cerr << " filename_in -- the input mesh (.off) (default: data/quint_tris.off) \n";
cerr << " filename_out -- the output mesh (.off) (default: result.off)" << endl;
return 1;
}
int d = (argc > 1) ? boost::lexical_cast<int>(argv[1]) : 1;
const char* in_file = (argc > 2) ? argv[2] : "data/quint_tris.off";
const char* out_file = (argc > 3) ? argv[3] : "result.off";
PolygonMesh pmesh;
std::ifstream in(in_file);
if(in.fail()) {
std::cerr << "Could not open input file " << in_file << std::endl;
return 1;
}
in >> pmesh;
Timer t;
t.start();
Subdivision_method_3::CatmullClark_subdivision(pmesh, params::number_of_iterations(d));
std::cerr << "Done (" << t.time() << " s)" << std::endl;
std::ofstream out(out_file);
out << pmesh;
return 0;
}

This example demonstrates the use of the Catmull-Clark subdivision method on a Surface_mesh. There is only one line deserving a detailed explanation:

Subdivision_method_3::CatmullClark_subdivision(pmesh, params::number_of_iterations(d));

Subdivision_method_3 specifies the namespace of the subdivision functions. CatmullClark_subdivision(P, params::number_of_iterations(d)) computes the Catmull-Clark subdivision surface of the polygon mesh pmesh after d iterations of the refinements. The polygon mesh pmesh is passed by reference, and is modified (i.e. subdivided) by the subdivision function.

This example shows how to subdivide a Surface_mesh with Subdivision_method_3. An application-defined polygon mesh might use a specialized kernel and/or a specialized internal container. There is one major restriction on the application-defined polygon mesh to work with Subdivision_method_3: The primitives (such as vertices, halfedges and faces) in the internal container are sequentially ordered (e.g. std::vector and std::list). This implies that the iterators traverse the primitives in the order of their creations/insertions.

Section Refinement Host gives detailed explanations on this two restrictions.

# Catmull-Clark Subdivision

Subdivision_method_3 is designed to allow customization of the subdivision methods. This section explains the implementation of the Catmull-Clark subdivision function in Subdivision_method_3. The implementation demonstrates the customization of the PQQ refinement to Catmull-Clark subdivision.

When a subdivision method is developed, a refinement pattern is chosen, and then a set of geometry masks is developed to position the new points. There are three key components to implement a subdivision method:

• a mesh data structure that can represent arbitrary 2-manifolds,
• a process that refines the mesh data structure,
• and the geometry masks that compute the new points.

E. Catmull and J. Clark chose the PQQ refinement for their subdivision method, and developed a set of geometry masks to generate (or more precisely, to approximate) the B-spline surface from the control mesh. Subdivision_method_3 provides a function that glues all three components of the Catmull-Clark subdivision method.

template <class PolygonMesh, class Mask, class NamedParameters>

PolygonMesh must be an instantiation of Polyhedron_3, Surface_mesh, or any other model of the concept MutableFaceGraph. It is a generic mesh data structure for arbitrary 2-manifolds. PQQ(), which refines the control mesh p, is a refinement host that uses a policy class Mask<PolygonMesh> as part of it geometry computation. During the refinement, PQQ() computes and assigns new points by cooperating with the mask. To implement Catmull-Clark subdivision, Mask, the geometry policy, has to realize the geometry masks of Catmull-Clark subdivision. The number of iterations as well as the vertex point map can be specified using the named parameter np.

To implement the geometry masks, we need to know how a refinement host communicates with its geometry masks. The PQQ refinement defines three stencils, and hence three geometry masks are required for Catmull-Clark subdivision. The following class defines the interfaces of the stencils for the PQQ refinement.

template <class PolygonMesh>
void face_node(boost::graph_traits<PolygonMesh>::face_descriptor face, Point_3& pt);
void edge_node(boost::graph_traits<PolygonMesh>::halfedge_descriptor edge, Point_3& pt);
void vertex_node(boost::graph_traits<PolygonMesh>::vertex_descriptor vertex, Point_3& pt);
};

Each class function in PQQMask_3 computes a new point based on the neighborhood of the primitive descriptor, and assigns the new point to Point_3& pt.

We realize each class function with the geometry masks of Catmull-Clark subdivision.

template <class PolygonMesh>
typedef boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
typedef boost::property_traits<VertexPointMap>::value_type Point_3;
Polygonmesh pmesh;
VertexPointMap vpm;
: pmesh(pmesh, vpm(vpm)
{}
void face_node(face_descriptor face, Point_3& pt) {
int n = 0;
Point_3 p(0,0,0);
for(halfedge_descriptor hd : halfedges_around_face(face,pmesh)){
p = p + get(vpm, (target(hd,pmesh)) - ORIGIN);
++n;
}
pt = ORIGIN + (p - ORIGIN)/FT(n);
}
};
void edge_node(halfedge_descriptor edge, Point_3& pt) {
Point_3 p1 = get(vpm,target(edge, pmesh));
Point_3 p2 = get(vpm,source(edge, pmesh));
Point_3 f1, f2;
face_node(face(edge,pmesh), f1);
face_node(face(opposite(edge,pmesh),pmesh), f2);
pt = Point_3((p1[0]+p2[0]+f1[0]+f2[0])/4,
(p1[1]+p2[1]+f1[1]+f2[1])/4,
(p1[2]+p2[2]+f1[2]+f2[2])/4 );
}
void vertex_node(vertex_descriptor vertex, Point_3& pt) {
Halfedge_around_target_circulator<PolygonMesh> vcir(vertex,pmesh);
typename boost::graph_traits<PolygonMesh>::degree_size_type n = degree(vertex,pmesh);
FT Q[] = {0.0, 0.0, 0.0}, R[] = {0.0, 0.0, 0.0};
Point_3 S = get(vpm,vertex);
Point_3 q;
for (int i = 0; i < n; i++, ++vcir) {
Point_3 p2 = get(vpm, target(opposite(*vcir,pmesh),pmesh));
R[0] += (S[0]+p2[0])/2;
R[1] += (S[1]+p2[1])/2;
R[2] += (S[2]+p2[2])/2;
face_node(face(*vcir,pmesh), q);
Q[0] += q[0];
Q[1] += q[1];
Q[2] += q[2];
}
R[0] /= n; R[1] /= n; R[2] /= n;
Q[0] /= n; Q[1] /= n; Q[2] /= n;
pt = Point_3((Q[0] + 2*R[0] + S[0]*(n-3))/n,
(Q[1] + 2*R[1] + S[1]*(n-3))/n,
(Q[2] + 2*R[2] + S[2]*(n-3))/n );
}
};

To invoke the Catmull-Clark subdivision method, we call PQQ() with the Catmull-Clark masks that we have just defined.

Loop, Doo-Sabin and $$\sqrt{3}$$ subdivisions are implemented using a similar process: pick a refinement host and implement the geometry policy. The key of developing your own subdivision method is implementing the right combination of the refinement host and the geometry policy. It is explained in the next two sections.

# Refinement Host

A refinement host is a template function of a polygon mesh class and a geometry mask class. It refines the input polygon mesh, and computes new points through the geometry masks. Subdivision_method_3 supports four refinement hosts: primal quadrilateral quadrisection (PQQ), primal triangle quadrisection (PTQ), dual quadrilateral quadrisection (DQQ) and $$\sqrt{3}$$ triangulation. Respectively, they are used by Catmull-Clark, Loop, Doo-Sabin and $$\sqrt{3}$$ subdivision.

namespace Subdivision_method_3 {
template <class PolygonMesh, class Mask, class NamedParameters>
template <class PolygonMesh, class Mask, class NamedParameters>
template <class PolygonMesh, class Mask, class NamedParameters>
template <class PolygonMesh, class Mask, class NamedParameters>
}

The mesh class must be a model of MutableFaceGraph and it must be a triangle mesh or a polygon mesh, and the mask is a policy class realizing the geometry masks of the subdivision method.

A refinement host refines the input polygon mesh, maintains the stencils (i.e., the mapping between the control mesh and the refined mesh), and calls the geometry masks to compute the new points. In Subdivision_method_3, refinements are implemented as a sequence of connectivity operations (mainly Euler operations). The order of the connectivity operations plays a key role when maintaining stencils. By matching the order of the source submeshes to the refined vertices, no flag in the primitives is required to register the stencils. It avoids the data dependency of the refinement host on the polygon mesh class. To make the ordering trick work, the polygon mesh class must have a sequential container, such as a vector or a linked-list, as the internal storage. A sequential container guarantees that the iterators of the polygon mesh always traverse the primitives in the order of their insertions. Non-sequential structures such as trees or maps do not provide the required ordering, and hence cannot be used with Subdivision_method_3.

Although Subdivision_method_3 does not require flags to support the refinements and the stencils, it still needs to know how to compute and store the geometry data (i.e. the points). The classes of Subdivision_method_3 have as optional template argument a vertex property map that provides a mapping between vertices and points.

Refinement hosts PQQ and DQQ work on a general polygon mesh, and PTQ and Sqrt3 work on a triangulated polygon mesh. The result of PTQ and Sqrt3 on a non-triangulated polygon mesh is undefined. Subdivision_method_3 does not verify the precondition of the mesh characteristics before the refinement.

For details of the refinement implementation, interested users should refer to [5].

# Geometry Policy

A geometry policy defines a set of geometry masks. Each geometry mask is realized as a member function that computes new points of the subdivision surface.

Each geometry mask receives a primitive descriptor (e.g. halfedge_descriptor) of the control mesh, and returns a Point_3 to the subdivided vertex. The function collects the vertex neighbors of the primitive descriptor (i.e. nodes on the stencil), and computes the new point based on the neighbors and the mask (i.e. the stencil weights).

This figure shows the geometry masks of Catmull-Clark subdivision. The weights shown here are unnormalized, and $$n$$ is the valence of the vertex. The new points are computed by the summation of the weighted points on their stencils. Following codes show an implementation of the geometry mask of the face-node. The complete listing of a Catmull-Clark geometry policy is in the Section Catmull-Clark Subdivision.

template <class PolygonMesh, class VertexPointMap>
typedef boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
typedef boost::property_traits<VertexPointMap>::value_type Point_3;
void face_node(face_descriptor face, Point_3& pt) {
int n = 0;
Point_3 p(0,0,0);
for(halfedge_descriptor hd : halfedges_around_face(face,pmesh)){
p = p + get(vpm, (target(hd,pmesh)) - ORIGIN);
++n;
}
pt = ORIGIN + (p - ORIGIN)/FT(n);
}
};

In this example, the computation is based on the assumption that the Point_3 is the CGAL::Point_3. It is an assumption, but not a restriction. You are allowed to use any point class as long as it is defined as the Point_3 in your polygon mesh. You may need to modify the geometry policy to support the computation and the assignment of the specialized point. This extension is not unusual in graphics applications. For example, you might want to subdivide the texture coordinates for your subdivision surface.

The refinement host of Catmull-Clark subdivision requires three geometry masks for polygon meshes without open boundaries: a vertex-node mask, an edge-node mask, and a face-node mask. To support polygon meshes with boundaries, a border-node mask is also required. The border-node mask for Catmull-Clark subdivision is listed below, where ept returns the new point splitting edge and vpt returns the new point on the vertex pointed by edge.

void border_node(halfedge_descriptor edge, Point_3& ept, Point_3& vpt) {
Point_3 ep1 = get(vpm, target(edge, pmesh));
Point_3 ep2 = get(vpm, target(opposite(edge, pmesh), pmesh));
ept = Point_3((ep1[0]+ep2[0])/2, (ep1[1]+ep2[1])/2, (ep1[2]+ep2[2])/2);
Halfedge_around_target_circulator<Poly> vcir(edge, pmesh);
Point_3 vp1 = get(vpm,target(opposite(*vcir, pmesh ), pmesh));
Point_3 vp0 = get(vpm, target(*vcir, pmesh));
--vcir;
Point_3 vp_1 = get(vpm, target(opposite(*vcir, pmesh), pmesh));
vpt = Point_3((vp_1[0] + 6*vp0[0] + vp1[0])/8,
(vp_1[1] + 6*vp0[1] + vp1[1])/8,
(vp_1[2] + 6*vp0[2] + vp1[2])/8 );
}

The mask interfaces of all four refinement hosts are listed below. DQQMask_3 does not have the border-node stencil because the refinement host of the DQQ refinement does not support global boundaries in the current release. This might change in future releases.

template <class PolygonMesh>
void face_node(boost::graph_traits<PolygonMesh>::face_descriptor, Point_3&);
void edge_node(boost::graph_traits<PolygonMesh>::halfedge_descriptor, Point_3&);
void vertex_node(boost::graph_traits<PolygonMesh>::vertex_descriptor, Point_3&);
void border_node(boost::graph_traits<PolygonMesh>::halfedge_descriptor, Point_3&, Point_3&);
};
template <class PolygonMesh>
void edge_node(boost::graph_traits<PolygonMesh>::halfedge_descriptor, Point_3&);
void vertex_node(boost::graph_traits<PolygonMesh>::vertex_descriptor, Point_3&);
void border_node(boost::graph_traits<PolygonMesh>::halfedge_descriptor, Point_3&, Point_3_&);
};
template <class PolygonMesh>
public:
void corner_node(boost::graph_traits<PolygonMesh>::halfedge_descriptor edge, Point_3& pt);
};
template <class PolygonMesh>
public:
void vertex_node(boost::graph_traits<PolygonMesh>::vertex_descriptor vertex, Point_3& pt);
};

The source codes of CatmullClark_mask_3, Loop_mask_3, DooSabin_mask_3, and Sqrt3_mask_3 are the best sources of learning these stencil interfaces.

# The Four Subdivision Methods

Subdivision_method_3 supports Catmull-Clark, Loop, Doo-Sabin and $$\sqrt{3}$$ subdivisions by specializing their respective refinement hosts. They are designed to work on any model of a MutableFaceGraph such as Polyhedron_3 and Surface_mesh. If your application uses a polygon mesh with a specialized geometry kernel, you need to specialize the refinement host with a geometry policy based on that kernel.

namespace Subdivision_method_3 {
template <class PolygonMesh, class NamedParameters>
void CatmullClark_subdivision(PolygonMesh& pmesh, NamedParameters np) {
}
template <class PolygonMesh, class NamedParameters>
void Loop_subdivision(PolygonMesh& pmesh, NamedParameters np) {
}
template <class PolygonMesh, class NamedParameters>
void DooSabin_subdivision(PolygonMesh& pmesh, NamedParameters np) {
}
template <class PolygonMesh, class NamedParameters>
void Sqrt3_subdivision(PolygonMesh& pmesh, NamedParameters np) {
}
}

# Example: Customizing a Subdivision Method

Subdivision_method_3 supports four practical subdivision methods on a Polyhedron_3 with points with Cartesian coordinates. More subdivision methods can be supported through the specialization of refinement hosts with custom geometry masks. The following example develops a subdivision method generating an improved Loop subdivision surface.

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/boost/graph/graph_traits_Surface_mesh.h>
#include <CGAL/Timer.h>
#include <boost/foreach.hpp>
#include <boost/lexical_cast.hpp>
#include <iostream>
#include <fstream>
using namespace std;
using namespace CGAL;
namespace params = CGAL::parameters;
// ======================================================================
template <class Poly>
typedef Poly PolygonMesh;
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::property_map<PolygonMesh, vertex_point_t>::type Vertex_pmap;
typedef typename boost::property_traits<Vertex_pmap>::value_type Point;
typedef typename boost::property_traits<Vertex_pmap>::reference Point_ref;
PolygonMesh& pmesh;
Vertex_pmap vpm;
public:
: pmesh(pmesh), vpm(get(CGAL::vertex_point, pmesh))
{}
void edge_node(halfedge_descriptor hd, Point& pt) {
Point_ref p1 = get(vpm, target(hd,pmesh));
Point_ref p2 = get(vpm, target(opposite(hd,pmesh),pmesh));
Point_ref f1 = get(vpm, target(next(hd,pmesh),pmesh));
Point_ref f2 = get(vpm, target(next(opposite(hd,pmesh),pmesh),pmesh));
pt = Point((3*(p1[0]+p2[0])+f1[0]+f2[0])/8,
(3*(p1[1]+p2[1])+f1[1]+f2[1])/8,
(3*(p1[2]+p2[2])+f1[2]+f2[2])/8 );
}
void vertex_node(vertex_descriptor vd, Point& pt) {
double R[] = {0.0, 0.0, 0.0};
Point_ref S = get(vpm,vd);
std::size_t n = 0;
BOOST_FOREACH(halfedge_descriptor hd, halfedges_around_target(vd, pmesh)){
++n;
Point_ref p = get(vpm, target(opposite(hd,pmesh),pmesh));
R[0] += p[0]; R[1] += p[1]; R[2] += p[2];
}
if (n == 6) {
pt = Point((10*S[0]+R[0])/16, (10*S[1]+R[1])/16, (10*S[2]+R[2])/16);
} else if (n == 3) {
double B = (5.0/8.0 - std::sqrt(3+2*std::cos(6.283/n))/64.0)/n;
double A = 1-n*B;
pt = Point((A*S[0]+B*R[0]), (A*S[1]+B*R[1]), (A*S[2]+B*R[2]));
} else {
double B = 3.0/8.0/n;
double A = 1-n*B;
pt = Point((A*S[0]+B*R[0]), (A*S[1]+B*R[1]), (A*S[2]+B*R[2]));
}
}
void border_node(halfedge_descriptor hd, Point& ept, Point& vpt) {
Point_ref ep1 = get(vpm, target(hd,pmesh));
Point_ref ep2 = get(vpm, target(opposite(hd,pmesh),pmesh));
ept = Point((ep1[0]+ep2[0])/2, (ep1[1]+ep2[1])/2, (ep1[2]+ep2[2])/2);
Point_ref vp1 = get(vpm, target(opposite(*vcir,pmesh),pmesh));
Point_ref vp0 = get(vpm, target(*vcir,pmesh));
--vcir;
Point_ref vp_1 = get(vpm,target(opposite(*vcir,pmesh),pmesh));
vpt = Point((vp_1[0] + 6*vp0[0] + vp1[0])/8,
(vp_1[1] + 6*vp0[1] + vp1[1])/8,
(vp_1[2] + 6*vp0[2] + vp1[2])/8 );
}
};
int main(int argc, char **argv) {
if (argc > 4) {
cerr << "Usage: Customized_subdivision [d] [filename_in] [filename_out] \n";
cerr << " d -- the depth of the subdivision (default: 1) \n";
cerr << " filename_in -- the input mesh (.off) (default: data/quint_tris.off) \n";
cerr << " filename_out -- the output mesh (.off) (default: result.off)" << endl;
return 1;
}
int d = (argc > 1) ? boost::lexical_cast<int>(argv[1]) : 1;
const char* in_file = (argc > 2) ? argv[2] : "data/quint_tris.off";
const char* out_file = (argc > 3) ? argv[3] : "result.off";
PolygonMesh pmesh;
std::ifstream in(in_file);
if(in.fail()) {
std::cerr << "Could not open input file " << in_file << std::endl;
return 1;
}
in >> pmesh;
Timer t;
t.start();
The points generated by the geometry mask are semantically required to converge to a smooth surface. This is the requirement imposed by the theory of the subdivision surface. Subdivision_method_3 does not enforce this requirement, nor will it verify the smoothness of the subdivided mesh. Subdivision_method_3 guarantees the topological properties of the subdivided mesh. A genus- $$n$$ 2-manifold is assured to be subdivided into a genus- $$n$$ 2-manifold. But when specialized with ill-designed geometry masks, Subdivision_method_3 may generate a surface that is odd, not smooth, or that does not even exist.
This package was initially developed by Le-Jeng Andy Shiue. For CGAL 4.11 it was generalized to work on any model of a MutableFaceGraph by Andreas Fabri and Mael Rouxel-Labbé.