Chapter 12
3D Polyhedral Surfaces

Lutz Kettner

Table of Contents

12.1 Introduction
12.2 Definition
12.3 Example Programs
   12.3.1   First Example Using Defaults
   12.3.2   Example with Geometry in Vertices
   12.3.3   Example for Affine Transformation
   12.3.4   Example Computing Plane Equations
   12.3.5   Example with a Vector Instead of a List Representation
   12.3.6   Example with Circulator Writing Object File Format (OFF)
   12.3.7   Example Using Euler Operators to Build a Cube
12.4 File I/O
12.5 Extending Vertices, Halfedges, and Facets
12.6 Advanced Example Programs
   12.6.1   Example Creating a Subdivision Surface
   12.6.2   Example Using the Incremental Builder and Modifier Mechanism

12.1   Introduction

Polyhedral surfaces in three dimensions are composed of vertices, edges, facets and an incidence relationship on them. The organization beneath is a halfedge data structure, which restricts the class of representable surfaces to orientable 2-manifolds - with and without boundary. If the surface is closed we call it a polyhedron, for example, see the following model of a hammerhead:

Hammerhead

The polyhedral surface is realized as a container class that manages vertices, halfedges, facets with their incidences, and that maintains the combinatorial integrity of them. It is based on the highly flexible design of the halfedge data structure, see the introduction in Chapter 13 and [Ket99]. However, the polyhedral surface can be used and understood without knowing the underlying design. Some of the examples in this chapter introduce also gradually into first applications of this flexibility.

12.2   Definition

A polyhedral surface CGAL::Polyhedron_3<PolyhedronTraits_3> in three dimensions consists of vertices V, edges E, facets F and an incidence relation on them. Each edge is represented by two halfedges with opposite orientations. The incidences stored with a halfedge are illustrated in the following figure:

Halfedge Diagram

Vertices represent points in space. Edges are straight line segments between two endpoints. Facets are planar polygons without holes. Facets are defined by the circular sequence of halfedges along their boundary. The polyhedral surface itself can have holes (with at least two facets surrounding it since a single facet cannot have a hole). The halfedges along the boundary of a hole are called border halfedges and have no incident facet. An edge is a border edge if one of its halfedges is a border halfedge. A surface is closed if it contains no border halfedges. A closed surface is a boundary representation for polyhedra in three dimensions. The convention is that the halfedges are oriented counterclockwise around facets as seen from the outside of the polyhedron. An implication is that the halfedges are oriented clockwise around the vertices. The notion of the solid side of a facet as defined by the halfedge orientation extends to polyhedral surfaces with border edges although they do not define a closed object. If normal vectors are considered for the facets, normals point outwards (following the right-hand rule).

The strict definition can be found in [Ket99]. One implication of this definition is that the polyhedral surface is always an orientable and oriented 2-manifold with border edges, i.e., the neighborhood of each point on the polyhedral surface is either homeomorphic to a disc or to a half disc, except for vertices where many holes and surfaces with boundary can join. Another implication is that the smallest representable surface avoiding self intersections is a triangle (for polyhedral surfaces with border edges) or a tetrahedron (for polyhedra). Boundary representations of orientable 2-manifolds are closed under Euler operations. They are extended with operations that create or close holes in the surface.

Other intersections besides the incidence relation are not allowed. However, this is not automatically verified in the operations, since self intersections are not easy to check efficiently. CGAL::Polyhedron_3<PolyhedronTraits_3> does only maintain the combinatorial integrity of the polyhedral surface (using Euler operations) and does not consider the coordinates of the points or any other geometric information.

CGAL::Polyhedron_3<PolyhedronTraits_3> can represent polyhedral surfaces as well as polyhedra. The interface is designed in such a way that it is easy to ignore border edges and work only with polyhedra.

12.3   Example Programs

The polyhedral surface is based on the highly flexible design of the halfedge data structure. Examples for this flexibility can be found in Section 12.5 and in Section 13.3. This design is not a prerequisite to understand the following examples. See also the Section 12.6 below for some advanced examples.

12.3.1   First Example Using Defaults

