CGAL + GEOS + Boost.Geometry, in depth: numeric models, validity, overlay and choosing a production geometry library

Predicates, constructions and why geometry code breaks

Every geometry algorithm splits into two kinds of operations. Predicates return a discrete answer: is point r left of, right of or on the line pq; is a point inside a circle. Constructions produce new coordinates: the intersection point of two segments, the centroid, a buffered outline. Predicates drive the combinatorial structure of the result, so a single wrong sign can make an algorithm loop forever or emit a polygon whose edges cross. Constructions only need to be close, unless you later feed the constructed points back into predicates, which overlay operations do all the time.

Floating-point arithmetic gets predicates wrong exactly when inputs are nearly degenerate, and real data is full of near-degeneracies: shared borders digitised twice, points snapped to a grid, coordinates that went through a projection and back. Each library is, at heart, a different answer to the question of how to stay consistent under those conditions.

Three libraries, three numeric models

CGAL (Computational Geometry Algorithms Library) follows the exact geometric computation paradigm. Every algorithm is parameterised by a kernel that defines number types, predicates and constructions. The two kernels most people use are Exact_predicates_inexact_constructions_kernel (Epick) and Exact_predicates_exact_constructions_kernel (Epeck). Epick evaluates each predicate in interval arithmetic first and falls back to exact arithmetic only when the interval straddles zero, so the common case runs at near double speed and every sign is right. Epeck additionally represents constructed points lazily as exact expression trees, so an intersection point fed back into later predicates is still exact. CGAL's breadth is unmatched: 2D and 3D triangulations, arrangements, Minkowski sums, Nef polyhedra, mesh generation and polygon-mesh processing.

GEOS (Geometry Engine, Open Source) is a C++ port of the Java JTS Topology Suite. Its data model is the OGC Simple Features model: Point, LineString, Polygon, their Multi variants and GeometryCollection, read and written as WKT or WKB. Coordinates are doubles. Robustness comes from algorithms designed for floating point: the overlay engine introduced in GEOS 3.9 (OverlayNG) nodes the input with snap-rounding, and every geometry can carry a precision model that snaps output to a fixed grid. GEOS is what PostGIS, Shapely, GeoPandas, QGIS and GDAL/OGR call underneath. Its C API is the supported, stable interface; the C++ API may change between releases.

Boost.Geometry is header-only generic C++. Algorithms work on any type that satisfies its concepts once you register it with an adapter macro, so you can run bg::intersects directly on your own point struct without copying. It implements the simple-features model plus an R-tree in boost::geometry::index and, unusually, coordinate systems beyond the plane: cs::spherical_equatorial and cs::geographic with ellipsoidal distance and area strategies. MySQL's spatial support is built on it.

The same overlay written three ways

The quickest way to feel the differences is to write the same job three times: intersect two squares, A from (0,0) to (4,4) and B from (2,2) to (6,6), and print the area of the overlap, which is 4. First CGAL. Boolean set operations require simple polygons with a counter-clockwise outer boundary, and the documentation recommends a kernel with exact constructions, because the intersection vertices become inputs to later predicates.

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Boolean_set_operations_2.h>
#include <iostream>
#include <list>

using K   = CGAL::Exact_predicates_exact_constructions_kernel;
using Pt  = K::Point_2;
using Pgn = CGAL::Polygon_2<K>;
using Pwh = CGAL::Polygon_with_holes_2<K>;

int main() {
  Pgn a, b;
  for (auto [x, y] : {std::pair{0,0}, {4,0}, {4,4}, {0,4}}) a.push_back(Pt(x, y));
  for (auto [x, y] : {std::pair{2,2}, {6,2}, {6,6}, {2,6}}) b.push_back(Pt(x, y));
  if (!a.is_simple() || !b.is_simple()) return 1;          // precondition, not a hint
  if (a.is_clockwise_oriented()) a.reverse_orientation();
  if (b.is_clockwise_oriented()) b.reverse_orientation();

  std::list<Pwh> out;
  CGAL::intersection(a, b, std::back_inserter(out));
  K::FT area = 0;
  for (const Pwh& p : out) {
    area += p.outer_boundary().area();
    for (auto h = p.holes_begin(); h != p.holes_end(); ++h) area += h->area();  // holes are CW: negative
  }
  std::cout << CGAL::to_double(area) << "\n";              // 4
}

