Category Archives: CGAL

Cork – A High Performance Library for Geometric Boolean/CSG Operations

Gilbert Bernstein is currently a Ph.D. student at Stanford and had published some remarkable papers on computational geometry.  I was first drawn to his work by his 2009 paper on Fast, Exact, Linear Booleans as my interest in 3D printing led me to create some tooling of my own.  The various libraries I found online for performing Constructive Solid Geometry (CSG) operations were certainly good but overall, very slow.  CGAL is one library I had worked with and I found that the time required for operations on even moderately complex meshes was quite long.  CGAL’s numeric precision and stability is impeccable, the 3D CSG operations in CGAL are based on 3D Nef Polyhedra but I found myself waiting quite a while for results.

I exchanged a couple emails with Gilbert and he pointed me to a new library he had published, Cork.  One challenge with the models he used in his Fast, Exact paper is that the internal representation of 3D meshes was not all that compatible with other toolsets.  Though the boolean operations were fast, using those algorithms imposed a conversion overhead and eliminates the ability to take other algorithms developed on standard 3D mesh representations and use them directly on the internal data structures.  Cork is fast but uses a ‘standard’ internal representation of 3D triangulated meshes, a win-win proposition.

I’ve always been one to tinker with code, so I forked Gilbert’s code to play with it.  I spent a fair amount of time working with the code and I don’t believe I found any defects but I did find a few ways to tune it and bump up the performance.  I also took a swag at parallelizing sections of the code to further reduce wall clock time required for operation, though with limited success.  I believe the main problem I ran into is related to cache invalidation within the x86 CPU.  I managed to split several of the most computationally intensive sections into multiple thread of execution – but the performance almost always dropped as a result.  I am not completely finished working on threading the library, I may write a post later on what I believe I have seen and how to approach parallelizing algorithms like Cork on current generation CPUs.

Building the Library

My fork of Cork can be found here:  At present, it only builds on MS Windows with MS Visual Studio, I use VS 2013 Community Edition.  There are dependencies on the Boost C++ libraries, the MPIR library, and Intel’s Threading Building Blocks library.  There are multiple build targets, both for Win32 and x64 as well as for Float, Double and SSE arithmetic based builds.  In general, the Win32 Float-SSE builds will be the quickest but will occasionally fail due to numeric over or underflow.  The Double x64 builds are 10 to 20% slower but seem solid as a rock numerically.  An ‘Environment.props’ file exists at the top level of the project and contains a set of macros pointing to dependencies.

The library builds as a DLL.  The external interface is straightforward, the only header to include is ‘cork.h’, it will include a few other files.  In a later post I will discuss using the library in a bit more detail but a good example of how to use it may be found in the ‘RegressionTest.cpp’ file in the ‘RegressionTest’ project.  At present, the library can only read ‘OFF’ file formats.

There is no reason why the library should not build on Linux platforms, the dependencies are all cross platform and the code itself is pretty vanilla C++.  There may be some issues compiling the SSE code in gcc, but the non-SSE builds should be syntactically fine.

Sample Meshes

I have a collection of sample OFF file meshes in the Github repository:  The regression test program loads a directory and performs all four boolean operations between each pair of meshes in the directory – writing the result mesh to a separate directory.

These sample meshes range from small and simple to large and very complex.  For the 32bit builds, the library is limited to one million points and some of the samples when meshed together will exceed that limit.  Also, for Float precision builds, there will be numeric over or underflows whereas for x64 Double precision builds, all operations across all meshes should complete successfully.

When Errors (Inevitably) Occur

I have tried to catch error conditions and return those in result objects with some descriptive text, the library should not crash.  The code is very sensitive to non-manifold meshes.  The algorithms assume both meshes are two manifold.  Given the way the optimizations work, a mesh may be self intersecting but if the self intersection is in a part of the model that does not intersect with the second model, the operation may run to completion perfectly.  A different mesh my intersect spatially with the self intersection and trigger an error.

Meshes randomly chosen from the internet are (in my experience) typically not two manifold.  I spent a fair amount of time cleaning up the meshes in the sample repository.  If you have a couple meshes and they do not what to union – use a separate program like MeshLab to look over the meshes and double check that they are both in fact 2 manifold.


If you are interested in CSG and need a fast boolean library, give Cork a shot.  If you run into crashes or errors, let me know – the code looks stable for my datasets but they are far from exhaustive.







Writing a CGAL Mesh to an STL File


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 ( 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

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

SubdomainIndex m_index;
std::string m_label;

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.

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



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() );}

class TriangulationPointList
const Triangulation m_triangulation;


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.

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 ))

// 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 );



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.