The first example instantiates a polyhedron using a kernel as traits class. It creates a tetrahedron and stores the reference to one of its halfedges in a Halfedge_handle. Handles, also know as trivial iterators, are used to keep references to halfedges, vertices, or facets for future use. There is also a Halfedge_iterator type for enumerating halfedges. Such an iterator type can be used wherever a handle is required. Respective Halfedge_const_handle and Halfedge_const_iterator for a constant polyhedron and similar handles and iterators with Vertex_ and Facet_ prefix are provided too.

The example continues with a test if the halfedge actually refers to a tetrahedron. This test checks the connected component referred to by the halfedge h and not the polyhedral surface as a whole. This examples works only on the combinatorial level of a polyhedral surface. The next example adds the geometry.

File: examples/Polyhedron/polyhedron_prog_simple.cpp
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_3.h>

typedef CGAL::Simple_cartesian<double>     Kernel;
typedef CGAL::Polyhedron_3<Kernel>         Polyhedron;
typedef Polyhedron::Halfedge_handle        Halfedge_handle;

int main() {
    Polyhedron P;
    Halfedge_handle h = P.make_tetrahedron();
    if ( P.is_tetrahedron(h))
        return 0;
    return 1;
}

12.3.2   Example with Geometry in Vertices

We add geometry to the our construction of a tetrahedron. Four points are passed as arguments to the construction. This example demonstrates in addition the use of the vertex iterator and the access to the point in the vertices. Note the natural access notation v->point(). Similarly, all information stored in a vertex, halfedge, and facet can be accessed with a member function given a handle or iterator. For example, given a halfedge handle h we can write h->next() to get a halfedge handle to the next halfedge, h->opposite() for the opposite halfedge, h->vertex() for the incident vertex at the tip of h, and so on. The output of the program will be ``1 0 0\n0 1 0\n0 0 1\n0 0 0\n''.

File: examples/Polyhedron/polyhedron_prog_tetra.cpp
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_3.h>
#include <iostream>

typedef CGAL::Simple_cartesian<double>     Kernel;
typedef Kernel::Point_3                    Point_3;
typedef CGAL::Polyhedron_3<Kernel>         Polyhedron;
typedef Polyhedron::Vertex_iterator        Vertex_iterator;

int main() {
    Point_3 p( 1.0, 0.0, 0.0);
    Point_3 q( 0.0, 1.0, 0.0);
    Point_3 r( 0.0, 0.0, 1.0);
    Point_3 s( 0.0, 0.0, 0.0);

    Polyhedron P;
    P.make_tetrahedron( p, q, r, s);
    CGAL::set_ascii_mode( std::cout);
    for ( Vertex_iterator v = P.vertices_begin(); v != P.vertices_end(); ++v)
        std::cout << v->point() << std::endl;
    return 0;
}

The polyhedron offers a point iterator for convenience. The above for loop simplifies to a single statement by using std::copy and the ostream iterator adaptor.

std::copy( P.points_begin(), P.points_end(), 
           std::ostream_iterator<Point_3>(std::cout,"\n"));

12.3.3   Example for Affine Transformation

An affine transformation A can act as a functor transforming points and a point iterator is conveniently defined for polyhedral surfaces. So, assuming we want only the point coordinates of a polyhedron P transformed, std::transform does the job in a single line.

std::transform( P.points_begin(), P.points_end(), P.points_begin(), A);

12.3.4   Example Computing Plane Equations

The polyhedral surface has already provisions to store a plane equation for each facet. However, it does not provide a function to compute plane equations.

This example computes the plane equations of a polyhedral surface. The actual computation is implemented in the compute_plane_equations function. Depending on the arithmetic (exact/inexact) and the shape of the facets (convex/non-convex) different methods are useful. We assume here strictly convex facets and exact arithmetic. In our example a homogeneous representation with int coordinates is sufficient. The four plane equations of the tetrahedron are the output of the program.

File: examples/Polyhedron/polyhedron_prog_planes.cpp
#include <CGAL/Homogeneous.h>
#include <CGAL/Polyhedron_3.h>
#include <iostream>
#include <algorithm>

struct Plane_equation {
    template <class Facet>
    typename Facet::Plane_3 operator()( Facet& f) {
        typename Facet::Halfedge_handle h = f.halfedge();
        typedef typename Facet::Plane_3  Plane;
        return Plane( h->vertex()->point(),
                      h->next()->vertex()->point(),
                      h->next()->next()->vertex()->point());
    }
};

