From 11af875f0a75e5e465a889880399fb9556eaba54 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 10 Jul 2026 11:43:15 +0000 Subject: [PATCH 1/2] phase 1 of issue #107 --- src_rust/src/coverage.rs | 171 ++++++++++++++++++++++++++++----- src_rust/src/coverage/tests.rs | 103 ++++++++++++++++++-- src_rust/src/sphere/tests.rs | 31 ++++++ 3 files changed, 274 insertions(+), 31 deletions(-) diff --git a/src_rust/src/coverage.rs b/src_rust/src/coverage.rs index 5555799..e0e3a7a 100644 --- a/src_rust/src/coverage.rs +++ b/src_rust/src/coverage.rs @@ -30,6 +30,7 @@ use std::cmp::Ordering; use std::collections::BinaryHeap; use std::f64::consts::PI; +use std::sync::OnceLock; use rayon::prelude::*; use smallvec::SmallVec; @@ -486,13 +487,29 @@ fn arc_crossing_parity( /// symbolic [`arcs_cross_sos`]. See [`arc_crossing_parity`]. #[inline] fn edge_crosses_probe(p: &Vec3, q: &Vec3, ip: PointId, iq: PointId, n_pq: &Vec3, e: &Edge) -> bool { + edge_crosses_probe_d(p, q, ip, iq, n_pq, e, dot(&e.n_ab, p), dot(&e.n_ab, q)) +} + +/// [`edge_crosses_probe`] with the probe endpoints' edge-plane determinants +/// `d_ab_p`, `d_ab_q` supplied by the caller — the fill-leg walk +/// ([`fill_leg_parity`]) computes them anyway for the #107 overlap gate. +#[inline] +#[allow(clippy::too_many_arguments)] +fn edge_crosses_probe_d( + p: &Vec3, + q: &Vec3, + ip: PointId, + iq: PointId, + n_pq: &Vec3, + e: &Edge, + d_ab_p: f64, + d_ab_q: f64, +) -> bool { // The four straddle determinants of the standard arcs-cross test. Near-zero // in any of them is the cell-centre-on-edge degeneracy (#11) where the plain // sign test is unstable; defer those to SoS. let d_pq_a = dot(n_pq, &e.a); let d_pq_b = dot(n_pq, &e.b); - let d_ab_p = dot(&e.n_ab, p); - let d_ab_q = dot(&e.n_ab, q); if d_pq_a.abs() < ORIENT_EPS || d_pq_b.abs() < ORIENT_EPS || d_ab_p.abs() < ORIENT_EPS @@ -504,6 +521,45 @@ fn edge_crosses_probe(p: &Vec3, q: &Vec3, ip: PointId, iq: PointId, n_pq: &Vec3, (d_pq_a > 0.0) != (d_pq_b > 0.0) && (d_ab_p > 0.0) != (d_ab_q > 0.0) } +/// Crossing parity of a **fill leg** `p → q`, or `None` when the leg is +/// **degenerate**: it overlaps some relevant edge's great circle, i.e. both +/// endpoints lie within [`ORIENT_EPS`] of that edge's plane (issue #107). +/// +/// A leg that runs *along* an edge's great circle over a finite span — a +/// probe chord on a HEALPix centre meridian shared with a long polygon edge — +/// cannot be chained: the intersection points of the two near-coincident +/// circles are determined by 1-ulp residues, and while every sidedness sign is +/// individually exact ([`arcs_cross_sos`]), the verdicts of the overlapped +/// edge and its transversal sibling at a shared vertex can be jointly wrong by +/// one crossing, silently inverting the subtree (the #107 reproducer). The +/// caller must classify such a cell **directly** instead +/// ([`DirectClassifier`]). The gate costs two dots per edge that +/// [`edge_crosses_probe_d`] needs anyway, so the clean path is unchanged. +#[cfg_attr(not(test), allow(dead_code))] // wired into the descent in phase 2 (#107) +fn fill_leg_parity( + p: &Vec3, + q: &Vec3, + ip: PointId, + iq: PointId, + relevant: &[usize], + edges: &[Edge], +) -> Option { + let n_pq = cross(p, q); + let mut crossings = 0u32; + for &i in relevant { + let e = &edges[i]; + let d_ab_p = dot(&e.n_ab, p); + let d_ab_q = dot(&e.n_ab, q); + if d_ab_p.abs() < ORIENT_EPS && d_ab_q.abs() < ORIENT_EPS { + return None; // leg overlaps this edge's great circle (#107) + } + if edge_crosses_probe_d(p, q, ip, iq, &n_pq, e, d_ab_p, d_ab_q) { + crossings += 1; + } + } + Some(crossings & 1 == 1) +} + /// Does polygon edge `e` cross **or exactly touch** the cell edge `c1 → c2` /// (normal `n_c = c1 × c2`)? The straddle test's closed-set form (#103). /// @@ -663,6 +719,87 @@ fn seed_fill(x: &Vec3, unit_normals: &[Vec3], rings: &[Vec]) -> Option bool { + let crossings = edges + .iter() + .filter(|e| arcs_cross_sos(r, x, &e.a, &e.b, ir, ix, e.ia, e.ib)) + .count(); + r_fill ^ (crossings % 2 == 1) +} + +/// Direct even-odd classification for **degenerate** descent cells (#107). +/// +/// A cell whose incoming fill leg overlaps a polygon edge's great circle +/// ([`fill_leg_parity`] returns `None`) is never chained; it is classified +/// directly instead — the same construction [`base_fills`] uses for seeds +/// whose centre lies on the boundary, generalized to any cell centre: pick an +/// unambiguous reference (a base centre off every edge plane, with a +/// trustworthy winding verdict) and run one [`arcs_cross_sos`] pass over all +/// edges ([`chain_fill_from`]). The invariant is one sentence: *degenerate +/// cells are classified directly, never chained.* +/// +/// Donor discovery is lazy (`OnceLock`): clean geometry never gates, so it +/// pays neither the per-edge normal normalization nor the seed PIPs. +#[cfg_attr(not(test), allow(dead_code))] // wired into the descent in phase 2 (#107) +struct DirectClassifier<'a> { + edges: &'a [Edge], + rings: &'a [Vec], + donors: OnceLock>, +} + +#[cfg_attr(not(test), allow(dead_code))] // wired into the descent in phase 2 (#107) +impl<'a> DirectClassifier<'a> { + fn new(edges: &'a [Edge], rings: &'a [Vec]) -> Self { + Self { + edges, + rings, + donors: OnceLock::new(), + } + } + + /// Fill at cell centre `x` (stable SoS id `ix`), classified directly from + /// the nearest unambiguous base-centre reference. Nearest keeps the + /// reference chord short; exact antipodality cannot occur (a sub-base cell + /// centre never coincides bit-exactly with a base centre or its antipode), + /// so the chord is always a valid minor arc for [`arcs_cross_sos`]. When + /// every base centre is ambiguous (boundary planes through all 12 — + /// pathological), keep the raw winding verdict, as [`base_fills`] does. + fn classify(&self, x: &Vec3, ix: PointId) -> bool { + let donors = self.donors.get_or_init(|| { + let units: Vec = self.edges.iter().map(|e| normalize(&e.n_ab)).collect(); + (0..12u64) + .filter_map(|b| { + let c = cell_center_vec(0, b); + seed_fill(&c, &units, self.rings).map(|f| (c, center_id(0, b), f)) + }) + .collect() + }); + match donors + .iter() + .max_by(|u, v| dot(&u.0, x).total_cmp(&dot(&v.0, x))) + { + Some(&(ref c, ir, f)) => chain_fill_from(c, ir, f, x, ix, self.edges), + None => parity_filled_robust(x, self.rings), + } + } +} + /// Are base cells `a` and `b` centred on antipodal points? The 12 HEALPix /// base centres pair up as north-cap/south-cap (0,10) (1,11) (2,8) (3,9) and /// equatorial (4,6) (5,7); a chain probe between an antipodal pair would be a @@ -735,28 +872,14 @@ fn base_fills(edges: &[Edge], rings: &[Vec], cap: &Cap, complement: bool) let d = (0..12u64) .find(|&d| known[d as usize] && !antipodal_base(b, d)) .expect("a non-antipodal known donor base centre"); - // Chain arcs between base centres reach ~122°, where the bare - // two-straddle fast path of `edge_crosses_probe` also fires on the - // *far* intersection of the two great circles (an edge antipodal to - // the chord) and silently inverts the seed. The full symbolic - // predicate is unconditional here — 11 arcs at most, off every hot - // path. - let crossings = edges - .iter() - .filter(|e| { - arcs_cross_sos( - ¢ers[d as usize], - ¢ers[b as usize], - &e.a, - &e.b, - center_id(0, d), - center_id(0, b), - e.ia, - e.ib, - ) - }) - .count(); - fill[b as usize] = fill[d as usize] ^ (crossings % 2 == 1); + fill[b as usize] = chain_fill_from( + ¢ers[d as usize], + center_id(0, d), + fill[d as usize], + ¢ers[b as usize], + center_id(0, b), + edges, + ); known[b as usize] = true; } fill diff --git a/src_rust/src/coverage/tests.rs b/src_rust/src/coverage/tests.rs index 1f2b2fd..e5bbcdb 100644 --- a/src_rust/src/coverage/tests.rs +++ b/src_rust/src/coverage/tests.rs @@ -589,15 +589,104 @@ fn test_base_fills_chain_no_antipodal_phantom() { } #[test] -#[ignore = "known limitation, out of #103 scope: a polygon edge collinear \ -with the probe lattice over tens of degrees on a hemisphere+ ring can still \ -desynchronize the descent's vertex-graze bookkeeping between the collinear \ -edge and its transversal sibling; predates this PR (the parity oracle made \ -it visible) — see PR #106 'Questions for review'"] +fn test_direct_classifier_reproducer_centres() { + // #107 phase 1: the direct classifier at the reproducer's diagnosed + // centres. Ring: hemisphere+ with a 40° edge along the lon-45 meridian, + // collinear with base 0's centre lattice. + // + // NOTE — this test *corrects* the #107 issue text. centre(3,1) at + // (66.44, 45) is truly INSIDE: walking (45, 50) [exterior wedge: right of + // edge 0, right of edge 1] → (60, 50) crosses edge 1 exactly once, and + // (60, 50) → centre(3,1) crosses nothing. The issue called it outside on + // the strength of `parity_filled_robust` — but the winding sum is + // antisymmetric under x → −x (the projections onto the plane ⊥x are + // identical for both poles; only the sign term flips), so it computes + // k(x) − k(−x), not k(x), and −centre(3,1) = (−66.44, 225) is *also* + // interior ⇒ the backend reads 0 (outside) for a truly interior point. + // The classifier, chained from an off-plane donor whose antipode is + // exterior (where the winding verdict IS correct), returns the truth. + let lats = vec![10.0, 50.0, -10.0, -70.0, -10.0]; + let lons = vec![45.0, 45.0, 170.0, 225.0, 280.0]; + let rings = build_rings(&[lats], &[lons], true); + let edges = build_edges(&rings, 6); + let cls = DirectClassifier::new(&edges, &rings); + // Seed base 0 centre (41.81, 45): exactly on the collinear edge's + // segment; the SoS perturbation resolves it east of the edge — the + // exterior wedge (consistent with base_fills' donor chain). + assert!( + !cls.classify(&cell_center_vec(0, 0), center_id(0, 0)), + "base 0 centre (on the collinear edge) must classify outside" + ); + // centre(3,1): truly inside (see NOTE) even though the winding backend + // reads it outside — pin both facts. + let c31 = cell_center_vec(1, 3); + assert!( + !parity_filled_robust(&c31, &rings), + "winding backend reads centre(3,1) outside (the antipodal-lens defect)" + ); + assert!( + cls.classify(&c31, center_id(1, 3)), + "centre(3,1) must classify inside (truth via transversal crossing parity)" + ); + // An off-plane centre whose antipode is exterior (the regime where the + // winding verdict is trustworthy): direct verdict must agree with it. + let c1 = cell_center_vec(0, 1); + assert_eq!( + cls.classify(&c1, center_id(0, 1)), + parity_filled_robust(&c1, &rings), + "direct verdict must match the winding oracle off the planes" + ); +} + +#[test] +fn test_fill_leg_parity_gates_collinear_leg() { + // #107 phase 1: the overlap gate. The diagnosed leg — base 0 centre + // (41.81, 45) → centre(3,1) (66.44, 45), both on the lon-45 near-plane + // shared with the ring's 40° collinear edge — must be flagged degenerate + // (None); a clean transversal leg must keep the chained parity and agree + // with arc_crossing_parity. + let lats = vec![10.0, 50.0, -10.0, -70.0, -10.0]; + let lons = vec![45.0, 45.0, 170.0, 225.0, 280.0]; + let rings = build_rings(&[lats], &[lons], true); + let edges = build_edges(&rings, 6); + let relevant: Vec = (0..edges.len()).collect(); + let (c00, c31) = (cell_center_vec(0, 0), cell_center_vec(1, 3)); + assert_eq!( + fill_leg_parity( + &c00, + &c31, + center_id(0, 0), + center_id(1, 3), + &relevant, + &edges + ), + None, + "collinear leg must be gated" + ); + // Clean leg: two off-plane points crossing the ring's lon-170 edge region. + let (p, q) = ( + latlon_to_unit_vec(0.0, 100.0), + latlon_to_unit_vec(0.0, 120.0), + ); + let (ip, iq) = (center_id(2, 5), center_id(2, 6)); + assert_eq!( + fill_leg_parity(&p, &q, ip, iq, &relevant, &edges), + Some(arc_crossing_parity(&p, &q, ip, iq, &relevant, &edges)), + "clean leg must keep the chained parity" + ); +} + +#[test] +#[ignore = "blocked on the #107 phase-1 finding: the debug parity oracle's \ +reference (parity_filled_robust) computes k(x) − k(−x), which is wrong \ +wherever a point and its antipode are both interior — this hemisphere+ \ +ring's lon-45 lens — so no descent-side fix can satisfy the oracle here; \ +needs the winding-backend decision recorded on issue #107"] fn test_descent_hemisphere_ring_collinear_edge_oracle() { // Runs the full descent under the debug parity oracle on the phase-2 - // review's adversarial ring. Un-ignore when the long-collinear-overlap - // consistency work lands. + // review's adversarial ring. Un-ignore when the winding backend is + // repaired for antipodally-interior points (see issue #107 / the phase-1 + // PR's "Questions for review"). let lats = vec![10.0, 50.0, -10.0, -70.0, -10.0]; let lons = vec![45.0, 45.0, 170.0, 225.0, 280.0]; let cov = polygon_to_morton_coverage(&lats, &lons, 6, true); diff --git a/src_rust/src/sphere/tests.rs b/src_rust/src/sphere/tests.rs index eddaafa..7973479 100644 --- a/src_rust/src/sphere/tests.rs +++ b/src_rust/src/sphere/tests.rs @@ -634,6 +634,37 @@ fn test_robust_pip_hemisphere_plus_band() { } } +#[test] +#[ignore = "known defect found in #107 phase 1: the short-way subtended-angle \ +sum is antisymmetric under x → −x (a and b project to the same vectors on the \ +plane ⊥x for both poles; only the sign term flips), so it computes \ +k(x) − k(−x) instead of the winding k(x) — every interior point whose \ +antipode is also interior reads outside. Sub-hemisphere interiors cannot \ +contain an antipodal pair, so only the hemisphere+ regime (#22) is affected. \ +Un-ignore when the winding backend is repaired (decision on issue #107)"] +fn test_point_in_ring_hemisphere_plus_antipodal_interior() { + // #22's flagship shape: the lat −10 band, CCW interior = everything north + // of −10 (hemisphere+). Every point with lat ∈ (−10, +10) has its + // antipode interior as well; the winding sum cancels to ~0 there and the + // backend calls a truly interior point outside. (The existing band test + // above only probes latitudes ≥ 13 or ≤ −11, and its cross-check oracle + // is another short-way winding sum — wrong in exactly the same region — + // which is how this slipped through.) + let band: Vec = (0..36) + .map(|k| latlon_to_unit_vec(-10.0, k as f64 * 10.0)) + .collect(); + for lon in [0.0, 33.0, 120.0, 260.0] { + for lat in [-5.0, 0.0, 5.0] { + assert!( + point_in_ring_robust(&latlon_to_unit_vec(lat, lon), &band), + "({lat},{lon}) is interior (north of the −10 band) but reads \ + outside: its antipode is interior too, so the antisymmetric \ + winding sum cancels" + ); + } + } +} + #[test] fn test_robust_pip_issue_11_meridian_box() { // #11 regression: the over-coverage flood was triggered by a polygon edge From 509bf651799337d63a29de23269ec419b54bf483 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 10 Jul 2026 12:06:35 +0000 Subject: [PATCH 2/2] address review on phase 1 of issue #107 --- src_rust/src/coverage.rs | 46 ++++++++++++++++++++++++++-------- src_rust/src/coverage/tests.rs | 5 +++- src_rust/src/sphere/tests.rs | 7 +++--- 3 files changed, 43 insertions(+), 15 deletions(-) diff --git a/src_rust/src/coverage.rs b/src_rust/src/coverage.rs index e0e3a7a..67c7524 100644 --- a/src_rust/src/coverage.rs +++ b/src_rust/src/coverage.rs @@ -748,11 +748,22 @@ fn chain_fill_from( /// A cell whose incoming fill leg overlaps a polygon edge's great circle /// ([`fill_leg_parity`] returns `None`) is never chained; it is classified /// directly instead — the same construction [`base_fills`] uses for seeds -/// whose centre lies on the boundary, generalized to any cell centre: pick an -/// unambiguous reference (a base centre off every edge plane, with a -/// trustworthy winding verdict) and run one [`arcs_cross_sos`] pass over all -/// edges ([`chain_fill_from`]). The invariant is one sentence: *degenerate -/// cells are classified directly, never chained.* +/// whose centre lies on the boundary, generalized to any cell centre: pick a +/// reference base centre clearly off every edge plane and run one +/// [`arcs_cross_sos`] pass over all edges ([`chain_fill_from`]). The +/// invariant is one sentence: *degenerate cells are classified directly, +/// never chained.* +/// +/// # Precondition on the reference verdict +/// +/// A donor's fill is the winding backend's verdict, so it is only as good as +/// [`parity_filled_robust`] at that centre — which is wrong wherever a point +/// *and its antipode* are both interior (the short-way winding sum computes +/// `k(x) − k(−x)`; the #107 phase-1 finding, pinned by the ignored +/// `test_point_in_ring_hemisphere_plus_antipodal_interior`). Being off the +/// edge planes does **not** rule that out: correctness here currently rests +/// on the donor being lens-free — exactly the caveat [`base_fills`]' seeds +/// already carry — until the winding backend is repaired. /// /// Donor discovery is lazy (`OnceLock`): clean geometry never gates, so it /// pays neither the per-edge normal normalization nor the seed PIPs. @@ -773,12 +784,22 @@ impl<'a> DirectClassifier<'a> { } } + /// A donor must sit **clearly** off every edge plane, not merely beyond + /// the [`ORIENT_EPS`] ambiguity cut: a donor barely past that cut would + /// put the reference chord in the same near-coincident regime the #107 + /// gate exists to escape. Degenerate lattice configurations sit at + /// ~1e-17 and clean base centres at O(0.1–1) relative to a polygon edge + /// plane, so 1e-9 (10³ × [`ORIENT_EPS`]) separates the two populations + /// with orders of magnitude to spare on both sides. + const DONOR_PLANE_MARGIN: f64 = 1e-9; + /// Fill at cell centre `x` (stable SoS id `ix`), classified directly from - /// the nearest unambiguous base-centre reference. Nearest keeps the - /// reference chord short; exact antipodality cannot occur (a sub-base cell - /// centre never coincides bit-exactly with a base centre or its antipode), - /// so the chord is always a valid minor arc for [`arcs_cross_sos`]. When - /// every base centre is ambiguous (boundary planes through all 12 — + /// the nearest clearly-off-plane base-centre reference + /// ([`Self::DONOR_PLANE_MARGIN`]). Nearest keeps the reference chord + /// short; exact antipodality cannot occur (a sub-base cell centre never + /// coincides bit-exactly with a base centre or its antipode), so the + /// chord is always a valid minor arc for [`arcs_cross_sos`]. When no + /// base centre clears the margin (boundary planes through all 12 — /// pathological), keep the raw winding verdict, as [`base_fills`] does. fn classify(&self, x: &Vec3, ix: PointId) -> bool { let donors = self.donors.get_or_init(|| { @@ -786,7 +807,10 @@ impl<'a> DirectClassifier<'a> { (0..12u64) .filter_map(|b| { let c = cell_center_vec(0, b); - seed_fill(&c, &units, self.rings).map(|f| (c, center_id(0, b), f)) + units + .iter() + .all(|n| dot(n, &c).abs() >= Self::DONOR_PLANE_MARGIN) + .then(|| (c, center_id(0, b), parity_filled_robust(&c, self.rings))) }) .collect() }); diff --git a/src_rust/src/coverage/tests.rs b/src_rust/src/coverage/tests.rs index e5bbcdb..997cf3b 100644 --- a/src_rust/src/coverage/tests.rs +++ b/src_rust/src/coverage/tests.rs @@ -664,11 +664,14 @@ fn test_fill_leg_parity_gates_collinear_leg() { "collinear leg must be gated" ); // Clean leg: two off-plane points crossing the ring's lon-170 edge region. + // These are not cell centres, so give them the one-shot derived-point ids + // (CORNER_ID_A/B): SoS identities only need to be pairwise distinct within + // a call, and nothing chains from this leg, so no stable id is required. let (p, q) = ( latlon_to_unit_vec(0.0, 100.0), latlon_to_unit_vec(0.0, 120.0), ); - let (ip, iq) = (center_id(2, 5), center_id(2, 6)); + let (ip, iq) = (CORNER_ID_A, CORNER_ID_B); assert_eq!( fill_leg_parity(&p, &q, ip, iq, &relevant, &edges), Some(arc_crossing_parity(&p, &q, ip, iq, &relevant, &edges)), diff --git a/src_rust/src/sphere/tests.rs b/src_rust/src/sphere/tests.rs index 7973479..272173f 100644 --- a/src_rust/src/sphere/tests.rs +++ b/src_rust/src/sphere/tests.rs @@ -647,9 +647,10 @@ fn test_point_in_ring_hemisphere_plus_antipodal_interior() { // of −10 (hemisphere+). Every point with lat ∈ (−10, +10) has its // antipode interior as well; the winding sum cancels to ~0 there and the // backend calls a truly interior point outside. (The existing band test - // above only probes latitudes ≥ 13 or ≤ −11, and its cross-check oracle - // is another short-way winding sum — wrong in exactly the same region — - // which is how this slipped through.) + // above does probe this co-band — lat 0 and −9 — but only against + // `winding_inside`, an identical short-way winding sum: both sides are + // wrong in exactly the same region, so the comparison passes and the + // defect slipped through. This test asserts ground truth instead.) let band: Vec = (0..36) .map(|k| latlon_to_unit_vec(-10.0, k as f64 * 10.0)) .collect();