CGAL 6.1.1 - 3D Surface Subdivision Methods

Loading...
Searching...
No Matches
User Manual

Author
  Le-Jeng Andy Shiue

# 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. In particular, the edge case where one wishes to subdivide a mesh while keeping its geometry fixed is offered, setting an optional parameter in Catmull-Clark subdivision and Loop subdivision.

# 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.

File Subdivision_method_3/CatmullClark_subdivision.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/boost/graph/graph_traits_Surface_mesh.h>
#include <CGAL/subdivision_method_3.h>
#include <CGAL/Timer.h>
#include <boost/lexical_cast.hpp>
#include <iostream>
#include <fstream>
typedef CGAL::Simple_cartesian<double> Kernel;
typedef CGAL::Surface_mesh<Kernel::Point_3> PolygonMesh;
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 std::string in_file = (argc > 2) ? argv[2] : CGAL::data_file_path("meshes/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;
}
CGAL::Surface_mesh
CGAL::Subdivision_method_3::CatmullClark_subdivision
void CatmullClark_subdivision(PolygonMesh &pmesh, const NamedParameters &np=parameters::default_values())
applies Catmull-Clark subdivision several times on the control mesh pmesh.
Definition: subdivision_methods_3.h:132
CGAL
CGAL::parameters
Kernel
std
STL namespace.
CGAL::Simple_cartesian
subdivision_method_3.h
Convenience header file including the headers for all the free functions of this package.

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(pmesh, 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.

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>
void PQQ(PolygonMesh& p, Mask mask, NamedParameters np)

`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>
class PQQMask_3 {
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);
};
CGAL::Point_3
PQQMask_3
Required member functions for the PQQMask_3 concept. This policy concept of geometric computations is...
Definition: PQQMask_3.h:22

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>
class CatmullClark_mask_3 {
typedef boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
typedef boost::property_traits<VertexPointMap>::value_type Point_3;
Polygonmesh pmesh;
VertexPointMap vpm;
CatmullClark_mask_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 );
}
};
CGAL::CatmullClark_mask_3
The geometry mask of Catmull-Clark subdivision.
Definition: subdivision_masks_3.h:153
CGAL::Halfedge_around_target_circulator
CGAL::halfedges_around_face
Iterator_range< Halfedge_around_face_iterator< Graph > > halfedges_around_face(typename boost::graph_traits< Graph >::halfedge_descriptor h, const Graph &g)
CGAL::ORIGIN
const CGAL::Origin ORIGIN
CGAL::opposite
Oriented_side opposite(const Oriented_side &o)

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

PQQ(pmesh, CatmullClark_mask_3(pmesh), params::number_of_iterations(depth));

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>
void PQQ(PolygonMesh& pmesh, Mask mask, NamedParameters np);
template <class PolygonMesh, class Mask, class NamedParameters>
void PTQ(PolygonMesh& pmesh, Mask mask, NamedParameters np);
template <class PolygonMesh, class Mask, class NamedParameters>
void DQQ(PolygonMesh& pmesh, Mask mask, NamedParameters np)
template <class PolygonMesh, class Mask, class NamedParameters>
void Sqrt3(PolygonMesh& pmesh, Mask mask, NamedParameters np)
}
CGAL::Subdivision_method_3::PTQ
void PTQ(PolygonMesh &pmesh, Mask mask, const NamedParameters &np)
applies the PTQ refinement several times on the control mesh pmesh.
Definition: subdivision_hosts_3.h:138
CGAL::Subdivision_method_3::PQQ
void PQQ(PolygonMesh &pmesh, Mask mask, const NamedParameters &np)
applies the PQQ refinement several times on the control mesh pmesh.
Definition: subdivision_hosts_3.h:82
CGAL::Subdivision_method_3::DQQ
void DQQ(PolygonMesh &pmesh, Mask mask, const NamedParameters &np)
applies the DQQ refinement several times on the control mesh pmesh.
Definition: subdivision_hosts_3.h:197
CGAL::Subdivision_method_3::Sqrt3
void Sqrt3(PolygonMesh &pmesh, Mask mask, const NamedParameters &np)
applies the refinement several times on the control mesh pmesh.
Definition: subdivision_hosts_3.h:258

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.

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>
class CatmullClark_mask_3 {
typedef boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
typedef boost::property_traits<VertexPointMap>::value_type Point_3;
CatmullClark_mask_3(PolygonMesh &pmesh, VertexPointMap 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);
}
};

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>
class PQQMask_3 {
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>
class PTQMask_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>
class DQQMask_3 {
public:
void corner_node(boost::graph_traits<PolygonMesh>::halfedge_descriptor edge, Point_3& pt);
};
template <class PolygonMesh>
class Sqrt3Mask_3 {
public:
void vertex_node(boost::graph_traits<PolygonMesh>::vertex_descriptor vertex, Point_3& pt);
};
DQQMask_3
Required member functions for the DQQMask_3 concept. This policy concept of geometric computations is...
Definition: DQQMask_3.h:21
PTQMask_3
Required member functions for the PTQMask_3 concept. This policy concept of geometric computations is...
Definition: PTQMask_3.h:19
Sqrt3Mask_3
Required member functions for the Sqrt3Mask_3 concept. This policy concept of geometric computations ...
Definition: Sqrt3Mask_3.h:20

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) {
PQQ(pmesh, CatmullClark_mask_3<PolygonMesh>(pmesh), np);
}
template <class PolygonMesh, class NamedParameters>
void Loop_subdivision(PolygonMesh& pmesh, NamedParameters np) {
PTQ(pmesh, Loop_mask_3<PolygonMesh>(pmesh) , np);
}
template <class PolygonMesh, class NamedParameters>
void DooSabin_subdivision(PolygonMesh& pmesh, NamedParameters np) {
DQQ(pmesh, DooSabin_mask_3<PolygonMesh>(pmesh), np);
}
template <class PolygonMesh, class NamedParameters>
void Sqrt3_subdivision(PolygonMesh& pmesh, NamedParameters np) {
Sqrt3(pmesh, Sqrt3_mask_3<PolygonMesh>(pmesh), np);
}
}
CGAL::DooSabin_mask_3
The geometry mask of Doo-Sabin subdivision.
Definition: subdivision_masks_3.h:440
CGAL::Loop_mask_3
The geometry mask of Loop subdivision.
Definition: subdivision_masks_3.h:278
CGAL::Sqrt3_mask_3
The geometry mask of Sqrt(3) subdivision.
Definition: subdivision_masks_3.h:534
CGAL::Subdivision_method_3::Sqrt3_subdivision
void Sqrt3_subdivision(PolygonMesh &pmesh, const NamedParameters &np=parameters::default_values())
applies -subdivision several times on the control mesh pmesh.
Definition: subdivision_methods_3.h:332
CGAL::Subdivision_method_3::DooSabin_subdivision
void DooSabin_subdivision(PolygonMesh &pmesh, const NamedParameters &np=parameters::default_values())
applies DooSabin subdivision several times on the control mesh pmesh.
Definition: subdivision_methods_3.h:270
CGAL::Subdivision_method_3::Loop_subdivision
void Loop_subdivision(PolygonMesh &pmesh, const NamedParameters &np=parameters::default_values())
applies Loop subdivision several times on the control mesh pmesh.
Definition: subdivision_methods_3.h:205

# Example: Customizing a Subdivision Method

`Subdivision_method_3` supports four practical subdivision methods on a polygon mesh 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.

File Subdivision_method_3/Customized_subdivision.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/boost/graph/graph_traits_Surface_mesh.h>
#include <CGAL/subdivision_method_3.h>
#include <CGAL/Timer.h>
#include <boost/lexical_cast.hpp>
#include <iostream>
#include <fstream>
typedef CGAL::Simple_cartesian<double> Kernel;
typedef CGAL::Surface_mesh<Kernel::Point_3> PolygonMesh;
using namespace std;
using namespace CGAL;
namespace params = CGAL::parameters;
// ======================================================================
template <class Poly>
class WLoop_mask_3 {
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:
WLoop_mask_3(PolygonMesh& pmesh)
: 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;
for(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);
Halfedge_around_target_circulator<Poly> vcir(hd,pmesh);
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 std::string in_file = (argc > 2) ? argv[2] : CGAL::data_file_path("meshes/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::PTQ(pmesh, WLoop_mask_3<PolygonMesh>(pmesh), params::number_of_iterations(d));
std::cerr << "Done (" << t.time() << " s)" << std::endl;
std::ofstream out(out_file);
out << pmesh;
return 0;
}
CGAL::halfedges_around_target
Iterator_range< Halfedge_around_target_iterator< Graph > > halfedges_around_target(typename boost::graph_traits< Graph >::halfedge_descriptor h, const Graph &g)

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.

# Implementation History

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é.