typedef CGAL::Homogeneous<int>      Kernel;
typedef Kernel::Point_3             Point_3;
typedef Kernel::Plane_3             Plane_3;
typedef CGAL::Polyhedron_3<Kernel>  Polyhedron;

int main() {
    Point_3 p( 1, 0, 0);
    Point_3 q( 0, 1, 0);
    Point_3 r( 0, 0, 1);
    Point_3 s( 0, 0, 0);
    Polyhedron P;
    P.make_tetrahedron( p, q, r, s);
    std::transform( P.facets_begin(), P.facets_end(), P.planes_begin(),
                    Plane_equation());
    CGAL::set_pretty_mode( std::cout);
    std::copy( P.planes_begin(), P.planes_end(),
               std::ostream_iterator<Plane_3>( std::cout, "\n"));
    return 0;
}

12.3.5   Example with a Vector Instead of a List Representation

The polyhedron class template has actually four parameters, where three of them have default values. Using the default values explicitly in our examples above for three parameter - ignoring the fourth parameter, which would be a standard allocator for container class - the definition of a polyhedron looks like:

typedef CGAL::Polyhedron_3< Traits, 
                            CGAL::Polyhedron_items_3, 
                            CGAL::HalfedgeDS_default>      Polyhedron;

The CGAL::Polyhedron_items_3 class contains the types used for vertices, edges, and facets. The CGAL::HalfedgeDS_default class defines the halfedge data structure used, which is a list-based representation in this case. An alternative is a vector-based representation. Using a vector provides random access for the elements in the polyhedral surface and is more space efficient, but elements cannot be deleted arbitrarily. Using a list allows arbitrary deletions, but provides only bidirectional iterators and is less space efficient. The following example creates again a tetrahedron with given points, but in a vector-based representation.

The vector-based representation resizes automatically if the reserved capacity is not sufficient for the new items created. Upon resizing all handles, iterators, and circulators become invalid. Their correct update in the halfedge data structure is costly, thus it is advisable to reserve enough space in advance as indicated with the alternative constructor in the comment.


begin of advanced section  advanced  begin of advanced section
Note that the polyhedron and not the underlying halfedge data structure triggers the resize operation, since the resize operation requires some preconditions, such as valid incidences, to be fulfilled that only the polyhedron can guarantee.
end of advanced section  advanced  end of advanced section

File: examples/Polyhedron/polyhedron_prog_vector.cpp
#include <CGAL/Cartesian.h>
#include <CGAL/HalfedgeDS_vector.h>
#include <CGAL/Polyhedron_3.h>
#include <iostream>

typedef CGAL::Cartesian<double>                        Kernel;
typedef Kernel::Point_3                                Point_3;
typedef CGAL::Polyhedron_3< Kernel,
                            CGAL::Polyhedron_items_3,
                            CGAL::HalfedgeDS_vector>   Polyhedron;

int main() {
    Point_3 p( 1.0, 0.0, 0.0);
    Point_3 q( 0.0, 1.0, 0.0);
    Point_3 r( 0.0, 0.0, 1.0);
    Point_3 s( 0.0, 0.0, 0.0);

    Polyhedron P;    // alternative constructor: Polyhedron P(4,12,4);
    P.make_tetrahedron( p, q, r, s);
    CGAL::set_ascii_mode( std::cout);
    std::copy( P.points_begin(), P.points_end(),
	       std::ostream_iterator<Point_3>( std::cout, "\n"));
    return 0;
}

12.3.6   Example with Circulator Writing Object File Format (OFF)

We create a tetrahedron and write it to std::cout using the Object File Format (OFF) [Phi96]. This example makes use of STL algorithms (std::copy, std::distance), STL std::ostream_iterator, and CGAL circulators. The polyhedral surface provides convenient circulators for the counterclockwise circular sequence of halfedges around a facet and the clockwise circular sequence of halfedges around a vertex.

However, the computation of the vertex index in the inner loop of the facet output is not advisable with the std::distance function, since it takes linear time for non random-access iterators, which leads to quadratic runtime. For better runtime the vertex index needs to be stored separately and computed once before writing the facets. It can be stored, for example, in the vertex itself or in a hash-structure. See also the following Section 12.4 for file I/O.

