Category Archives: Computational Geometry

New Revision of Cork Computational Geometry Library – runs on Linux !

I have had enough free time lately to return to Cork and have made a couple key improvements to the build :

  • Builds and runs on Linux !
  • Generally 10% faster !
  • Moved to CMake build system
  • Script available to build 3rd Party dependencies
  • Updated to C++20
  • Updated to most recent Boost, TBB and MPIR libraries
  • Started vectorization with AVX2 SIMD instruction set
  • A few improvements to the regression test app
  • Added a few more unit tests
  • Faster OFF file output

Combined this makes for a much smoother ‘getting started’ experience. I will publish a Packer script that can be used to create a Ubuntu Mate 20.04 VM in VirtualBox or Proxmox for development.

The Github repository is here : At present I am working in the v0.9.0 branch.

I plan to move forward and bring the 3rd party dependencies up to date and build out more unit tests while working on performance improvements. I believe there are a number of places in the code that will benefit from AVX2 vectorization.

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.







Metal Parts from 3D Prints


Although 3D printing technology is advancing rapidly and home 3D printing is becoming both increasingly accessible and reliable, it will likely be a while before metal printing catches up with plastic Fused Filament Fabrication technology.  That said, there is a way to fabricate metal parts from some sets of 3D FFF designs.  In this blog post, I will describe a technique I have been using with some success to produce high quality metal parts from 3D prints.

Metal Clays and ‘Lost HIPS’ Molding

Metal clay is essentially a very low-tech approach to powder metallurgy.  Metal clays are a combination of atomized metal powder and organic, water soluble binders.  When soft, metal clay can be worked like a regular ceramic clay, dried to a hard yet brittle state and finally sintered in a conventional kiln to produce a solid metal piece.  The first metal clays were silver compounds but today metals such as bronze, brass, copper and steel are also readily available.

High Impact Polystyrene (HIPS) is a standard FFF filament, though definitely less popular than PLA or ABS.  My experience with HIPS as a general printing filament is quite good, it is easy to print with and can be printed very successfully on a Elmer’s glue coated, heated borosilicate glass plate.  An advantage of HIPS is that it is soluble in Limonene, a solvent derived from citrus fruit rinds.  As organic solvents go, limonene is about as safe as you can find, it is used medicinally for heartburn and GERD.  There is anecdotal evidence of it making a good margarita mixer…

Putting 3D printing with HIPS, metal clay, limonene as a solvent and finally kiln sintering together, we come up with the ‘Lost HIPS’ technique for creating metal parts from 3D prints.

Step 1 : Create a 3D Mold

The beauty and power of CAD/CAM is the ability to define, manipulate, visualize and refine 3D parts numerically prior to actually creating the physical part.  For our purposes, it is straightforward to take a 3D object definition in an STL file and create a mold for that part by performing a boolean difference operation between the 3D part and a rectangular cubiod (i.e. block).  For this post, I used the ‘Sun Medallion’ design I found on Thingiverse.  I then used OpenScad to create the cuboid and perform the binary difference to create a mold of the medallion.  I tweaked the mold along the way to strengthen the connections of the arms to the central solar disk but the design is still quite obviously that of Hank Dietz.

When printing a mold, try to find a good middle ground between a mold that is physically strong enough to work with but contains a minimum of HIPS material.  In ‘Lost HIPS’, all the mold material has to be dissolved by the limonene – so less is definitely more.


Sun Medallion mold printed in HIPS on a MakerGear M2 Printer.

HIPS is soluable in acetone as well, which means it can also be vapor polished in the same way as ABS.  I find vapor polishing to be helpful in smoothing the surface of the mold and sealing up any small holes or creases that may be left in the mold after printing – particularly when printing with thin layers.

Step 2: Fill the Mold with Metal Clay

At present, I am using FastFire Bronze Clay as it is relatively cheap (~$200/kg) and easy to work with though I have found it very sensitive to sintering temperatures.  I have also worked with PMC+ and it is easier to work with and very forgiving with respect to sintering temperatures but it is expensive (~$1500/kg).

When filling the mold, I have had the best luck painting the mold with water containing just a tiny amount of dish soap.  The water will cause the clay to form a thinner slurry next to the mold (much like ‘slip‘) and the detergent acts as a surfactant to help insure the slurry covers the entire base of the mold.  NB – do not use much water/detergent solution in the mold, as making the clay runny has lead to poor results for me.  I just paint the surface of the mold with a brush and that is it.  I usually put down a first layer of clay with an emphasis on insuring all the corners, nooks and crannies are filled and then fill the rest of the mold.  I use an old credit card to scrape off excess clay.

Once the mold is filled, I let it stand for a day to dry and sand the whole thing with a 200 grit sanding sponge to remove any excess clay.  It doesn’t take much sanding to get to a point where the finer features of the mold are visible again.  Finally, I wanted to make the sun medallion into necklaces for my daughters, so I added a loop to the back of the piece.  To make the loop, I used three pieces of HIPS together and placed a bit more clay over the HIPS and onto the back of the medallion.


