Spherical RelateNG kernel

docs/plans/2026-06-15-relateng-spherical-kernel-design.md

@@ -0,0 +1,146 @@
+# RelateNG on the Spherical manifold — implementation design
+
+2026-06-15. Builds on the spike `2026-06-12-relateng-spherical-spike.md` /
+`.jl`. This document records the validated design decisions; the bite-sized
+task list lives in `2026-06-15-relateng-spherical-kernel.md`.
+
+## Goal
+
+Make `relate(RelateNG(; manifold = Spherical()), a, b)` work end-to-end —
+Phases 1+2+3 of the spike: the spherical kernel, the acceleration, and the
+exact paths — verified by a standalone spherical kernel-conformance suite and
+end-to-end tests.
+
+## Decisions (from the 2026-06-15 brainstorm)
+
+1. **Scope: full implementation**, not a stub scaffold. All 13 `rk_*`
+   functions for `::Spherical`, including the un-prototyped angle-ordering
+   cluster and the `Rational{BigInt}` exact paths.
+2. **Conformance: a separate standalone suite**
+   (`test/methods/relateng/kernel_conformance_spherical.jl`). The planar
+   `kernel_conformance.jl` is left untouched. The spherical suite uses
+   exactly-representable integer-xyz inputs in general position (so sign
+   predicates stay exact) and carries its own `Rational{BigInt}` differential
+   reference for the crossing-apex property.
+3. **Boundary: full end-to-end** (Phases 1+2+3). Out of scope: Phase 4's
+   broad differential-testing program (JTS-XML port, densified-geodesic
+   oracle); we keep validation to the conformance suite plus targeted
+   end-to-end smoke tests.
+4. **Spherical kernel point type is `UnitSphericalPoint{Float64}`** (not a
+   plain `NTuple{3,Float64}`): isbits, `GI.x/y/z`-accessible, and the kernel's
+   cross/dot math runs on it directly.
+
+## Architecture
+
+The engine is **already manifold-generic in its control flow**: `RelateNG`,
+`RelateGeometry`, `TopologyComputer`, both point locators,
+`RayCrossingCounter` and `edge_intersector` all take `m::Manifold`, and the
+accelerator selection already dispatches `Planar` vs generic `Manifold`. Two
+things are missing:
+
+1. **Point element type `P` is pinned to `Tuple{Float64,Float64}`** at
+   construction (`topology_computer.jl:53`, the `Set`/`Dict`/struct fields in
+   `relate_geometry.jl`, `point_locator.jl`, `indexed_point_in_area.jl`). The
+   structs are *already* parameterized on `P`; only the construction fixes it.
+   The enabling refactor derives `P` from the manifold and converts lon/lat →
+   xyz once at ingest.
+
+2. **The spherical kernel itself** (`kernel_spherical.jl`) — the 13 `rk_*`
+   methods on `::Spherical`.
+
+Correctness vs acceleration are cleanly separable: at `point_locator.jl:525`
+non-`Planar` already falls through to `rk_point_in_ring`, and
+`edge_intersector` already falls back to `NestedLoop()` for non-`Planar`. So
+once the kernel and `P`-threading exist, `relate(Spherical(), …)` is
+*correct* (if unaccelerated). Acceleration (3D edge index, lon-interval
+point locator, STR, prepared mode) is layered on after.
+
+### Point-type plumbing
+
+- `_kernel_point_type(::Planar) = NTuple{2,Float64}`,
+  `_kernel_point_type(::Spherical) = UnitSphericalPoint{Float64}`.
+- `_to_kernel_point(m, geopoint)`: identity `(x,y)` on `Planar`; lon/lat → xyz
+  via `UnitSphereFromGeographic`, **renormalized** (`normalize`) and
+  signed-zero-normalized per component, on `Spherical`.
+- `_node_point` / `_node_points` / `_canonical_segment` / `crossing_node` in
+  `kernel.jl` generalize from `(x,y)` to N components.
+- Thread the derived `P` through `RelateGeometry` →
+  `TopologyComputer`/`RelateSegmentString`/`RelatePointLocator`.
+
+### The spherical kernel
+
+Every predicate reduces to a sign of `det(u,v,w) = (u×v)·w`:
+
+- `rk_orient` → `sign((a×b)·c)`; exact via
+  `ExactPredicates.orient(tup(a),tup(b),tup(c),(0,0,0))`.
+- `rk_point_on_segment` → coplanarity + arc-span dot tests.
+- `rk_classify_intersection` → full `SegSegClass`: `SS_PROPER` (candidate
+  direction `d = n_a×n_b` strictly interior to both arcs), `SS_COLLINEAR`
+  (same great circle, overlapping spans), `SS_TOUCH`/endpoint flags.
+- `rk_point_in_ring` → meridian-arc crossing parity to a pole reference;
+  pole-insideness derived from the ring's signed area (S2 convention,
+  interior on the left of the directed ring).
+- angle cluster (`rk_quadrant`, `rk_compare_edge_dir`, `rk_crossing_dirs_ccw`,
+  `rk_is_crossing`, `rk_is_interior_segment`) → tangent-plane hemisphere
+  split with a reference direction `r` (a pole; fall back when the apex is
+  that pole).
+- `rk_nodes_coincide` → `Rational{BigInt}` on xyz components.
+- `rk_interaction_bounds` → `arc_extent` (spike-proven); area elements
+  extended by the six ±eᵢ axis-point test.
+- antipodal degeneracy → informative `ArgumentError` naming `AntipodalEdgeSplit`.
+
+### Acceleration
+
+- 3D `Extent{(:X,:Y,:Z)}` flow through `_relate_edge_index` unchanged
+  (`NaturalIndex` is dimension-generic). **Done**, but the per-segment extent
+  table had to be made manifold-aware too: it built a 2D endpoint box, which
+  misses a long arc's bulge and prunes real crossings — `_segment_extent`
+  now uses `arc_extent` on the sphere.
+- `edge_intersector.jl`: `_select_edge_set_accelerator(::Spherical, …)` →
+  tree accelerator; `_segment_envs_disjoint(::Spherical, …)` → 3D arc-extent
+  disjoint test. **Done.**
+- Prepared geometries: `_build_prepared_edge_index(::Spherical, …)` indexes A
+  in 3D; the dimension-generic `NaturalIndex` needs no separate STR ordering,
+  so prepared spherical relate matches unprepared. **Done.**
+- **SKIPPED (time-boxed):** the lon-interval indexed spherical point locator
+  (`SortedPackedIntervalRTree` over longitude intervals; antimeridian crossers
+  split). This is a pure *optimization* over the already-correct fall-through
+  at `point_locator.jl`'s `locate_on_polygonal` — Spherical takes the direct
+  `rk_point_in_ring` ring loop (the `loc.m isa Planar` guard skips the indexed
+  arm), which is correct but O(n) per query. Revisit if spherical point-in-area
+  becomes a hot path.
+
+### Antipodal edges
+
+Kernel throws on exactly-antipodal consecutive vertices. The opt-in
+`AntipodalEdgeSplit` correction (in `transformations/correction/`) inserts the
+lon/lat midpoint; this is the final, independently-deferrable task.
+
+### Ring/direction containment (`_ring_contains_dir`)
+
+Whether a direction `N` is interior to a ring (S2 interior-on-left) is decided
+by the ring's **winding** about `N`: project each edge onto the plane ⊥ `N` and
+sum the signed turn. Each edge contributes a turn of magnitude < π, so the sum
+has no atan2 branch ambiguity — a per-triangle signed-solid-angle sum
+(Van Oosterom–Strackee) does *not* have this property and reports spurious
+interiors for reference directions far from the loop's winding axis (e.g. an
+equatorial axis against a polar-cap ring). A single interior-on-left enclosure
+sweeps +2π, so `N` is interior iff the total exceeds π.
+
+This is exact for rings smaller than a hemisphere — the common geographic case,
+and all rings the conformance/end-to-end suites build. **Limitation:** a ring
+whose interior is *larger* than a hemisphere (interior = the bigger region)
+would under-report a non-encircled interior axis (winding 0 but interior).
+`_ring_contains_pole` (used for `rk_point_in_ring`'s pole reference) and the
+area-bounds axis extension (`_widen_area_axes`) both rest on this assumption.
+The axis extension only over-prunes when wrong, so it stays conservative-safe;
+super-hemisphere rings in `rk_point_in_ring` are explicitly out of scope here.
+
+## Out of scope
+
+Phase 4 differential-testing program (JTS-XML spherical port,
+densified-geodesic oracle). The conformance suite + end-to-end smoke tests are
+the proof for this work.
+
+Super-hemisphere rings (interior larger than a hemisphere) in
+`_ring_contains_dir` — see the limitation above.

src/GeometryOps.jl

@@ -94,6 +94,7 @@ include("methods/geom_relations/relateng/relate_predicates.jl")
 # topology-layer files (Task 6+) that will call the kernel functions.
 include("methods/geom_relations/relateng/kernel.jl")
 include("methods/geom_relations/relateng/kernel_planar.jl")
+include("methods/geom_relations/relateng/kernel_spherical.jl")
 # Node sections: after the kernel (uses `NodeKey`), before the point locator
 # (AdjacentEdgeLocator builds NodeSections).
 include("methods/geom_relations/relateng/node_sections.jl")
@@ -136,6 +137,7 @@ include("transformations/forcedims.jl")
 include("transformations/correction/geometry_correction.jl")
 include("transformations/correction/closed_ring.jl")
 include("transformations/correction/intersecting_polygons.jl")
+include("transformations/correction/antipodal_edge_split.jl")
 
 # Import all names from GeoInterface and Extents, so users can do `GO.extent` or `GO.trait`.
 for name in names(GeoInterface)

src/methods/geom_relations/relateng/edge_intersector.jl

@@ -232,6 +232,19 @@ function _select_edge_set_accelerator(::Planar, ssa_list, ssb_list)
         return DoubleSTRtree()
     end
 end
+# The same threshold heuristic on the sphere: the tree path is valid because
+# the segment extents are 3D great-circle arc extents (`_segment_extent`), and
+# `NaturalIndex` / the dual DFS / `Extents.intersects` are dimension-generic.
+function _select_edge_set_accelerator(::Spherical, ssa_list, ssb_list)
+    na = _total_segment_count(ssa_list)
+    nb = _total_segment_count(ssb_list)
+    if na < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS &&
+            nb < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS
+        return NestedLoop()
+    else
+        return DoubleSTRtree()
+    end
+end
 _select_edge_set_accelerator(::Manifold, ssa_list, ssb_list) = NestedLoop()
 
 _total_segment_count(ss_list) =
@@ -270,6 +283,11 @@ _segment_envs_disjoint(::Planar, a0, a1, b0, b1) =
     max(a0[1], a1[1]) < min(b0[1], b1[1]) ||
     min(a0[2], a1[2]) > max(b0[2], b1[2]) ||
     max(a0[2], a1[2]) < min(b0[2], b1[2])
+# On the sphere the per-segment extent is the bulge-aware 3D great-circle arc
+# extent, which is exact in xyz and has no antimeridian pathology — so this
+# prune is valid (and tighter than no pruning).
+_segment_envs_disjoint(::Spherical, a0, a1, b0, b1) =
+    !Extents.intersects(arc_extent(a0, a1), arc_extent(b0, b1))
 _segment_envs_disjoint(::Manifold, a0, a1, b0, b1) = false
 
 # The spatial index built over per-segment extents for the tree-accelerated
@@ -289,8 +307,8 @@ function process_edge_intersections!(tc::TopologyComputer,
         ssb_list::AbstractVector{<:RelateSegmentString},
         ::IntersectionAccelerator;
         m::Manifold = _manifold(tc), exact = _exact(tc))
-    extents_a, owners_a = _segment_extent_table(ssa_list)
-    extents_b, owners_b = _segment_extent_table(ssb_list)
+    extents_a, owners_a = _segment_extent_table(m, ssa_list)
+    extents_b, owners_b = _segment_extent_table(m, ssb_list)
     (isempty(extents_a) || isempty(extents_b)) && return nothing
     tree_a = _relate_edge_index(extents_a)
     tree_b = _relate_edge_index(extents_b)
@@ -387,7 +405,7 @@ function process_self_intersections!(tc::TopologyComputer,
         ss_list::AbstractVector{<:RelateSegmentString},
         ::IntersectionAccelerator;
         m::Manifold = _manifold(tc), exact = _exact(tc))
-    extents, owners = _segment_extent_table(ss_list)
+    extents, owners = _segment_extent_table(m, ss_list)
     isempty(extents) && return nothing
     tree = _relate_edge_index(extents)
     SpatialTreeInterface.dual_depth_first_search(Extents.intersects, tree, tree) do ia, ib
@@ -405,25 +423,31 @@ end
 
 # Flat per-segment extent list for a segment-string list, with the offset
 # table mapping each flat index back to (string index, segment index).
-function _segment_extent_table(ss_list)
-    extents = Extents.Extent{(:X, :Y), NTuple{2, NTuple{2, Float64}}}[]
+# Per-segment extent in the manifold's coordinate space: a planar 2D box from
+# the endpoints, or the bulge-aware 3D great-circle arc extent on the sphere
+# (the endpoint box would miss a long arc's bulge and wrongly prune crossings).
+_segment_extent(::Planar, p, q) = Extents.Extent(X = minmax(p[1], q[1]), Y = minmax(p[2], q[2]))
+_segment_extent(::Spherical, p, q) = arc_extent(p, q)
+_segment_extent_type(::Planar) = Extents.Extent{(:X, :Y), NTuple{2, NTuple{2, Float64}}}
+_segment_extent_type(::Spherical) = Extents.Extent{(:X, :Y, :Z), NTuple{3, NTuple{2, Float64}}}
+
+function _segment_extent_table(m::Manifold, ss_list)
+    extents = _segment_extent_type(m)[]
     owners = NTuple{2, Int}[]
     nseg = _total_segment_count(ss_list)
     sizehint!(extents, nseg)
     sizehint!(owners, nseg)
     for (si, ss) in enumerate(ss_list)
-        _push_segment_extents!(extents, owners, si, ss.pts)
+        _push_segment_extents!(m, extents, owners, si, ss.pts)
     end
     return extents, owners
 end
 
 # Function barrier: dispatch once per segment string, so the per-segment
 # loop stays statically typed even if `ss_list` has a non-concrete eltype.
-function _push_segment_extents!(extents::Vector, owners::Vector, si::Int, pts::Vector)
+function _push_segment_extents!(m::Manifold, extents::Vector, owners::Vector, si::Int, pts::Vector)
     for k in 1:(length(pts) - 1)
-        p = pts[k]
-        q = pts[k + 1]
-        push!(extents, Extents.Extent(X = minmax(p[1], q[1]), Y = minmax(p[2], q[2])))
+        push!(extents, _segment_extent(m, pts[k], pts[k + 1]))
         push!(owners, (si, k))
     end
     return nothing