File: examples/Polyhedron/polyhedron_prog_off.cpp
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_3.h>
#include <iostream>

typedef CGAL::Simple_cartesian<double>               Kernel;
typedef Kernel::Point_3                              Point_3;
typedef CGAL::Polyhedron_3<Kernel>                   Polyhedron;
typedef Polyhedron::Facet_iterator                   Facet_iterator;
typedef Polyhedron::Halfedge_around_facet_circulator Halfedge_facet_circulator;

int main() {
    Point_3 p( 0.0, 0.0, 0.0);
    Point_3 q( 1.0, 0.0, 0.0);
    Point_3 r( 0.0, 1.0, 0.0);
    Point_3 s( 0.0, 0.0, 1.0);

    Polyhedron P;
    P.make_tetrahedron( p, q, r, s);

    // Write polyhedron in Object File Format (OFF).
    CGAL::set_ascii_mode( std::cout);
    std::cout << "OFF" << std::endl << P.size_of_vertices() << ' '
              << P.size_of_facets() << " 0" << std::endl;
    std::copy( P.points_begin(), P.points_end(),
               std::ostream_iterator<Point_3>( std::cout, "\n"));
    for (  Facet_iterator i = P.facets_begin(); i != P.facets_end(); ++i) {
        Halfedge_facet_circulator j = i->facet_begin();
        // Facets in polyhedral surfaces are at least triangles.
        CGAL_assertion( CGAL::circulator_size(j) >= 3);
        std::cout << CGAL::circulator_size(j) << ' ';
        do {
            std::cout << ' ' << std::distance(P.vertices_begin(), j->vertex());
        } while ( ++j != i->facet_begin());
        std::cout << std::endl;
    }
    return 0;
}

12.3.7   Example Using Euler Operators to Build a Cube

Euler operators are the natural way of modifying polyhedral surfaces. We provide a set of operations for polyhedra: split_facet(), join_facet(), split_vertex(), join_vertex(), split_loop(), and join_loop(). We add further convenient operators, such as split_edge(). However, they could be implemented using the six operators above. Furthermore, we provide more operators to work with polyhedral surfaces with border edges, for example, creating and deleting holes. We refer to the references manual for the definition and illustrative figures of the Euler operators.

The following example implements a function that appends a unit cube to a polyhedral surface. To keep track of the different steps during the creation of the cube a sequence of sketches might help with labels for the different handles that occur in the program code. The following Figure shows six selected steps from the creation sequence. These steps are also marked in the program code.

Steps in making a cube.

File: examples/Polyhedron/polyhedron_prog_cube.cpp
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_3.h>
#include <iostream>

template <class Poly>
typename Poly::Halfedge_handle make_cube_3( Poly& P) {
    // appends a cube of size [0,1]^3 to the polyhedron P.
    CGAL_precondition( P.is_valid());
    typedef typename Poly::Point_3         Point;
    typedef typename Poly::Halfedge_handle Halfedge_handle;
    Halfedge_handle h = P.make_tetrahedron( Point( 1, 0, 0),
                                            Point( 0, 0, 1),
                                            Point( 0, 0, 0),
                                            Point( 0, 1, 0));
    Halfedge_handle g = h->next()->opposite()->next();             // Fig. (a)
    P.split_edge( h->next());
    P.split_edge( g->next());
    P.split_edge( g);                                              // Fig. (b)
    h->next()->vertex()->point()     = Point( 1, 0, 1);
    g->next()->vertex()->point()     = Point( 0, 1, 1);
    g->opposite()->vertex()->point() = Point( 1, 1, 0);            // Fig. (c)
    Halfedge_handle f = P.split_facet( g->next(),
                                       g->next()->next()->next()); // Fig. (d)
    Halfedge_handle e = P.split_edge( f);
    e->vertex()->point() = Point( 1, 1, 1);                        // Fig. (e)
    P.split_facet( e, f->next()->next());                          // Fig. (f)
    CGAL_postcondition( P.is_valid());
    return h;
}

typedef CGAL::Simple_cartesian<double>     Kernel;
typedef CGAL::Polyhedron_3<Kernel>         Polyhedron;
typedef Polyhedron::Halfedge_handle        Halfedge_handle;

