3D Surface Subdivision Methods

Subdivision methods are simple yet powerful ways to generate smooth surfaces from arbitrary polyhedral 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*, designed to work on the class
*Polyhedron_3*, 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.

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[CC78], Loop, Doo-Sabin
and $$*sqrt(3)* subdivision (left to right in the
figure). We name these patterns by their topological
characteristics instead of the associated subdivision methods.
PQQ indicates the *P*rimal *Q*uadtrateral *Q*uadrisection.
PTQ indicates the *P*rimal *T*riangle *Q*uadrisection.
DQQ indicates the *D*ual *Q*uadtrateral *Q*uadrisection.
$$*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/facet.
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 facet-node stencil,
which defines an input facet. 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 facet-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 facet-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. [WW02] explains the
details on designing masks for a quality subdivision surface.

File:examples/Subdivision_method_3/CatmullClark_subdivision.cpp

#include <CGAL/Cartesian.h> #include <CGAL/Subdivision_method_3.h> #include <iostream> #include <CGAL/Polyhedron_3.h> #include <CGAL/IO/Polyhedron_iostream.h> typedef CGAL::Cartesian<double> Kernel; typedef CGAL::Polyhedron_3<Kernel> Polyhedron; using namespace std; using namespace CGAL; int main(int argc, char** argv) { if (argc != 2) { cout << "Usage: CatmullClark_subdivision d < filename" << endl; cout << " d: the depth of the subdivision (0 < d < 10)" << endl; cout << " filename: the input mesh (.off)" << endl; return 0; } int d = argv[1][0] - '0'; Polyhedron P; cin >> P; // read the .off Subdivision_method_3::CatmullClark_subdivision(P,d); cout << P; // write the .off return 0; }

This example demonstrates the use of the Catmull-Clark subdivision method
on a *Polyhedron_3*. The polyhedron is restricted in the Cartesian
space, where most subdivision applications are designed to work.
There is only one line deserving a detailed explanation:

Subdivision_method_3::CatmullClark_subdivision(P,d);

This example shows how to subdivide a simple *Polyhedron_3*
with *Subdivision_method_3*.
An application-defined polyhedron might use a specialized kernel and/or
a specialized internal container. There are two major restrictions on the
application-defined polyhedron to work with
*Subdivision_method_3*.

*Point_3*is type-defined by the kernel. Without*Point_3*and the associated operations being defined,*Subdivision_method_3*can not know how to compute and store the new vertex points.- The primitives (such as vertices, halfedges and facets)
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 37.5 gives detailed explanations on those two restrictions.

When a subdivision method is developed, a refinement pattern is chosen, and then a set of geometry masks are 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 picked 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 Polyhedron_3, template <typename> class Mask> void PQQ(Polyhedron_3& p, Mask<Polyhedron_3> mask, int depth)

*Polyhedron_3* 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<Polyhedron_3>* 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.
*depth* specifies the iterations of the refinement
on the control mesh.

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 Polyhedron_3> class PQQ_stencil_3 { void facet_node(Facet_handle facet, Point_3& pt); void edge_node(Halfedge_handle edge, Point_3& pt); void vertex_node(Vertex_handle vertex, Point_3& pt); };

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

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