GEOS through its re-entrant C API. Every call takes a context handle, and every object you receive is yours to destroy. Real code wraps these in RAII types.

#include <geos_c.h>
#include <stdio.h>

int main(void) {
  GEOSContextHandle_t ctx = GEOS_init_r();
  GEOSWKTReader *rd = GEOSWKTReader_create_r(ctx);
  GEOSGeometry *a = GEOSWKTReader_read_r(ctx, rd, "POLYGON((0 0,4 0,4 4,0 4,0 0))");
  GEOSGeometry *b = GEOSWKTReader_read_r(ctx, rd, "POLYGON((2 2,6 2,6 6,2 6,2 2))");

  GEOSGeometry *x = GEOSIntersection_r(ctx, a, b);   /* NULL on error: check it */
  double area = 0.0;
  if (x && GEOSArea_r(ctx, x, &area)) printf("%g\n", area);   /* 4 */

  GEOSGeom_destroy_r(ctx, x);
  GEOSGeom_destroy_r(ctx, b);
  GEOSGeom_destroy_r(ctx, a);
  GEOSWKTReader_destroy_r(ctx, rd);
  GEOS_finish_r(ctx);
  return 0;
}

Boost.Geometry. Note the template arguments: the default model::polygon is clockwise and closed, and bg::correct fixes orientation and closure in place. Skip it on counter-clockwise input and you can get negative areas and wrong overlay results without any exception.

#include <boost/geometry.hpp>
#include <iostream>
namespace bg = boost::geometry;

using Pt   = bg::model::d2::point_xy<double>;
using Poly = bg::model::polygon<Pt>;              // clockwise = true, closed = true
using MPoly = bg::model::multi_polygon<Poly>;

int main() {
  Poly a, b;
  bg::read_wkt("POLYGON((0 0,4 0,4 4,0 4,0 0))", a);   // counter-clockwise as written
  bg::read_wkt("POLYGON((2 2,6 2,6 6,2 6,2 2))", b);
  bg::correct(a); bg::correct(b);                       // now clockwise, as the type promises
  std::string why;
  if (!bg::is_valid(a, why)) { std::cerr << why << "\n"; return 1; }
  MPoly out;
  bg::intersection(a, b, out);
  std::cout << bg::area(out) << "\n";                   // 4
}

In Python, Shapely exposes GEOS with vectorised functions, which is usually the fastest route to a correct prototype: shapely.intersection(a, b).area gives the same 4.

Worked example: an invalid bowtie polygon

Overlay algorithms in all three libraries assume valid input, and real data often is not. Take the "bowtie" POLYGON((0 0, 4 4, 4 0, 0 4, 0 0)): its boundary crosses itself at (2, 2). Work it through:

1. The shoelace formula sums signed triangle areas: 0, then −16, then +16, then 0. The total is exactly 0, so a naive area() reports that a shape covering 8 square units covers nothing.
2. A validity check catches it. GEOS's GEOSisValid_r returns 0 and GEOSisValidReason_r names a self-intersection at the crossing point. Boost.Geometry's bg::is_valid(g, message) reports self-intersections the same way, and in CGAL Polygon_2::is_simple() returns false.
3. Repair splits the ring at the crossing. GEOS's GEOSMakeValid_r (and shapely.make_valid) return a MultiPolygon of two triangles, (0,0)-(2,2)-(0,4) and (2,2)-(4,0)-(4,4). Each has area 4, for a total of 8.
4. Decide whether repair is acceptable. If a parcel boundary self-intersects, splitting it may be fine for rendering and wrong for a land registry. Log the reason, count repairs per source, and route anything you cannot repair automatically to review.

Since GEOS 3.10 there are two repair algorithms, selectable through GEOSMakeValidWithParams_r: the original linework method, which re-nodes all edges and rebuilds faces, and a structure method that treats shells and holes separately and usually matches intent better for overlapping-hole errors. Older GEOS builds only provide the parameterless call.

