Direct classification for degenerate descent cells (issue #107)

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<bool> {
+    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,111 @@ fn seed_fill(x: &Vec3, unit_normals: &[Vec3], rings: &[Vec<Vec3>]) -> Option<boo
     Some(parity_filled_robust(x, rings))
 }
 
+/// Fill at `x` (stable SoS id `ix`) classified **directly** from an
+/// already-classified reference point `r`: `r`'s fill XOR the crossing parity
+/// of the single chord `r → x` over **all** edges.  The full symbolic
+/// predicate is unconditional: reference chords are long (up to ~122° between
+/// base centres), where the bare two-straddle fast path of
+/// [`edge_crosses_probe`] can also fire on the *far* intersection of the two
+/// great circles (an edge antipodal to the chord) and silently invert the
+/// verdict.  Off every hot path — this runs only for on-boundary seeds
+/// ([`base_fills`]) and #107-gated degenerate cells ([`DirectClassifier`]).
+fn chain_fill_from(
+    r: &Vec3,
+    ir: PointId,
+    r_fill: bool,
+    x: &Vec3,
+    ix: PointId,
+    edges: &[Edge],
+) -> 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 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.
+#[cfg_attr(not(test), allow(dead_code))] // wired into the descent in phase 2 (#107)
+struct DirectClassifier<'a> {
+    edges: &'a [Edge],
+    rings: &'a [Vec<Vec3>],
+    donors: OnceLock<Vec<(Vec3, PointId, bool)>>,
+}
+
+#[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<Vec3>]) -> Self {
+        Self {
+            edges,
+            rings,
+            donors: OnceLock::new(),
+        }
+    }
+
+    /// 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 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(|| {
+            let units: Vec<Vec3> = self.edges.iter().map(|e| normalize(&e.n_ab)).collect();
+            (0..12u64)
+                .filter_map(|b| {
+                    let c = cell_center_vec(0, b);
+                    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()
+        });
+        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 +896,14 @@ fn base_fills(edges: &[Edge], rings: &[Vec<Vec3>], 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(
-                    &centers[d as usize],
-                    &centers[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(
+            &centers[d as usize],
+            center_id(0, d),
+            fill[d as usize],
+            &centers[b as usize],
+            center_id(0, b),
+            edges,
+        );
         known[b as usize] = true;
     }
     fill

src_rust/src/coverage/tests.rs

@@ -589,15 +589,107 @@ 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<usize> = (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.
+    // 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) = (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)),
+        "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);

src_rust/src/sphere/tests.rs

@@ -634,6 +634,38 @@ 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 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<Vec3> = (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