Sun Medallion mold after drying, sanding to remove excess clay and the addition of the necklace loop.

Step 3: Lose the HIPS

Once the clay is dry, place the mold into a container of limonene and let the solvent do its work.  I use a glass container with a flourinated plastic lid that I found at Bed, Bath and Beyond (don’t forget your coupons).  Limonene is a solvent and will attack non flourinated plastics, though plastic gas cans and many consumer plastics are flourinated these days.  The more HIPS in the mold, the longer it will take to remove the material so expect anywhere from overnight to a couple days to get all the HIPS removed.  Fortunately, the metal clay does not appear to be nearly as sensitive to limonene as does the HIPS, so a couple days in a limonene bath does not appear to effect the clay.

Once the bulk of the HIPS is gone, I soak the piece in fresh limonene for a couple hours to get rid of the rest of the mold material and then dunk it in acetone for a minute or two.  The acetone serves two purposes.  First, it removes any gooey HIPS /limonene emulsion from the surface of the piece and second, it is a drying agent so after just a few minutes in air the piece is dry and can be worked a bit before sintering.


The bronze metal clay Sun Medallion after HIPS removal.  Note a bit of stringy HIPS material on the mesh holder and some HIPS left around the crease between the central disk and the sun arms.  This extra HIPS on the piece will burn off in the kiln.

Step 4: Make Repairs to the Clay Piece before Sintering

In its current state, the dried clay can be worked just as green clay can be worked.  I will typically sand off visible printing artifacts (i.e steps between layers), fill any voids with fresh clay and file off any excess material from the piece.  At this point as it is much easier to add/remove the clay material compared to post sintering.  I also find it helpful to use the water/detergent solution again to paint the surface a couple times to get a smoother finish.  The metal clay will absolutely reproduce every detail in the printed mold, so it you want a smoother aesthetic look – now is the time to take off the rough edges.

 Step 5: Burn out HIPS and Binder then Sinter

Once you are happy with the appearance of the piece, it is time to sinter.  I have a Paragon Caldera kiln which I love.  I did not get the digitally controlled version which I would suggest strongly for anyone looking to purchase a kiln.  I find the difference between a beautiful finished piece and an under-fired or over-fired piece to be just tens of degrees F.  Thus I end up having to watch my kiln closely as it finishes its ramp to insure it gets into the right temperature range and holds that range long enough to fully sinter the piece.

Pretty much any material other than silver needs to be fired in an anoxic (i.e oxygen free) environment.  For firing metal clays, someone far more clever than I figured out that one could easily create a locally oxygen free environment by burying the piece in carbon granules during firing.  This process works spectacularly well.  I will not go into the details here, there are plenty of references online.

I use a ceramic container for firing.  Firing in a stainless steel vessel leaves lots of black oxide in my kiln whereas the ceramic fiber pot leaves no residue whatsoever.  Having sintered with both, I also expect the pot to outlast a stainless vessel  as well.

I typically rest the piece on a piece of fiber kiln paper and then put the piece on the paper into the container filled with an inch or two of acid washed carbon granules.  I do not cover the piece but ramp my kiln to 400F and leave it there for an hour to burn off any remaining HIPS and the binder in the metal clay.


The cleaned up Sun Medallion on a piece of kiln paper.


The bronze clay Sun Medallion and kiln paper on a bed of carbon in the firing vessel.

After burning off any organic compounds left on the piece, I put another piece of kiln paper over the top of the piece and fill the container with carbon granules to within an inch of the top.  I put the lid on the container and ramp my kiln to 1450F and leave it there for an hour.  I then turn the kiln off and crack the lid to cool the piece quickly.

Step 6: Cleaning the Piece

It can take several hours for everything to cool to a point where it can be touched.  In particular, the vessel and carbon will hold heat well.

Once everything has cooled, I remove the piece from the carbon granules and clean it.  I use a Dremel tool with a wire brush to take the black scale off the surface of the piece and then use a bath of Picklean to remove the rest of the oxidation.  It may take a couple Picklean baths to really get the piece cleaned up but it is the only way I have found to get all the little details in the piece bright and shiny.


The finished product

Other Examples

I have created a number of other designs as well, below is a Tudor Rose extracted from a Thingiverse design.  What is interesting about this design is that there are regions in which upper layers overlap lower layers and if your printer does a decent job of bridging and you can force the metal clay into the mold, you can get a fairly intricate design which would be hard to fabricate using other means – like straight up stamping.


A Tudor Rose in Silver and Bronze.  In this example, the bronze piece has been slightly overfired and lost some of the detail of the original.  In contrast, the silver has retained much of the original detail to the point where the individual printing layers are clearly visible.


Though there are a number of steps involved with this process, for folks with more engineering talent than artistic talent this provides a way to create some gorgeous pieces simply by ‘turning the crank’.  After a few practice runs, I have found the ‘Lost HIPS’ process to be fairly straightforward.

Next I will probably try to fabricate some structural pieces using steel clay.  I have tinkered with steel clay once early on and I expect to have similar success with that material as well.

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.