Add RelateNG: exact DE-9IM predicate engine ported from JTS, on planar and spherical manifolds by asinghvi17 · Pull Request #416 · JuliaGeo/GeometryOps.jl

asinghvi17 · 2026-06-11T15:57:29Z

Formerly the head of a six-PR stack — #422, #438, #439, #440, and #441 are collapsed into this PR so the port can be reviewed as a whole. Each closed PR remains readable as a focused diff of its layer (review map below).

Merge order note: the branch merges spherical-arc-extent (#434) to consume the shared spherical substrate, so the diff shows #434's content until that lands on main. Merge #434 first.

Summary

Implements the JTS RelateNG DE-9IM predicate engine in GeometryOps, on both Planar and Spherical manifolds, per the designs in docs/plans/ (included on the branch).

Engine (planar port). Exact, construction-free evaluation: every geometric decision goes through an rk_* kernel backed by exact orientation predicates; proper crossings are symbolic nodes keyed by the canonicalized segment pair, with coincidence decided in exact Rational{BigInt} arithmetic — no floating-point intersection coordinate is ever constructed for node identity. The topology layer is a faithful, file-by-file diffable port (TopologyComputer, RelateGeometry, RelatePointLocator, node/edge topology, predicate framework), with predicate-type specialization, tree-accelerated edge enumeration with early exit, and prepared mode (prepare). Public surface is opt-in (design D4): relate, DE9IM, RelateNG, prepare, boundary-node rules — existing defaults untouched.

Spherical kernel. The same engine runs on the sphere by swapping the kernel: points convert once to unit-sphere kernel points, orientation and arc crossings are decided exactly (3D orientation + exact-rational crossing classification, following S2's predicate structure), boundary/interior ring classification stays exact, and exactly-antipodal edges are rejected at ingest with a pointer to the AntipodalEdgeSplit correction. Acceleration carries over — segment extents are bulge-aware great-circle arc extents, and the same trees prune the dual traversals in 3D.

Spherical interior semantics. A lon/lat ring on the sphere bounds two finite regions, so the port defines which one a polygon means — and follows the ecosystem consensus (R s2/sf, spherely, BigQuery, S2's unoriented ingest): by default the interior is the enclosed (smaller-than-hemisphere) region, independent of winding, so shapefile-wound real data (e.g. Natural Earth's CW shells) means what its producers intend. Spherical(; oriented = true) opts into winding-authoritative interiors (S2 InitOriented: interior on the left, shells CCW/holes CW), which makes regions larger than a hemisphere representable. One orientation bit — an S2-style turning-angle curvature with exact predicate signs — is the single authority shared by edge topology, point location, and Extents.extent(Spherical(), ...).

Spherical validity contract. A ring that is simple in planar lon/lat can self-intersect when its edges are reinterpreted as geodesics (Natural Earth's 110m Sudan does — a 0.5 km² needle lobe that no planar validity tool can see), and undetected it corrupts containment. Following the S2-ecosystem contract: prepare validates ring self-crossings by default on Spherical (off by default on Planar, where standard tooling already covers validity and prepared-geometry precedent is no validation — validate = true opts in), erroring with the crossing edge pair and pointing to the new CrossingEdgeSplit correction, which splits rings at geodesic crossings and keeps both lobes as shells (even-odd, matching S2Builder's undirected repair). Unprepared calls don't pay a validation toll; instead, enclosed-region containment parity is anchored at a definitionally-exterior point (antipode of the ring's vertex mass), so unvalidated invalid input degrades even-odd-gracefully like the planar path instead of inverting globally — spherical vs planar intersects now agree on all 15,051 Natural Earth 110m country pairs, Sudan included.

Performance. Engine internals are type-erased over input geometry types and a PrecompileTools workload covers the common predicates, fixing a Julia 1.12 first-call compile blowup (CI matrix jobs: ~113 min → 13–27 min; crosses(poly, poly) first call 9.5 s → 0.01 s). Unprepared point-in-area drops 862 → 63 µs on Natural Earth Canada (single-pass extent kernel + JTS-parity envelope short-circuits; 22 µs with extent-stamped inputs, faster than plain GEOS), and prepared spherical point queries drop 17.6 ms → 2.2 µs via a longitude-interval indexed locator (the 1-D-stab analog of JTS's y-interval trick, with a meridian reference arc and a build-time pole anchor). The spherical exact kernel gains certified Float64 fast paths (the proper-crossing branch reduced algebraically to four orientation signs, plus a certified arc-span triage) — audited sign-identical against the pure-exact path on 1.65M arc pairs, bringing spherical classification below planar cost per clearly-separated pair (GADM Canada × USA spherical edge-pair collection: 3.12 s → 674 ms). Real-data and TTFX benchmark suites are committed under benchmarks/ with recorded representative output.

Shared substrate. The kernel consumes the shared geometry layer instead of forking it: interaction bounds compose spherical_arc_extent and the crossing-parity region extents from #434 over kernel-converted points; ring containment is the promoted UnitSpherical.spherical_ring_contains with injectable predicates (the kernel injects its exact ones; the extent path keeps the tolerance-based defaults). Index structures are shared too: FlexibleRTrees.RTree now carries the collection it indexes (tree.data, with precomputed-extents support and a zero-copy Unsorted path), the relate edge index is a natural-order RTree of segment owners, and the prepared point-in-area interval index is RTree(STR(), segs; extents = y_intervals) — STR in one dimension is JTS's midpoint sort. No private index structures remain in the port.

Also on this branch. Spherical area's Eriksson fan triangulation is restored to its native signed Van Oosterom–Strackee form (atan2 instead of atan∘abs): unsigned fan triangles overcounted concave polygons badly on real data (Norway +174%, Chile +432% vs a BigFloat reference), and the unsigned form also mis-branched for fan triangles with negative denominator. Magnitude is bit-identical where the denominator is positive, so the small-polygon accuracy that motivated Eriksson over Girard is retained (verified against BigFloat down to 1e-5° features). Found by area-conservation checks during OverlayNG design spikes on this foundation.

Review map

The collapsed PRs each cover one layer, in stack order:

PR	Layer
#416 (original commits)	Planar RelateNG engine + topology layer
#422	Spherical `rk_` kernel, kernel-point threading, spherical acceleration, S2 crossings conformance
#438	Substrate unification: extent delegation, `spherical_ring_contains` promotion/graft, `rk_bounds_*` pruning
#439	Cleanups: twin-method merges, `DoubleNaturalTree` naming, lazy point-in-area auto-index removal (Java `isPrepared`-only rule restored), dead code
#440	Payload-carrying `RTree`; edge-index side-tables (`PreparedEdgeIndex`, `owners`) dissolved
#441	`SortedPackedIntervalRTree` re-expressed as 1-D `RTree(STR(), ...)`

Validation

All 64 methods of JTS RelateNGTest plus every component JUnit suite, ported and passing (engine file: 1,868 assertions).
Full vendored JTS XML conformance suites (20 files, 6,537 ops) with a 2-entry documented skiplist — both are cases where GEOS ≥ 3.13 RelateNG agrees with GeometryOps against the legacy XML expectations.
Seeded differential fuzzing vs GEOS native RelateNG: zero divergences (50,000-case deep runs on the planar engine; 500 cases per gate run).
Spherical: 51,943-assertion kernel conformance suite, S2 edge-crossing conformance ported from s2geometry, and 186 end-to-end spherical relates including prepared mode, both interior modes, winding-independence, bowtie even-odd degradation, and the antipodal ingest guard.
Spherical vs planar intersects agree on all 15,051 Natural Earth 110m country pairs (previously 161 disagreements, all from Sudan's geodesically self-intersecting 110m ring — resolved by the definitional-exterior parity anchor; S2 forced past validation inverts the same ring, verified against s2geometry directly).
GADM full-resolution census: 43 countries / 25.7M vertices, zero spherical or planar validity defects (the Sudan needle is a low-resolution simplification artifact) and zero cross-engine disagreements on real land-border pairs; benchmarks at 3.9M-vertex Canada committed (benchmarks/relateng_gadm.jl) — prepared point-in-area 1.8 µs planar / 2.3 µs spherical.
Allocation budgets and type-stability checks.