int main() {
    Polyhedron P;
    Halfedge_handle h = make_cube_3( P);
    return (P.is_tetrahedron(h) ? 1 : 0);
}

12.4   File I/O

Simple file I/O for polyhedral surfaces is already provided in the library. The file I/O considers so far only the topology of the surface and its point coordinates. It ignores a possible plane equation or any user-added attributes, such as color.

The default file format supported in CGAL for output as well as for input is the Object File Format, OFF, with file extension .off, which is also understood by Geomview [Phi96]. For OFF an ASCII and a binary format exist. The format can be selected with the CGAL modifiers for streams, set_ascii_mode and set_binary_mode respectively. The modifier set_pretty_mode can be used to allow for (a few) structuring comments in the output. Otherwise, the output would be free of comments. The default for writing is ASCII without comments. Both, ASCII and binary format, can be read independent of the stream setting. Since this file format is the default format, iostream operators are provided for it.

#include <CGAL/IO/Polyhedron_iostream.h>

template <class PolyhedronTraits_3>
ostream& ostream& out << CGAL::Polyhedron_3<PolyhedronTraits_3> P

template <class PolyhedronTraits_3>
istream& istream& in >> CGAL::Polyhedron_3<PolyhedronTraits_3>& P

Additional formats supported for writing are OpenInventor (.iv) [Wer94], VRML 1.0 and 2.0 (.wrl) [BPP95, VRM96, HW96], and Wavefront Advanced Visualizer object format (.obj). Another convenient output function writes a polyhedral surface to a Geomview process spawned from the CGAL program. These output functions are provided as stream operators, now acting on the stream type of the respective format.

#include <CGAL/IO/Polyhedron_inventor_ostream.h>


#include <CGAL/IO/Polyhedron_VRML_1_ostream.h>


#include <CGAL/IO/Polyhedron_VRML_2_ostream.h>


#include <CGAL/IO/Polyhedron_geomview_ostream.h>

template <class PolyhedronTraits_3>
Inventor_ostream& Inventor_ostream& out << CGAL::Polyhedron_3<PolyhedronTraits_3> P

template <class PolyhedronTraits_3>
VRML_1_ostream& VRML_1_ostream& out << CGAL::Polyhedron_3<PolyhedronTraits_3> P

template <class PolyhedronTraits_3>
VRML_2_ostream& VRML_2_ostream& out << CGAL::Polyhedron_3<PolyhedronTraits_3> P

template <class PolyhedronTraits_3>
Geomview_stream& Geomview_stream& out << CGAL::Polyhedron_3<PolyhedronTraits_3> P

All these file formats have in common that they represent a surface as a set of facets. Each facet is a list of indices pointing into a set of vertices. Vertices are represented as coordinate triples. The file I/O for polyhedral surfaces CGAL::Polyhedron_3 imposes certain restrictions on these formats. They must represent a permissible polyhedral surface, e.g., a 2-manifold and no isolated vertices, see Section 12.1.

Some example programs around the different file formats are provided in the distribution under examples/Polyhedron_IO/ and demo/Polyhedron_IO/. We show an example converting OFF input into VRML 1.0 output.

// examples/Polyhedron_IO/polyhedron2vrml.cpp
// ----------------------------------------

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

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

int main() {
    Polyhedron P;
    std::cin >> P;
    CGAL::VRML_1_ostream out( std::cout);
    out << P;
    return ( std::cin && std::cout) ? 0 : 1;
}

12.5   Extending Vertices, Halfedges, and Facets

In Section 12.3.5 we have seen how to change the default list representation

typedef CGAL::Polyhedron_3< Traits, 
                            CGAL::Polyhedron_items_3, 
                            CGAL::HalfedgeDS_default>      Polyhedron;

to a vector based representation of the underlying halfedge data structure. Now we want to look a bit closer at the second template argument, Polyhedron_items_3, that specifies what kind of vertex, halfedge, and facet is used. The implementation of Polyhedron_items_3 looks a bit involved with nested wrapper class templates. But ignoring this technicality, what remains are three local typedefs that define the Vertex, the Halfedge, and the Face for the polyhedral surface. Note that we use here Face instead of facet. Face is the term used for the halfedge data structure. Only the top layer of the polyhedral surface gives alias names renaming face to facet.

