BooleanOps

Boolean operations and mesh intersection

Boolean operations

Geogram has functions to compute accurate intersections between triangulated meshes. The boolean_op demo program loads two meshes, applies a boolean operation to them, and saves the result. It is used as follows:

  $ boolean_op operation=union|intersection|difference filename1 filename2 resultfilename

The corresponding functions are implemented in geogram/mesh/mesh_surface_intersection.h

• mesh_union()
• mesh_intersection()
• mesh_difference()

Tetrahedral meshing with intersecting surfaces

The surface mesh intersection algorithm can also be used to tetrahedralize the interior of a collection of intersecting closed meshes. This operation is a bit different from the classical boolean operations, since one can optimize how to detect and discard internal boundary (in contrast with boolean operations that require a more complicated classification for all the triangles).

The figure shows (from left to right):

• a union of mutually intersecting closed meshes
• a cross-section, revealing many internal boundaries
• the computed intersection, with all internal boundaries removed
• the resulting tetrahedral mesh

The bundled vorpalite utility automatically does this pre-processing if the input mesh has intersections:

  $ vorpalite profile=tet inputfilename outputfilename

Internally, it uses the following functions:

• mesh_repair(M, MESH_REPAIR_DEFAULT, epsilon);
• fill_holes(M, max_hole_area);
• MeshSurfaceIntersection intersection(M);
• intersection.intersect();
• intersection.remove_internal_shells();
• mesh_repair(M, MESH_REPAIR_DEFAULT, epsilon);

They are called iteratively until no intersection is found. This is needed for two reasons:

• geogram does not have (yet) a snap-rounding algorithm that ensures that the exact points snapped to double precision do not yield additional intersections.
• tetgen sometimes moves the points and generates additional intersections

A simple intersect example program is also provided, used as follows:

  $ intersect inputmesh.obj <option1=true|false> .. <optionn=true|false>

It demonstrates several options of the MeshSurfaceIntersection class, described below:

option	description	default
verbose	display lots of messages	on
Delaunay	triangulate intersected facets with Delaunay	on
detect_intersecting_neighbors	also test neigboring triangles for intersection	on
normalize	normalize coord in [0,1] before intersection	on
remove_internal_shells	only keep external skin of union	off
dry_run (debug)	do not store intersection, only compute it	off
monster_threshold=nnn (debug)	save facets that have more than nnn intersections	off

Internals

The algorithm is implemented in MeshSurfaceIntersection::intersect(). It has the following steps:

1. compute the list of candidate triangles, using an AABB tree. See also tutorial on AABBs
2. compute the exact intersections between each pair of candidate triangles, see geogram/triangle_intersection.h
3. re-triangulate the triangles that have intersections. This uses MeshInTriangle, a subclass of the Constrained 2D Delaunay class that creates the new vertices in exact precision and that uses special predicates Internally, this uses vec2E, vec3E (vectors with arbitrary-precision coordinates stored as expansion_nt), and vec2HE, vec3HE (same thing in homogeneous coordinates, which makes it possible to divide coordinates without dividing (one multiplies the homogeneous factor w instead).
4. (optional) to implement boolean operations, triangles are sorted around non-manifold edges to indentify the formed 3D region. This is done by the RadialSort class.

Robustness

Some meshes, like this one from Thingi10K have a huge number degenerate intersections.

This specific mesh is a bit special, it has dramatic modeling errors that created a 3D array of triangles intersecting in three directions, generating N^3 intersections. Although it corresponds to nothing in real life, it is a good benchmark to test the robustness of the algorithms.

The test base contains also more realistic (though challenging) cases, here, like this one with many degenerate coplanar triangles:

The co-planar regions are coherently triangulated, thanks to the uniqueness of the Delaunay triangulation. The result can be filled with tetrahedra (in yellow on the second image). The key is the implementation of the in_circle predicate, here.

The optional geogramplus expansion package

Intersection points are represented using exact arithmetics, either with expansion_nt (Shewchuk's arithmetic expansions), or with the optional faster exact_nt number types implemented in the geogramplus expansion package (marketed by the Tessael company).

The geogramplus expansion package provides the exact_nt number type (arbitratry-precision floating point numbers), as well as geometric predicates based on it. It is both much faster (10x to 20x) and supports a much larger range of numbers (expansion_nt may underflow or overflow when multiplying very small or very large numbers, which may happen when computing surface intersections with very degenerate configurations). Geogramplus/exact_nt offers stronger guarantees, with a much wider exponent range:

	number type	speed	exponent range
geogram	expansion_nt	1x	[-1022,1023]
geogram+	exact_nt	10x-20x	[-2M, 2M]

Geogramplus provides:

• exact_nt: exact floating point number type based on GMP. Supports +,-,*,sign, conversion from double, conversion to interval_nt.
• vector types based on exact_nt: vec2Ex (2d), vec2Ex (3d), vec2HEx (2d homogeneous), vec3HEx (3d homogeneous) (through instanciations of vec2g<>, vec3g<>, vec2Hg<>, vec3Hg<>)
• predicates: orient_2d, orient_3d, orient_2d_projected, incircle_2d_SOS with airthmetic filters