template <class Polyhedron_3> class CatmullClark_mask_3 { void facet_node(Facet_handle facet, Point_3& pt) { Halfedge_around_facet_circulator hcir = facet->facet_begin(); int n = 0; Point_3 p(0,0,0); do { p = p + (hcir->vertex()->point() - ORIGIN); ++n; } while (++hcir != facet->facet_begin()); pt = ORIGIN + (p - ORIGIN)/FT(n); } void edge_node(Halfedge_handle edge, Point_3& pt) { Point_3 p1 = edge->vertex()->point(); Point_3 p2 = edge->opposite()->vertex()->point(); Point_3 f1, f2; facet_node(edge->facet(), f1); facet_node(edge->opposite()->facet(), 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_handle vertex, Point_3& pt) { Halfedge_around_vertex_circulator vcir = vertex->vertex_begin(); int n = circulator_size(vcir); FT Q[] = {0.0, 0.0, 0.0}, R[] = {0.0, 0.0, 0.0}; Point_3& S = vertex->point(); Point_3 q; for (int i = 0; i < n; i++, ++vcir) { Point_3& p2 = vcir->opposite()->vertex()->point(); R[0] += (S[0]+p2[0])/2; R[1] += (S[1]+p2[1])/2; R[2] += (S[2]+p2[2])/2; facet_node(vcir->facet(), 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 ); } };

This example shows the default implementation of Catmull-Clark
masks in *Subdivision_method_3*.
This default implementation assumes the *types*
(such as *Point_3* and *Facet_handle*) are defined
within *Polyhedron_3*. *CatmullClark_mask_3*
is designed to work on a *Polyhedron_3* with the *Cartesian*
kernel. You may need to rewrite the geometry computation
to match the kernel geometry of your application.

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

PQQ(p, CatmullClark_mask_3<Polyhedron_3>(), depth);

Loop, Doo-Sabin and $$*sqrt(3)* subdivisions are implemented
in the 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.

namespace Subdivision_method_3 { template <class Polyhedron_3, template <typename> class Mask> void PQQ(Polyhedron_3& p, Mask<Polyhedron_3> mask, int step); template <class Polyhedron_3, template <typename> class Mask> void PTQ(Polyhedron_3& p, Mask<Polyhedron_3> mask, int step); template <class Polyhedron_3, template <typename> class Mask> void DQQ(Polyhedron_3& p, Mask<Polyhedron_3> mask, int step) template <class Polyhedron_3, template <typename> class Mask> void Sqrt3(Polyhedron_3& p, Mask<Polyhedron_3> mask, int step) }

The polyhedron class is a specialization of
*Polyhedron_3*, and the mask is a policy
class realizing the geometry masks of the subdivision
method.

A refinement host refines the input polyhedron, 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 polyhedron class.
To make the ordering trick work, the polyhedron 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
polyhedron 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
can not 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). *Subdivision_method_3*
expects that the typename *Point_3* is
defined in the geometry kernel of the polyhedron
(i.e. the *Polyhedron_3::Traits::Kernel*).
A point of the type *Point_3* is returned by the geometry
policy and is then assigned to the new vertex.
The geometry policy is explained in next section.

Refinement hosts *PQQ* and *DQQ* work on a general
polyhedron, and *PTQ* and *Sqrt3* work on a triangulated
polyhedron. The result of *PTQ* and *Sqrt3* on a non-triangulated
polyhedron 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 [SP05].

Each geometry mask receives a primitive handle
(e.g. *Halfedge_handle*) of the control mesh,
and returns a *Point_3* to the subdivided vertex.
The function collects the vertex neighbors of the primitive handle
(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 facet-node. The complete listing
of a Catmull-Clark geometry policy is in the Section 37.4.

template <class Polyhedron_3> class CatmullClark_mask_3 { void facet_node(Facet_handle facet, Point_3& pt) { Halfedge_around_facet_circulator hcir = facet->facet_begin(); int n = 0; Point_3 p(0,0,0); do { p = p + (hcir->vertex()->point() - ORIGIN); ++n; } while (++hcir != facet->facet_begin()); 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 polyhedron.
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
polyhedrons without open boundaries: a vertex-node
mask, an edge-node mask, and a facet-node mask.
To support polyhedrons 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_handle edge, Point_3& ept, Point_3& vpt) { Point_3& ep1 = edge->vertex()->point(); Point_3& ep2 = edge->opposite()->vertex()->point(); ept = Point_3((ep1[0]+ep2[0])/2, (ep1[1]+ep2[1])/2, (ep1[2]+ep2[2])/2); Halfedge_around_vertex_circulator vcir = edge->vertex_begin(); Point_3& vp1 = vcir->opposite()->vertex()->point(); Point_3& vp0 = vcir->vertex()->point(); Point_3& vp_1 = (--vcir)->opposite()->vertex()->point(); 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.
*DQQ_stencil_3* and *Sqrt3_stencil_3*
do not have the border-node stencil because the refinement hosts of
DQQ and $$*sqrt(3)* refinements do not support global boundaries in the
current release. This might be changed in the future releases.

template <class Polyhedron_3> class PQQ_stencil_3 { void facet_node(Facet_handle, Point_3&); void edge_node(Halfedge_handle, Point_3&); void vertex_node(Vertex_handle, Point_3&); void border_node(Halfedge_handle, Point_3&, Point_3&); }; template <class Polyhedron_3> class PTQ_stencil_3 { void edge_node(Halfedge_handle, Point_3&); void vertex_node(Vertex_handle, Point_3&); void border_node(Halfedge_handle, Point_3&, Point_&); }; template <class Polyhedron_3> class DQQ_stencil_3 { public: void corner_node(Halfedge_handle edge, Point_3& pt); }; template <class Polyhedron_3> class Sqrt3_stencil_3 { public: void vertex_node(Vertex_handle 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.

namespace Subdivision_method_3 { template <class Polyhedron_3> void CatmullClark_subdivision(Polyhedron_3& p, int step = 1) { PQQ(p, CatmullClark_mask_3<Polyhedron_3>(), step); } template <class Polyhedron_3> void Loop_subdivision(Polyhedron_3& p, int step = 1) { PTQ(p, Loop_mask_3<Polyhedron_3>() , step); } template <class Polyhedron_3> void DooSabin_subdivision(Polyhedron_3& p, int step = 1) { DQQ(p, DooSabin_mask_3<Polyhedron_3>(), step); } template <class Polyhedron_3> void Sqrt3_subdivision(Polyhedron_3& p, int step = 1) { Sqrt3(p, Sqrt3_mask_3<Polyhedron_3>(), step); } }

The following example demonstrates the use of Doo-Sabin subdivision on a polyhedral mesh.

File:examples/Subdivision_method_3/DooSabin_subdivision.cpp

#include <CGAL/Cartesian.h> #include <CGAL/Subdivision_method_3.h> #include <iostream> #include <CGAL/Polyhedron_3.h> #include <CGAL/IO/Polyhedron_iostream.h> typedef CGAL::Cartesian<double> Kernel; typedef CGAL::Polyhedron_3<Kernel> Polyhedron; using namespace std; using namespace CGAL; int main(int argc, char **argv) { if (argc != 2) { cout << "Usage: DooSabin_subdivision d < filename" << endl; cout << " d: the depth of the subdivision (0 < d < 10)" << endl; cout << " filename: the input mesh (.off)" << endl; return 0; } int d = argv[1][0] - '0'; Polyhedron P; cin >> P; // read the .off Subdivision_method_3::DooSabin_subdivision(P,d); cout << P; // write the .off return 0; }

File:examples/Subdivision_method_3/Customized_subdivision.cpp

#include <CGAL/Cartesian.h> #include <CGAL/Subdivision_method_3.h> #include <cstdio> #include <iostream> #include <CGAL/Polyhedron_3.h> #include <CGAL/IO/Polyhedron_iostream.h> typedef CGAL::Cartesian<double> Kernel; typedef CGAL::Polyhedron_3<Kernel> Polyhedron; using namespace std; using namespace CGAL; // ====================================================================== template <class Poly> class WLoop_mask_3 { typedef Poly Polyhedron; typedef typename Polyhedron::Vertex_iterator Vertex_iterator; typedef typename Polyhedron::Halfedge_iterator Halfedge_iterator; typedef typename Polyhedron::Facet_iterator Facet_iterator; typedef typename Polyhedron::Halfedge_around_facet_circulator Halfedge_around_facet_circulator; typedef typename Polyhedron::Halfedge_around_vertex_circulator Halfedge_around_vertex_circulator; typedef typename Polyhedron::Traits Traits; typedef typename Traits::Kernel Kernel; typedef typename Kernel::FT FT; typedef typename Kernel::Point_3 Point; typedef typename Kernel::Vector_3 Vector; public: void edge_node(Halfedge_iterator eitr, Point& pt) { Point& p1 = eitr->vertex()->point(); Point& p2 = eitr->opposite()->vertex()->point(); Point& f1 = eitr->next()->vertex()->point(); Point& f2 = eitr->opposite()->next()->vertex()->point(); 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_iterator vitr, Point& pt) { float R[] = {0.0, 0.0, 0.0}; Point& S = vitr->point(); Halfedge_around_vertex_circulator vcir = vitr->vertex_begin(); int n = circulator_size(vcir); for (int i = 0; i < n; i++, ++vcir) { Point& p = vcir->opposite()->vertex()->point(); 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_iterator eitr, Point& ept, Point& vpt) { Point& ep1 = eitr->vertex()->point(); Point& ep2 = eitr->opposite()->vertex()->point(); ept = Point((ep1[0]+ep2[0])/2, (ep1[1]+ep2[1])/2, (ep1[2]+ep2[2])/2); Halfedge_around_vertex_circulator vcir = eitr->vertex_begin(); Point& vp1 = vcir->opposite()->vertex()->point(); Point& vp0 = vcir->vertex()->point(); Point& vp_1 = (--vcir)->opposite()->vertex()->point(); 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 != 2) { cout << "Usage: Customized_subdivision d < filename" << endl; cout << " d: the depth of the subdivision (0 < d < 10)" << endl; cout << " filename: the input mesh (.off)" << endl; return 0; } int d = argv[1][0] - '0'; Polyhedron P; cin >> P; // read the .off Subdivision_method_3::PTQ(P, WLoop_mask_3<Polyhedron>(), d); cout << P; // write the .off return 0; }

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 not even exist.