class Polyhedron_items_3 {
public:
    template < class Refs, class Traits>
    struct Vertex_wrapper {
        typedef typename Traits::Point_3 Point;
        typedef CGAL::HalfedgeDS_vertex_base<Refs, CGAL::Tag_true, Point> Vertex;
    };
    template < class Refs, class Traits>
    struct Halfedge_wrapper {
        typedef CGAL::HalfedgeDS_halfedge_base<Refs>                      Halfedge;
    };
    template < class Refs, class Traits>
    struct Face_wrapper {
        typedef typename Traits::Plane_3 Plane;
        typedef CGAL::HalfedgeDS_face_base<Refs, CGAL::Tag_true, Plane>   Face;
    };
};

If we look up in the reference manual the definitions of the three classes used in the typedefs, we will see the confirmation that the default polyhedron uses all supported incidences, a point in the vertex class, and a plane equation in the face class. Note how the wrapper class provides two template parameters, Refs, which we discuss a bit later, and Traits, which is the geometric traits class used by the polyhedral surface and which provides us here with the types for the point and the plane equation.

Using this example code we can write our own items class. Instead, we illustrate an easier way if we only want to exchange one class. We use a simpler face without the plane equation but with a color attribute added. To simplify the creation of a vertex, halfedge, or face class, it is always recommended to derive from one of the given base classes. Even if the base class would contain no data it would provide convenient type definitions. So, we derive from the base class, repeat the mandatory constructors if necessary - which is not the case for faces but would be for vertices - and add the color attribute.

template <class Refs>
struct My_face : public CGAL::HalfedgeDS_face_base<Refs> {
    CGAL::Color color;
};

The new items class is derived from the old items class and the wrapper containing the face typedef gets overridden. Note that the name of the wrapper and its template parameters are fixed. They cannot be changed even if, as in this example, a template parameter is not used.

struct My_items : public CGAL::Polyhedron_items_3 {
    template <class Refs, class Traits>
    struct Face_wrapper {
        typedef My_face<Refs> Face;
    };
};

When we use our new items class with the polyhedral surface, our new face class is used in the halfedge data structure and the color attribute is available in the type Polyhedron::Facet. However, Polyhedron::Facet is not the same type as our local face typedef for My_face, but it is derived therefrom. Thus, everything that we put in the local face type except constructors is then available in the Polyhedron::Facet type. For more details, see the Chapter 13 on the halfedge data structure design.

Pulling all pieces together, the full example program illustrates how easy the color attribute can be accessed once it is defined.

File: examples/Polyhedron/polyhedron_prog_color.cpp
#include <CGAL/Simple_cartesian.h>
#include <CGAL/IO/Color.h>
#include <CGAL/Polyhedron_3.h>

// A face type with a color member variable.
template <class Refs>
struct My_face : public CGAL::HalfedgeDS_face_base<Refs> {
    CGAL::Color color;
};

// An items type using my face.
struct My_items : public CGAL::Polyhedron_items_3 {
    template <class Refs, class Traits>
    struct Face_wrapper {
        typedef My_face<Refs> Face;
    };
};

typedef CGAL::Simple_cartesian<double>        Kernel;
typedef CGAL::Polyhedron_3<Kernel, My_items>  Polyhedron;
typedef Polyhedron::Halfedge_handle           Halfedge_handle;

int main() {
    Polyhedron P;
    Halfedge_handle h = P.make_tetrahedron();
    h->facet()->color = CGAL::RED;
    return 0;
}

We come back to the first template parameter, Refs, of the wrapper classes. This parameter provides us with local types that allow us to make further references between vertices, halfedges, and facets, which have not already been prepared for in the current design. These local types are Vertex_handle, Halfedge_handle, Face_handle, and there respective ..._const_handle. We add now a new vertex reference to a face class as follows. Encapsulation and access functions could be added for a more thorough design, but we omit that here for the sake of brevity. The integration of the face class with the items class works as illustrated above.

template <class Refs>
struct My_face : public CGAL::HalfedgeDS_face_base<Refs> {
    typedef typename Refs::Vertex_handle Vertex_handle;
    Vertex_handle vertex_ref;
};

