Introduction

The CGAL library for computational geometry is truly a work of art.  It focuses on precision and accuracy above all else and yet manages to stay very flexible through perhaps the most comprehensive use of C++ metaprogramming that I have ever encountered.  CGAL deals efficiently and elegantly with rounding errors in IEEE floating point operations by escalating from IEEE floating point to exact numeric computation when rounding errors may occur.  The library is the product of 15+ years of development by some of the best computational geometry developers on the planet.

Using CGAL

CGAL is not the most accessible library to just pick up and use.  The template based generic programming paradigm can be difficult to wrap your head around at first but there are just enough examples to jump-start HelloWorld style apps.  Beyond that, the data structures are optimized for computational geometry not for obviousness.  Traversing the data structures requires a bit of study and thought to accomplish a task.

One of the more powerful aspects of the CGAL library is the Delaunay 3D Mesh Generator.  This generator can take a variety of geometric elements, such as polyhedrons, and generate a 3D triangulated mesh.  This operation is key to 3D printing as it is that triangulated mesh which is then sliced to create the layer-by-layer extrusion paths.  The OpenSCAD 3D CAD Modeller (http://www.openscad.org/) uses CGAL for binary polyhedral operations and mesh generation.

Writing a CGAL Mesh in STL file Format

There are not a lot of persistence formats supported in the CGAL library itself.  For 3D printing, the primary file format is arguably the STL (Standard Tessellation Language) format.  An STL file contains a list of triangular faces defined by 3 vertices and a normal to the facet.

Though the STL format is straightforward, it took a bit of poking around and experimentation to figure out how to traverse the mesh and order the vertices to insure that the mesh is manifold.  Getting the various template arguments right was also an occasional issue.

The code below takes a CGAL  Mesh_complex_3_in_triangulation_3 instance, a list of subdomains within the mesh complex and an open stream.  The template function writes the listed subdomains to the stream.  Each subdomain is written as a distinct solid in the STL file. It compiles under VS2010 and newer g++ releases.

#ifndef __MESH_TO_STL_H__
#define __MESH_TO_STL_H__</code>

#include &lt;CGAL/bounding_box.h&gt;
#include &lt;CGAL/number_utils.h&gt;

#include #include

template
struct SubdomainRecord
{
SubdomainRecord( const SubdomainIndex index,
const std::string label )
: m_index( index ),
m_label( label )
{}

SubdomainIndex m_index;
std::string m_label;
};

template
class SubdomainList : public std::list&lt;SubdomainRecord&gt;
{};

//
// The TriangulationPointIterator and TriangulationPointList template classes
// ease the task of iterating over the points associated with vertices in the mesh.
//

template
class TriangulationPointIterator : public std::iterator&lt;std::forward_iterator_tag, typename Triangulation::Point&gt;
{
typename Triangulation::Finite_vertices_iterator m_currentLoc;

public:

TriangulationPointIterator()
{}

TriangulationPointIterator( typename Triangulation::Finite_vertices_iterator&amp; vertIterator )
: m_currentLoc( vertIterator )
{}

TriangulationPointIterator(const TriangulationPointIterator&amp; mit)
: m_currentLoc( mit.m_currentLoc )
{}

TriangulationPointIterator&amp; operator++() {++m_currentLoc;return *this;}
TriangulationPointIterator operator++(int) {TriangulationPointIterator tmp(*this); operator++(); return tmp;}
bool operator==(const TriangulationPointIterator&amp; rhs) { return( m_currentLoc == rhs.m_currentLoc ); }
bool operator!=(const TriangulationPointIterator&amp; rhs) { return( m_currentLoc != rhs.m_currentLoc ); }
typename Triangulation::Point&amp; operator*() {return( m_currentLoc-&gt;point() );}
};

template
class TriangulationPointList
{
const Triangulation m_triangulation;

public:

TriangulationPointList( const Triangulation&amp; triangulation )
: m_triangulation( triangulation )
{}

TriangulationPointIterator begin()
{
typename Triangulation::Finite_vertices_iterator beginningVertex = m_triangulation.finite_vertices_begin();

return( TriangulationPointIterator( beginningVertex ));
}

TriangulationPointIterator end()
{
typename Triangulation::Finite_vertices_iterator endingVertex = m_triangulation.finite_vertices_end();

return( TriangulationPointIterator( endingVertex ));
}

};

//
// This function writes the ASCII version of an STL file. Writing the binary version should be
// a straightforward modification of this code.
//

template
std::ostream&amp;
output_boundary_of_c3t3_to_stl( const C3T3&amp; c3t3,
const SubdomainList&amp; subdomainsToWrite,
std::ostream&amp; outputStream )
{
typedef typename C3T3::Triangulation Triangulation;
typedef typename Triangulation::Vertex_handle VertexHandle;

// This is an ugly path to the Kernel type but this works and is all compile time anyway

typedef typename Triangulation::Geom_traits::Compute_squared_radius_3::To_exact::Source_kernel Kernel;

// Get the bounding box for the mesh so we can offset it into the all positive quadrant

TriangulationPointList pointList( c3t3.triangulation() );

typename Kernel::Iso_cuboid_3 boundingBox = CGAL::bounding_box( pointList.begin(), pointList.end() );

typename Kernel::Vector_3 offset( 1 - boundingBox.xmin(), 1 - boundingBox.ymin(), 1 - boundingBox.zmin() );

// Iterate over the facets in the mesh

std::array&lt;VertexHandle,3&gt; vertices;

for( SubdomainList::const_iterator itrSubdomain = subdomainsToWrite.begin(); itrSubdomain != subdomainsToWrite.end(); itrSubdomain++ )
{
// Write the solid prologue to the stream

outputStream &lt;&lt; "solid " &lt;&lt; itrSubdomain-&gt;m_label &lt;&lt; std::endl;
outputStream &lt;&lt; std::scientific; for( typename C3T3::Facets_in_complex_iterator itrFacet = c3t3.facets_in_complex_begin(), end = c3t3.facets_in_complex_end(); itrFacet != end; ++itrFacet) { // Get the subdomain index for the cell and the opposite cell typename C3T3::Subdomain_index cell_sd = c3t3.subdomain_index( itrFacet-&gt;first );
typename C3T3::Subdomain_index opp_sd = c3t3.subdomain_index( itrFacet-&gt;first-&gt;neighbor( itrFacet-&gt;second ));

// Both cells must be in the subdomain we are writing

if(( cell_sd != itrSubdomain-&gt;m_index ) &amp;&amp; ( opp_sd != itrSubdomain-&gt;m_index ))
{
continue;
}

// Get the vertices of the facet

for( int j=0, i = 0; i &lt; 4; ++i ) { if( i != itrFacet-&gt;second )
{
vertices[j++] = (*itrFacet).first-&gt;vertex(i);
}
}

// If the facet is not oriented properly, swap the first two vertices to flip it

if(( cell_sd == itrSubdomain-&gt;m_index ) != ( itrFacet-&gt;second%2 == 1 ))
{
std::swap( vertices[0], vertices[1] );
}

// Get the unit normal so we can write it

const typename Kernel::Vector_3 unit_normal = CGAL::unit_normal( vertices[0]-&gt;point(), vertices[1]-&gt;point(), vertices[2]-&gt;point() );

// Write the facet record to the file

outputStream &lt;&lt; "facet normal " &lt;&lt; unit_normal &lt;&lt; std::endl;
outputStream &lt;&lt; "outer loop" &lt;&lt; std::endl;
outputStream &lt;&lt; "vertex " &lt;&lt; vertices[0]-&gt;point() + offset &lt;&lt; std::endl;
outputStream &lt;&lt; "vertex " &lt;&lt; vertices[1]-&gt;point() + offset &lt;&lt; std::endl;
outputStream &lt;&lt; "vertex " &lt;&lt; vertices[2]-&gt;point() + offset &lt;&lt; std::endl;
outputStream &lt;&lt; "endloop" &lt;&lt; std::endl;
outputStream &lt;&lt; "endfacet" &lt;&lt; std::endl &lt;&lt; std::endl;
}

// Write the epilog for the solid

outputStream &lt;&lt; "endsolid " &lt;&lt; itrSubdomain-&gt;m_label &lt;&lt; std::endl;
}

// Return the stream and we are done

return( outputStream );
}

#endif // __MESH_TO_STL_H__

The code above includes a pair of helper template classes to ease iterating over mesh points for determining the bounding box for the mesh. The STL format requires that all points be positive but it doesn’t care about units.

The following code snippet demonstrates how to call the template function. The template parameter is inferred from the function arguments.

 SubdomainList<C3t3::Subdomain_index> subdomainsToWrite;

 subdomainsToWrite.push_back( SubdomainRecord<C3t3::Subdomain_index>( 0, std::string( "elephant" ) ));


 std::ofstream outputStream( "elephant.stl" );

 output_boundary_of_c3t3_to_stl( c3t3, subdomainsToWrite, outputStream );

 outputStream.close();

Conclusion

In later posts, I will follow up with CGAL examples of using Nef Polyhedra and performing the kinds of binary operations needed for CSG (Constructive Solid Geometry) applications.  Having the facility to mesh the polyhedra and persist the mesh in a file format that can then be consumed by a slicer and printed is a valuable stepping stone.