Spatial indexes and prepared geometry

Real workloads are many-against-many: which of 5 million building footprints intersect each of 40,000 flood polygons. Exact overlay on every pair is hopeless, so every library pairs a spatial index for the filter step with exact predicates for the refine step.

• GEOS provides an STR-packed R-tree (GEOSSTRtree_create_r, GEOSSTRtree_query_r) and prepared geometries (GEOSPrepare_r, GEOSPreparedIntersects_r) which cache an internal index on one large polygon so that thousands of point-in-polygon tests against it stop re-scanning its edges. In Shapely, STRtree.query(geoms, predicate="intersects") runs filter and refine in one vectorised call.
• Boost.Geometry's bgi::rtree<Value, bgi::rstar<16>> supports bulk loading through its range constructor, spatial predicates such as bgi::intersects(box) and k-nearest queries with bgi::nearest(pt, k).
• CGAL leans toward computational-geometry structures: AABB_tree for segment and triangle soups, kd-trees for neighbour search, and arrangements when you need the full planar subdivision rather than pairwise answers.

The ordering in the pipeline diagram below is the single most useful habit: validate before you index, filter with bounding boxes, and only construct new geometry for pairs that survive the exact predicate. How R-trees answer window and nearest-neighbour queries explains why bulk-loaded trees beat incrementally built ones for static layers.

Failure modes

• Inexact constructions fed back into predicates. With CGAL's Epick kernel the predicates are exact, but intersection points are rounded to doubles. A Boolean operation that reuses those points can trip a precondition or produce a non-simple result. Use Epeck for overlays, or confine Epick to algorithms such as Delaunay triangulation that construct nothing that later predicates consume (see Delaunay triangulation).
• TopologyException in GEOS. Invalid input, or near-coincident edges after a projection, can make overlay fail. Validate first. OverlayNG retries internally with snapping heuristics, and setting a precision grid (GEOSGeom_setPrecision_r or shapely.set_precision) often removes the slivers that trigger it.
• Orientation and closure mismatches in Boost.Geometry. The type declares a winding, the data has another, and nothing checks unless you call bg::correct or bg::is_valid.
• Degrees treated as metres. GEOS is purely planar: the area of a polygon in longitude/latitude comes back in square degrees. Project to a suitable CRS first, or use Boost.Geometry's geographic strategies when you need ellipsoidal answers.
• Ownership and thread safety in GEOS's C API. Forgetting to destroy results leaks memory; sharing one context handle between threads races on its error state. Use one context per thread.
• Licence surprises. Much of CGAL is GPLv3+. Linking it into proprietary distributed software requires either compliance or a commercial licence from GeometryFactory. GEOS is LGPL-2.1, so dynamic linking is the norm. Boost.Geometry uses the permissive Boost Software License.

Trade-offs and how to choose

A rule of thumb that holds up in practice: if your data is GIS vector data and your users think in WKT, use GEOS (directly, or via PostGIS or Shapely). If you are writing a C++ service and want spatial predicates on your own structs with no licence friction, use Boost.Geometry. If you need guaranteed combinatorial correctness, 3D, or algorithms the others lack (constrained triangulations, mesh Booleans, arrangements), use CGAL and choose the kernel deliberately. Mixing is common: a pipeline might use GEOS for ingestion and repair, then CGAL for one exact step, converting through WKT or raw coordinates at the boundary.

What to do next

1. Collect 50 real geometries from your own data, including the ugliest ones, and run each library's validity check on all of them. Record the reasons, not just counts.
2. Write the three-square example above in the library you are leaning toward, then add the bowtie and confirm your code reports and repairs it instead of trusting area.
3. Put validate, then normalise orientation, then index, then exact predicate, then construct, in that order, into your pipeline and add a metric for each rejection.
4. Pin the library version in your build, and keep a regression folder of every geometry that ever crashed or produced a wrong result.
5. Check the licence against how you distribute your software before you write significant code against CGAL.
6. Go deeper on the primitives: convex polygon intersection and line-sweep algorithms are what these libraries' overlay engines are built from.