More advanced examples can be found in the Section 13.3 illustrating further the design of the halfedge data structure.

12.6   Advanced Example Programs

12.6.1   Example Creating a Subdivision Surface

This program reads a polyhedral surface from the standard input and writes a refined polyhedral surface to the standard output. Input and output are in the Object File Format, OFF, with the common file extension .off, which is also understood by Geomview [Phi96].

The refinement is a single step of the sqrt(3)-scheme for creating a subdivision surface [Kob00]. Each step subdivides a facet into triangles around a new center vertex, smoothes the position of the old vertices, and flips the old edges. The program is organized along this outline. In each of these parts, the program efficiently uses the knowledge that the newly created vertices, edges, and facets have been added to the end of the sequences. The program needs additional processing memory only for the smoothing step of the old vertices.

subdivision examples

The above figure shows three example objects, each subdivided four times. The initial object for the left sequence is the closed surface of three unit cubes glued together to a corner. The example program shown here can handle only closed surfaces, but the extended example examples/Polyhedron/polyhedron_prog_subdiv_with_boundary.cpp handles surfaces with boundary. So, the middle sequence starts with the same surface where one of the facets has been removed. The boundary subdivides to a nice circle. The third sequence creates a sharp edge using a trick in the object presentation. The sharp edge is actually a hole whose vertex coordinates pinch the hole shut to form an edge. The example directory examples/Polyhedron/ contains the OFF files used here.

File: examples/Polyhedron/polyhedron_prog_subdiv.cpp
#include <CGAL/Cartesian.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/IO/Polyhedron_iostream.h>
#include <iostream>
#include <algorithm>
#include <vector>
#include <cmath>

typedef CGAL::Cartesian<double>                              Kernel;
typedef Kernel::Vector_3                                     Vector;
typedef Kernel::Point_3                                      Point;
typedef CGAL::Polyhedron_3<Kernel>                           Polyhedron;

typedef Polyhedron::Vertex                                   Vertex;
typedef Polyhedron::Vertex_iterator                          Vertex_iterator;
typedef Polyhedron::Halfedge_handle                          Halfedge_handle;
typedef Polyhedron::Edge_iterator                            Edge_iterator;
typedef Polyhedron::Facet_iterator                           Facet_iterator;
typedef Polyhedron::Halfedge_around_vertex_const_circulator  HV_circulator;
typedef Polyhedron::Halfedge_around_facet_circulator         HF_circulator;

void create_center_vertex( Polyhedron& P, Facet_iterator f) {
    Vector vec( 0.0, 0.0, 0.0);
    std::size_t order = 0;
    HF_circulator h = f->facet_begin();
    do {
        vec = vec + ( h->vertex()->point() - CGAL::ORIGIN);
        ++ order;
    } while ( ++h != f->facet_begin());
    CGAL_assertion( order >= 3); // guaranteed by definition of polyhedron
    Point center =  CGAL::ORIGIN + (vec / order);
    Halfedge_handle new_center = P.create_center_vertex( f->halfedge());
    new_center->vertex()->point() = center;
}

struct Smooth_old_vertex {
    Point operator()( const Vertex& v) const {
        CGAL_precondition((CGAL::circulator_size( v.vertex_begin()) & 1) == 0);
        std::size_t degree = CGAL::circulator_size( v.vertex_begin()) / 2;
        double alpha = ( 4.0 - 2.0 * std::cos( 2.0 * CGAL_PI / degree)) / 9.0;
        Vector vec = (v.point() - CGAL::ORIGIN) * ( 1.0 - alpha);
        HV_circulator h = v.vertex_begin();
        do {
            vec = vec + ( h->opposite()->vertex()->point() - CGAL::ORIGIN)
                       * alpha / degree;
            ++ h;
            CGAL_assertion( h != v.vertex_begin()); // even degree guaranteed
            ++ h;
        } while ( h != v.vertex_begin());
        return (CGAL::ORIGIN + vec);
    }
};

void flip_edge( Polyhedron& P, Halfedge_handle e) {
    Halfedge_handle h = e->next();
    P.join_facet( e);
    P.split_facet( h, h->next()->next());
}

