diff --git a/docs/plans/2026-07-16-overlayng-noding-substrate.md b/docs/plans/2026-07-16-overlayng-noding-substrate.md new file mode 100644 index 0000000000..2d21d1a697 --- /dev/null +++ b/docs/plans/2026-07-16-overlayng-noding-substrate.md @@ -0,0 +1,377 @@ +# OverlayNG port — exact-arrangement design, and the phase 1 (noding substrate) implementation spec + +Status: design settled 2026-07-16, after a four-agent survey of JTS `overlayng`/`noding`, +s2geometry's boolean-operation stack, and this repo's foundation, followed by four validation +spikes (S1–S4, verdicts in the appendix). Phase 1 is specified here at implementation depth; +phases 2–3 are outlined so phase-1 authors know their consumer. This document is written to be +consumed by implementation subagents: where it names an existing function, treat the name as a +strong pointer but verify the signature against the code before use. + +Reference checkouts (read-only): `/Users/anshul/temp/GO_jts/jts` (JTS, see +`modules/core/src/main/java/org/locationtech/jts/operation/overlayng/` and `…/noding/`), +`/Users/anshul/temp/GO_jts/s2geometry`. This repo: branch `relateng` (tip ≥ `69e416484`). + +--- + +## 0. The one governing decision + +Overlay (`intersection` / `union` / `difference` / `symdifference`) is computed as an **exact +arrangement with symbolic nodes, rounded once at emission** — the extension of the RelateNG +port's design D2 to constructive output. + +- **No constructed coordinate ever enters a decision.** Node identity, ordering, angular sort, + labeling, and result assembly are all decided by exact kernel predicates over input vertices + and symbolic crossing keys. Float64 values appear only as *filters* (certified bounds, + escalate on uncertainty) and at final *emission*. +- **There is no snapping, no snap-rounding, no precision model, no retry ladder — permanently.** + JTS's `SnappingNoder`/`SnapRoundingNoder`/`OverlayNGRobust` machinery exists to patch an + inexact-noding substrate (constructed Float64 points feeding later float decisions). This + substrate has no such defect, so that machinery is dead weight; do not port it, do not + reintroduce tolerances. If a fixed-precision *output* mode is ever wanted, it is a + post-process on the exact result, never part of the computation. +- The pipeline keeps JTS OverlayNG's shape — **node → half-edge graph → label → build → emit** + — so the engine port stays file-by-file diffable against Java, exactly like the RelateNG port. + +Motivation: the current Foster–Hormann clipping's recurring bug class (order-dependent wrong +results and `TracingError`s from inexact constructed intersections + ent/exit tracing — +issues #114, #281, #193, #335, #191, #337, #417, #129). The exact arrangement makes that class +unrepresentable, and the spikes showed it is *faster* than GEOS at the noding stage, not slower. + +Evidence base (details in appendix): S1 measured the symbolic noder at **0.03×–0.17× of +LibGEOS's entire intersection op** with zero exact-predicate fallbacks on real data; S2 +validated the graph/label/build structural claims end-to-end against LibGEOS bit-for-bit on +planar cases; S3 audited ~930k noded edges and found **zero rounding artifacts on valid +input**, and proved a certified double-double emission path (100% certified, 273× faster than +rational); S4 (landed, `8d938e832`) brought spherical exact classification to planar parity. + +## 1. Architecture and phasing + +File layout (names mirror JTS where a counterpart exists, for diffability): + +``` +src/methods/clipping/overlayng/ +├── noding/ # ← PHASE 1 (this spec) +│ ├── noded_arrangement.jl # types: NodedArrangement, NodedEdge, node table +│ ├── collect.jl # stage 1: candidate enumeration + classification +│ ├── node_identity.jl # stage 3: two-tier node grouping +│ ├── split.jl # stage 4: noded-edge emission (topology only) +│ └── emit.jl # coordinate realization (certified fast path + exact) +├── edge_source.jl # ← PHASE 2: EdgeSourceInfo/depth_delta consumption +├── overlay_label.jl # OverlayLabel (plain struct, per JTS) +├── overlay_graph.jl # half-edge graph, EdgeMerger by node pair +├── overlay_labeller.jl # multi-pass label propagation +├── maximal_edge_ring.jl # result ring linking (+ minimal-ring split) +├── polygon_builder.jl # shells/holes via rk_point_in_ring +├── line_builder.jl # result lines (skip JTS's dead merged path) +├── intersection_point_builder.jl +└── overlay_ng.jl # ← PHASE 3: isResultOfOp, driver, ops, public opt-in API +``` + +Kernel additions (phase 1) go in the existing kernel files +(`src/methods/geom_relations/relateng/kernel*.jl`), not in overlayng — they are manifold +predicates, not overlay logic. GeometryOps is a single module, so placement is file +organization only; graduating `noding/` to a top-level `src/noding/` happens if and when a +second consumer (buffer) materializes — not before. + +Branch/PR mechanics (decided): branch `overlayng` off `relateng`, one stacked PR per phase. +Phase 1 is internal-only — every new name is un-exported and `_`-prefixed or clearly internal; +no public API or docs-page changes (avoids the docstring `@ref` trap entirely; if any docstring +is written, `@ref` may target **exported names only** — this killed CI twice on relateng). + +- **Phase 1 — noding substrate** (§2): geometries → `NodedArrangement` + emission. ~800–1,000 + SLOC + tests. +- **Phase 2 — engine core** (§3): half-edge graph, labels, builders. ~4–4.5k SLOC, faithful + JTS port with the amendments in §3. +- **Phase 3 — ops + API** (§4): the four ops, mixed dimensions, `OverlayNG{M}` opt-in + algorithm, differential validation harnesses, PrecompileTools workload, benchmarks. + +## 2. Phase 1: the noding substrate + +### 2.1 Contract: `NodedArrangement` + +```julia +struct NodedEdge + string_idx :: Int32 # index into segstrings + seg_idx :: Int32 # segment within the parent string + node_lo :: Int32 # node id at sub-segment start (in traversal order of the parent) + node_hi :: Int32 # node id at sub-segment end +end + +struct NodedArrangement{P} # P = kernel point type; exactly two instantiations + segstrings :: Vector{RelateSegmentString{P}} # ingested inputs (NOT text — JTS "SegmentString") + nodes :: NodeTable{P} # see §2.4 + seg_nodes :: … # per-segment ordered interior node-id lists + edges :: Vector{NodedEdge} +end +``` + +Constructor-guaranteed invariants (each one is a test): + +1. Every proper crossing between any two indexed segments appears as one shared node id on + both parents. +2. Node ids are geometrically unique — coincident symbolic keys are merged (§2.4). +3. Per-segment interior node lists are exactly ordered along the parent (§2.5) and deduped; + no `NodedEdge` is zero-length (S1 finding: one geometric node is routinely reported by + several candidate pairs — dedup before splitting is mandatory). +4. No Float64-constructed coordinate influenced any of 1–3. Floats appear only inside + certified filters whose uncertain cases escalated to exact predicates. +5. A `NodedEdge` carries no geometry: its shape is a lookup into its parent segment. All + source metadata (owner, ring id, dimension) is reached through `string_idx` — nothing is + copied, so nothing can desynchronize. (Stronger than JTS, which copies a `data` payload + onto split substrings.) + +Naming note: the field is `segstrings`, never `strings` — these are JTS SegmentStrings +(polyline chains of kernel points), and the bare name collides confusingly with `Base.String`. + +### 2.2 Ingest + +Reuse `RelateSegmentString{P}` unchanged (fields incl. `is_a`, `dim`, `id`, `ring_id`, +`pts::Vector{P}`; built by `_rss_create_ring`/`_rss_create_line`, repeated points removed with +kernel `==` at ingest). Do **not** add fields: `is_hole` is derivable (`ring_id > 1` for +polygon rings — verify the convention in `relate_geometry.jl`), and `depth_delta` belongs to +phase 2's edge-source step, not to the string. Kernel-point conversion happens once here (the +S1 lesson that extraction dominates sparse pairs → keep arrangement construction separable +from ingest so a future prepared overlay converts once; phase 1 just keeps the two steps as +separate functions). + +Validity contract: inputs are contractually valid, same as relate — spherical ring +self-crossings are caught by the existing `prepare(...; validate)` path (which uses the same +enumeration machinery), and invalid input is the user's problem to `fix` (`CrossingEdgeSplit` +/ `AntipodalEdgeSplit`). The noder itself performs no validation and no self-noding of A +against A. Input spikes (a–b–a) need no special handling here: they surface as collinear +coincident edges and phase 2's merger resolves them (JTS `DIM_COLLAPSE` semantics). + +### 2.3 Stage 1 — collect (`collect.jl`) + +- Candidate enumeration: per side, `_relate_edge_index(m, ss_list)` (the existing payload + `RTree(Unsorted(), owners; extents…)` over per-segment extents — planar boxes / + `spherical_arc_extent` 3-D boxes), then + `SpatialTreeInterface.dual_depth_first_search(Extents.intersects, tree_a, tree_b)` for the + A×B join. Accept optional caller-supplied prebuilt trees (a `PreparedRelate` carries exactly + this index) — an optional argument, **no** new prepare type, no auto-preparing. +- Per candidate pair: `rk_classify_intersection(m, a0,a1,b0,b1; exact=True())` (post-S4 this + has certified Float64 fast paths; do not re-filter above it). + - `SS_PROPER` → `crossing_node(...)` recorded on **both** parent segments. + - `SS_TOUCH` / `SS_COLLINEAR` → `vertex_node` copies recorded on each segment in whose + *interior* the flagged vertex lies (the classification's `*_on_*` flags say which). S1 + verified across ~150k real classifications that every touch/collinear intersection is an + input vertex; keep the defensive assertion anyway — it is cheap and the claim is + load-bearing. +- Single-threaded (decided). The collect loop is embarrassingly parallel if ever needed; S1 + says it is not where time goes. + +### 2.4 Stage 3 — node identity (`node_identity.jl`): two tiers, no canonical key + +Node ids are assigned by a two-tier scheme (this supersedes the earlier +`rk_canonical_node_key` idea from spike S2 — computing a canonical rational key costs +rational arithmetic *per node* (same order as emission, 6–18 µs) to optimize a path S1 +measured firing **zero** times on real data; the float sweep costs ~1 ms total): + +- **Tier 1 (free):** a `Dict{NodeKey,Int32}` — egal-equal keys merge (kernel points normalize + signed zeros, so vertex-node bit-equality is already canonical). +- **Tier 2 (rare):** geometric coincidence across *different* keys (two distinct segment + pairs crossing at one point; a crossing landing exactly on a third string's vertex). + Compute a throwaway Float64 approximation of each crossing node with a certified error + radius; sort by one coordinate; sweep for overlapping intervals; confirm candidates with + `rk_nodes_coincide(...; exact=True())`; merge confirmed pairs with a union-find over + **candidates only** (bounded by candidate count — zero on all real data measured, a handful + on constructed degree-6 cases). + +The approximate position + radius used here is the same computation as the emission fast path +(§2.6) minus the final certification — share the code. + +### 2.5 Stage 2 — ordering: new kernel predicate `rk_compare_along_segment` + +Added to the `rk_` contract in `kernel.jl` with per-manifold implementations, documented in +the same style as its siblings: + +``` +rk_compare_along_segment(m, s0, s1, na, nb; exact) -> -1 / 0 / +1 + Order of two nodes (NodeKeys) along the oriented segment (s0, s1). +``` + +- Float stage: approximate along-segment parameter (planar) / signed discriminant against the + crossing directions (spherical), each with a certified error bound; adjacent pairs whose + float gap exceeds the summed bounds are decided in float. +- Exact fallback (lazy — computed only for pairs the filter cannot separate): planar — + compare exact `Rational{BigInt}` parameters derived from `_exact_crossing_point`; spherical + — triple-product sign tests against `_sph_crossing_dir` directions. Returning 0 is legal + and means the nodes coincide — but by construction stage 3 ran first, so equal-order nodes + have already merged; assert this. +- S1 measured the float filter resolving 100% of 124,500 adjacent comparisons in the dense + regime (240–590× vs exact-always) and the whole stage a no-op on real data (no real segment + acquired ≥2 interior nodes). The implementation must be zero-cost for the <2-node case. + +### 2.6 Emission (`emit.jl`) — the only lossy step, memoized in the node table + +`node_point(arr, id)::P_out` realizes a node's output coordinate on demand and caches it: + +- Vertex nodes: the input vertex, bit-exact pass-through. +- Planar crossing nodes: **certified double-double fast path** (S3-proven): TwoSum on endpoint + differences → compensated 2×2 determinants → dd division → dd recombination; the result is + *accepted as correctly rounded* iff `|residual| + dd_error_bound < ½·ulp` of the candidate + (with the determinant-conditioning term, so near-parallel pairs fail the certificate). + Fallback: `_exact_crossing_point` (`Rational{BigInt}`) rounded to Float64. S3: 64,982/64,982 + real crossings certified, 0 disagreements with the rational answer, 0.024 µs vs 6.44 µs. +- Spherical crossing nodes: dd-certified `na×nb` direction, normalized; lon/lat conversion via + Float64 `atan2` (documented: the conversion itself is not certified — worst observed + deviation 1.4e-14° ≈ 1.5 nm; full certification through trig is not worth it because **no + decision ever consumes emitted coordinates**). Fallback on certificate failure: + `_sph_crossing_dir` exact direction → high-precision normalize → lon/lat (the + `CrossingEdgeSplit` precedent, `src/transformations/correction/crossing_edge_split.jl`). + +Substrate refactor bundled here (S2 finding g): make `_exact_crossing_point` / +`_sph_crossing_dir` the single shared exact-crossing authority with one call shape, consumed +by (a) this fallback, (b) §2.5's exact ordering keys, (c) `rk_nodes_coincide` — today their +call shapes differ slightly. Grow them in place in the kernel files; no wrappers. + +### 2.7 Substrate fix — material-interior authority (for phase 2, landed in phase 1) + +`_ring_interior_on_left(m, pts, is_hole; exact)` answers with the **denoted-region** side — +for a hole, the cavity. Overlay labeling needs the polygon's **material-interior** side (the +flip for holes). Add the sibling next to the authority (in `relate_geometry.jl` / +`kernel_spherical.jl`, wherever `_ring_interior_on_left` lives): + +```julia +_ring_material_interior_on_left(m, pts, is_hole; exact) = + is_hole ? !_ring_interior_on_left(m, pts, is_hole; exact) + : _ring_interior_on_left(m, pts, is_hole; exact) +``` + +so the flip exists exactly once and relate/overlay/extents agree by construction. Phase 2 +derives JTS's per-edge `depth_delta::Int8` from it (`material_interior_on_left ? -1 : +1`, +matching JTS `Edge.locationLeft/Right`: positive delta ⇒ Left=EXTERIOR, Right=INTERIOR). +Getting this wrong produced a wrong-area hole intersection instantly in spike S2 — write the +regression test from that case. + +### 2.8 Explicitly out of scope for phase 1 + +Half-edge graph, `OverlayLabel`, merger, builders (phase 2). Any op, export, or docs page +(phase 3). Snapping/precision machinery (never). Self-noding invalid inputs (contract, §2.2). +Threading. Prepared-overlay types (optional prebuilt-tree *argument* only). + +### 2.9 Tests and gates + +New test files under `test/methods/clipping/overlayng/` (one file per julia process when run +locally; they will be wired into the runtests tree): + +- **Invariant tests** for §2.1 items 1–5 on constructed cases: two crossing quads; the + degree-6 node (three segments through one exact point — Tier 2 must merge and the assertion + in §2.5 must hold); crossing-exactly-on-a-third-string's-vertex; collinear shared boundary + (GADM-style vertex-identical border → zero phantom crossings, zero interior splits); a–b–a + spike input. +- **Ordering cross-check**: dense synthetic (comb/zigzag, ~10³ crossings on few segments): + float-filtered order == exact-always order, elementwise; both manifolds. +- **Emission certificate audit**: for every crossing node in the test corpus, certified + result == rational result (asserted, not sampled); spherical direction deviation bound + checked against the exact direction. +- **Rounded-arrangement audit** (S3's census, small): emit all nodes for a small NE-110m + subset (a few shifted-self pairs), re-classify **crossing-incident** edges on the rounded + coordinates (the S3 invariant: only crossing nodes move, so only their incident edges need + auditing), assert zero introduced proper crossings. +- **Classification census invariants** on the NE subset: every touch/collinear carries ≥1 + vertex flag; planar and spherical produce the identical proper/touch/collinear multiset on + identical data (S1 observed this; encode it). +- Slow/optional (env-gated like the existing realdata benchmarks, not CI): the full NE-110m + sweep and one GADM pair. + +Performance budgets (regression bars from S1/S3, generous ×1.5 for production overhead — +record in the PR, enforce by judgment not CI): NE110 planar sweep noding ~1 ms; GADM CAN×USA +~82 ms planar / ~674 ms spherical; synthetic 250×250 ~19 ms; planar emission ≤0.05 µs/node +fast path. + +### 2.10 Conventions for implementers (repo-specific, learned the hard way) + +- Floats filter, exact predicates decide; every filter must have a written error-bound + justification in a comment at the constant's definition (see S4's landed style in + `kernel_spherical.jl`). +- Engine types stay type-erased over input geometry types (parameterize on kernel point `P` + only — the Julia 1.12 first-call compile blowup was caused by geometry-typed engine + internals; do not reintroduce). +- Match JTS names where a counterpart exists; do not invent jargon. Comments state + constraints, not narration. +- No adapters: if two pieces don't fit, grow the existing abstraction at its source (the §2.6 + and §2.7 refactors are the sanctioned examples). +- Docstring `@ref` targets exported names only. Never `git add -A`. One test file per julia + process locally; redirect julia output to log files. + +## 3. Phase 2 outline (engine core) — for context; port faithfully from JTS with these amendments + +JTS sources: `edgegraph/HalfEdge.java` (~500 lines, port required), then the overlayng package +(`Edge`/`EdgeMerger`/`EdgeKey` → `OverlayGraph`/`OverlayEdge` → `OverlayLabel` (plain fields, +NOT bit-packed) → `OverlayLabeller` (5 passes) → `MaximalEdgeRing`/`OverlayEdgeRing`/ +`PolygonBuilder` → `LineBuilder`/`IntersectionPointBuilder` → point/mixed paths). Skip: +`ElevationModel`, `FastOverlayFilter` (dead code), strict-mode branches, `LineBuilder`'s +merged path, `RingClipper`/`LineLimiter` (replaced by construct-free whole-ring extent pruning ++ one PIP per pruned ring). + +Spike-S2-validated amendments (each was tested end-to-end): + +1. Half-edge direction points are the parent segment's **far endpoint** (an original vertex); + all angular sorting goes through `rk_compare_edge_dir` with a `NodeKey` apex. At a Tier-2 + coincidence-merged node, directions foreign to the representative key take the kernel's + exact-rational slow path — it exists and is correct; do not special-case. +2. `EdgeMerger` keys on the **unordered node-id pair** (uniqueness of the segment/minor-arc + between two nodes; antipodal ingest guard covers the sphere) — never on coordinates, and on + node ids (post-Tier-2), never raw `NodeKey`s. +3. `depth_delta` from `_ring_material_interior_on_left` (§2.7), so `oriented` semantics flow + through with no further plumbing. +4. Porting trap: `linkResultAreaMaxRingAtNode` must be called **ungated** for every in-result + edge (JTS's unfulfilled `// TODO: skip already-linked` is deliberate); gating empties + degree-2 nodes. +5. `PolygonBuilder` containment via `rk_point_in_ring` / the indexed locators — never planar + even-odd on emitted coordinates. +6. Spherical empty-vs-full: a boundaryless nonempty result is disambiguated by locating one + input vertex under the op's semantics (topology requirement of the closed manifold, not a + robustness measure). +7. Labeling's disconnected-edge PIP runs against original inputs via the existing indexed + locators; symbolic endpoints locate exactly (rare path). + +Known future workload for phase 2/3 from the spikes: full minimal-ring split + +hole-nesting (`buildMinimalRings`) — the S3 prototype's union imprecision on 355-island +France and antimeridian Russia came from faking exactly this. + +## 4. Phase 3 outline (ops + API) + +`isResultOfOp` + driver; the four ops incl. new `symdifference`; mixed-dimension results per +`OverlayUtil` rules; `OverlayNG{M} <: Algorithm{M}` **opt-in** (`intersection(OverlayNG(...), +a, b)`); Foster–Hormann remains the default until differential validation completes. +Validation: JTS overlay JUnit/XML suites; differential fuzz vs LibGEOS (OverlayNG is GEOS's +default engine ≥3.9); spherical differential vs compiled s2 `S2BooleanOperation` (the +Sudan-experiment C++ harness pattern); NE/GADM sweeps with area conservation +(`area(A∪B)+area(A∩B) ≈ area(A)+area(B)` — the signed-area fix `69e416484` makes this a +machine-precision gate). PrecompileTools workload extension; `benchmarks/` additions. + +## Appendix — spike verdicts (2026-07-15/16, condensed) + +Spike code lives in the session scratchpad +(`/private/tmp/claude-501/-Users-anshul-temp-GO-jts/36b4dc87-418b-4e8c-8a1f-b2b36a1ca475/scratchpad/spikes/{s1_noder,s2_skeleton,s3_census,s4_prefilter,s5_area}/`) +— **ephemeral**: consult if still present, but this document stands alone. + +- **S1 (noder throughput)**: symbolic noder on the relateng substrate ran at 0.03×–0.17× of + LibGEOS's full intersection op (GADM CAN×USA 5.9M edges: 82 ms vs 1.41 s). Zero exact + fallbacks on real data in ordering and dedup; float ordering filter mandatory in dense + regimes (240–590×). Touch/collinear = input vertex: zero counterexamples in ~150k. GADM + shared borders are vertex-identical → zero phantom crossings. Emission was the per-node cost + center (6/18 µs rational) → S3 fixed. Extraction dominates sparse pairs → keep ingest + separable. +- **S2 (walking skeleton)**: all structural claims validated; planar results area- and + vertex-exact vs LibGEOS; crossing emission bit-exact on an integer grid; degree-6 and + spherical tangent angular ordering correct. Produced amendments §3.1–5 and the + material-interior fix (§2.7). +- **S3 (rounding census)**: ~930k noded edges audited across NE/GADM/graticule/synthetic/ + spherical: **0** rounding-introduced crossings, **0** float collisions, order inversions only + on antimeridian-wrapping Russia; area conservation ≤8e-16 (valid outputs); worst LibGEOS + differential 1.7e-9 on a 1e-4 sliver (their snap-rounding). Only artifacts sat inside + already-invalid input (NE110 Sudan), where our result matched makeValid-then-overlay while + raw LibGEOS threw. Certified-dd emission: 100.0000% certified on 64,982 crossings, 0 + soundness violations, 273×; spherical float direction 704× at ≤1.4e-14°. +- **S4 (kernel prefilter — landed `8d938e832`)**: spherical exact classification reduced to + four `rk_orient` signs (proper-crossing branch, exact path only) + certified span triage + + repeated-vertex orient short-circuit; 1.65M-pair audit, 0 disagreements; spherical classify + now ≈0.08 µs/pair (below planar); GADM CAN×USA spherical collect 3.12 s → 674 ms. +- **S5 (signed area — landed `69e416484`)**: `area(Spherical())` Eriksson fan restored to the + native signed Van Oosterom–Strackee `atan2` form; concave overcounts (Norway +174%, Chile + +432%) fixed; small-polygon accuracy bit-identical where denom>0; latent denom<0 misbranch + confirmed and fixed. Relevant here because it makes area conservation usable as a + machine-precision overlay gate. diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index b9ae685c3e..f55a042f1e 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -123,6 +123,14 @@ include("methods/geom_relations/relateng/topology_computer.jl") include("methods/geom_relations/relateng/edge_intersector.jl") # The RelateNG engine: drives all of the above through the phased evaluation. include("methods/geom_relations/relateng/relate_ng.jl") + +# OverlayNG noding substrate (phase 1): geometries → exactly-noded arrangement. +include("methods/clipping/overlayng/noding/noded_arrangement.jl") +include("methods/clipping/overlayng/noding/collect.jl") +include("methods/clipping/overlayng/noding/node_identity.jl") +include("methods/clipping/overlayng/noding/split.jl") +include("methods/clipping/overlayng/noding/emit.jl") + include("methods/orientation.jl") include("methods/polygonize.jl") include("methods/minimum_bounding_circle.jl") diff --git a/src/methods/clipping/overlayng/noding/collect.jl b/src/methods/clipping/overlayng/noding/collect.jl new file mode 100644 index 0000000000..07bdbe3675 --- /dev/null +++ b/src/methods/clipping/overlayng/noding/collect.jl @@ -0,0 +1,75 @@ +# # Stage 1 — collect (design §2.3) +# +# Enumerate candidate segment pairs (A × B) through the reused RelateNG edge +# index and classify each with the exact kernel (`rk_classify_intersection`), +# recording symbolic node keys onto the parent segments. No intersection +# coordinate is ever constructed here. +# +# `seg_nodes` accumulates, per parent segment `(string_idx, seg_idx)`, the ids of +# nodes lying strictly in that segment's interior. `string_idx` is global into +# the arrangement's `segstrings`: A strings occupy `1:na`, B strings `na+1:end`. + +@inline function _record_interior!(seg_nodes, string_idx::Int32, seg_idx::Int32, nid::Int32) + push!(get!(() -> Int32[], seg_nodes, (string_idx, seg_idx)), nid) + return nothing +end + +function _collect_crossings!(m::Manifold, table::NodeTable{P}, seg_nodes, + ssa::AbstractVector{RelateSegmentString{P}}, + ssb::AbstractVector{RelateSegmentString{P}}, na::Int32; + exact = True(), tree_a = nothing, tree_b = nothing) where {P} + ta = tree_a === nothing ? _relate_edge_index(m, ssa) : tree_a + tb = tree_b === nothing ? _relate_edge_index(m, ssb) : tree_b + (ta === nothing || tb === nothing) && return nothing + SpatialTreeInterface.dual_depth_first_search(Extents.intersects, ta, tb) do ia, ib + (sa, ka) = ta.data[ia] + (sb, kb) = tb.data[ib] + _classify_pair!(m, table, seg_nodes, ssa, ssb, na, sa, ka, sb, kb; exact) + return nothing + end + return nothing +end + +# Function barrier: statically-typed classification of one candidate pair +# (the do-block above is a dynamic closure over the tree traversal). +function _classify_pair!(m::Manifold, table::NodeTable{P}, seg_nodes, + ssa, ssb, na::Int32, sa::Int, ka::Int, sb::Int, kb::Int; exact) where {P} + a0 = ssa[sa].pts[ka]; a1 = ssa[sa].pts[ka + 1] + b0 = ssb[sb].pts[kb]; b1 = ssb[sb].pts[kb + 1] + gsa = Int32(sa) # A string global index + gsb = na + Int32(sb) # B string global index + ksa = Int32(ka); ksb = Int32(kb) + + cls = rk_classify_intersection(m, a0, a1, b0, b1; exact) + kind = cls.kind + if kind == SS_DISJOINT + return nothing + elseif kind == SS_PROPER + #-- a proper crossing is strictly interior to BOTH segments + nid = _intern_node!(table, crossing_node(a0, a1, b0, b1)) + _record_interior!(seg_nodes, gsa, ksa, nid) + _record_interior!(seg_nodes, gsb, ksb, nid) + return nothing + end + #-- SS_TOUCH / SS_COLLINEAR: every intersection point is an input vertex, + #-- reported via the incidence flags (design §2.3, S1 census). The claim is + #-- load-bearing; assert it. + @assert (cls.a0_on_b || cls.a1_on_b || cls.b0_on_a || cls.b1_on_a) """ + touch/collinear intersection with no vertex incidence flag — \ + the touch=input-vertex invariant (design §2.3) is violated""" + #-- a?/b? are vertices of A/B; record on the OTHER segment only where the + #-- vertex lies strictly in that segment's interior (not at its endpoints). + if cls.a0_on_b && a0 != b0 && a0 != b1 + _record_interior!(seg_nodes, gsb, ksb, _intern_node!(table, vertex_node(a0))) + end + if cls.a1_on_b && a1 != b0 && a1 != b1 + _record_interior!(seg_nodes, gsb, ksb, _intern_node!(table, vertex_node(a1))) + end + if cls.b0_on_a && b0 != a0 && b0 != a1 + _record_interior!(seg_nodes, gsa, ksa, _intern_node!(table, vertex_node(b0))) + end + if cls.b1_on_a && b1 != a0 && b1 != a1 + _record_interior!(seg_nodes, gsa, ksa, _intern_node!(table, vertex_node(b1))) + end + return nothing +end diff --git a/src/methods/clipping/overlayng/noding/emit.jl b/src/methods/clipping/overlayng/noding/emit.jl new file mode 100644 index 0000000000..3a1b0eb579 --- /dev/null +++ b/src/methods/clipping/overlayng/noding/emit.jl @@ -0,0 +1,161 @@ +# # Emission (design §2.6) — the only lossy step +# +# `node_point(arr, id)` realizes a node's OUTPUT coordinate on demand and caches +# it. Nothing in the arrangement's decisions ever consumed a constructed +# coordinate; this is where the exact symbolic result is rounded to Float64 for +# output. Both manifolds emit `(x, y)` / `(lon, lat)` as `Tuple{Float64,Float64}`. +# +# - Vertex nodes: the input vertex, bit-exact pass-through. +# - Planar crossings: a **certified double-double** fast path (spike S3, 100% +# certified on 64,982 real crossings, 0 disagreements with the rational answer, +# 273×) — TwoSum on endpoint differences, compensated 2×2 determinants, +# dd division, dd recombination; the coordinate is accepted iff its residual +# plus the dd error bound is below ½ ulp, with the determinant-conditioning +# term so near-parallel pairs fail the certificate. Fallback: the exact +# `Rational{BigInt}` crossing point, rounded. +# - Spherical crossings: the Float64 crossing direction `±(na×nb)`, accepted when +# the arcs clear a near-tangency conditioning gate (spike S3 measured the float +# direction at ≤1.4e-14° ≈ 1.5 nm; the lon/lat trig itself is uncertified by +# design — no decision ever consumes an emitted coordinate). Fallback: the +# exact `_sph_crossing_dir`, normalized and converted. + +# ## Error-free transforms and double-double primitives (spike S3, productionized) + +@inline function _twosum(a::Float64, b::Float64) + s = a + b; bb = s - a + return (s, (a - (s - bb)) + (b - bb)) +end +@inline function _twoproduct(a::Float64, b::Float64) + p = a * b + return (p, fma(a, b, -p)) +end +@inline _diff_dd(a::Float64, b::Float64) = _twosum(a, -b) # (hi, lo) == a - b exactly + +@inline function _ddmul(ah, al, bh, bl) # (ah+al)*(bh+bl) + (ph, pl) = _twoproduct(ah, bh) + return _twosum(ph, pl + (ah * bl + al * bh)) +end +@inline function _ddsub(ah, al, bh, bl) # (ah+al) - (bh+bl) + (sh, se) = _twosum(ah, -bh) + return _twosum(sh, (se + al) - bl) +end +@inline function _ddadd(ah, al, bh, bl) # (ah+al) + (bh+bl) + (sh, se) = _twosum(ah, bh) + return _twosum(sh, (se + al) + bl) +end +# det = a*d - b*c, each operand a double-double 2-tuple +@inline function _det2_ddfull(a, b, c, d) + (mh, ml) = _ddmul(a[1], a[2], d[1], d[2]) + (nh, nl) = _ddmul(b[1], b[2], c[1], c[2]) + return _ddsub(mh, ml, nh, nl) +end +# dd / dd -> dd (Dekker) +@inline function _div_dd(ah, al, bh, bl) + q1 = ah / bh + (ph, pl) = _twoproduct(q1, bh) + (sh, sl) = _twosum(ah, -ph) + r = ((sh - pl) + sl) + al - q1 * bl + return _twosum(q1, r / bh) +end + +# Certified correctly-rounded emit of one coordinate: `xf = fl(hi+lo)`, exact +# residual `rem` from TwoSum, accept iff `|rem| + dderr < ½ ulp(xf)`. +@inline function _certify_coord(hi, lo, dderr) + (xf, rem) = _twosum(hi, lo) + return (xf, abs(rem) + dderr < 0.5 * eps(xf)) +end + +# Fast certified planar crossing of (a0,a1) × (b0,b1). Returns (x, y, certified). +function _certified_crossing(a0, a1, b0, b1) + ax0, ay0 = Float64(GI.x(a0)), Float64(GI.y(a0)) + ax1, ay1 = Float64(GI.x(a1)), Float64(GI.y(a1)) + bx0, by0 = Float64(GI.x(b0)), Float64(GI.y(b0)) + bx1, by1 = Float64(GI.x(b1)), Float64(GI.y(b1)) + #-- exact endpoint differences as double-doubles + da_x = _diff_dd(ax1, ax0); da_y = _diff_dd(ay1, ay0) + db_x = _diff_dd(bx1, bx0); db_y = _diff_dd(by1, by0) + c0_x = _diff_dd(bx0, ax0); c0_y = _diff_dd(by0, ay0) + (denh, denl) = _det2_ddfull(da_x, da_y, db_x, db_y) # da × db + (tnh, tnl) = _det2_ddfull(c0_x, c0_y, db_x, db_y) # (b0-a0) × db + (th, tl) = _div_dd(tnh, tnl, denh, denl) + (txh, txl) = _ddmul(th, tl, da_x[1], da_x[2]); (xh, xl) = _ddadd(ax0, 0.0, txh, txl) + (tyh, tyl) = _ddmul(th, tl, da_y[1], da_y[2]); (yh, yl) = _ddadd(ay0, 0.0, tyh, tyl) + #-- dd error bounds amplified by determinant conditioning: near-parallel ⇒ + #-- small |denom| ⇒ large condK ⇒ certificate correctly fails (spike S3) + u2 = eps(Float64)^2 + scale = abs(da_x[1]) + abs(da_y[1]) + abs(db_x[1]) + abs(db_y[1]) + condK = (abs(da_x[1] * db_y[1]) + abs(da_y[1] * db_x[1])) / max(abs(denh), floatmin(Float64)) + tmag = abs(th) + ex = 64 * u2 * (abs(xh) + tmag * abs(da_x[1]) * condK + scale) + ey = 64 * u2 * (abs(yh) + tmag * abs(da_y[1]) * condK + scale) + (xf, cx) = _certify_coord(xh, xl, ex) + (yf, cy) = _certify_coord(yh, yl, ey) + return (xf, yf, cx & cy) +end + +# ## Per-manifold node coordinate realization (dispatched on the kernel point `P`) + +# Planar: vertex pass-through; crossing via the certified dd path with a rational +# fallback (identical to the fallback the S3 audit compared against). +function _emit_node_coord(k::NodeKey{Tuple{Float64, Float64}}) + k.is_crossing || return k.pt + (x, y, cert) = _certified_crossing(k.pt, k.a1, k.b0, k.b1) + cert && return (x, y) + rx, ry = _exact_crossing_point(k) + return (Float64(rx), Float64(ry)) +end + +# Near-tangency gate for the spherical float direction: |na×nb|² ≥ tol²·|na|²·|nb|² +# means the arcs' planes meet at ≥ ~1e-9 rad, so the float direction's relative +# error (≈ eps / sin θ) is bounded well below the ≤1.4e-14° the design accepts. +# Below the gate the crossing is near-tangent and falls to the exact direction. +const _SPH_TANGENT_GATE = 1e-9 + +function _emit_node_coord(k::NodeKey{<:UnitSphericalPoint}) + k.is_crossing || return _usp_to_lonlat(k.pt) + #-- float na, nb, d = na × nb, with the conditioning gate + A0 = _vec3(False(), k.pt); A1 = _vec3(False(), k.a1) + B0 = _vec3(False(), k.b0); B1 = _vec3(False(), k.b1) + na = _cross3(A0, A1); nb = _cross3(B0, B1) + d = _cross3(na, nb) + d2 = _dot3(d, d); na2 = _dot3(na, na); nb2 = _dot3(nb, nb) + if d2 >= _SPH_TANGENT_GATE^2 * na2 * nb2 + dir = _sph_crossing_dir(False(), k) # picks the interior candidate + return _dir_to_lonlat(dir) + end + #-- near-tangent fallback: exact direction (Rational), then normalize + convert + return _dir_to_lonlat(_sph_crossing_dir(True(), k)) +end + +@inline function _dir_to_lonlat(d) + s = sqrt(Float64(d[1])^2 + Float64(d[2])^2 + Float64(d[3])^2) + return _usp_to_lonlat(UnitSphericalPoint(Float64(d[1]) / s, Float64(d[2]) / s, Float64(d[3]) / s)) +end + +@inline function _usp_to_lonlat(u) + ll = GeographicFromUnitSphere()(u) + return (Float64(ll[1]), Float64(ll[2])) +end + +# ## Public-internal accessor + +""" + node_point(arr::NodedArrangement, id) -> Tuple{Float64,Float64} + +The realized output coordinate of node `id` (planar `(x, y)` / spherical +`(lon, lat)`), memoized in the node table (design §2.6). The only place a +constructed coordinate enters the substrate. +""" +function node_point(arr::NodedArrangement, id::Integer) + t = arr.nodes + i = Int(id) + @inbounds t.realized[i] && return t.coords[i] + c = _emit_node_coord(t.keys[i]) + @inbounds t.coords[i] = c + @inbounds t.realized[i] = true + return c +end + +# The two output coordinates of a noded edge (convenience for callers/tests). +edge_endpoints(arr::NodedArrangement, e::NodedEdge) = + (node_point(arr, e.node_lo), node_point(arr, e.node_hi)) diff --git a/src/methods/clipping/overlayng/noding/node_identity.jl b/src/methods/clipping/overlayng/noding/node_identity.jl new file mode 100644 index 0000000000..656d7aa2d1 --- /dev/null +++ b/src/methods/clipping/overlayng/noding/node_identity.jl @@ -0,0 +1,161 @@ +# # Stage 3 — node identity (design §2.4) +# +# Two tiers, no canonical key. Tier 1 (egal `NodeKey` equality) already ran +# during collect via `_intern_node!`. Tier 2 finds *geometric* coincidence across +# distinct keys — two segment pairs crossing at one point, or a crossing landing +# on a third string's vertex — which egal equality misses. A throwaway Float64 +# proximity sweep flags candidate pairs; the exact `rk_nodes_coincide` confirms; +# a union-find over confirmed pairs (bounded by candidate count — zero on all +# real data measured, a handful on constructed degree-≥3 nodes) merges them. +# +# The approximate positions come from the same kernel filter code as the emission +# fast path (`_approx_node_point` / `_exact_node_dir(False(), …)`), minus the +# certification — design §2.4/§2.6, shared code. + +@inline function _uf_find(parent::Vector{Int32}, i::Int32) + while parent[i] != i + parent[i] = parent[parent[i]] # path halving + i = parent[i] + end + return i +end + +@inline function _uf_union!(parent::Vector{Int32}, i::Int32, j::Int32) + ri = _uf_find(parent, i); rj = _uf_find(parent, j) + ri != rj && (parent[ri] = rj) + return nothing +end + +function _merge_coincident_nodes!(m::Manifold, table::NodeTable{P}, seg_nodes; exact) where {P} + n = num_nodes(table) + if n < 2 + #-- no cross-key coincidence possible, but one node is routinely reported + #-- by several candidate pairs (S1; e.g. an a–b–a spike), so the interior + #-- lists must still be deduped before ordering/splitting (invariant 3) + for (seg, nids) in seg_nodes + seg_nodes[seg] = unique(nids) + end + return nothing + end + parent = collect(Int32(1):Int32(n)) + #-- tier-2 confirmation sweep (manifold-specific proximity geometry) + nmerges = _coincidence_sweep!(m, table, parent; exact) + if nmerges == 0 + #-- the overwhelmingly common case (S1: zero coincidences on real data): + #-- no ids collapse, so skip the O(n) key-rehash compaction entirely. + #-- Only the interior lists need deduping (a node reported by several pairs); + #-- the coordinate cache is sized later by `_ensure_coord_cache!`. + for (seg, nids) in seg_nodes + seg_nodes[seg] = unique(nids) + end + return nothing + end + + #-- compact: provisional id -> final id, representative = first-seen member + remap = Vector{Int32}(undef, n) + final_keys = NodeKey{P}[] + root_final = Dict{Int32, Int32}() + for i in 1:n + r = _uf_find(parent, Int32(i)) + fid = get(root_final, r, Int32(0)) + if fid == 0 + push!(final_keys, table.keys[i]) + fid = Int32(length(final_keys)) + root_final[r] = fid + end + remap[i] = fid + end + + #-- re-point the interner so later endpoint interning resolves to final ids, + #-- then swap in the compacted key table and a fresh coordinate cache + for (k, oldid) in table.ids + table.ids[k] = remap[oldid] + end + table.keys = final_keys + #-- coords/realized are sized once after splitting (`_ensure_coord_cache!`) + + #-- rewrite the interior node lists through the remap, deduping ids a merge + #-- collapsed together (order is re-established in `split.jl`) + for (seg, nids) in seg_nodes + seen = Set{Int32}() + out = Int32[] + for id in nids + fid = remap[id] + (fid in seen) && continue + push!(seen, fid); push!(out, fid) + end + seg_nodes[seg] = out + end + return nothing +end + +# Planar proximity sweep: a single x-sort, forward scan while the x-gap is within +# the summed proximity radii, exact confirm, union. A true coincidence has +# near-identical float positions (≪ the crossing radius), so the sweep never +# misses one; over-broad radii only cost extra exact confirms. +function _coincidence_sweep!(m::Planar, table::NodeTable{P}, parent; exact) where {P} + n = num_nodes(table) + xs = Vector{Float64}(undef, n) + ys = Vector{Float64}(undef, n) + rad = Vector{Float64}(undef, n) + for i in 1:n + k = table.keys[i] + x, y, _ = _approx_node_point(k) + xs[i] = x; ys[i] = y + #-- crossings can land off their float approximation; vertices are exact. + #-- 1e-8·|coord| is a generous proximity gate (the exact test confirms) + rad[i] = k.is_crossing ? 1e-8 * max(1.0, abs(x), abs(y)) : 0.0 + end + order = sortperm(xs) + nmerges = 0 + @inbounds for a in 1:n + i = order[a] + for b in (a + 1):n + j = order[b] + rr = rad[i] + rad[j] + xs[j] - xs[i] > rr && break + abs(ys[i] - ys[j]) > rr && continue + _uf_find(parent, Int32(i)) == _uf_find(parent, Int32(j)) && continue + if rk_nodes_coincide(m, table.keys[i], table.keys[j]; exact) + _uf_union!(parent, Int32(i), Int32(j)); nmerges += 1 + end + end + end + return nmerges +end + +# Spherical proximity sweep: same shape over the (float, normalized) crossing +# directions, gated by chordal proximity along the first coordinate. +function _coincidence_sweep!(m::Spherical, table::NodeTable{P}, parent; exact) where {P} + n = num_nodes(table) + dirs = Vector{NTuple{3, Float64}}(undef, n) + for i in 1:n + d = _exact_node_dir(False(), table.keys[i]) + s = sqrt(d[1]^2 + d[2]^2 + d[3]^2) + dirs[i] = (d[1] / s, d[2] / s, d[3] / s) + end + xs = [d[1] for d in dirs] + #-- 1e-11 chord ≈ 6e-4 m on the unit sphere: generous vs the ≤1.4e-14° float + #-- direction error, so no true coincidence is missed + rr = 1e-11 + order = sortperm(xs) + nmerges = 0 + @inbounds for a in 1:n + i = order[a] + for b in (a + 1):n + j = order[b] + xs[j] - xs[i] > rr && break + di = dirs[i]; dj = dirs[j] + dotp = di[1] * dj[1] + di[2] * dj[2] + di[3] * dj[3] + cx = di[2] * dj[3] - di[3] * dj[2] + cy = di[3] * dj[1] - di[1] * dj[3] + cz = di[1] * dj[2] - di[2] * dj[1] + (dotp > 0 && sqrt(cx^2 + cy^2 + cz^2) <= rr) || continue + _uf_find(parent, Int32(i)) == _uf_find(parent, Int32(j)) && continue + if rk_nodes_coincide(m, table.keys[i], table.keys[j]; exact) + _uf_union!(parent, Int32(i), Int32(j)); nmerges += 1 + end + end + end + return nmerges +end diff --git a/src/methods/clipping/overlayng/noding/noded_arrangement.jl b/src/methods/clipping/overlayng/noding/noded_arrangement.jl new file mode 100644 index 0000000000..9582d49292 --- /dev/null +++ b/src/methods/clipping/overlayng/noding/noded_arrangement.jl @@ -0,0 +1,151 @@ +# # OverlayNG noding substrate — the `NodedArrangement` +# +# Phase 1 of the OverlayNG port (design doc `2026-07-16-overlayng-noding-substrate.md`). +# Nodes two input geometries into a shared, exactly-noded edge arrangement with +# **symbolic** node identity: node identity, ordering, and coincidence are decided +# by exact kernel predicates over the input vertices and the symbolic crossing +# keys (`NodeKey`); Float64 appears only as certified filters and at +# emission. There is no snapping, no tolerance in any decision (design §0). +# +# The pipeline is four stages over the reused RelateNG substrate: +# 1. collect (`collect.jl`) — candidate enumeration + exact classification +# 2. identity (`node_identity.jl`) — two-tier node grouping +# 3. order+split (`split.jl`) — along-segment ordering, noded-edge emission +# 4. emit (`emit.jl`) — the sole lossy step, coordinate realization +# +# Everything here is internal to GeometryOps — nothing is exported. + +# A noded sub-segment: a piece of one parent segment between two nodes. It +# carries no geometry (design §2.1 invariant 5) — its shape is a lookup into its +# parent segment, and all source metadata (owner, ring id, dimension) is reached +# through `string_idx` into the arrangement's `segstrings`, so nothing can +# desynchronize. +struct NodedEdge + string_idx :: Int32 # index into `NodedArrangement.segstrings` + seg_idx :: Int32 # segment within that parent string (`pts[seg_idx:seg_idx+1]`) + node_lo :: Int32 # node id at the sub-segment start (in the parent's traversal order) + node_hi :: Int32 # node id at the sub-segment end +end + +#= +The node identity table (design §2.4). `ids` is the tier-1 egal interner +(`NodeKey` bit-equality is canonical — kernel points normalize signed zeros), +mapping every known key to its node id; after tier-2 merging it maps every +provisional key to its *final* id, so keys interned later (segment endpoints in +`split.jl`) resolve to the merged id. `keys[id]` is the group representative. +`coords`/`realized` memoize emitted output coordinates (design §2.6), realized +lazily by `node_point` and grown as endpoint nodes are interned. + +Mutable so tier-2 merging can compact `keys` and re-point `ids` in place. +=# +mutable struct NodeTable{P} + ids :: Dict{NodeKey{P}, Int32} + keys :: Vector{NodeKey{P}} + coords :: Vector{Tuple{Float64, Float64}} + realized :: Vector{Bool} +end + +NodeTable{P}() where {P} = + NodeTable{P}(Dict{NodeKey{P}, Int32}(), NodeKey{P}[], Tuple{Float64, Float64}[], Bool[]) + +num_nodes(t::NodeTable) = length(t.keys) + +# Intern a key, returning its node id (tier-1 egal merge). The output-coordinate +# cache is an emission concern (`_ensure_coord_cache!`), not grown here. +function _intern_node!(t::NodeTable{P}, key::NodeKey{P}) where {P} + id = get(t.ids, key, Int32(0)) + id != 0 && return id + push!(t.keys, key) + id = Int32(length(t.keys)) + t.ids[key] = id + return id +end + +# Size the (lazily-realized) output-coordinate cache to the final node count. +# Called once after splitting, so noding itself never touches it. +function _ensure_coord_cache!(t::NodeTable) + n = length(t.keys) + t.coords = Vector{Tuple{Float64, Float64}}(undef, n) + t.realized = fill(false, n) + return nothing +end + +""" + NodedArrangement{P} + +The exactly-noded arrangement of two input geometries (design §2.1). `P` is the +manifold's kernel point type — exactly two instantiations, +`Tuple{Float64,Float64}` (planar) and `UnitSphericalPoint{Float64}` (spherical) — +so the engine is type-erased over the input geometry types. + +Fields: +- `segstrings`: the ingested inputs as `RelateSegmentString`s (A side + first, then B side); `NodedEdge.string_idx` indexes here. +- `nodes`: the symbolic node table (`NodeTable`). +- `seg_nodes`: per-parent-segment ordered interior node-id lists, keyed by + `(string_idx, seg_idx)`; absent for unsplit segments. +- `edges`: every noded sub-segment of every parent segment. + +Construct with `NodedArrangement(m, a, b)` (raw geometries) or +`NodedArrangement(m, ssa, ssb)` (pre-ingested segment strings). +""" +struct NodedArrangement{P} + segstrings :: Vector{RelateSegmentString{P}} + nodes :: NodeTable{P} + seg_nodes :: Dict{Tuple{Int32, Int32}, Vector{Int32}} + edges :: Vector{NodedEdge} +end + +num_nodes(arr::NodedArrangement) = num_nodes(arr.nodes) +num_edges(arr::NodedArrangement) = length(arr.edges) + +# Whether a segment string is a polygon hole. Derived (not stored — design §2.2): +# shells are `ring_id == 0`, holes `ring_id >= 1`, lines `ring_id == -1` +# (`_extract_ring_to_segment_string!` in relate_geometry.jl). +_ss_is_hole(ss::RelateSegmentString) = ss.ring_id > 0 + +# ## Ingest (design §2.2) +# +# Reuse `RelateGeometry` / `extract_segment_strings` unchanged: kernel-point +# conversion, repeated-point removal, and ring orientation all happen here, once. +# Kept separate from arrangement construction so a future prepared overlay can +# convert once and re-arrange many times (S1: extraction dominates sparse pairs). +_overlay_segstrings(m::Manifold, geom, is_a::Bool; exact = True()) = + extract_segment_strings(RelateGeometry(m, geom; exact), is_a, nothing) + +# ## Construction + +# From raw geometries: ingest each side, then arrange. +function NodedArrangement(m::Manifold, a, b; exact = True(), tree_a = nothing, tree_b = nothing) + ssa = _overlay_segstrings(m, a, true; exact) + ssb = _overlay_segstrings(m, b, false; exact) + return NodedArrangement(m, ssa, ssb; exact, tree_a, tree_b) +end + +# From pre-ingested segment strings. `tree_a`/`tree_b` accept caller-supplied +# prebuilt segment indices (a `PreparedRelate` carries exactly `_relate_edge_index` +# output) — an optional argument only, no new prepare type (design §2.3). +function NodedArrangement(m::Manifold, + ssa::AbstractVector{RelateSegmentString{P}}, + ssb::AbstractVector{RelateSegmentString{P}}; + exact = True(), tree_a = nothing, tree_b = nothing) where {P} + na = length(ssa) + segstrings = Vector{RelateSegmentString{P}}(undef, na + length(ssb)) + @inbounds for i in 1:na + segstrings[i] = ssa[i] + end + @inbounds for i in eachindex(ssb) + segstrings[na + i] = ssb[i] + end + + table = NodeTable{P}() + seg_nodes = Dict{Tuple{Int32, Int32}, Vector{Int32}}() + # stage 1 + _collect_crossings!(m, table, seg_nodes, ssa, ssb, Int32(na); exact, tree_a, tree_b) + # stage 3 + _merge_coincident_nodes!(m, table, seg_nodes; exact) + # stages 2 + 4 + edges = _split_edges!(m, table, seg_nodes, segstrings; exact) + _ensure_coord_cache!(table) + return NodedArrangement{P}(segstrings, table, seg_nodes, edges) +end diff --git a/src/methods/clipping/overlayng/noding/split.jl b/src/methods/clipping/overlayng/noding/split.jl new file mode 100644 index 0000000000..a3687eda26 --- /dev/null +++ b/src/methods/clipping/overlayng/noding/split.jl @@ -0,0 +1,64 @@ +# # Stages 2 + 4 — order and split (design §2.5, §2.1) +# +# For every parent segment of every string, order its interior nodes along the +# segment (stage 2, `rk_compare_along_segment`), cap the chain with the segment's +# endpoint vertex nodes, and emit one `NodedEdge` per non-degenerate link +# (stage 4). Interior node lists are already deduped by `node_identity.jl`; the +# consecutive-dedup here guards the endpoint↔interior joins. +# +# Ordering is zero-cost for the <2-node case (the overwhelming majority — S1: no +# real segment acquired ≥2 interior nodes), which is the only case that touches +# the kernel ordering predicate. + +function _split_edges!(m::Manifold, table::NodeTable{P}, seg_nodes, + segstrings::Vector{RelateSegmentString{P}}; exact) where {P} + edges = NodedEdge[] + chain = Int32[] + for gsi in eachindex(segstrings) + ss = segstrings[gsi] + pts = ss.pts + gsi32 = Int32(gsi) + for k in 1:(length(pts) - 1) + interior = get(seg_nodes, (gsi32, Int32(k)), nothing) + if interior !== nothing && length(interior) >= 2 + _order_along_segment!(m, interior, pts[k], pts[k + 1], table; exact) + end + #-- build the node chain: start vertex, ordered interior, end vertex + lo = _intern_node!(table, vertex_node(pts[k])) + hi = _intern_node!(table, vertex_node(pts[k + 1])) + empty!(chain) + push!(chain, lo) + if interior !== nothing + append!(chain, interior) + end + push!(chain, hi) + #-- emit one edge per link, skipping any zero-length (coincident) join + for c in 1:(length(chain) - 1) + a = chain[c]; b = chain[c + 1] + a == b && continue + push!(edges, NodedEdge(gsi32, Int32(k), a, b)) + end + end + end + return edges +end + +# Sort a segment's interior node ids along the oriented segment (s0, s1). The +# kernel predicate is the sole authority; by construction (stage 3 ran first) no +# two distinct ids coincide, so a `0` comparison is an invariant violation. +function _order_along_segment!(m::Manifold, ids::Vector{Int32}, s0, s1, + table::NodeTable; exact) + sort!(ids; lt = (a, b) -> + rk_compare_along_segment(m, s0, s1, table.keys[a], table.keys[b]; exact) < 0) + #-- assert no residual coincidence among ordered neighbours (design §2.5) + @assert begin + ok = true + for c in 2:length(ids) + if rk_compare_along_segment(m, s0, s1, table.keys[ids[c - 1]], table.keys[ids[c]]; exact) == 0 + ok = false; break + end + end + ok + end "two distinct node ids coincide along a segment after node-identity merging (design §2.5)" + return nothing +end diff --git a/src/methods/geom_relations/relateng/kernel.jl b/src/methods/geom_relations/relateng/kernel.jl index 53e651d178..dfd2cea85f 100644 --- a/src/methods/geom_relations/relateng/kernel.jl +++ b/src/methods/geom_relations/relateng/kernel.jl @@ -105,6 +105,19 @@ constructed apex coordinate; foreign directions (incident edges of a D3 coincidence-merged node, from other segment pairs crossing at the same point) are compared exactly around the rational apex (slow path). + rk_compare_along_segment(m, s0, s1, na, nb; exact)::Int + +Order of two nodes `na`, `nb` (`NodeKey`s, both lying on the oriented segment +`(s0, s1)`) along that segment: negative / zero / positive as `na` precedes / +coincides with / follows `nb` in the direction `s0 → s1`. Node ordering for +the noding substrate (OverlayNG phase 1, design §2.5). A certified Float64 +filter decides the pair when the approximate along-segment gap exceeds its +error bound; only otherwise is the exact key computed — planar via +`Rational{BigInt}` parameters from `_exact_crossing_point`, spherical via +triple-product signs against `_sph_crossing_dir` directions. Returning `0` +means the nodes coincide, which cannot happen for distinct node ids after +node-identity merging (design §2.4) has run — callers assert this. + rk_crossing_dirs_ccw(m, a0, a1, b0, b1; exact) CCW cyclic order of the four half-edge directions incident to the proper diff --git a/src/methods/geom_relations/relateng/kernel_planar.jl b/src/methods/geom_relations/relateng/kernel_planar.jl index 72d17086ba..9ba38cbf0c 100644 --- a/src/methods/geom_relations/relateng/kernel_planar.jl +++ b/src/methods/geom_relations/relateng/kernel_planar.jl @@ -250,7 +250,7 @@ function rk_compare_edge_dir(m::Planar, node::NodeKey, p, q; exact) around the exact rational apex instead (slow path; only reachable on the rare self-noding merge path). =# - apex = _exact_crossing_point(node.pt, node.a1, node.b0, node.b1) + apex = _exact_crossing_point(node) return _compare_angle_exact(apex, p, q) end @@ -324,8 +324,13 @@ function _exact_crossing_point(a0, a1, b0, b1) return (ax0 + t * dax, ay0 + t * day) end +# NodeKey call shape of the exact-crossing authority (design §2.6): one shape +# consumed by node-ordering keys, coincidence, and emission fallback. A crossing +# node keys its defining pair in `(pt, a1)`/`(b0, b1)`. +_exact_crossing_point(k::NodeKey) = _exact_crossing_point(k.pt, k.a1, k.b0, k.b1) + _exact_node_point(k::NodeKey) = k.is_crossing ? - _exact_crossing_point(k.pt, k.a1, k.b0, k.b1) : + _exact_crossing_point(k) : (Rational{BigInt}(GI.x(k.pt)), Rational{BigInt}(GI.y(k.pt))) function rk_nodes_coincide(::Planar, k1::NodeKey, k2::NodeKey; exact) @@ -334,6 +339,71 @@ function rk_nodes_coincide(::Planar, k1::NodeKey, k2::NodeKey; exact) return _exact_node_point(k1) == _exact_node_point(k2) end +# ## Node ordering along a segment (design §2.5) + +# Float64 crossing solve without the rational lift: the *filter* seed shared by +# node ordering (`rk_compare_along_segment`), coincidence proximity +# (node_identity), and the emission fast path (design §2.4/§2.6). Returns the +# approximate crossing coordinate together with `condK`, the determinant +# conditioning `(|dax·dby| + |day·dbx|) / |denom|` — near-parallel pairs (small +# `|denom|`) get a large `condK`, which every consumer folds into its error +# bound so the certified filter escalates instead of trusting the approximation. +@inline function _approx_crossing_point(a0, a1, b0, b1) + ax0, ay0 = Float64(GI.x(a0)), Float64(GI.y(a0)) + ax1, ay1 = Float64(GI.x(a1)), Float64(GI.y(a1)) + bx0, by0 = Float64(GI.x(b0)), Float64(GI.y(b0)) + bx1, by1 = Float64(GI.x(b1)), Float64(GI.y(b1)) + dax, day = ax1 - ax0, ay1 - ay0 + dbx, dby = bx1 - bx0, by1 - by0 + denom = dax * dby - day * dbx # nonzero for a proper crossing + t = ((bx0 - ax0) * dby - (by0 - ay0) * dbx) / denom + condK = (abs(dax * dby) + abs(day * dbx)) / max(abs(denom), floatmin(Float64)) + return (ax0 + t * dax, ay0 + t * day, condK) +end + +# Approximate coordinate of a node key (crossing: the float solve; vertex: the +# exact input coordinate) with a conservative absolute position error radius. +@inline function _approx_node_point(k::NodeKey) + if !k.is_crossing + return (Float64(GI.x(k.pt)), Float64(GI.y(k.pt)), 0.0) + end + x, y, condK = _approx_crossing_point(k.pt, k.a1, k.b0, k.b1) + #-- position error ~ (1 + conditioning) ulp of the coordinate magnitude; the + #-- 8× cushion covers the handful of rounding ops in the float solve above + err = 8 * eps(Float64) * (1.0 + condK) * (abs(x) + abs(y) + 1.0) + return (x, y, err) +end + +# Order two nodes along the oriented segment (s0, s1). The along parameter is +# the projection onto the segment direction `d = s1 - s0`; the float gap is +# trusted only when it exceeds the summed projected position errors plus the +# dot-product rounding (the certified filter). Otherwise the exact rational +# parameters decide (design §2.5). Zero means the nodes coincide. +function rk_compare_along_segment(m::Planar, s0, s1, na::NodeKey, nb::NodeKey; exact) + s0x, s0y = Float64(GI.x(s0)), Float64(GI.y(s0)) + dx = Float64(GI.x(s1)) - s0x + dy = Float64(GI.y(s1)) - s0y + scale = abs(dx) + abs(dy) + xa, ya, ea = _approx_node_point(na) + xb, yb, eb = _approx_node_point(nb) + ta = (xa - s0x) * dx + (ya - s0y) * dy + tb = (xb - s0x) * dx + (yb - s0y) * dy + #-- projected position error (ea, eb) times |d| ≈ scale, plus the ~4 ulp of + #-- the two subtract/multiply/add chains that form each parameter + tol = (ea + eb) * scale + 8 * eps(Float64) * (abs(ta) + abs(tb) + scale * scale) + gap = tb - ta + abs(gap) > tol && return gap < 0 ? 1 : -1 + #-- exact fallback (lazy): rational along-parameters, exact since Float64 are + #-- dyadic. Reached only for pairs the filter cannot separate. + R = Rational{BigInt} + dxr = R(GI.x(s1)) - R(GI.x(s0)); dyr = R(GI.y(s1)) - R(GI.y(s0)) + s0xr = R(GI.x(s0)); s0yr = R(GI.y(s0)) + pa = _exact_node_point(na); pb = _exact_node_point(nb) + tar = (pa[1] - s0xr) * dxr + (pa[2] - s0yr) * dyr + tbr = (pb[1] - s0xr) * dxr + (pb[2] - s0yr) * dyr + return tar < tbr ? -1 : (tar > tbr ? 1 : 0) +end + #= `_compare_angle` with an exact rational apex: the slow path of `rk_compare_edge_dir` for direction points incident to a D3 diff --git a/src/methods/geom_relations/relateng/kernel_spherical.jl b/src/methods/geom_relations/relateng/kernel_spherical.jl index 5a40c6258d..2d8224d514 100644 --- a/src/methods/geom_relations/relateng/kernel_spherical.jl +++ b/src/methods/geom_relations/relateng/kernel_spherical.jl @@ -419,6 +419,36 @@ function rk_nodes_coincide(::Spherical, k1::NodeKey, k2::NodeKey; exact) return _iszero3(_cross3(d1, d2)) && _dot3(d1, d2) > 0 end +# ## Node ordering along an arc (design §2.5) +# +# Two nodes on the minor arc (s0, s1) are ordered by the sign of the +# discriminant `(da × db) · N`, `N = s0 × s1` the arc's plane normal and `da`, +# `db` the nodes' on-sphere directions (`_exact_node_dir`): `da` precedes `db` +# along `s0 → s1` iff the discriminant is positive (both directions lie strictly +# interior to the minor arc, so they sit in the half where `N` points out of the +# turning plane). The float filter uses the `False()` (Float64) directions; it +# is trusted only when `|disc|` clears a conditioning-inflated bound (a +# near-tangent crossing has a large-error direction and must escalate). The exact +# fallback recomputes the discriminant over `Rational{BigInt}` directions. +function rk_compare_along_segment(m::Spherical, s0, s1, na::NodeKey, nb::NodeKey; exact) + S0 = _vec3(False(), s0); S1 = _vec3(False(), s1) + N = _cross3(S0, S1) + da = _exact_node_dir(False(), na); db = _exact_node_dir(False(), nb) + disc = _dot3(_cross3(da, db), N) + #-- |disc| ≈ |da||db||N| sin(Δ); the directions carry ≲ few·ulp relative + #-- error (unit-ish vectors), amplified for crossing nodes by their arc + #-- geometry. 64 ulp of the product magnitude is a conservative escalation + #-- trigger — coincident/near-coincident pairs fall to the exact path. + mag = sqrt(_dot3(da, da) * _dot3(db, db) * _dot3(N, N)) + tol = 64 * eps(Float64) * mag + abs(disc) > tol && return disc > 0 ? -1 : 1 + #-- exact fallback (lazy): rational directions, exact for Float64 inputs. + Se0 = _vec3(True(), s0); Se1 = _vec3(True(), s1); Ne = _cross3(Se0, Se1) + ea = _exact_node_dir(True(), na); eb = _exact_node_dir(True(), nb) + o = _dot3(_cross3(ea, eb), Ne) + return o > 0 ? -1 : (o < 0 ? 1 : 0) +end + # ## Ring orientation #= diff --git a/src/methods/geom_relations/relateng/relate_geometry.jl b/src/methods/geom_relations/relateng/relate_geometry.jl index c3d9c67db5..464fd98949 100644 --- a/src/methods/geom_relations/relateng/relate_geometry.jl +++ b/src/methods/geom_relations/relateng/relate_geometry.jl @@ -602,6 +602,21 @@ the stored winding is authoritative. _ring_interior_on_left(m, pts::Vector, is_hole::Bool; exact) = _ring_is_ccw(m, pts; exact) +#= +The polygon's **material-interior** side of a ring, as opposed to the ring's +denoted region (`_ring_interior_on_left`). For a shell the denoted region *is* +the material interior; for a hole the denoted region is the cavity, so the +material interior is on the opposite side — one flip, defined exactly once so +relate, overlay, and extents agree by construction (OverlayNG design §2.7). +Overlay labeling derives JTS `Edge.depth_delta` from this +(`material_interior_on_left ? -1 : +1`, matching JTS `locationLeft/Right`: +a positive delta means Left = EXTERIOR, Right = INTERIOR). Manifold-generic: +it delegates to the manifold-dispatched `_ring_interior_on_left`. +=# +_ring_material_interior_on_left(m, pts::Vector, is_hole::Bool; exact) = + is_hole ? !_ring_interior_on_left(m, pts, is_hole; exact) : + _ring_interior_on_left(m, pts, is_hole; exact) + #= Port of JTS `Orientation.isCCW(CoordinateSequence)` with the orientation index routed through the kernel (`rk_orient`): whether the closed `ring` diff --git a/src/methods/geom_relations/relateng/topology_computer.jl b/src/methods/geom_relations/relateng/topology_computer.jl index 6a1b38d11a..20b80c3a11 100644 --- a/src/methods/geom_relations/relateng/topology_computer.jl +++ b/src/methods/geom_relations/relateng/topology_computer.jl @@ -725,6 +725,6 @@ end # BigFloat's default precision is deterministic, though not strictly # correctly rounded). function _crossing_locate_point(key::NodeKey) - xr, yr = _exact_crossing_point(key.pt, key.a1, key.b0, key.b1) + xr, yr = _exact_crossing_point(key) return (Float64(BigFloat(xr)), Float64(BigFloat(yr))) end diff --git a/test/methods/clipping/overlayng/noding.jl b/test/methods/clipping/overlayng/noding.jl new file mode 100644 index 0000000000..ff27ed4c21 --- /dev/null +++ b/test/methods/clipping/overlayng/noding.jl @@ -0,0 +1,312 @@ +# Tests for the OverlayNG phase-1 noding substrate (design §2.9): the +# `NodedArrangement` invariants, along-segment ordering vs the exact authority, +# the certified emission fast paths, and the rounded-arrangement / classification +# censuses on a small Natural Earth subset. + +using Test +import GeometryOps as GO +import GeometryOps: Planar, Spherical, True, False +import GeometryOps.UnitSpherical: UnitSphericalPoint, UnitSphereFromGeographic +import GeoInterface as GI +import Extents +using LinearAlgebra: cross, dot, norm + +# --------------------------------------------------------------------------- +# helpers +# --------------------------------------------------------------------------- + +_crossing_ids(arr) = [Int32(i) for i in 1:GO.num_nodes(arr) if arr.nodes.keys[i].is_crossing] + +# every proper crossing appears as one shared node id on exactly one A-segment and +# one B-segment interior list (invariant 1); node ids are unique (invariant 2); +# no NodedEdge is zero-length (invariant 3). +function check_invariants(arr, na_strings) + @test length(unique(arr.nodes.keys)) == GO.num_nodes(arr) # 2 + for e in arr.edges + @test e.node_lo != e.node_hi # 3 + end + for cid in _crossing_ids(arr) + a_hits = 0; b_hits = 0 + for ((si, _), ids) in arr.seg_nodes + if cid in ids + si <= na_strings ? (a_hits += 1) : (b_hits += 1) + end + end + @test a_hits >= 1 && b_hits >= 1 # 1 + end +end + +# exact-always along-parameter order of a segment's interior node ids +function exact_order(arr, s0, s1, ids) + R = Rational{BigInt} + dxr = R(GI.x(s1)) - R(GI.x(s0)); dyr = R(GI.y(s1)) - R(GI.y(s0)) + param(id) = let p = GO._exact_node_point(arr.nodes.keys[id]) + (p[1] - R(GI.x(s0))) * dxr + (p[2] - R(GI.y(s0))) * dyr + end + return sort(ids; by = param) +end + +# exact-always along-parameter order on the sphere +function exact_order_sph(arr, s0, s1, ids) + Ne = GO._cross3(GO._vec3(True(), s0), GO._vec3(True(), s1)) + dir(id) = GO._exact_node_dir(True(), arr.nodes.keys[id]) + return sort(ids; lt = (i, j) -> GO._dot3(GO._cross3(dir(i), dir(j)), Ne) > 0) +end + +shift_geom(g, dx, dy) = GO.apply(GI.PointTrait(), g) do p + (GI.x(p) + dx, GI.y(p) + dy) +end + +# --------------------------------------------------------------------------- +# 1. Invariants on constructed cases (§2.1) +# --------------------------------------------------------------------------- + +@testset "two crossing quads" begin + A = GI.Polygon([[(0.0, 0.0), (4.0, 0.0), (4.0, 4.0), (0.0, 4.0), (0.0, 0.0)]]) + B = GI.Polygon([[(2.0, 2.0), (6.0, 2.0), (6.0, 6.0), (2.0, 6.0), (2.0, 2.0)]]) + for m in (Planar(), Spherical()) + arr = GO.NodedArrangement(m, A, B; exact = True()) + na = count(ss -> ss.is_a, arr.segstrings) + @test length(_crossing_ids(arr)) == 2 + check_invariants(arr, na) + if m isa Planar + cpts = sort([GO.node_point(arr, i) for i in _crossing_ids(arr)]) + @test cpts == [(2.0, 4.0), (4.0, 2.0)] # bit-exact on the integer grid + end + end +end + +@testset "degree-6 node (tier-2 merge of distinct keys)" begin + # two A lines and one B line all through the origin -> two distinct crossing + # keys, crossing_node(lineA1,B) and crossing_node(lineA2,B), coincident at the + # origin: tier 2 must merge them into one node. + A = GI.MultiLineString([[(-1.0, -1.0), (1.0, 1.0)], [(-1.0, 1.0), (1.0, -1.0)]]) + B = GI.LineString([(-1.0, 0.0), (1.0, 0.0)]) + for m in (Planar(), Spherical()) + arr = GO.NodedArrangement(m, A, B; exact = True()) + cids = _crossing_ids(arr) + @test length(cids) == 1 # merged into one node + p = GO.node_point(arr, cids[1]) + @test isapprox(p[1], 0.0; atol = 1e-12) && isapprox(p[2], 0.0; atol = 1e-12) + # the merged node is incident to both A lines and B (three parent strings) + na = count(ss -> ss.is_a, arr.segstrings) + check_invariants(arr, na) + end +end + +@testset "crossing exactly on a third string's vertex" begin + # A horizontal line crosses B-line-1 (vertical) at the origin, which is also + # B-line-2's endpoint vertex: the crossing key and the vertex key coincide + # and tier 2 merges them. + A = GI.LineString([(-2.0, 0.0), (2.0, 0.0)]) + B = GI.MultiLineString([[(0.0, -2.0), (0.0, 2.0)], [(0.0, 0.0), (1.0, 1.0)]]) + for m in (Planar(), Spherical()) + arr = GO.NodedArrangement(m, A, B; exact = True()) + # the origin is a single node shared by A, B-line-1 and B-line-2's vertex + origin_ids = [i for i in 1:GO.num_nodes(arr) + if isapprox(GO.node_point(arr, i)[1], 0.0; atol = 1e-12) && + isapprox(GO.node_point(arr, i)[2], 0.0; atol = 1e-12)] + @test length(origin_ids) == 1 + end +end + +@testset "collinear shared boundary — zero phantom crossings" begin + # edge-adjacent squares sharing the vertex-identical edge x = 2 + A = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]]) + B = GI.Polygon([[(2.0, 0.0), (4.0, 0.0), (4.0, 2.0), (2.0, 2.0), (2.0, 0.0)]]) + for m in (Planar(), Spherical()) + arr = GO.NodedArrangement(m, A, B; exact = True()) + @test length(_crossing_ids(arr)) == 0 # zero phantom crossings + @test sum(length, values(arr.seg_nodes); init = 0) == 0 # zero interior splits + end +end + +@testset "a-b-a spike input" begin + # B retraces (0,0)->(1,1)->(0,0); A crosses the retraced segment at one point, + # reported by both candidate pairs but the same canonical crossing key. + A = GI.LineString([(-1.0, 0.5), (2.0, 0.5)]) + B = GI.LineString([(0.0, 0.0), (1.0, 1.0), (0.0, 0.0)]) + for m in (Planar(), Spherical()) + arr = GO.NodedArrangement(m, A, B; exact = True()) # must not throw + cids = _crossing_ids(arr) + @test length(cids) == 1 # one merged node + m isa Planar && @test GO.node_point(arr, cids[1]) == (0.5, 0.5) + end +end + +# --------------------------------------------------------------------------- +# 2. Ordering cross-check: float-filtered order == exact-always order (§2.5) +# --------------------------------------------------------------------------- + +@testset "dense comb ordering matches exact (planar)" begin + A = GI.LineString([(0.0, 0.0), (201.0, 0.0)]) + B = GI.MultiLineString([[(Float64(i) + 0.3, -1.0), (Float64(i) + 0.3, 1.0)] for i in 1:200]) + arr = GO.NodedArrangement(Planar(), A, B; exact = True()) + # A is string 1, its single segment carries all 200 interior crossings + ids = arr.seg_nodes[(Int32(1), Int32(1))] + @test length(ids) == 200 + s0 = arr.segstrings[1].pts[1]; s1 = arr.segstrings[1].pts[2] + @test ids == exact_order(arr, s0, s1, ids) # elementwise + # strictly increasing along the segment + for c in 2:length(ids) + @test GO.rk_compare_along_segment(Planar(), s0, s1, arr.nodes.keys[ids[c-1]], arr.nodes.keys[ids[c]]; exact = True()) < 0 + end +end + +@testset "dense comb ordering matches exact (spherical)" begin + A = GI.LineString([(0.0, 0.0), (60.0, 0.0)]) + B = GI.MultiLineString([[(Float64(i) * 0.25 + 0.1, -1.0), (Float64(i) * 0.25 + 0.1, 1.0)] for i in 1:200]) + arr = GO.NodedArrangement(Spherical(), A, B; exact = True()) + ids = arr.seg_nodes[(Int32(1), Int32(1))] + @test length(ids) == 200 + s0 = arr.segstrings[1].pts[1]; s1 = arr.segstrings[1].pts[2] + @test ids == exact_order_sph(arr, s0, s1, ids) +end + +# --------------------------------------------------------------------------- +# 3. Emission certificate audit (§2.6) +# --------------------------------------------------------------------------- + +@testset "planar emission: certified == rational, every node" begin + # dense generic-slope grid + shifted-self coastline-like crossings + Ag = GI.MultiLineString([[(Float64(k) * 4.0, 0.0), (Float64(k) * 4.0 + 0.31, 1000.0)] for k in 1:60]) + Bg = GI.MultiLineString([[(0.0, Float64(j) * 4.0), (1000.0, Float64(j) * 4.0 + 0.29)] for j in 1:60]) + arr = GO.NodedArrangement(Planar(), Ag, Bg; exact = True()) + ncert = 0; ntot = 0 + for i in _crossing_ids(arr) + k = arr.nodes.keys[i] + (x, y, cert) = GO._certified_crossing(k.pt, k.a1, k.b0, k.b1) + rx, ry = GO._exact_crossing_point(k) + rat = (Float64(rx), Float64(ry)) + ntot += 1 + if cert + ncert += 1 + @test (x, y) == rat # certified must equal rational + end + @test GO.node_point(arr, i) == rat # node_point rounds to rational either way + end + @test ntot > 1000 + @test ncert == ntot # 100% certified on clean data (S3) +end + +@testset "spherical emission: direction within bound of exact" begin + Ag = GI.MultiLineString([[(Float64(k) * 0.09 + 0.05, 0.0), (Float64(k) * 0.09 + 0.05 + 0.031, 20.0)] for k in 1:60]) + Bg = GI.MultiLineString([[(0.0, Float64(j) * 0.09 + 0.05), (20.0, Float64(j) * 0.09 + 0.05 + 0.029)] for j in 1:60]) + arr = GO.NodedArrangement(Spherical(), Ag, Bg; exact = True()) + maxdev = 0.0 + for i in _crossing_ids(arr) + k = arr.nodes.keys[i] + emitted = GO.node_point(arr, i) + exact = GO._dir_to_lonlat(GO._sph_crossing_dir(True(), k)) + maxdev = max(maxdev, abs(emitted[1] - exact[1]), abs(emitted[2] - exact[2])) + end + @test length(_crossing_ids(arr)) > 1000 + @test maxdev <= 1e-8 # measured ≤1.4e-14° (S3) +end + +# --------------------------------------------------------------------------- +# 4 & 5. Natural Earth subset: rounded-arrangement audit + classification census +# --------------------------------------------------------------------------- + +# direct A×B classification census (proper/touch/collinear counts, flag check) +function classify_census(m, ssa, ssb) + ta = GO._relate_edge_index(m, ssa); tb = GO._relate_edge_index(m, ssb) + (ta === nothing || tb === nothing) && return (0, 0, 0, true) + nprop = 0; ntouch = 0; ncol = 0; flags_ok = true + GO.SpatialTreeInterface.dual_depth_first_search(Extents.intersects, ta, tb) do ia, ib + (sa, ka) = ta.data[ia]; (sb, kb) = tb.data[ib] + a0 = ssa[sa].pts[ka]; a1 = ssa[sa].pts[ka+1] + b0 = ssb[sb].pts[kb]; b1 = ssb[sb].pts[kb+1] + c = GO.rk_classify_intersection(m, a0, a1, b0, b1; exact = True()) + if c.kind == GO.SS_PROPER + nprop += 1 + elseif c.kind == GO.SS_TOUCH + ntouch += 1 + (c.a0_on_b || c.a1_on_b || c.b0_on_a || c.b1_on_a) || (flags_ok = false) + elseif c.kind == GO.SS_COLLINEAR + ncol += 1 + (c.a0_on_b || c.a1_on_b || c.b0_on_a || c.b1_on_a) || (flags_ok = false) + end + return nothing + end + return (nprop, ntouch, ncol, flags_ok) +end + +# rounded-arrangement audit: the crossing-incident edges of the emitted geometry +# must not properly cross each other (only crossing nodes move, so re-classifying +# their incident edges suffices — S3). Returns the count of introduced crossings. +function rounded_crossings(arr) + incident = Dict{Int32, Vector{GO.NodedEdge}}() + for e in arr.edges + if arr.nodes.keys[e.node_lo].is_crossing + push!(get!(() -> GO.NodedEdge[], incident, e.node_lo), e) + end + if arr.nodes.keys[e.node_hi].is_crossing + push!(get!(() -> GO.NodedEdge[], incident, e.node_hi), e) + end + end + introduced = 0 + for (_, es) in incident + for i in 1:length(es), j in (i+1):length(es) + ea = es[i]; eb = es[j] + # only audit A-vs-B incident pairs (opposite sides can spuriously cross) + (arr.segstrings[ea.string_idx].is_a == arr.segstrings[eb.string_idx].is_a) && continue + pa0 = GO.node_point(arr, ea.node_lo); pa1 = GO.node_point(arr, ea.node_hi) + pb0 = GO.node_point(arr, eb.node_lo); pb1 = GO.node_point(arr, eb.node_hi) + GO.rk_classify_intersection(Planar(), pa0, pa1, pb0, pb1; exact = True()).kind == GO.SS_PROPER && + (introduced += 1) + end + end + return introduced +end + +ne_ok = false +ne_names = String[]; ne_geoms = Any[] +try + import NaturalEarth, GeoJSON + fc = NaturalEarth.naturalearth("admin_0_countries", 110) + for f in fc + g = GeoJSON.geometry(f) + (g === nothing || GI.npoint(g) == 0) && continue + nm = try; string(f.NAME); catch; "?"; end + push!(ne_names, nm); push!(ne_geoms, GO.tuples(g)) + end + global ne_ok = length(ne_geoms) > 0 +catch err + @info "Natural Earth subset skipped (data unavailable)" err +end + +@testset "Natural Earth subset (rounded-arrangement + census)" begin + if !ne_ok + @test_skip "Natural Earth data unavailable" + else + import LibGEOS as LG + picks = String["Brazil", "France", "Egypt", "India", "Australia"] + tested = 0 + for nm in picks + idx = findfirst(==(nm), ne_names) + idx === nothing && continue + A = ne_geoms[idx] + LG.isValid(GI.convert(LG, A)) || continue + B = shift_geom(A, 0.5, 0.0) + tested += 1 + # census: planar & spherical proper-crossing multiset identical, flags ok + ssa_p = GO._overlay_segstrings(Planar(), A, true); ssb_p = GO._overlay_segstrings(Planar(), B, false) + ssa_s = GO._overlay_segstrings(Spherical(), A, true); ssb_s = GO._overlay_segstrings(Spherical(), B, false) + (pp, pt, pc, pf) = classify_census(Planar(), ssa_p, ssb_p) + (sp, st, sc, sf) = classify_census(Spherical(), ssa_s, ssb_s) + @test pf && sf # every touch/collinear carries a flag (§2.3) + #-- NB: proper counts differ across manifolds on real coastlines + #-- (great-circle arcs bow, so they cross a shifted copy differently + #-- than straight segments) — genuine geometry, validated independently + #-- by the rounded-arrangement audit below on each manifold. + @test pp > 0 && sp > 0 + # rounded-arrangement audit on both manifolds + arr_p = GO.NodedArrangement(Planar(), ssa_p, ssb_p; exact = True()) + arr_s = GO.NodedArrangement(Spherical(), ssa_s, ssb_s; exact = True()) + @test rounded_crossings(arr_p) == 0 + @test rounded_crossings(arr_s) == 0 + end + @test tested >= 2 + end +end diff --git a/test/methods/relateng/kernel_along_segment.jl b/test/methods/relateng/kernel_along_segment.jl new file mode 100644 index 0000000000..49997bfe79 --- /dev/null +++ b/test/methods/relateng/kernel_along_segment.jl @@ -0,0 +1,116 @@ +# Tests for the OverlayNG phase-1 kernel additions (design §2.5–§2.7): +# `rk_compare_along_segment` (planar + spherical, float filter vs exact +# fallback), the unified exact-crossing authority NodeKey call shape, and +# `_ring_material_interior_on_left`. + +using Test +import GeometryOps as GO +import GeometryOps: Planar, Spherical, True, False +import GeometryOps.UnitSpherical: UnitSphericalPoint, UnitSphereFromGeographic +import GeoInterface as GI +using LinearAlgebra: cross, dot +using Random + +# crossing node of horizontal segment (0,0)-(L,0) with the vertical through x. +_planar_xnode(L, x) = GO.crossing_node((0.0, 0.0), (L, 0.0), (x, -1.0), (x, 1.0)) + +@testset "rk_compare_along_segment (planar)" begin + m = Planar() + s0, s1 = (0.0, 0.0), (10.0, 0.0) + k3 = _planar_xnode(10.0, 3.0) + k7 = _planar_xnode(10.0, 7.0) + v5 = GO.vertex_node((5.0, 0.0)) + + @test GO.rk_compare_along_segment(m, s0, s1, k3, k7; exact = True()) == -1 + @test GO.rk_compare_along_segment(m, s0, s1, k7, k3; exact = True()) == 1 + @test GO.rk_compare_along_segment(m, s0, s1, k3, v5; exact = True()) == -1 + @test GO.rk_compare_along_segment(m, s0, s1, v5, k7; exact = True()) == -1 + #-- coincidence returns 0 (same node) + @test GO.rk_compare_along_segment(m, s0, s1, k3, k3; exact = True()) == 0 + + #-- reversed segment direction flips every order + @test GO.rk_compare_along_segment(m, s1, s0, k3, k7; exact = True()) == 1 +end + +# exact-always along-parameter (rational) for a node key on segment (s0,s1) +function _exact_param(s0, s1, k) + R = Rational{BigInt} + dxr = R(s1[1]) - R(s0[1]); dyr = R(s1[2]) - R(s0[2]) + p = GO._exact_node_point(k) + return (p[1] - R(s0[1])) * dxr + (p[2] - R(s0[2])) * dyr +end + +@testset "planar order matches exact-always (dense)" begin + m = Planar() + rng = MersenneTwister(42) + s0, s1 = (0.0, 0.0), (1.0, 0.0) + for _ in 1:200 + xs = sort!(rand(rng, 6) .* 0.9 .+ 0.05) + # crossing nodes at those xs, with random (well-conditioned) slopes + ks = [GO.crossing_node(s0, s1, (x, -rand(rng) - 0.2), (x, rand(rng) + 0.2)) for x in xs] + # comparator sort + order_filter = sortperm(collect(1:length(ks)); + lt = (i, j) -> GO.rk_compare_along_segment(m, s0, s1, ks[i], ks[j]; exact = True()) < 0) + order_exact = sortperm([_exact_param(s0, s1, k) for k in ks]) + @test order_filter == order_exact + # since xs are strictly increasing, both must equal 1:length + @test order_filter == collect(1:length(ks)) + end +end + +@testset "planar near-parallel escalates soundly" begin + m = Planar() + s0, s1 = (0.0, 0.0), (1.0, 0.0) + # two crossings extremely close together, from near-parallel secondaries: + # the filter must return the same sign as the exact parameters. + x1 = 0.5 + x2 = nextfloat(nextfloat(0.5)) + ka = GO.crossing_node(s0, s1, (x1, -1e-9), (x1, 1e-9)) + kb = GO.crossing_node(s0, s1, (x2, -1e-9), (x2, 1e-9)) + got = GO.rk_compare_along_segment(m, s0, s1, ka, kb; exact = True()) + want = _exact_param(s0, s1, ka) < _exact_param(s0, s1, kb) ? -1 : 1 + @test got == want +end + +@testset "rk_compare_along_segment (spherical)" begin + ms = Spherical() + usp(lon, lat) = UnitSphereFromGeographic()((lon, lat)) + s0, s1 = usp(0.0, 0.0), usp(10.0, 0.0) + k3 = GO.crossing_node(s0, s1, usp(3.0, -1.0), usp(3.0, 1.0)) + k7 = GO.crossing_node(s0, s1, usp(7.0, -1.0), usp(7.0, 1.0)) + v5 = GO.vertex_node(usp(5.0, 0.0)) + + @test GO.rk_compare_along_segment(ms, s0, s1, k3, k7; exact = True()) == -1 + @test GO.rk_compare_along_segment(ms, s0, s1, k7, k3; exact = True()) == 1 + @test GO.rk_compare_along_segment(ms, s0, s1, k3, v5; exact = True()) == -1 + @test GO.rk_compare_along_segment(ms, s0, s1, v5, k7; exact = True()) == -1 + @test GO.rk_compare_along_segment(ms, s0, s1, k3, k3; exact = True()) == 0 + @test GO.rk_compare_along_segment(ms, s1, s0, k3, k7; exact = True()) == 1 + + #-- dense monotone check: crossings at increasing longitudes order correctly + lons = collect(1.0:0.5:9.0) + ks = [GO.crossing_node(s0, s1, usp(x, -1.0), usp(x, 1.0)) for x in lons] + order = sortperm(collect(1:length(ks)); + lt = (i, j) -> GO.rk_compare_along_segment(ms, s0, s1, ks[i], ks[j]; exact = True()) < 0) + @test order == collect(1:length(ks)) +end + +@testset "unified exact-crossing authority" begin + k = _planar_xnode(10.0, 4.0) + @test GO._exact_crossing_point(k) == GO._exact_crossing_point((0.0, 0.0), (10.0, 0.0), (4.0, -1.0), (4.0, 1.0)) + @test GO._exact_node_point(k) == GO._exact_crossing_point(k) + v = GO.vertex_node((2.0, 3.0)) + @test GO._exact_node_point(v) == (Rational{BigInt}(2), Rational{BigInt}(3)) + @test Float64.(GO._exact_crossing_point(k)) == (4.0, 0.0) +end + +@testset "_ring_material_interior_on_left flip" begin + for m in (Planar(), Spherical()) + # a CW square (shell) and its use as a hole + ptsll = [(0.0, 0.0), (0.0, 4.0), (4.0, 4.0), (4.0, 0.0), (0.0, 0.0)] + pts = m isa Planar ? ptsll : [UnitSphereFromGeographic()(p) for p in ptsll] + il = GO._ring_interior_on_left(m, pts, false; exact = True()) + @test GO._ring_material_interior_on_left(m, pts, false; exact = True()) == il + @test GO._ring_material_interior_on_left(m, pts, true; exact = True()) == !il + end +end diff --git a/test/methods/relateng/runtests.jl b/test/methods/relateng/runtests.jl index ab5269128c..f7b3da3377 100644 --- a/test/methods/relateng/runtests.jl +++ b/test/methods/relateng/runtests.jl @@ -4,6 +4,7 @@ using SafeTestsets @safetestset "Predicates" begin include("predicates.jl") end @safetestset "Kernel" begin include("kernel.jl") end @safetestset "Kernel conformance" begin include("kernel_conformance.jl") end +@safetestset "Kernel along-segment ordering" begin include("kernel_along_segment.jl") end @safetestset "Spherical kernel conformance" begin include("kernel_conformance_spherical.jl") end @safetestset "S2 edge-crossing conformance" begin include("s2_crossings_conformance.jl") end @safetestset "Point locator" begin include("point_locator.jl") end diff --git a/test/runtests.jl b/test/runtests.jl index 55f9cc8a9e..ec62ff8c8b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -44,6 +44,7 @@ end @safetestset "Intersection Point" begin include("methods/clipping/intersection_points.jl") end @safetestset "Polygon Clipping" begin include("methods/clipping/polygon_clipping.jl") end @safetestset "Sutherland-Hodgman" begin include("methods/clipping/sutherland_hodgman.jl") end +@safetestset "OverlayNG noding" begin include("methods/clipping/overlayng/noding.jl") end # Transformations @safetestset "Embed Extent" begin include("transformations/extent.jl") end @safetestset "Reproject" begin include("transformations/reproject.jl") end