void subdiv( Polyhedron& P) {
    if ( P.size_of_facets() == 0)
        return;
    // We use that new vertices/halfedges/facets are appended at the end.
    std::size_t nv = P.size_of_vertices();
    Vertex_iterator last_v = P.vertices_end();
    -- last_v;  // the last of the old vertices
    Edge_iterator last_e = P.edges_end();
    -- last_e;  // the last of the old edges
    Facet_iterator last_f = P.facets_end();
    -- last_f;  // the last of the old facets

    Facet_iterator f = P.facets_begin();    // create new center vertices
    do {
        create_center_vertex( P, f);
    } while ( f++ != last_f);

    std::vector<Point> pts;                    // smooth the old vertices
    pts.reserve( nv);  // get intermediate space for the new points
    ++ last_v; // make it the past-the-end position again
    std::transform( P.vertices_begin(), last_v, std::back_inserter( pts),
                    Smooth_old_vertex());
    std::copy( pts.begin(), pts.end(), P.points_begin());

    Edge_iterator e = P.edges_begin();              // flip the old edges
    ++ last_e; // make it the past-the-end position again
    while ( e != last_e) {
        Halfedge_handle h = e;
        ++e; // careful, incr. before flip since flip destroys current edge
        flip_edge( P, h);
    };
    CGAL_postcondition( P.is_valid());
}

int main() {
    Polyhedron P;
    std::cin >> P;
    P.normalize_border();
    if ( P.size_of_border_edges() != 0) {
        std::cerr << "The input object has border edges. Cannot subdivide."
                  << std::endl;
        std::exit(1);
    }
    subdiv( P);
    std::cout << P;
    return 0;
}

12.6.2   Example Using the Incremental Builder and Modifier Mechanism

A utility class CGAL::Polyhedron_incremental_builder_3 helps in creating polyhedral surfaces from a list of points followed by a list of facets that are represented as indices into the point list. This is particularly useful for implementing file reader for common file formats. It is used here to create a triangle.

A modifier mechanism allows to access the internal representation of the polyhedral surface, i.e., the halfedge data structure, in a controlled manner. A modifier is basically a callback mechanism using a function object. When called, the function object receives the internal halfedge data structure as a parameter and can modify it. On return, the polyhedron can check the halfedge data structure for validity. Such a modifier object must always return with a halfedge data structure that is a valid polyhedral surface. The validity check is implemented as an expensive postcondition at the end of the delegate() member function, i.e., it is not called by default, only when expensive checks are activated.

In this example, Build_triangle is such a function object derived from CGAL::Modifier_base<HalfedgeDS>. The delegate() member function of the polyhedron accepts this function object and calls its operator() with a reference to its internally used halfedge data structure. Thus, this member function in Build_triangle can create the triangle in the halfedge data structure.

File: examples/Polyhedron/polyhedron_prog_incr_builder.cpp
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_incremental_builder_3.h>
#include <CGAL/Polyhedron_3.h>

// A modifier creating a triangle with the incremental builder.
template <class HDS>
class Build_triangle : public CGAL::Modifier_base<HDS> {
public:
    Build_triangle() {}
    void operator()( HDS& hds) {
        // Postcondition: `hds' is a valid polyhedral surface.
        CGAL::Polyhedron_incremental_builder_3<HDS> B( hds, true);
        B.begin_surface( 3, 1, 6);
        typedef typename HDS::Vertex   Vertex;
        typedef typename Vertex::Point Point;
        B.add_vertex( Point( 0, 0, 0));
        B.add_vertex( Point( 1, 0, 0));
        B.add_vertex( Point( 0, 1, 0));
        B.begin_facet();
        B.add_vertex_to_facet( 0);
        B.add_vertex_to_facet( 1);
        B.add_vertex_to_facet( 2);
        B.end_facet();
        B.end_surface();
    }
};

typedef CGAL::Simple_cartesian<double>     Kernel;
typedef CGAL::Polyhedron_3<Kernel>         Polyhedron;
typedef Polyhedron::HalfedgeDS             HalfedgeDS;

int main() {
    Polyhedron P;
    Build_triangle<HalfedgeDS> triangle;
    P.delegate( triangle);
    CGAL_assertion( P.is_triangle( P.halfedges_begin()));
    return 0;
}