From 0df34e9d424506173754a10807956c2516b673cb Mon Sep 17 00:00:00 2001 From: sth-v Date: Fri, 25 Sep 2026 16:16:45 +0300 Subject: [PATCH 1/3] Restore complete curved overlap regions at CAD tolerance --- .../2026-09-25-case17-curved-overlap.md | 104 ++++ examples/ssx/case_17.py | 52 ++ mmcore/numeric/intersection/csx/_bez_csx4.py | 209 +++++--- mmcore/numeric/intersection/ssx/_bez_ssx5.py | 25 +- .../numeric/intersection/ssx/_ssx5_overlap.py | 277 +++++++++- .../ssx/_ssx_spherical_overlap.py | 496 ++++++++++++++++++ tests/test_csx_cad_span_contacts.py | 204 +++++++ tests/test_csx_overlap_normal_gap.py | 28 +- tests/test_ssx_curved_overlap_corners.py | 135 +++++ tests/test_ssx_spherical_overlap.py | 287 ++++++++++ tests/test_ssx_spherical_reduction.py | 246 +++++++++ 11 files changed, 1975 insertions(+), 88 deletions(-) create mode 100644 docs/superpowers/reviews/2026-09-25-case17-curved-overlap.md create mode 100644 examples/ssx/case_17.py create mode 100644 mmcore/numeric/intersection/ssx/_ssx_spherical_overlap.py create mode 100644 tests/test_csx_cad_span_contacts.py create mode 100644 tests/test_ssx_curved_overlap_corners.py create mode 100644 tests/test_ssx_spherical_overlap.py create mode 100644 tests/test_ssx_spherical_reduction.py diff --git a/docs/superpowers/reviews/2026-09-25-case17-curved-overlap.md b/docs/superpowers/reviews/2026-09-25-case17-curved-overlap.md new file mode 100644 index 00000000..c305ee85 --- /dev/null +++ b/docs/superpowers/reviews/2026-09-25-case17-curved-overlap.md @@ -0,0 +1,104 @@ +# Case 17: complete curved overlap boundaries + +`examples/ssx/case_17.py` contains two partially overlapping spherical +NURBS patches. At `atol=1e-3`, their common boundary is a closed spherical +quadrilateral: two arcs from each surface. The failing result contained +only the two arcs from the second surface, stitched into an open L, and +no overlap region. + +## Causes and corrections + +1. CSX rejected the first surface's two coincident edges because their + signed normal gaps changed sign. The deviations were about `1e-7`, far + below the requested `1e-3`. Default CAD overlap classification now + accepts such spans and retains paired contact proposals for SSX to + refine. Explicit exact-mode behavior remains separate. Normal-gap + roundoff uses Cartesian geometry, so homogeneous weight scaling does + not change which contacts survive. +2. Curved-rim discovery found all four edges independently, but stopped + each endpoint at its own tolerance fringe. Two descriptions of one + corner were `0.00365` apart; the loop builder discarded the open chain. + Incident domain edges now resolve one shared paired corner. Parallel + edges, different parameter faces, and remote contact spans are not + joined by enlarging a spatial radius. +3. Seventeen fixed rim samples were insufficient for these curved arcs. + Adaptive subdivision now bounds both lifted parameter chords against + the reported spatial chord over the whole interval, using positive + homogeneous Bernstein weights. Shared poles retain their distinct + parameter representatives. Initial samples, corrected corners, and + adaptive samples consistently use paired surface midpoints; mixing + those with owner-only points had rejected valid gaps close to `atol`. +4. A narrow internal CSX overlap fell entirely between the global sample + locations even though its two boundary contacts were already known. + Those contacts now seed a candidate interval that undergoes the same + membership and domain-end checks. The captured cut changed from + 22.54 seconds / 67,104 work units to 0.233 seconds / 1,126 units. + +Correcting the boundary still left the ordinary one-dimensional search +subdividing a two-dimensional coincident area. A separate geometric +reduction now handles supported rational quadratic spherical octants. +The fitted template is only a proposal: a whole-patch rational difference +bound must place each actual chart within the allocated CAD tolerance of +its ideal chart. A complete analytic hemisphere census is then matched +against the actual assembled boundary, including all corners, cyclic +ordering, monotone parameter intervals, and whole-chord accuracy. +The distance bound also accounts for error amplification at nearly +opposed domain boundaries; an uncertain footprint uses general SSX. + +This reduction uses the original surfaces for returned parameters and +geometry. It makes no early empty-result claim. Other charts continue +through general SSX and receive the curved-rim and CSX fixes above. A +surface with identical boundary curves and one matching interior point, +but an interior bulge larger than tolerance, cannot select the spherical +reduction. Small within-tolerance control and weight perturbations can. + +Shared spherical poles need geometric treatment as well: different +longitude parameters at a pole are not separate cusp curves or crossing +branches. A Bernstein normal-hull check excludes rank defects outside a +small pole cap. The entire cap must fit within a reported CAD corner, or +be separated from the other actual patch. Only then is the ordinary C1 +enumeration unnecessary. The validated simple region boundary also +resolves C3; other reductions retain their existing singularity checks. + +The first integrated NURBS run returned one region, four referenced rim +branches, no isolated points or singularities, and `complete=True` in +about 0.31 seconds. The failing version took about 81 seconds on the same +machine. These are single-run diagnostic timings, not a general speed claim. + +## Tests + +The spherical reference comes from independently intersecting six inward +hemispheres. Case 17 has perimeter approximately `26.16026806` and area +`44.90430351`; both collapsed poles lie outside its common region. Tests +check every boundary component, area and perimeter, both source images, +intermediate chord points, and the absence of leftover branches or points. + +Generated cases cover six-sided intersections, swapping and reversing +charts, scaling, a narrow overlap with a shared pole, homogeneous gauge, +separated patches, and lower-dimensional contacts. Additional controls +cover hidden interior bulges, weight-only distortions, ordinary nonspherical +patches with modest work allowances, and interrupted output publication. + +The former test requiring a line at signed distance `+/-1e-7` to be rejected +as a CAD overlap at `1e-3` was inappropriate. It now requires the overlap +and retained contact proposals. Above-tolerance crossings and explicit +exact-mode roots remain checked. The new narrow-strip test has an +independent nonlinear parameter correspondence; its interval is only +`0.0005` in curve parameters but spans fifty model tolerances. + +Contact proposals stored inside overlap records share the result allowance +with isolated roots. A denied CAD-classification step reports the work +stop and retains already found geometry. Completed paired rims also survive +work denial during spherical-region assembly, without invalid region indices. + +## Historical comparison + +Using the exact supplied rounded fixture and the same native extensions, +the open L also reproduced at `76735f9`; the pre-PR-48 `552d972` version +returned no geometry. The working version of this particular example has +therefore not been identified. The corner and normal-side assumptions +predate the latest merge; this investigation does not attribute the whole +failure to PR #48. + +Final validation is recorded in the repair pull request. Focused runs and +the timing above do not replace the full-suite and GitHub checks. diff --git a/examples/ssx/case_17.py b/examples/ssx/case_17.py new file mode 100644 index 00000000..2f50f303 --- /dev/null +++ b/examples/ssx/case_17.py @@ -0,0 +1,52 @@ +import numpy as np +from mmcore.nurbs._nurbs_eval import NURBSSurfaceTuple +from mmcore.numeric.intersection.ssx import nurbs_ssx + +s1 = NURBSSurfaceTuple( + order_u=3, + order_v=3, + knot_u=np.array([0., 0., 0., 1., 1., 1.]), + knot_v=np.array([0., 0., 0., 1., 1., 1.]), + control_points=np.array([[[ 27.72140823, -10.76294763, -2.32132869], + [ 31.15136161, -9.17865474, 3.77780409], + [ 30.12306049, -15.88894641, 6.09913285]], + + [[ 21.50433486, -8.77901557, 0.65961025], + [ 24.93428823, -7.19472268, 6.75874302], + [ 30.12306049, -15.88894641, 6.09913285]], + + [[ 20.47603374, -15.48930724, 2.980939 ], + [ 23.90598711, -13.90501435, 9.08007178], + [ 30.12306049, -15.88894641, 6.09913285]]]), + weights=np.array([[1. , 0.70710678, 1. ], + [0.70710678, 0.5 , 0.70710678], + [1. , 0.70710678, 1. ]]) +) +import numpy as np +from mmcore.nurbs._nurbs_eval import NURBSSurfaceTuple + + +s2 = NURBSSurfaceTuple( + order_u=3, + order_v=3, + knot_u=np.array([0., 0., 0., 1., 1., 1.]), + knot_v=np.array([0., 0., 0., 1., 1., 1.]), + control_points=np.array([[[ 26.69310698, -10.29870217, -0. ], + [ 26.69310698, -10.29870217, 7.17453725], + [ 26.69310698, -17.47323942, 7.17453725]], + + [[ 19.51856973, -10.29870217, -0. ], + [ 19.51856973, -10.29870217, 7.17453725], + [ 26.69310698, -17.47323942, 7.17453725]], + + [[ 19.51856973, -17.47323942, -0. ], + [ 19.51856973, -17.47323942, 7.17453725], + [ 26.69310698, -17.47323942, 7.17453725]]]), + weights=np.array([[1. , 0.70710678, 1. ], + [0.70710678, 0.5 , 0.70710678], + [1. , 0.70710678, 1. ]]) +) + +if __name__ == "__main__": + result = nurbs_ssx(s1, s2, atol=1e-3) + print(result) diff --git a/mmcore/numeric/intersection/csx/_bez_csx4.py b/mmcore/numeric/intersection/csx/_bez_csx4.py index 2d46963a..fa850b9f 100644 --- a/mmcore/numeric/intersection/csx/_bez_csx4.py +++ b/mmcore/numeric/intersection/csx/_bez_csx4.py @@ -730,22 +730,24 @@ def _certify_affine_csx_overlap(C, S, a, b, rational): def _tolerance_csx_overlap_certificate(C, S, atol, rational, ptol_t, strict_context, n_samples=65, - on_dense=None): + on_dense=None, on_contact_work=None, + max_contacts=None, parameter_range=(0., 1.)): """Modeling-tolerance curve-on-surface span classification. Endpoint evidence, domain-pinned projections and sampled interior - witnesses classify within-atol spans. A sampled normal-side sign flip - refuses promotion, but finite sampling cannot guarantee separation of - every exact crossing. The legacy 'exact' label here describes the + witnesses classify within-atol spans. Normal-side sign changes do not + invalidate CAD coincidence; their paired contact proposals are retained + for consumers that need to refine individual crossings. Finite sampling + does not promise a census of every algebraic root. The 'exact' label describes the sampled residual tier, not an exact polynomial identity certificate. Exact-topology callers bypass this helper and use source identities. Returns overlap dictionaries (possibly multiple clipped spans), or None when no modeling span can be classified. """ - diag = float(np.linalg.norm( - np.max(S.reshape(-1, S.shape[-1]), axis=0) - - np.min(S.reshape(-1, S.shape[-1]), axis=0))) or 1.0 + cartesian = S[..., :-1] / S[..., -1:] if rational else S + points = np.asarray(cartesian).reshape(-1, cartesian.shape[-1]) + diag = float(np.linalg.norm(points.max(axis=0) - points.min(axis=0))) or 1.0 tiny = max(4096.0 * float(np.finfo(np.float64).eps), 1e-12) * max(1.0, diag) edge_pad = 1e-9 @@ -763,7 +765,8 @@ def _member(t, seed_uv): return p, u_cl, v_cl, float(dist) # --- (1) arming: coarse scan for ANY on-surface stretch ------------- - coarse = np.linspace(0.0, 1.0, 17) + parameter_low, parameter_high = map(float, parameter_range) + coarse = np.linspace(parameter_low, parameter_high, 17) seed = (0.5, 0.5) hits = [] for t in coarse: @@ -782,7 +785,7 @@ def _member(t, seed_uv): on_dense() # --- (2) dense membership + witnesses -------------------------------- - ts = np.linspace(0.0, 1.0, int(n_samples)) + ts = np.linspace(parameter_low, parameter_high, int(n_samples)) inside = np.zeros(len(ts), dtype=bool) res = np.zeros(len(ts)) uvs = np.zeros((len(ts), 2)) @@ -830,6 +833,7 @@ def _refine(t_in, t_out, seed_uv): return t_in, seed_uv overlaps = [] + contact_count = 0 for k0, k1 in runs: if k1 == k0: # Single-sample runs are the CORNER-CONTACT signature (the @@ -863,30 +867,54 @@ def _pinned(t_end, outward): # Probe one grid step OUTWARD: a domain-clipped span clamps to # an edge there; a genuine interior fade-out (the offset-twin # signature, which must stay refused) does not. - t_probe = min(1.0, max(0.0, t_end + outward / (len(ts) - 1.0))) + probe_step = (parameter_high - parameter_low) / (len(ts) - 1.) + t_probe = min(1.0, max(0.0, t_end + outward * probe_step)) _p, u, v, _d = _member(t_probe, (u, v)) return (min(u, 1.0 - u) <= 1e-6 or min(v, 1.0 - v) <= 1e-6) if not (_pinned(t_lo, -1.0) and _pinned(t_hi, +1.0)): continue - # (3) flip guard on the normal gaps INSIDE the span: root-like - # samples (abs(signed) <= tiny) are bridged; consecutive gap samples with - # opposite normal-side signs = crossing structure -> refuse. - # END-ADJACENT flips are exempt: a genuine touch AT a pinned span - # end (real-world-inexact data sits above the roundoff floor, so - # bridging cannot cover it) is the span's own endpoint root — the - # theorem terminates the overlap there anyway. INTERIOR flips - # still refuse (the never-merge-crossings invariant). - # A stopped projection can retain a tangential residual larger - # than this floor even on a planar overlap. Its normal component - # can still be pure rounding noise, so the total distance must - # not make that component's sign meaningful (or BLAS rounding - # alone can reject a full overlap). + # A within-atol span remains CAD-coincident when rounding puts + # samples on both sides of the target surface. Preserve bracketed + # contacts separately so SSX can refine them on its owning face; + # an overlap's unpaired UV ranges cannot stand in for those seeds. + # Use normal gaps only: tangential projection lag must not amplify + # an otherwise meaningless roundoff normal sign. idx = [k for k in range(k0, k1 + 1) if abs(signed[k]) > tiny] - flip = any(signed[a] * signed[b] < 0.0 - for a, b in zip(idx, idx[1:]) - if a != k0 and b != k1) - if flip: - continue + contacts = [] + contact_results_exhausted = False + for a, b in zip(idx, idx[1:]): + if signed[a] * signed[b] >= 0.0: + continue + if max_contacts is not None and contact_count >= max_contacts: + contact_results_exhausted = True + break + if on_contact_work is not None and not on_contact_work(1): + break + fraction = signed[a] / (signed[a] - signed[b]) + seed_t = float(ts[a] + fraction * (ts[b] - ts[a])) + seed_uv = uvs[a] + fraction * (uvs[b] - uvs[a]) + t, u, v, _residual, root_ok = _polish_csx_root( + C, S, seed_t, *seed_uv, rational, strict_context) + if not (root_ok and ts[a] <= t <= ts[b] + and 0. <= u <= 1. and 0. <= v <= 1.): + if on_contact_work is not None and not on_contact_work(1): + break + t = seed_t + point, u, v, distance = _member(t, seed_uv) + else: + point = eval_curve(C, t, rational=rational) + distance = float(np.linalg.norm( + point - eval_surface(S, u, v, rational=rational))) + if (np.isfinite([t, u, v, distance]).all() + and np.isfinite(point).all() and distance <= atol): + # This is a source-validated proximity proposal. The SSX + # consumer applies its fixed-face root correction before + # registering a continuation crossing, even if this first + # polishing step already reached numerical equality. + contacts.append(dict(t=float(t), u=float(u), v=float(v), + point=np.asarray(point), + certification='tolerance', d_min=distance)) + contact_count += 1 span_res = float(res[k0:k1 + 1].max()) cert = "exact" if span_res <= tiny else "tolerance" u_rng = (float(uvs[k0:k1 + 1, 0].min()), @@ -902,6 +930,10 @@ def _pinned(t_end, outward): "certification": cert, "residual_max": span_res, }) + if contacts: + overlaps[-1]['boundary_contacts'] = contacts + if contact_results_exhausted: + overlaps[-1]['_contact_results_exhausted'] = True return overlaps or None @@ -2162,8 +2194,9 @@ def bez_csx( Maximum total cells shared by boundary analysis, nested CCX calls, and Phase 2 (safety limit). max_results : int - Maximum isolated roots materialized before returning a partial - result. Positive-dimensional sets must be classified separately; + Maximum contact records materialized, including isolated roots + retained on overlaps, before returning a partial result. + Positive-dimensional sets must be classified separately; this cap prevents an unrecognized set from turning dedup quadratic. tolerance_tier : bool Allow curve/surface contacts and overlaps within ``atol``. @@ -2409,6 +2442,12 @@ def _append_fiber(candidate): t_exclude = [] # t-intervals to cut from the curve + def _contact_count(): + # Moving a root onto an overlap preserves its materialization + # cost; it must not create new allowance for the remaining search. + return len(isolated) + sum(len(span.get('boundary_contacts', ())) + for span in overlaps) + def _absorb_boundary_contacts(overlap, low, high): # A tolerance span can surround an isolated boundary contact. # Keep its already-polished parameters for consumers such as SSX @@ -2416,7 +2455,7 @@ def _absorb_boundary_contacts(overlap, low, high): # surrounding tolerance neighborhood as an overlap. These entries # were admitted under max_results; no new roots are generated here. nonlocal isolated - contacts, remaining = [], [] + contacts, remaining = list(overlap.get('boundary_contacts', ())), [] for entry in isolated: if low <= entry['t'] <= high: contacts.append({key: entry[key] @@ -2442,7 +2481,7 @@ def _absorb_boundary_contacts(overlap, low, high): # boundary. Newton refinement converges to the real root if one # exists nearby; otherwise the candidate is rejected. for bz in csx_boundary_zeros: - if len(isolated) >= max_results: + if _contact_count() >= max_results: budget_exhausted = True obligations_complete = False truncation_cause = truncation_cause or 'results' @@ -2624,46 +2663,76 @@ def _absorb_boundary_contacts(overlap, low, high): for o in overlaps] + unresolved_roots, } - # L59 (USER DECISION: theorem-first tier): certify tolerance/exact - # overlap spans whenever cheap endpoint evidence arms — a curve t-end - # lying ON the surface (an axis-0 boundary zero) or a valley pair — - # INDEPENDENTLY of the exact-affine identity. The real-data class has - # only ONE boundary zero (the other span end pins on the uv-domain - # edge in the projected sense: measured 6.7e-4 clearance from the edge - # line — no exact 3-D root exists there for the boundary phase to - # find). Certified spans bypass the Phase-2 grind entirely. - # Arming: a curve-end zero on-surface, a valley pair, OR a pair with - # NO boundary zeros at all (measured on user data: a coincident - # stretch that enters AND exits through patch edges with sub-atol - # clearance produces zero exact boundary roots — the certificate must - # still get its 17-sample look; transversal SSX-nested calls always - # carry boundary zeros, so their cost profile is untouched). - if (tolerance_tier and not boundary_exhausted and not overlaps - and (_valley_pair_seen - or not csx_boundary_zeros - or any(bz.axis == 0 for bz in csx_boundary_zeros))): - # Split pricing (§11.5 drift-gate lesson): the 17-projection - # arming scan is billed always; the dense pass (65 witnesses + - # refines) only when the scan HITS — nested cut-face calls - # routinely arm-and-miss and must not pay the full tier. - _tier_price = 17 - _dense_price = 65 + 80 - if cells.remaining > _tier_price + _dense_price: - cells.spend(_tier_price) - _tol_overlaps = _tolerance_csx_overlap_certificate( + # Classify CAD-coincident spans independently of affine identity. + # A whole-curve grid alone can miss a short domain-clipped interval, + # even when the boundary pass has already found both of its ends. + if tolerance_tier and not boundary_exhausted and not overlaps: + tier_price = 17 + dense_price = 65 + 80 + seed_contacts = sorted(isolated, key=lambda point: point['t']) + + def _contact_work(amount): + nonlocal budget_exhausted, truncation_cause + if amount > cells.remaining: + budget_exhausted = True + truncation_cause = truncation_cause or 'cells' + return False + cells.spend(amount) + return True + + def _sample_span(parameter_range): + nonlocal budget_exhausted, truncation_cause + nonlocal non_affine_overlap_span + if cells.remaining <= tier_price + dense_price: + budget_exhausted = True + truncation_cause = truncation_cause or 'cells' + return + cells.spend(tier_price) + candidates = _tolerance_csx_overlap_certificate( C, S, atol, rational, ptol_t, strict_root_tol, + parameter_range=parameter_range, on_dense=lambda: cells.spend( - min(_dense_price, max(0, cells.remaining)))) - if _tol_overlaps: - for _o in _tol_overlaps: - overlaps.append(_o) - _t_lo, _t_hi = _o["t_range"] - t_exclude.append((_t_lo - ptol_t, _t_hi + ptol_t)) - _absorb_boundary_contacts(_o, _t_lo - atol, _t_hi + atol) - # the span(s) are certified: the L42 uncertified-span - # fallback no longer applies to them. + min(dense_price, max(0, cells.remaining))), + on_contact_work=_contact_work, + max_contacts=max(0, int(max_results) - _contact_count())) + for overlap in candidates or (): + if overlap.pop('_contact_results_exhausted', False): + budget_exhausted = True + truncation_cause = truncation_cause or 'results' + overlaps.append(overlap) + low, high = overlap['t_range'] + t_exclude.append((low - ptol_t, high + ptol_t)) + _absorb_boundary_contacts(overlap, low - atol, high + atol) non_affine_overlap_span = None + if (_valley_pair_seen or not csx_boundary_zeros + or any(point.axis == 0 for point in csx_boundary_zeros)): + _sample_span((0., 1.)) + + # Consecutive validated boundary contacts bound candidate spans. + # A cheap middle projection proposes coincidence; the same dense + # membership and domain-pinning checks must still validate it. + # This also covers entering and leaving through the same UV face. + for left, right in zip(seed_contacts, seed_contacts[1:]): + low, high = float(left['t']), float(right['t']) + if high - low < 4. * ptol_t: + continue + if any(a - ptol_t <= low and high <= b + ptol_t + for a, b in (overlap['t_range'] for overlap in overlaps)): + continue + if cells.remaining <= tier_price + dense_price: + budget_exhausted = True + truncation_cause = truncation_cause or 'cells' + break + cells.spend(1) + middle = .5 * (low + high) + point = eval_curve(C, middle, rational=rational) + _u, _v, distance = _project_point_on_surface( + point, S, .5 * (left['u'] + right['u']), + .5 * (left['v'] + right['v']), atol, rational) + if np.isfinite(distance) and distance <= atol: + _sample_span((low, high)) + # =================================================================== # PHASE 2: Isolated intersection search on remaining curve intervals # =================================================================== @@ -2732,7 +2801,7 @@ def _interval_hits_span(a, b): non_affine_overlap_span is not None and any(not _interval_hits_span(a, b) for a, b in t_intervals[_ii:] if (b - a) >= ptol_t)) - if cells.remaining <= 0 or len(isolated) >= max_results: + if cells.remaining <= 0 or _contact_count() >= max_results: budget_exhausted = True truncation_cause = truncation_cause or ( "cells" if cells.remaining <= 0 else "results") @@ -2757,7 +2826,7 @@ def _interval_hits_span(a, b): t_lo, t_hi, atol, rational, ptol_t, ptol_u, ptol_v, known_points=isolated, max_depth=max_depth, max_cells=phase2_cell_limit, - max_results=max_results - len(isolated), + max_results=max_results - _contact_count(), exact_topology=exact_topology, unresolved_root_boxes=search_unresolved_boxes, obligation_diagnostics=phase2_obligations, diff --git a/mmcore/numeric/intersection/ssx/_bez_ssx5.py b/mmcore/numeric/intersection/ssx/_bez_ssx5.py index f48b76b4..0182a04f 100644 --- a/mmcore/numeric/intersection/ssx/_bez_ssx5.py +++ b/mmcore/numeric/intersection/ssx/_bez_ssx5.py @@ -309,11 +309,15 @@ def _distance_net_lower_bound(S1, S2, F, rational): def _try_isoline_intersection(first, second, atol, budget): - """Use planar or one-variable reductions when the surface charts allow it.""" + """Use planar, spherical, or one-variable reductions when available.""" from mmcore.numeric.intersection.ssx._ssx_isolines import try_isoline_ssx, try_plane_coincidence from mmcore.numeric.intersection.ssx._ssx_matched_isolines import try_matched_isoline_ssx from mmcore.numeric.intersection.ssx._ssx5_overlap import try_planar_intersection + from mmcore.numeric.intersection.ssx._ssx_spherical_overlap import try_spherical_overlap result = try_planar_intersection(first, second, atol, budget) + if result is not None or budget.exhausted: + return result + result = try_spherical_overlap(first, second, atol, budget) if result is not None or budget.exhausted: return result result = try_plane_coincidence(first, second, atol, budget) @@ -6953,7 +6957,7 @@ def _discover_c1_singularities(S1_h, S2_h, atol, ptol4, budget, existing=()): def _link_c1_and_discover_c3(S1_h_top, S2_h_top, all_branches, all_singularities, c1_singularities, - atol, ptol4_global, budget): + atol, ptol4_global, budget, *, check_c3=True): """Apply the same singularity postprocessing to every geometry path.""" # Attach C1 links after branch filters; the singularities themselves # have already been published by their discovery path. @@ -6981,6 +6985,9 @@ def _link_c1_and_discover_c3(S1_h_top, S2_h_top, all_branches, links.append((bi, vertex)) singularity.branch_links = links + if not check_c3: + return + # --- C3 pass (paper §5.4): 3D self-intersections of the SSI image --- # Runs AFTER tracing (branch geometry drives the candidate search). # Ledger L8: the old trigger — "some traced cell failed the per-cell @@ -7289,14 +7296,20 @@ def _result(branches=None, points=None, singularities=None, *bez_surface_param_tolerance(S2_h_top, atol, rational=True), ], dtype=np.float64), 1e-9) singularities = isolines.setdefault("singularities", []) - c1 = _discover_c1_singularities( - S1_h_top, S2_h_top, atol, ptol4, budget, - existing=singularities) + # A validated spherical area can also resolve its source + # regularity outside point-confined pole fibres. Their distinct + # parameter representatives belong to one geometric boundary + # point, not extra cusp curves or rim self-intersections. + c1 = ([] if getattr(isolines, '_c1_resolved', False) else + _discover_c1_singularities( + S1_h_top, S2_h_top, atol, ptol4, budget, + existing=singularities)) singularities.extend(g for g in c1 if not any(g is known for known in singularities)) _link_c1_and_discover_c3( S1_h_top, S2_h_top, isolines.get("branches", []), - singularities, c1, atol, ptol4, budget) + singularities, c1, atol, ptol4, budget, + check_c3=not getattr(isolines, '_c3_resolved', False)) return _result(**isolines) # The current distance-net constructor forms a pairwise Gram tensor over diff --git a/mmcore/numeric/intersection/ssx/_ssx5_overlap.py b/mmcore/numeric/intersection/ssx/_ssx5_overlap.py index eefd2b34..71a2b4cb 100644 --- a/mmcore/numeric/intersection/ssx/_ssx5_overlap.py +++ b/mmcore/numeric/intersection/ssx/_ssx5_overlap.py @@ -14,12 +14,13 @@ from __future__ import annotations from dataclasses import dataclass, field +from math import comb import numpy as np from numpy.typing import NDArray from mmcore.numeric._bezier_common import ( - eval_curve, eval_surface, eval_surface_d1, + eval_curve, eval_surface, eval_surface_d1, restrict_net_axis_v, ) from mmcore.numeric.intersection.ssx._ssx_substrate import SSXBranch @@ -147,6 +148,53 @@ def value(u, v): return float(np.max(np.linalg.norm(elevated-line, axis=1))/np.min(c[:, 3])) +def _surface_roundoff(surface, xyz=()): + """Arithmetic scale for closed CAD-distance comparisons.""" + net = np.asarray(surface, dtype=float) + coordinates = net[..., :3]/net[..., 3, None] + return 128.*np.finfo(float).eps*max( + 1., float(np.max(abs(coordinates))), + float(np.max(np.abs(xyz))) if np.size(xyz) else 0.) + + +def _surface_chord_error(surface, start, end, xyz): + """Bound an arbitrary-degree rational chart over one affine UV chord. + + Restriction to the UV rectangle followed by its diagonal gives a + degree p+q homogeneous curve. Subtracting the reported XYZ chord times + its weight gives a Bernstein numerator bound over the whole segment. + """ + net = np.asarray(surface, dtype=float).copy() + if not np.isfinite(net).all() or np.min(net[..., 3]) <= 0.: + return np.inf + net /= np.max(net[..., 3]) + origin = np.asarray(xyz[0], dtype=float) + net[..., :3] -= origin * net[..., 3, None] + for axis in range(2): + a, b = float(start[axis]), float(end[axis]) + if not (0. <= a <= 1. and 0. <= b <= 1.): + return np.inf + net = restrict_net_axis_v(net, axis, min(a, b), max(a, b), 0., 1.) + if b < a: + net = np.flip(net, axis=axis) + p, q = np.asarray(net.shape[:2])-1 + degree = p+q + curve = np.zeros((degree+1, 4)) + for i in range(p+1): + for j in range(q+1): + curve[i+j] += (comb(p, i)*comb(q, j)/comb(degree, i+j))*net[i, j] + if not np.isfinite(curve).all() or np.min(curve[:, 3]) <= 0.: + return np.inf + fraction = np.arange(degree+2)/(degree+1.) + elevated = np.zeros((degree+2, 3)) + elevated[:-1] += (1.-fraction[:-1, None])*curve[:, :3] + elevated[1:] += fraction[1:, None]*curve[:, :3] + chord = np.zeros_like(elevated) + chord[1:] = (fraction[1:]*curve[:, 3])[:, None]*(np.asarray(xyz[1])-origin) + bound = float(np.max(np.linalg.norm(elevated-chord, axis=1))/np.min(curve[:, 3])) + return bound+_surface_roundoff(surface, xyz) + + def _planar_convex_rims(S1_h, S2_h, atol, charge, context, *, include_boundary_contacts=False): """Clip convex bilinear chart images before inverting their shared rims. @@ -422,6 +470,9 @@ def _sample_edge_rims(S_own_h, S_other_h, owner, axis, side, atol, Returns a list of rims, each a dict with synchronized arrays: ``stuv`` (N,4) global parameters, ``xyz`` (N,3), ``resid_max``. """ + membership_tol = atol+_surface_roundoff(S_own_h)+_surface_roundoff(S_other_h) + if not np.isfinite(membership_tol): + return [] ws = np.linspace(0.0, 1.0, coarse) onsurf = np.zeros(coarse, dtype=bool) inv = np.zeros((coarse, 2), dtype=np.float64) @@ -432,8 +483,8 @@ def _sample_edge_rims(S_own_h, S_other_h, owner, axis, side, atol, u, v, resid = _invert_point(S_other_h, xyz, seed=seed) # A wandering warm start can hand the next sample a foreign basin # after an off-surface stretch; only propagate on-surface seeds. - seed = (u, v) if resid <= atol else None - onsurf[k] = resid <= atol + seed = (u, v) if resid <= membership_tol else None + onsurf[k] = np.isfinite(resid) and resid <= membership_tol inv[k] = (u, v) rims = [] @@ -464,13 +515,14 @@ def _refine_end(S_own_h, S_other_h, axis, side, atol, w_in, w_out, seed): w_out = min(1.0, max(0.0, w_out)) if w_in == w_out: return w_in + membership_tol = atol+_surface_roundoff(S_own_h)+_surface_roundoff(S_other_h) lo, hi = (w_in, w_out) for _ in range(24): mid = 0.5 * (lo + hi) uo, vo = _edge_own_uv(axis, side, mid) xyz = eval_surface(S_own_h, uo, vo, rational=True) _, _, resid = _invert_point(S_other_h, xyz, seed=tuple(seed)) - if resid <= atol: + if np.isfinite(resid) and resid <= membership_tol: lo = mid else: hi = mid @@ -481,6 +533,9 @@ def _resample_rim(S_own_h, S_other_h, owner, axis, side, atol, w_lo, w_hi, dense): if w_hi - w_lo <= 1e-12: return None + membership_tol = atol+_surface_roundoff(S_own_h)+_surface_roundoff(S_other_h) + if not np.isfinite(membership_tol): + return None ws = np.linspace(w_lo, w_hi, dense) stuv = np.zeros((dense, 4), dtype=np.float64) xyz = np.zeros((dense, 3), dtype=np.float64) @@ -490,12 +545,16 @@ def _resample_rim(S_own_h, S_other_h, owner, axis, side, atol, uo, vo = _edge_own_uv(axis, side, float(w)) p = eval_surface(S_own_h, uo, vo, rational=True) u, v, resid = _invert_point(S_other_h, p, seed=seed) - if resid > atol: + if not (np.isfinite(resid) and resid <= membership_tol): return None # the refined span must certify end-to-end seed = (u, v) resid_max = max(resid_max, resid) stuv[k] = _stuv_sample(owner, (uo, vo), (u, v)) - xyz[k] = p + # Use the same paired-point convention at discovery, corrected + # corners, and adaptive midpoints. Mixing owner points with averages + # can leave an irreducible bound at an old sample even when the + # complete two-surface gap is smaller than the requested tolerance. + xyz[k] = .5*(p+eval_surface(S_other_h, u, v, rational=True)) if float(np.linalg.norm(xyz[-1] - xyz[0])) < atol and dense > 2: return None # degenerate (corner-touch) span return {"stuv": stuv, "xyz": xyz, "resid_max": resid_max, @@ -506,6 +565,187 @@ def _resample_rim(S_own_h, S_other_h, owner, axis, side, atol, # Loop assembly # --------------------------------------------------------------------------- +def _joint_rim_corner(first, second, rim_a, end_a, rim_b, end_b, + atol, ptol4, charge): + """Correct a pair of rim ends on their actual source-domain edges. + + Opposite-chart clamping identifies the other incident edge. Native + source corners stay fixed. Other coordinates stay in the endpoint's + half of its own observed rim span, so Newton cannot jump to the rim's + other corner or to a disjoint contact span. + """ + if rim_a['owner'] == rim_b['owner']: + return None + rims, ends = ((rim_a, rim_b), (end_a, end_b)) if rim_a['owner'] == 1 else ((rim_b, rim_a), (end_b, end_a)) + endpoint = [np.asarray(rim['stuv'][end], dtype=float) for rim, end in zip(rims, ends)] + tolerance = np.minimum(.25, np.maximum(np.asarray(ptol4), 64.*np.finfo(float).eps)) + for owner in (0, 1): + axis = 2*owner+rims[owner]['axis'] + if abs(endpoint[1-owner][axis]-rims[owner]['side']) > tolerance[axis]: + return None + point = np.r_[endpoint[0][:2], endpoint[1][2:]] + fixed = np.zeros(4, dtype=bool) + limits = np.array([[0., 1.]]*4) + for owner, (rim, end) in enumerate(zip(rims, ends)): + axis = 2*owner+rim['axis'] + free = 2*owner+1-rim['axis'] + point[axis] = rim['side'] + fixed[axis] = True + low, high = sorted(rim['stuv'][[0, -1], free]) + middle = .5*(low+high) + limits[free] = (low, middle) if end == 0 else (middle, high) + if point[free] in (0., 1.): + fixed[free] = True + free_axes = np.flatnonzero(~fixed) + surfaces = (first, second) + + def evaluate(parameters): + values = [eval_surface_d1(surface, *parameters[2*i:2*i+2], rational=True) + for i, surface in enumerate(surfaces)] + residual = values[0][0]-values[1][0] + jacobian = np.column_stack((values[0][1], values[0][2], + -values[1][1], -values[1][2]))[:, free_axes] + return values, residual, jacobian + + converged = not len(free_axes) + for _ in range(32): + if not charge(1): + raise _OverlapWorkStopped + values, residual, jacobian = evaluate(point) + if not len(free_axes): + break + try: + singular = np.linalg.svd(jacobian, compute_uv=False) + if singular[-1] <= 64.*np.finfo(float).eps*singular[0]: + return None # parallel/coincident edges do not define a corner + step = np.linalg.lstsq(jacobian, -residual, rcond=None)[0] + except np.linalg.LinAlgError: + return None + projected = np.clip(point[free_axes]+step, 0., 1.)-point[free_axes] + if np.max(abs(projected)) <= 128.*np.finfo(float).eps: + converged = True + break + previous = float(residual@residual) + advanced = False + for exponent in range(12): + candidate = point.copy() + candidate[free_axes] = np.clip(point[free_axes]+step*2.**(-exponent), + limits[free_axes, 0], limits[free_axes, 1]) + other = [eval_surface(surface, *candidate[2*i:2*i+2], rational=True) + for i, surface in enumerate(surfaces)] + error = other[0]-other[1] + if float(error@error) < previous: + point = candidate + advanced = True + break + if not advanced: + break + if not converged: + return None + values, residual, _ = evaluate(point) + membership_tol = atol+_surface_roundoff(first)+_surface_roundoff(second) + if (not np.isfinite(point).all() or not np.isfinite(residual).all() + or not np.isfinite(membership_tol) + or np.linalg.norm(residual) > membership_tol): + return None + xyz = .5*(values[0][0]+values[1][0]) + score = sum(float(np.linalg.norm(point[2*i:2*i+2]-endpoint[i][2*i:2*i+2])) + for i in (0, 1)) + return point, xyz, score + + +def _refine_curved_rim_corners(rims, first, second, atol, ptol4, charge): + """Share jointly corrected corners before endpoint graph construction.""" + proposals = [] + for i, rim in enumerate(rims): + for j in range(i+1, len(rims)): + other = rims[j] + if rim['owner'] == other['owner']: + continue + for end_i in (0, -1): + for end_j in (0, -1): + result = _joint_rim_corner(first, second, rim, end_i, other, + end_j, atol, ptol4, charge) + if result is not None: + point, xyz, score = result + proposals.append((score, i, end_i, j, end_j, point, xyz)) + used = set() + for _, i, end_i, j, end_j, point, xyz in sorted(proposals, key=lambda p: p[0]): + if (i, end_i) in used or (j, end_j) in used: + continue + for index, end in ((i, end_i), (j, end_j)): + rim = rims[index] + rim['stuv'][end] = point + rim['xyz'][end] = xyz + used.add((index, end)) + return rims + + +def _adaptive_curved_rim(rim, first, second, atol, charge): + """Refine a paired rim until both lifted parameter chords meet atol.""" + owner = rim['owner']-1 + own, other = (first, second) if owner == 0 else (second, first) + roundoff = _surface_roundoff(first, rim['xyz'])+_surface_roundoff(second, rim['xyz']) + if not np.isfinite(roundoff) or roundoff >= atol: + return None # this arithmetic cannot resolve the requested accuracy + free = 2*owner+1-rim['axis'] + start, end = np.asarray(rim['stuv'][0]), np.asarray(rim['stuv'][-1]) + # Correcting an endpoint retracts its tolerance fringe. Remove old + # samples outside the new interval before refining the retained path. + low, high = sorted((start[free], end[free])) + samples = [(q.copy(), x.copy()) for q, x in zip(rim['stuv'], rim['xyz']) + if low <= q[free] <= high] + if len(samples) < 2: + return None + result = [samples[0]] + pending = list(reversed(list(zip(samples[:-1], samples[1:])))) + residual_max = 0. + while pending: + a, b = pending.pop() + if not charge(1): + raise _OverlapWorkStopped + errors = [_surface_chord_error(surface, a[0][2*i:2*i+2], b[0][2*i:2*i+2], + np.array([a[1], b[1]])) + for i, surface in enumerate((first, second))] + # The sum also bounds the distance between the two lifted source + # chords. Individual source errors alone could add up to 2*atol. + # Each bound already includes its arithmetic pad. A second pad + # covers evaluation/comparison rounding at the CLOSED atol boundary; + # otherwise a constant gap equal to atol subdivides forever. + if not np.isfinite(errors).all(): + return None + if sum(errors) <= atol+2.*roundoff: + result.append(b) + continue + value = .5*(a[0][free]+b[0][free]) + if value == a[0][free] or value == b[0][free]: + return None + own_uv = _edge_own_uv(rim['axis'], rim['side'], value) + xyz = eval_surface(own, *own_uv, rational=True) + opposite = slice(2*(1-owner), 2*(1-owner)+2) + seed = .5*(a[0][opposite]+b[0][opposite]) + u, v, residual = _invert_point(other, xyz, seed=seed) + if not (np.isfinite(residual) and residual <= atol+roundoff): + return None + point = _stuv_sample(owner+1, own_uv, (u, v)) + projected = eval_surface(other, u, v, rational=True) + middle = (np.asarray(point), .5*(xyz+projected)) + residual_max = max(residual_max, residual) + pending.extend(((middle, b), (a, middle))) + path = np.asarray([q for q, _ in result]) + xyz = np.asarray([x for _, x in result]) + for parameters, point in zip(path, xyz): + source_points = [] + for side, surface in enumerate((first, second)): + source_point = eval_surface( + surface, *parameters[2*side:2*side+2], rational=True) + source_points.append(source_point) + residual_max = max(residual_max, float(np.linalg.norm(source_point-point))) + residual_max = max(residual_max, + float(np.linalg.norm(source_points[0]-source_points[1]))) + return dict(rim, stuv=path, xyz=xyz, resid_max=residual_max) + + def _dist_point_polyline(p, poly): a, b = poly[:-1], poly[1:] ab = b - a @@ -622,7 +862,8 @@ def _loop_paths(rims, loop): for idx, (ri, rev) in enumerate(loop): stuv = rims[ri]["stuv"][::-1] if rev else rims[ri]["stuv"] xyz = rims[ri]["xyz"][::-1] if rev else rims[ri]["xyz"] - if idx > 0: + if (idx > 0 and np.array_equal(stuv_parts[-1][-1], stuv[0]) + and np.array_equal(xyz_parts[-1][-1], xyz[0])): stuv, xyz = stuv[1:], xyz[1:] stuv_parts.append(stuv) xyz_parts.append(xyz) @@ -714,6 +955,7 @@ def assemble_overlap_regions( overlap_boxes=(), charge=None, _planar_data=None, + _completed_rims=None, ): """Assemble certified SSXOverlapRegion entities from rim evidence. @@ -759,8 +1001,25 @@ def _charge(n): own = S1_h if owner == 1 else S2_h other = S2_h if owner == 1 else S1_h rims.extend(_sample_edge_rims(own, other, owner, axis, side, atol)) - rims = _dedup_rims(rims, atol) - rims = _peel_dangling_rims(rims, atol) + try: + rims = _refine_curved_rim_corners(rims, S1_h, S2_h, atol, ptol4, _charge) + rims = _dedup_rims(rims, atol) + rims = _peel_dangling_rims(rims, atol) + refined = [] + for rim in rims: + validated = _adaptive_curved_rim(rim, S1_h, S2_h, atol, _charge) + if validated is None: + return empty + refined.append(validated) + if _completed_rims is not None: + # Private progress sink: only whole-chord-validated + # paired paths survive a later work-budget denial. + _completed_rims.append(validated) + except _OverlapWorkStopped: + return empty + if any(rim is None for rim in refined): + return empty + rims = refined if not rims: return empty diff --git a/mmcore/numeric/intersection/ssx/_ssx_spherical_overlap.py b/mmcore/numeric/intersection/ssx/_ssx_spherical_overlap.py new file mode 100644 index 00000000..01682c1e --- /dev/null +++ b/mmcore/numeric/intersection/ssx/_ssx_spherical_overlap.py @@ -0,0 +1,496 @@ +"""A bounded CAD reduction for two standard rational spherical octants. + +The control template proposes an ideal chart only. A homogeneous Bernstein +bound compares every point of each actual patch to that chart. An analytic +hemisphere census supplies the complete footprint; an existing numerical +assembler supplies the actual paired rim geometry. Unsupported inputs stay +with ordinary SSX. No empty-intersection result is inferred here. +""" +from __future__ import annotations + +from dataclasses import dataclass +from itertools import permutations, product + +import numpy as np + +from mmcore.numeric._bezier_common import eval_surface, geometry_collapsed +from mmcore.numeric.bern import bernstein_product_conv, de_casteljau_split_nd +from mmcore.numeric.intersection.ssx._ssx5_overlap import ( + assemble_overlap_regions, _surface_chord_error, +) +from mmcore.numeric.intersection.ssx._ssx_substrate import SSXBranch + +_PATTERN = np.array([ + [[1., 0., 0.], [1., 0., 1.], [0., 0., 1.]], + [[1., 1., 0.], [1., 1., 1.], [0., 0., 1.]], + [[0., 1., 0.], [0., 1., 1.], [0., 0., 1.]], +]) +_WEIGHTS = np.outer([1., np.sqrt(.5), 1.], [1., np.sqrt(.5), 1.]) +_DESIGN = np.column_stack((np.ones(9), _PATTERN.reshape(-1, 3))) +_PRODUCT = bernstein_product_conv(2) +_EPS = np.finfo(float).eps + + +@dataclass +class _Chart: + center: np.ndarray + radius: float + frame: np.ndarray + transpose: bool + reverse_u: bool + reverse_v: bool + control_error: float + + def canonical_net(self, net): + result = net.swapaxes(0, 1) if self.transpose else net + if self.reverse_u: + result = result[::-1] + if self.reverse_v: + result = result[:, ::-1] + return result + + def original_net(self, net): + result = net[::-1] if self.reverse_u else net + if self.reverse_v: + result = result[:, ::-1] + return result.swapaxes(0, 1) if self.transpose else result + + def canonical_parameters(self, parameters): + result = np.asarray(parameters, dtype=float).copy() + if self.transpose: + result = result[:, ::-1].copy() + if self.reverse_u: + result[:, 0] = 1.-result[:, 0] + if self.reverse_v: + result[:, 1] = 1.-result[:, 1] + return result + + +def _propose_chart(surface, budget): + best = None + for transpose, reverse_u, reverse_v in product((False, True), repeat=3): + if not budget.charge_cells(1, 'spherical_fit'): + return None + candidate = _Chart(None, 0., None, transpose, reverse_u, reverse_v, np.inf) + net = candidate.canonical_net(surface) + controls = (net[..., :3]/net[..., 3, None]).reshape(-1, 3) + origin = controls.mean(axis=0) + try: + fitted = np.linalg.lstsq(_DESIGN, controls-origin, rcond=None)[0] + left, values, right = np.linalg.svd(fitted[1:].T) + except np.linalg.LinAlgError: + continue + radius = float(values.mean()) + if not np.isfinite(radius) or radius <= 0.: + continue + frame, center = left@right, origin+fitted[0] + expected = center+radius*(_PATTERN@frame.T) + error = float(np.linalg.norm(expected.reshape(-1, 3)-controls, axis=1).max()) + if best is None or error < best.control_error: + candidate.center, candidate.radius, candidate.frame = center, radius, frame + candidate.control_error = error + best = candidate + return best + + +def _ideal_net(center, radius, frame): + points = center+radius*(_PATTERN@frame.T) + return np.concatenate((points*_WEIGHTS[..., None], _WEIGHTS[..., None]), axis=2) + + +def _chart_error(actual, chart, center, radius, frame, budget): + """Whole same-parameter rational difference, not a control-point fit.""" + if not budget.charge_cells(112, 'spherical_bound'): + return np.inf, None + first = chart.canonical_net(actual).copy() + first /= first[..., 3].max() + # Translation and radius scaling keep product cancellation independent + # of the model's coordinate placement and of homogeneous gauge. + original_size = max(float(np.max(abs(first[..., :3]))), float(np.max(abs(center)))) + first[..., :3] = (first[..., :3]-center*first[..., 3, None])/radius + second = np.concatenate(((_PATTERN@frame.T)*_WEIGHTS[..., None], + _WEIGHTS[..., None]), axis=2) + numerator = ( + np.einsum('aij,bkl,ikd,jl->abd', _PRODUCT, _PRODUCT, first[..., :3], second[..., 3]) + - np.einsum('aij,bkl,ikd,jl->abd', _PRODUCT, _PRODUCT, second[..., :3], first[..., 3])) + denominator = np.einsum('aij,bkl,ik,jl->ab', _PRODUCT, _PRODUCT, + first[..., 3], second[..., 3]) + if (not np.isfinite(numerator).all() or not np.isfinite(denominator).all() + or np.min(denominator) <= 0.): + return np.inf, None + operands = max(1., float(np.max(abs(first[..., :3]))), + float(np.max(abs(second[..., :3]))), original_size/radius) + roundoff = 256.*_EPS*operands + error = radius*float(np.max((np.linalg.norm(numerator, axis=2)+roundoff) + / denominator)) + ideal = chart.original_net(_ideal_net(center, radius, frame)) + return error, ideal + + +def _footprint(frames, radius, atol, budget): + """All vertices and edges of a convex intersection of six hemispheres.""" + normals = np.concatenate([frame.T for frame in frames]) + vertices = [] + positional_error = 0. + for i in range(6): + for j in range(i+1, 6): + if not budget.charge_cells(1, 'spherical_footprint'): + return None + axis = np.cross(normals[i], normals[j]) + length = float(np.linalg.norm(axis)) + if length <= 256.*_EPS: + continue + error = 512.*_EPS/length + if radius*error > atol/16.: + return None # ill-conditioned footprint: ordinary SSX owns it + positional_error = max(positional_error, error) + for direction in (axis/length, -axis/length): + if np.min(normals@direction) < -error: + continue + if not any(np.linalg.norm(direction-old) <= 4.*error for old in vertices): + vertices.append(direction) + if len(vertices) < 3: + return None # never claim absence or curve-only contact here + vertices = np.asarray(vertices) + inside = vertices.sum(axis=0) + if np.linalg.norm(inside) <= 256.*_EPS: + return None + inside /= np.linalg.norm(inside) + # Positive area must have an interior in all six inward halfspaces. + if np.min(normals@inside) <= positional_error: + return None + east = vertices[0]-inside*np.dot(vertices[0], inside) + if np.linalg.norm(east) <= 256.*_EPS: + return None + east /= np.linalg.norm(east) + north = np.cross(inside, east) + vertices = vertices[np.argsort(np.arctan2(vertices@north, vertices@east))] + edges = [] + for first, second in zip(vertices, np.roll(vertices, -1, axis=0)): + cosine = float(np.clip(first@second, -1., 1.)) + angle = float(np.arctan2(np.linalg.norm(np.cross(first, second)), cosine)) + if angle <= 0. or angle >= np.pi or radius*angle <= atol/4.: + return None + owners = np.flatnonzero((abs(normals@first) <= 8.*positional_error) + & (abs(normals@second) <= 8.*positional_error)) + if not len(owners): + return None + tangent = second-cosine*first + tangent /= np.linalg.norm(tangent) + edges.append((first, tangent, angle, set(map(int, owners)))) + return edges + + +def _matches_edge(branch, reverse, edge, charts, ideals, center, radius, atol, budget): + parameters, xyz = map(np.asarray, branch.curve) + if reverse: + parameters, xyz = parameters[::-1], xyz[::-1] + if len(parameters) < 2 or parameters.shape != (len(xyz), 4): + return set() + first, tangent, angle, owners = edge + matches = set() + for owner, (chart, ideal) in enumerate(zip(charts, ideals)): + uv = parameters[:, 2*owner:2*owner+2] + canonical = chart.canonical_parameters(uv) + for fixed, value, normal_index in ((0, 0., 1), (0, 1., 0), (1, 0., 2)): + if 3*owner+normal_index not in owners: + continue + if np.any(abs(canonical[:, fixed]-value) > 128.*_EPS): + continue + difference = np.diff(canonical[:, 1-fixed]) + if not (np.all(difference >= -128.*_EPS) + or np.all(difference <= 128.*_EPS)): + continue + if not budget.charge_cells(len(parameters), 'spherical_footprint'): + return set() + directions = np.array([eval_surface(ideal, *q, rational=True)-center for q in uv])/radius + angles = np.arctan2(directions@tangent, directions@first) + sense = 1 if angles[-1] >= angles[0] else -1 + oriented = angles if sense == 1 else angles[::-1] + if (np.any(np.diff(oriented) < -512.*_EPS) + or radius*max(abs(oriented[0]), abs(oriented[-1]-angle)) > atol/4.): + continue + valid = True + for q0, q1, x0, x1 in zip(uv[:-1], uv[1:], xyz[:-1], xyz[1:]): + if not budget.charge_cells(1, 'spherical_footprint'): + return set() + # The ideal edge is a monotone circular arc. This bound + # covers its ENTIRE parameter interval in both directions + # against the corresponding XYZ chord, without grid sampling. + if _surface_chord_error(ideal, q0, q1, np.array([x0, x1])) > .75*atol: + valid = False + break + if valid: + matches.add(sense) + return matches + + +def _completed_rims_accounted_for(completed, assembled, sources, charts, ideals, + atol, budget): + """An unreferenced tolerance stub may only contract to a shared pole. + + Ordinary extra rims always decline the reduction. At a shared native + collapsed edge, free longitude parameters denote the same pole; retain + that identity rather than imposing an arbitrary pole parameter value. + """ + branches = assembled['rim_branches'] + extra = [rim for rim in completed if not any( + np.array_equal(rim['stuv'], branch.curve[0]) + and np.array_equal(rim['xyz'], branch.curve[1]) for branch in branches)] + if not extra: + return True + poles = [] + for source, chart, ideal in zip(sources, charts, ideals): + actual_edge = chart.canonical_net(source)[:, -1] + if not geometry_collapsed(actual_edge[:, :3]/actual_edge[:, 3, None]): + return False + ideal_edge = chart.canonical_net(ideal)[:, -1] + poles.append(ideal_edge[0, :3]/ideal_edge[0, 3]) + roundoff = 512.*_EPS*max(1., float(np.max(np.abs(poles)))) + if np.linalg.norm(poles[0]-poles[1]) > roundoff: + return False + endpoints = [np.asarray(branch.curve[1])[end] for branch in branches for end in (0, -1)] + represented = [point for point in endpoints if np.linalg.norm(point-poles[0]) <= atol/4.] + if not represented: + return False + target = min(represented, key=lambda point: np.linalg.norm(point-poles[0])) + for rim in extra: + if np.max(np.linalg.norm(rim['xyz']-target, axis=1))+roundoff > atol: + return False + for a, b in zip(rim['stuv'][:-1], rim['stuv'][1:]): + if not budget.charge_cells(2, 'spherical_footprint'): + return False + for owner, source in enumerate(sources): + if _surface_chord_error(source, a[2*owner:2*owner+2], + b[2*owner:2*owner+2], np.array([target, target])) > atol: + return False + return True + + +def _validate_footprint(assembled, edges, charts, ideals, center, radius, atol, budget, + completed, sources): + regions, branches = assembled['regions'], assembled['rim_branches'] + if not _completed_rims_accounted_for(completed, assembled, sources, charts, + ideals, atol, budget): + return False + if (len(regions) != 1 or len(regions[0].boundary) != 1 + or not regions[0].certification.get('orientation_consistent', False)): + return False + if not budget.charge_cells(2, 'spherical_footprint'): + return False + witness = np.asarray(regions[0].interior_stuv) + if witness.shape != (4,) or np.any(witness < 0.) or np.any(witness > 1.): + return False + for owner, ideal in enumerate(ideals): + direction = (eval_surface(ideal, *witness[2*owner:2*owner+2], rational=True)-center)/radius + if any(float(np.cross(first, tangent)@direction) <= 256.*_EPS + for first, tangent, _, _ in edges): + return False + loop = regions[0].boundary[0] + if len(loop) != len(edges) or len(branches) != len(edges): + return False + if len({index for index, _ in loop}) != len(branches): + return False + matches = {} + for i, (index, reverse) in enumerate(loop): + for j, edge in enumerate(edges): + matches[i, j] = _matches_edge(branches[index], reverse, edge, charts, + ideals, center, radius, atol, budget) + if budget.exhausted: + return False + for sense in (-1, 1): + for start in range(len(edges)): + if all(sense in matches[i, (start+sense*i) % len(edges)] for i in range(len(loop))): + return True + return False + + +class _SphericalResult(dict): + """Private postprocessing facts; these are not public result fields.""" + + _c1_resolved = False + _c3_resolved = False + + +def _resolve_spherical_c1(sources, charts, result, atol, budget): + """Exclude rank defects away from represented or disjoint pole caps. + + The whole-map position bound alone does not imply derivative regularity. + A signed Bernstein normal hull proves it on each truncated chart; the + omitted cap must fit inside a CAD point already on the returned boundary, + or be separated from the other actual surface by more than ``atol``. + """ + from mmcore.numeric.intersection.ssx._ssx5_singular import sigma_normal_net + + if not result.get('overlap_regions'): + return False + corners = np.array([np.asarray(branch.curve[1])[end] + for branch in result['branches'] for end in (0, -1)]) + if not len(corners): + return False + cartesian_sources = [net[..., :3]/net[..., 3, None] for net in sources] + for owner, (source, chart) in enumerate(zip(sources, charts)): + # Charge before constructing the rational normal numerator. The + # input is biquadratic, so this setup always produces an 8x8 net. + if not budget.charge_cells(256, 'spherical_c1'): + return False + net = chart.canonical_net(source).copy() + net /= net[..., 3].max() + coordinate_scale = max(1., float(np.max(abs(net[..., :3]))), + float(np.max(abs(chart.center)))) + spatial_roundoff = 512.*_EPS*coordinate_scale/float(net[..., 3].min()) + pole = net[0, -1, :3]/net[0, -1, 3] + net[..., :3] = (net[..., :3]-chart.center*net[..., 3, None])/chart.radius + normal = sigma_normal_net(net, rational=True) + direction = np.linalg.det(chart.frame)*chart.frame.sum(axis=1) + scalar = normal@direction + normal_roundoff = (4096.*_EPS*max(1., float(np.max(abs(normal)))) + * max(1., coordinate_scale/chart.radius)) + normalized_pole = (pole-chart.center)/chart.radius + # Dyadic restriction ends before a rounded cut could equal one. + # Failure to confine a cap simply leaves ordinary C1 enabled. + cap_radius = None + for depth in range(1, 53): + if not budget.charge_cells(1, 'spherical_c1'): + return False + cut = 1.-2.**-depth + _, cap = de_casteljau_split_nd(net, 1, cut) + cap_points = cap[..., :3]/cap[..., 3, None] + radius = (chart.radius*float(np.max(np.linalg.norm( + cap_points-normalized_pole, axis=2)))+spatial_roundoff) + if radius <= atol/4.: + cap_radius = radius + regular, _ = de_casteljau_split_nd(scalar[..., None], 1, cut) + if (not np.isfinite(regular).all() + or float(np.min(regular)) <= normal_roundoff): + return False + break + if cap_radius is None: + return False + if float(np.min(np.linalg.norm(corners-pole, axis=1)))+cap_radius <= atol: + continue + other_controls = cartesian_sources[1-owner] + separation_roundoff = 512.*_EPS*max( + coordinate_scale, float(np.max(abs(other_controls)))) + separated = False + for direction in charts[1-owner].frame.T: + # The fitted frame is orthogonal; normalizing keeps the + # support distance in model units despite SVD roundoff. + direction = direction/np.linalg.norm(direction) + distance = float(np.min((other_controls-pole)@direction)) + if distance-separation_roundoff > cap_radius+atol: + separated = True + break + if not separated: + return False + return True + + +def _publish(branches, regions, budget): + result = _SphericalResult(branches=[], points=[], singularities=[], overlap_regions=[]) + budget.extend_output(result['branches'], branches, 'spherical_overlap') + if len(result['branches']) == len(branches): + budget.extend_output(result['overlap_regions'], regions, 'spherical_overlap') + return result + + +def _partial(completed, budget): + branches = [SSXBranch(curve=(rim['stuv'], rim['xyz']), kind='overlap', overlap=True) + for rim in completed] + return _publish(branches, [], budget) + + +def try_spherical_overlap(first, second, atol, budget): + """Return an actual paired overlap region, or decline this reduction.""" + sources = tuple(np.asarray(net, dtype=float) for net in (first, second)) + if any(net.shape != (3, 3, 4) or not np.isfinite(net).all() + or np.any(net[..., 3] <= 0.) for net in sources): + return None + charts = tuple(_propose_chart(net, budget) for net in sources) + if any(chart is None for chart in charts) or budget.exhausted: + return None + # Rejection-only proposal screens keep ordinary biquadratic planes and + # graphs out of the more expensive rational product checks. Passing a + # screen never accepts a carrier: the whole-map bound is still required. + if (any(chart.control_error > atol for chart in charts) + or np.linalg.norm(charts[0].center-charts[1].center) > atol + or abs(charts[0].radius-charts[1].radius) > atol): + return None + center = .5*(charts[0].center+charts[1].center) + radius = .5*(charts[0].radius+charts[1].radius) + if not np.isfinite(radius) or radius <= atol: + return None + # Coalescing an almost identical octant is only another proposal. Its + # whole-patch bound must pass too; actual controls are never snapped. + frame_pairs = [(charts[0].frame, charts[0].frame[:, permutation]) + for permutation in permutations(range(3))] + frame_pairs.append((charts[0].frame, charts[1].frame)) + first_validated = _chart_error(sources[0], charts[0], center, radius, charts[0].frame, budget) + if budget.exhausted or first_validated[0] > atol/8.: + return None + candidates = [] + second_cartesian = charts[1].canonical_net(sources[1]) + second_cartesian = second_cartesian[..., :3]/second_cartesian[..., 3, None] + for frames in frame_pairs: + proposal = center+radius*(_PATTERN@frames[1].T) + if np.max(np.linalg.norm(proposal-second_cartesian, axis=2)) > atol: + continue + validated = [first_validated, + _chart_error(sources[1], charts[1], center, radius, frames[1], budget)] + if budget.exhausted: + return None + if validated[1][0] > atol/8.: + continue + edges = _footprint(frames, radius, atol, budget) + if budget.exhausted: + return None + if edges is not None: + # Whole-footprint stability, not just matching sampled corners. + # An actual common point is at most e from each ideal octant, + # hence every inward normal has n.d >= -e/(R-e). Moving its + # direction toward an interior direction a by lambda restores + # all six inequalities. Normalization adds at most lambda again. + # Reject poorly conditioned narrow footprints rather than let a + # small chart error conceal a long, separate boundary fringe. + interior = np.sum([edge[0] for edge in edges], axis=0) + interior /= np.linalg.norm(interior) + normals = np.concatenate([frame.T for frame in frames]) + margin = float(np.min(normals@interior))-1024.*_EPS + error = max(value for value, _ in validated) + if margin <= 0. or error >= radius: + continue + displacement = error/((radius-error)*margin) + if displacement >= 1. or error+2.*radius*displacement > atol/4.: + continue + candidates.append((edges, [ideal for _, ideal in validated])) + if not candidates: + return None + from mmcore.nurbs._nurbs_param_tol import bez_surface_param_tolerance + ptol = np.array([*bez_surface_param_tolerance(first, .5*atol, rational=True), + *bez_surface_param_tolerance(second, .5*atol, rational=True)]) + completed = [] + assembled = assemble_overlap_regions( + *sources, atol=.5*atol, ptol4=ptol, + charge=lambda amount: budget.charge_cells(amount, 'spherical_region'), + _completed_rims=completed) + if budget.exhausted: + return _partial(completed, budget) + if (not assembled['regions'] + or assembled['unmatched_branches'] + or assembled['unmatched_intersection_branches']): + return None + for edges, ideals in candidates: + if _validate_footprint(assembled, edges, charts, ideals, center, radius, atol, budget, + completed, sources): + result = _publish(assembled['rim_branches'], assembled['regions'], budget) + if result['overlap_regions']: + # A fully represented convex footprint has no 1D rim + # self-crossings. Shared corners and pole UV aliases are + # already part of its boundary, not C3 events. + result._c3_resolved = True + result._c1_resolved = _resolve_spherical_c1( + sources, charts, result, atol, budget) + return result + if budget.exhausted: + return _partial(completed, budget) + return None diff --git a/tests/test_csx_cad_span_contacts.py b/tests/test_csx_cad_span_contacts.py new file mode 100644 index 00000000..a4ce319d --- /dev/null +++ b/tests/test_csx_cad_span_contacts.py @@ -0,0 +1,204 @@ +"""CAD coincidence retains the numerical contact seeds needed by SSX.""" +from math import comb + +import numpy as np +import pytest + +import mmcore.numeric.intersection.csx._bez_csx4 as csx +from mmcore.numeric._bezier_common import eval_curve, eval_surface +from mmcore.numeric.intersection.ssx._bez_ssx5 import _cut_face_contacts + + +def _three_crossings(amplitude): + power = amplitude * np.polynomial.polynomial.polyfromroots([.2, .5, .8]) + coefficients = np.array([sum(power[k]*comb(i, k)/comb(3, k) + for k in range(i+1)) for i in range(4)]) + curve = np.column_stack((np.linspace(.1, .9, 4), np.full(4, .5), coefficients)) + plane = np.array([[[s, t, 0.] for t in (0., 1.)] for s in (0., 1.)]) + return curve, plane + + +def _validate_contacts(curve, surface, contacts): + for point in contacts: + assert all(0. <= point[key] <= 1. for key in ('t', 'u', 'v')) + curve_point = eval_curve(curve, point['t'], rational=False) + surface_point = eval_surface(surface, point['u'], point['v'], rational=False) + assert np.linalg.norm(curve_point-surface_point) <= 1e-3 + assert np.linalg.norm(point['point']-curve_point) <= 1e-3 + + +def test_full_cad_span_keeps_three_distinct_interior_contact_proposals(): + curve, surface = _three_crossings(.005) + result = csx.bez_csx(curve, surface, atol=1e-3, rational=False) + assert not result['budget_exhausted'] + assert len(result['overlaps']) == 1 + assert result['overlaps'][0]['t_range'] == (0., 1.) + assert result['overlaps'][0]['certification'] == 'tolerance' + contacts = list(_cut_face_contacts(result)) + for expected in (.2, .5, .8): + matches = [point for point in contacts if abs(point['t']-expected) <= 1e-3] + assert matches, (expected, contacts) + assert all(point['certification'] == 'tolerance' for point in matches) + _validate_contacts(curve, surface, contacts) + + +@pytest.mark.parametrize('amplitude, exact_mode', [(1., False), (.005, True)]) +def test_isolated_crossings_remain_when_no_cad_span_is_admitted(amplitude, exact_mode): + curve, surface = _three_crossings(amplitude) + result = csx.bez_csx(curve, surface, atol=1e-3, rational=False, + tolerance_tier=not exact_mode) + assert result['overlaps'] == [] + assert not result['budget_exhausted'] + assert sorted(point['t'] for point in result['isolated']) == pytest.approx([.2, .5, .8], abs=1e-3) + _validate_contacts(curve, surface, result['isolated']) + + +def test_contact_refinement_denial_keeps_the_valid_cad_span(monkeypatch): + curve, surface = _three_crossings(.005) + def forbidden(*args, **kwargs): + raise AssertionError('contact polishing ran after work denial') + monkeypatch.setattr(csx, '_polish_csx_root', forbidden) + spans = csx._tolerance_csx_overlap_certificate( + curve, surface, 1e-3, False, 1e-4, None, + on_contact_work=lambda amount: False) + assert len(spans) == 1 and spans[0]['t_range'] == (0., 1.) + assert spans[0].get('boundary_contacts', []) == [] + + +def test_contact_result_cap_is_reported_without_losing_the_span(): + curve, surface = _three_crossings(.005) + result = csx.bez_csx(curve, surface, atol=1e-3, rational=False, + max_results=3) + assert len(result['overlaps']) == 1 + assert result['overlaps'][0]['t_range'] == (0., 1.) + assert result['budget_exhausted'] + assert result['truncation_cause'] == 'results' + contacts = list(_cut_face_contacts(result)) + assert len(contacts) <= 3 + _validate_contacts(curve, surface, contacts) + + +@pytest.mark.parametrize('max_results', [3, 10]) +def test_overlap_contacts_share_the_cap_with_separate_isolated_roots(max_results): + # The curve lies within CAD tolerance over its first in-domain span + # [0,.3], crosses the plane at .15 there, then reenters the patch and + # has a separate isolated root at .8. All three span contacts remain + # charged when they move from the isolated list onto the overlap. + degree = 8 + def coefficients(power): + power = np.pad(power, (0, degree+1-len(power))) + return [sum(power[k]*comb(i, k)/comb(degree, k) + for k in range(i+1)) for i in range(degree+1)] + z = np.polynomial.polynomial.polymul( + [0.]*6+[1.], np.polynomial.polynomial.polyfromroots([.15, .8])) + curve = np.column_stack((coefficients([.04, 1.8, -2.]), + coefficients([.5]), coefficients(z))) + surface = np.array([[[.4*u, v, 0.] for v in (0., 1.)] + for u in (0., 1.)]) + result = csx.bez_csx(curve, surface, atol=1e-3, rational=False, + max_results=max_results) + assert len(result['overlaps']) == 1 + contacts = list(_cut_face_contacts(result)) + assert len(contacts) <= max_results + _validate_contacts(curve, surface, contacts) + assert any(abs(point['t']-.15) <= 1e-3 for point in contacts) + if max_results == 3: + assert result['budget_exhausted'] + assert result['truncation_cause'] == 'results' + else: + assert not result['budget_exhausted'] + assert any(abs(point['t']-.8) <= 1e-3 + for point in result['isolated']) + + +def test_denied_cad_span_classification_does_not_report_complete(): + curve = np.array([[0., .5, .0005], [1., .5, .0005]]) + surface = np.array([[[u, v, 0.] for v in (0., 1.)] + for u in (0., 1.)]) + result = csx.bez_csx(curve, surface, atol=1e-3, rational=False, + max_cells=200) + # Exact zero exclusion cannot finish a CAD span-classification task + # which was skipped for lack of work allowance. + assert result['budget_exhausted'] + assert result['truncation_cause'] == 'cells' + assert result['cells_processed'] <= 200 + _validate_contacts(curve, surface, list(_cut_face_contacts(result))) + + +def _narrow_nonaffine_strip(start, reverse=False, bulge=False): + end = start + .05 + curve = np.array([[0., 0., 0.], [100., 0., 0.]]) + if reverse: + curve = curve[::-1].copy() + surface = np.array([[[start+(end-start)*u, v-.5+y, z] + for v in (0., 1.)] + for u, y, z in zip((0., .5, 1.), (0., 0., .2), + (0., .025 if bulge else 0., 0.))]) + return curve, surface, end + + +@pytest.mark.parametrize('start,reverse', [(71.23, False), (13.71, False), (71.23, True)]) +def test_boundary_contacts_resolve_a_narrow_nonaffine_clipped_span(start, reverse): + curve, surface, end = _narrow_nonaffine_strip(start, reverse) + result = csx.bez_csx(curve, surface, atol=1e-3, rational=False) + assert not result['budget_exhausted'] + assert result['overlaps'], result + # The true interval is .0005 in curve parameters but 50*atol long + # in model space. Neither fixed global sampling grid visits it. + for u in np.linspace(0., 1., 101): + x = start + (end-start)*u + t = (100.-x)/100. if reverse else x/100. + assert any(low-1e-3/100. <= t <= high+1e-3/100. + for low, high in (span['t_range'] for span in result['overlaps'])) + for span in result['overlaps']: + for t in np.linspace(*span['t_range'], 101): + point = eval_curve(curve, t, rational=False) + u = float(np.clip((point[0]-start)/(end-start), 0., 1.)) + target = eval_surface(surface, u, .5-.2*u*u, rational=False) + assert np.linalg.norm(point-target) <= 1e-3 + + +def test_two_boundary_contacts_do_not_promote_an_out_of_tolerance_bulge(): + curve, surface, end = _narrow_nonaffine_strip(71.23, bulge=True) + result = csx.bez_csx(curve, surface, atol=1e-3, rational=False) + assert result['overlaps'] == [] + assert not result['budget_exhausted'] + assert len(result['isolated']) == 2 + _validate_contacts(curve, surface, result['isolated']) + + +@pytest.mark.parametrize('owner,axis,side', [ + (owner, axis, side) for owner in (0, 1) for axis in (0, 1) for side in (0, 1)]) +def test_rounded_spherical_boundary_spans_cover_the_independent_common_rim(owner, axis, side): + from examples.ssx.case_17 import s1, s2 + from test_ssx_spherical_overlap import _fit_octant_frame, _reference, _homogeneous + from examples.ssx.ssx5_analytic_audit import distances_to_polylines + + pair = (s1, s2) + fits = [_fit_octant_frame(surface) for surface in pair] + center = np.mean([fit[0] for fit in fits], axis=0) + radius = np.mean([fit[1] for fit in fits]) + reference = _reference(center, radius, [fit[2] for fit in fits]) + sources = [_homogeneous(surface) for surface in pair] + curve = sources[owner][2*side] if axis == 0 else sources[owner][:, 2*side] + result = csx.bez_csx(curve, sources[1-owner], atol=1e-3, rational=True) + assert not result['budget_exhausted'] + normal_index = (1-side) if axis == 0 else (2 if side == 0 else None) + wanted = [] + if normal_index is not None: + normal = fits[owner][2][:, normal_index] + wanted = [arc for arc in reference['arcs'] + if np.max(abs((arc-center) @ normal)) <= .01*1e-3] + if not wanted: + assert result['overlaps'] == [] and result['isolated'] == [] + return + assert len(result['overlaps']) == 1 + paths = [np.array([eval_curve(curve, t, rational=True) + for t in np.linspace(*span['t_range'], 1025)]) + for span in result['overlaps']] + for arc in wanted: + assert distances_to_polylines(arc, paths).max() <= 1e-3 + for point in _cut_face_contacts(result): + a = eval_curve(curve, point['t'], rational=True) + b = eval_surface(sources[1-owner], point['u'], point['v'], rational=True) + assert np.linalg.norm(a-b) <= 1e-3 diff --git a/tests/test_csx_overlap_normal_gap.py b/tests/test_csx_overlap_normal_gap.py index fc4d2039..11372f27 100644 --- a/tests/test_csx_overlap_normal_gap.py +++ b/tests/test_csx_overlap_normal_gap.py @@ -8,7 +8,8 @@ @pytest.mark.parametrize("normal_gap, overlap", [ (4.0 * np.finfo(float).eps, True), - (1e-7, False), + (1e-7, True), + (2e-3, False), ]) def test_tangential_projection_error_does_not_amplify_normal_signs( monkeypatch, normal_gap, overlap): @@ -31,6 +32,27 @@ def lagging_projection(point, surface, u, v, atol, rational): assert len(result) == 1 assert result[0]["t_range"] == (0., 1.) else: - # A resolved normal sign change remains a crossing even though - # the whole curve lies inside the modeling tolerance band. + # The in-tolerance middle interval has neither a curve-domain end + # nor a surface-domain end. It remains an isolated crossing. assert result is None + + +@pytest.mark.parametrize('curve_scale,surface_scale', [(1., 1.), (1e8, 1e8), (1e-8, 1e8)]) +def test_normal_crossing_contact_is_independent_of_homogeneous_gauge(curve_scale, surface_scale): + plane = np.array([[[0., 0., 0., 1.], [0., 1., 0., 1.]], + [[1., 0., 0., 1.], [1., 1., 0., 1.]]]) * surface_scale + curve = np.array([[.1, .5, -1e-7, 1.], [.9, .5, 1e-7, 1.]]) * curve_scale + result = csx.bez_csx(curve, plane, atol=1e-3, rational=True) + assert not result['budget_exhausted'] + assert len(result['overlaps']) == 1 + contacts = result['overlaps'][0]['boundary_contacts'] + crossing = [point for point in contacts if abs(point['t']-.5) <= 1e-3] + assert crossing + for point in crossing: + a = csx.eval_curve(curve, point['t'], rational=True) + b = csx.eval_surface(plane, point['u'], point['v'], rational=True) + assert np.linalg.norm(a-b) <= 1e-3 + exact = csx.bez_csx(curve, plane, atol=1e-3, rational=True, tolerance_tier=False) + assert exact['overlaps'] == [] + assert len(exact['isolated']) == 1 + assert abs(exact['isolated'][0]['t']-.5) <= 1e-3 diff --git a/tests/test_ssx_curved_overlap_corners.py b/tests/test_ssx_curved_overlap_corners.py new file mode 100644 index 00000000..c5cc82f5 --- /dev/null +++ b/tests/test_ssx_curved_overlap_corners.py @@ -0,0 +1,135 @@ +"""Shared curved-overlap corners and whole-chord CAD error bounds.""" +import numpy as np +import pytest + +from mmcore.numeric._bezier_common import eval_surface +from mmcore.numeric.intersection.ssx import _ssx5_overlap as overlap + + +ATOL = 1e-3 + + +def _plane(x0=0., x1=1., y0=0., y1=1., z=0.): + return np.array([[[x, y, z, 1.] for y in (y0, y1)] for x in (x0, x1)]) + + +def _rim(owner, axis, side, parameters, surface): + parameters = np.asarray(parameters, dtype=float) + xyz = np.array([eval_surface(surface, *q[2*(owner-1):2*owner], rational=True) + for q in parameters]) + return dict(owner=owner, axis=axis, side=side, stuv=parameters, + xyz=xyz, resid_max=0.) + + +def test_two_tolerance_fringes_become_one_paired_domain_corner(): + first, second = _plane(), _plane(.5, 1., -.5, .5) + a = _rim(1, 1, 0., [[.4994, 0., 0., .5], [1., 0., 1., .5]], first) + b = _rim(2, 0, 0., [[.5, 0., 0., .4994], [.5, .5, 0., 1.]], second) + corner = overlap._joint_rim_corner(first, second, a, 0, b, 0, + ATOL, np.full(4, ATOL), lambda n: True) + assert corner is not None + parameters, xyz, _ = corner + assert np.linalg.norm(xyz-[.5, 0., 0.]) <= ATOL + assert parameters[1] == parameters[2] == 0. # fixed domain faces + for surface, uv in ((first, parameters[:2]), (second, parameters[2:])): + assert np.linalg.norm(eval_surface(surface, *uv, rational=True)-xyz) <= ATOL + overlap._refine_curved_rim_corners([a, b], first, second, ATOL, + np.full(4, ATOL), lambda n: True) + # Shared object values establish the graph node; its accuracy is checked + # geometrically above, rather than by requiring particular solver digits. + np.testing.assert_array_equal(a['stuv'][0], b['stuv'][0]) + np.testing.assert_array_equal(a['xyz'][0], b['xyz'][0]) + + +def test_parallel_nearby_edges_do_not_invent_a_unique_corner(): + first, second = _plane(), _plane(z=.5*ATOL) + parameters = [[.2, 0., .2, 0.], [.8, 0., .8, 0.]] + a, b = _rim(1, 1, 0., parameters, first), _rim(2, 1, 0., parameters, second) + assert overlap._joint_rim_corner(first, second, a, 0, b, 0, + ATOL, np.full(4, ATOL), lambda n: True) is None + + +def test_nearby_opposite_partner_faces_keep_their_parameter_identity(): + first, second = _plane(), _plane(.5, .5+.5*ATOL, -.5, .5) + a = _rim(1, 1, 0., [[.5, 0., 0., .5], [.7, 0., 1., .5]], first) + b = _rim(2, 0, 1., [[.5+.5*ATOL, 0., 1., .5], [.5+.5*ATOL, .5, 1., 1.]], second) + # Even a poorly resolved partner axis cannot make its two domain faces + # the same face. XYZ-only endpoint merging would accept this pair. + assert overlap._joint_rim_corner(first, second, a, 0, b, 0, + ATOL, np.ones(4), lambda n: True) is None + + +def test_corner_correction_cannot_jump_to_a_remote_contact_span(): + first, second = _plane(), _plane(.8, 1., -.5, .5) + a = _rim(1, 1, 0., [[.1, 0., 0., .5], [.3, 0., 0., .5]], first) + b = _rim(2, 0, 0., [[.8, 0., 0., .49], [.8, .5, 0., 1.]], second) + assert overlap._joint_rim_corner(first, second, a, -1, b, 0, + ATOL, np.full(4, ATOL), lambda n: True) is None + + +def test_homogeneous_chord_bound_catches_an_inflection_bulge(): + surface = np.array([[[s, t, z, 1.] for t in (0., 1.)] + for s, z in zip(np.linspace(0., 1., 4), (-.125, .125, -.125, .125))]) + start, end = np.array([0., 0.]), np.array([1., 0.]) + xyz = np.array([eval_surface(surface, *q, rational=True) for q in (start, end)]) + middle = eval_surface(surface, .5, 0., rational=True) + assert np.linalg.norm(middle-xyz.mean(axis=0)) <= ATOL + bound = overlap._surface_chord_error(surface, start, end, xyz) + sampled = max(np.linalg.norm(eval_surface(surface, t, 0., rational=True) + - ((1.-t)*xyz[0]+t*xyz[1])) + for t in np.linspace(0., 1., 257)) + assert sampled > ATOL and sampled <= bound + + +def test_corner_work_denial_keeps_the_original_rim_endpoints(monkeypatch): + first, second = _plane(), _plane(.5, 1., -.5, .5) + a = _rim(1, 1, 0., [[.4994, 0., 0., .5], [1., 0., 1., .5]], first) + b = _rim(2, 0, 0., [[.5, 0., 0., .4994], [.5, .5, 0., 1.]], second) + before = [(rim['stuv'].copy(), rim['xyz'].copy()) for rim in (a, b)] + def forbidden(*args, **kwargs): + raise AssertionError('corner derivative evaluation after work denial') + monkeypatch.setattr(overlap, 'eval_surface_d1', forbidden) + with pytest.raises(overlap._OverlapWorkStopped): + overlap._refine_curved_rim_corners([a, b], first, second, ATOL, + np.full(4, ATOL), lambda n: False) + for rim, (parameters, xyz) in zip((a, b), before): + np.testing.assert_array_equal(rim['stuv'], parameters) + np.testing.assert_array_equal(rim['xyz'], xyz) + + +@pytest.mark.parametrize('gap_fraction', [.9, 1.]) +def test_within_tolerance_offset_keeps_the_complete_paired_region(gap_fraction): + # Degree elevation selects the general curved-rim path. Both entire + # planes are within atol; different sampling stages must not impose a + # stricter gap by mixing owner points with paired midpoint samples. + first = np.array([[[i/2, j/2, gap_fraction*ATOL, 1.] for j in range(3)] + for i in range(3)]) + second = first.copy() + second[..., 2] = 0. + result = overlap.assemble_overlap_regions( + first, second, atol=ATOL, ptol4=np.full(4, ATOL)) + assert len(result['regions']) == 1 + assert len(result['rim_branches']) == 4 + for branch in result['rim_branches']: + parameters, points = branch.curve + for fraction in (0., .25, .5, .75, 1.): + for q, point in zip( + (1.-fraction)*parameters[:-1]+fraction*parameters[1:], + (1.-fraction)*points[:-1]+fraction*points[1:]): + first_point = eval_surface(first, *q[:2], rational=True) + second_point = eval_surface(second, *q[2:], rational=True) + assert np.linalg.norm(first_point-point) <= ATOL + assert np.linalg.norm(second_point-point) <= ATOL + # The closed-boundary control has analytic gap exactly + # atol; evaluating that constant may differ by one ULP. + assert np.linalg.norm(first_point-second_point) <= ATOL+64*np.finfo(float).eps + + +def test_offset_above_tolerance_does_not_make_an_overlap_region(): + first = np.array([[[i/2, j/2, 1.01*ATOL, 1.] for j in range(3)] + for i in range(3)]) + second = first.copy() + second[..., 2] = 0. + result = overlap.assemble_overlap_regions( + first, second, atol=ATOL, ptol4=np.full(4, ATOL)) + assert result['regions'] == [] diff --git a/tests/test_ssx_spherical_overlap.py b/tests/test_ssx_spherical_overlap.py new file mode 100644 index 00000000..3f35b441 --- /dev/null +++ b/tests/test_ssx_spherical_overlap.py @@ -0,0 +1,287 @@ +"""Spherical CAD overlap checked against independent hemisphere geometry.""" +from itertools import combinations + +import numpy as np +import pytest + +from examples.ssx.ssx5_analytic_audit import distances_to_polylines +from mmcore.nurbs._nurbs_eval import NURBSSurfaceTuple, evaluate_nurbs_surface +from mmcore.numeric.intersection.ssx._ssx5_overlap import assemble_overlap_regions + + +ATOL = 1e-3 +_OCTANT_CONTROLS = np.array([ + [[1., 0., 0.], [1., 0., 1.], [0., 0., 1.]], + [[1., 1., 0.], [1., 1., 1.], [0., 0., 1.]], + [[0., 1., 0.], [0., 1., 1.], [0., 0., 1.]], +]) + + +def _octant(frame, radius=1., center=(0., 0., 0.)): + """Independent rational tensor product of two quarter-circle arcs.""" + control = np.asarray(center) + radius * (_OCTANT_CONTROLS @ np.asarray(frame).T) + weights = np.outer([1., np.sqrt(.5), 1.], [1., np.sqrt(.5), 1.]) + knots = np.array([0., 0., 0., 1., 1., 1.]) + return NURBSSurfaceTuple(3, 3, knots.copy(), knots.copy(), control, weights) + + +def _rotation(axis, angle): + direction = np.asarray(axis, dtype=float) + direction /= np.linalg.norm(direction) + x, y, z = direction + cross = np.array([[0., -z, y], [z, 0., -x], [-y, x, 0.]]) + return (np.cos(angle) * np.eye(3) + + (1.-np.cos(angle)) * np.outer(direction, direction) + + np.sin(angle) * cross) + + +def _reverse_u(surface): + return NURBSSurfaceTuple( + surface.order_u, surface.order_v, surface.knot_u.copy(), surface.knot_v.copy(), + surface.control_points[::-1].copy(), surface.weights[::-1].copy()) + + +def _homogeneous(surface): + return np.concatenate((surface.control_points * surface.weights[..., None], + surface.weights[..., None]), axis=-1) + + +def _point(surface, uv): + return np.asarray(evaluate_nurbs_surface(surface, *uv, d_order=0)['S']) + + +def _fit_octant_frame(surface): + # Used only to recover the known construction from rounded user data. + # No SSI algorithm, sampled intersection, or solver result enters this fit. + design = np.column_stack((np.ones(9), _OCTANT_CONTROLS.reshape(-1, 3))) + solution = np.linalg.lstsq(design, surface.control_points.reshape(-1, 3), + rcond=None)[0] + center, scaled_frame = solution[0], solution[1:].T + radius = np.linalg.norm(scaled_frame, axis=0).mean() + left, _, right = np.linalg.svd(scaled_frame / radius) + frame = left @ right + fitted = center + radius * (_OCTANT_CONTROLS @ frame.T) + assert np.linalg.norm(fitted-surface.control_points, axis=-1).max() < .01*ATOL + # Confirm that the supplied rational weights describe this sphere too; + # fitting only Euclidean controls would not establish the fixture oracle. + for u in np.linspace(0., 1., 9): + for v in np.linspace(0., 1., 9): + point = _point(surface, (u, v))-center + assert abs(np.linalg.norm(point)-radius) < .01*ATOL + assert np.min(point @ frame) >= -.01*ATOL + return center, radius, frame + + +def _polygon(frames): + # Each octant is the intersection of three inward hemispheres. Polygon + # vertices are admissible intersections of pairs of their great circles. + normals = np.concatenate([frame.T for frame in frames]) + vertices = [] + for first, second in combinations(normals, 2): + axis = np.cross(first, second) + norm = np.linalg.norm(axis) + if norm < 1e-12: + continue + for point in (axis/norm, -axis/norm): + if (np.min(normals @ point) >= -1e-10 + and not any(np.linalg.norm(point-old) < 1e-9 for old in vertices)): + vertices.append(point) + vertices = np.asarray(vertices) + assert len(vertices) >= 3, 'This reference requires positive overlap area' + interior = vertices.sum(axis=0) + interior /= np.linalg.norm(interior) + east = vertices[0] - interior * np.dot(vertices[0], interior) + east /= np.linalg.norm(east) + north = np.cross(interior, east) + vertices = vertices[np.argsort(np.arctan2(vertices @ north, vertices @ east))] + return vertices, normals, interior + + +def _reference(center, radius, frames): + vertices, normals, interior = _polygon(frames) + arcs, angles, owners = [], [], [] + max_angle = 2*np.arccos(1.-ATOL/(8*radius)) + for first, second in zip(vertices, np.roll(vertices, -1, axis=0)): + angle = np.arctan2(np.linalg.norm(np.cross(first, second)), np.dot(first, second)) + # The independent reference's chord sagitta is at most atol/8. + count = max(257, int(np.ceil(angle/max_angle))+1) + fraction = np.linspace(0., 1., count) + directions = (np.sin((1.-fraction)*angle)[:, None]*first + + np.sin(fraction*angle)[:, None]*second) / np.sin(angle) + arcs.append(np.asarray(center) + radius*directions) + angles.append(angle) + active = np.nonzero((abs(normals @ first) < 1e-9) + & (abs(normals @ second) < 1e-9))[0] + owners.append({int(index)//3 for index in active}) + return dict(center=np.asarray(center), radius=radius, vertices=vertices, + normals=normals, interior=interior, arcs=arcs, + owners=owners, length=radius*sum(angles)) + + +def _distance_to_rim(points, reference): + """Distance to finite great-circle arcs, not their sampled polylines.""" + relative = np.asarray(points)-reference['center'] + radius, vertices = reference['radius'], reference['vertices'] + distances = [] + for first, second in zip(vertices, np.roll(vertices, -1, axis=0)): + normal = np.cross(first, second) + normal /= np.linalg.norm(normal) + projected = relative - (relative @ normal)[:, None]*normal + norms = np.linalg.norm(projected, axis=1) + directions = projected / np.maximum(norms[:, None], np.finfo(float).tiny) + within = ((np.cross(first, directions) @ normal >= -1e-12) + & (np.cross(directions, second) @ normal >= -1e-12) + & (norms > 0.)) + endpoint_distance = np.minimum(np.linalg.norm(relative-radius*first, axis=1), + np.linalg.norm(relative-radius*second, axis=1)) + distances.append(np.where(within, + np.linalg.norm(relative-radius*directions, axis=1), endpoint_distance)) + return np.min(distances, axis=0) + + +def _spherical_area(directions, interior): + # Triangulation about an independent interior point, using solid angles. + area = 0. + for first, second in zip(directions, np.roll(directions, -1, axis=0)): + numerator = np.dot(interior, np.cross(first, second)) + denominator = 1.+np.dot(interior, first)+np.dot(first, second)+np.dot(second, interior) + area += 2*np.arctan2(numerator, denominator) + return abs(area) + + +def _assert_region(result, surfaces, reference, normal_agreement=1): + assert len(result['overlap_regions']) == 1 + region = result['overlap_regions'][0] + assert len(region.boundary) == 1 + assert region.normal_agreement == normal_agreement + assert {index for index, _ in region.boundary[0]} == set(range(len(result['branches']))) + paths, owner_edges = [], set() + for index, reverse in region.boundary[0]: + branch = result['branches'][index] + assert branch.kind == 'overlap' + parameters, xyz = map(np.asarray, branch.curve) + if reverse: + parameters, xyz = parameters[::-1], xyz[::-1] + assert len(xyz) >= 2 and parameters.shape == (len(xyz), 4) + assert np.isfinite(parameters).all() and np.isfinite(xyz).all() + assert np.all((-1e-12 <= parameters) & (parameters <= 1.+1e-12)) + # Corners and both source parameter curves must describe the same + # model-space geometry, including between published vertices. + for fraction in (0., .25, .5, .75, 1.): + for q, point in zip((1.-fraction)*parameters[:-1]+fraction*parameters[1:], + (1.-fraction)*xyz[:-1]+fraction*xyz[1:]): + a, b = _point(surfaces[0], q[:2]), _point(surfaces[1], q[2:]) + assert np.linalg.norm(a-point) <= ATOL + assert np.linalg.norm(b-point) <= ATOL + assert np.linalg.norm(a-b) <= ATOL + assert _distance_to_rim(xyz, reference).max() <= ATOL + assert _distance_to_rim(.5*(xyz[:-1]+xyz[1:]), reference).max() <= ATOL + for owner, surface in enumerate(surfaces): + uv = parameters[:, 2*owner:2*owner+2] + for axis in (0, 1): + for side in (0., 1.): + on_edge = uv.copy() + on_edge[:, axis] = side + if all(np.linalg.norm(_point(surface, edge)-point) <= ATOL + for edge, point in zip(on_edge, xyz)): + owner_edges.add(owner) + paths.append(xyz) + assert owner_edges == {0, 1}, 'The overlap rim must include both domain boundaries' + for first, second in zip(paths, paths[1:]+paths[:1]): + assert np.linalg.norm(first[-1]-second[0]) <= ATOL + for arc in reference['arcs']: + assert distances_to_polylines(arc, paths).max() <= ATOL + length = sum(np.linalg.norm(np.diff(path, axis=0), axis=1).sum() for path in paths) + assert abs(length-reference['length']) <= 8*np.pi*ATOL + loop = np.concatenate(paths) + directions = loop-reference['center'] + directions /= np.linalg.norm(directions, axis=1)[:, None] + expected_area = _spherical_area(reference['vertices'], reference['interior']) + actual_area = _spherical_area(directions, reference['interior']) + assert abs(actual_area-expected_area)*reference['radius']**2 <= 2*reference['length']*ATOL + witness = np.asarray(region.interior_stuv) + a, b = _point(surfaces[0], witness[:2]), _point(surfaces[1], witness[2:]) + assert np.linalg.norm(a-b) <= ATOL + assert np.min(reference['normals'] @ (a-reference['center'])) > ATOL + assert len(region.uv1_loops) == len(region.uv2_loops) == 1 + first_loop, second_loop = region.uv1_loops[0], region.uv2_loops[0] + assert np.shape(first_loop) == np.shape(second_loop) + for uv1, uv2 in zip(first_loop, second_loop): + assert np.linalg.norm(_point(surfaces[0], uv1)-_point(surfaces[1], uv2)) <= ATOL + for surface, uv in zip(surfaces, (first_loop, second_loop)): + # At a collapsed pole, different u values denote the same point. + assert np.linalg.norm(_point(surface, uv[0])-_point(surface, uv[-1])) <= ATOL + + +def _assemble(surfaces): + from mmcore.nurbs._nurbs_param_tol import bez_surface_param_tolerance + homogeneous = tuple(map(_homogeneous, surfaces)) + parameter_tolerance = np.concatenate([ + np.asarray(bez_surface_param_tolerance(net, ATOL, rational=True)) + for net in homogeneous]) + assembled = assemble_overlap_regions(*homogeneous, atol=ATOL, + ptol4=parameter_tolerance) + return {'overlap_regions': assembled['regions'], 'branches': assembled['rim_branches']} + + +@pytest.mark.parametrize('variant', ['identity', 'swap', 'reverse_u']) +def test_case17_public_overlap_has_the_whole_spherical_quadrilateral(variant): + from examples.ssx.case_17 import s1, s2 + from mmcore.numeric.intersection.ssx import nurbs_ssx + fits = [_fit_octant_frame(surface) for surface in (s1, s2)] + assert np.linalg.norm(fits[0][0]-fits[1][0]) < .01*ATOL + assert abs(fits[0][1]-fits[1][1]) < .01*ATOL + center, radius = np.mean([f[0] for f in fits], axis=0), np.mean([f[1] for f in fits]) + reference = _reference(center, radius, [f[2] for f in fits]) + assert len(reference['vertices']) == 4 # property of the independent fixture + assert [sum(owners == {owner} for owners in reference['owners']) + for owner in (0, 1)] == [2, 2] + # Neither collapsed pole is part of this particular common region. + assert all(np.min(reference['normals'] @ fit[2][:, 2])*radius < -ATOL + for fit in fits) + surfaces = [s1, s2] + if variant == 'swap': + surfaces.reverse() + elif variant == 'reverse_u': + surfaces[1] = _reverse_u(surfaces[1]) + result = nurbs_ssx(*surfaces, atol=ATOL) + assert result['complete'] is True, result['status'] + assert result['status']['reasons'] == [] + _assert_region(result, surfaces, reference, + normal_agreement=-1 if variant == 'reverse_u' else 1) + assert result['points'] == [] and result['singularities'] == [] + + +@pytest.mark.parametrize('variant', ['identity', 'swap', 'reverse_u']) +def test_rotated_octant_overlap_keeps_all_six_independent_rims(variant): + frames = [np.eye(3), _rotation([1., 2., 3.], .35)] + surfaces = [_octant(frame, center=(2., -3., .5)) for frame in frames] + reference = _reference([2., -3., .5], 1., frames) + assert len(reference['vertices']) == 6 + if variant == 'swap': + surfaces.reverse() + elif variant == 'reverse_u': + surfaces[1] = _reverse_u(surfaces[1]) + _assert_region(_assemble(surfaces), surfaces, reference, + normal_agreement=-1 if variant == 'reverse_u' else 1) + + +def test_narrow_spherical_overlap_keeps_distinct_nearby_corners_and_shared_pole(): + frames = [np.eye(3), _rotation([0., 0., 1.], np.pi/2-.02)] + surfaces = [_octant(frame) for frame in frames] + reference = _reference([0., 0., 0.], 1., frames) + # The equatorial corners are 20*atol apart, even though both arcs meet + # at a common collapsed pole. An enlarged join radius must not merge them. + assert min(np.linalg.norm(a-b) for a, b in combinations(reference['vertices'], 2)) > 19*ATOL + _assert_region(_assemble(surfaces), surfaces, reference) + + +@pytest.mark.parametrize('kind', ['radius_gap', 'shared_meridian', 'shared_pole']) +def test_lower_dimensional_or_separated_spherical_contacts_do_not_make_a_region(kind): + first = _octant(np.eye(3)) + if kind == 'radius_gap': + second = _octant(np.eye(3), radius=1.+4*ATOL) + else: + angle = np.pi/2 + (.02 if kind == 'shared_pole' else 0.) + second = _octant(_rotation([0., 0., 1.], angle)) + assert _assemble((first, second))['overlap_regions'] == [] diff --git a/tests/test_ssx_spherical_reduction.py b/tests/test_ssx_spherical_reduction.py new file mode 100644 index 00000000..f1c436f1 --- /dev/null +++ b/tests/test_ssx_spherical_reduction.py @@ -0,0 +1,246 @@ +"""Whole-chart spherical CAD reduction, including ineligible lookalikes.""" +import numpy as np +import pytest + +from mmcore.numeric._bezier_common import eval_surface +from mmcore.numeric._work_budget import SoftWorkBudget +from mmcore.numeric.intersection.ssx import _ssx_spherical_overlap as spherical +from mmcore.numeric.intersection.ssx import _ssx5_overlap as overlap +from test_ssx_spherical_overlap import ( + ATOL, _octant, _rotation, _reference, _homogeneous, _assert_region, +) + + +def _budget(**kwargs): + return SoftWorkBudget(max_cells=kwargs.pop('max_cells', 60000), + max_csx_calls=1000, max_output_items=kwargs.pop('max_output_items', 1024), + **kwargs) + + +def _pair(radius=1., center=(0., 0., 0.)): + frames = (np.eye(3), _rotation([1., 2., 3.], .35)) + return tuple(_octant(frame, radius, center) for frame in frames), frames + + +def test_case17_direct_reduction_returns_all_four_actual_rims(): + from examples.ssx.case_17 import s1, s2 + from test_ssx_spherical_overlap import _fit_octant_frame + surfaces = (s1, s2) + fits = [_fit_octant_frame(surface) for surface in surfaces] + center = .5*(fits[0][0]+fits[1][0]) + radius = .5*(fits[0][1]+fits[1][1]) + result = spherical.try_spherical_overlap(*map(_homogeneous, surfaces), ATOL, _budget()) + assert result is not None + assert result._c1_resolved and result._c3_resolved + assert '_c1_resolved' not in result and '_c3_resolved' not in result + _assert_region(result, surfaces, _reference(center, radius, [fit[2] for fit in fits])) + + +@pytest.mark.parametrize('transform', ['identity', 'reverse', 'transpose', 'gauge']) +def test_rotated_scaled_charts_preserve_complete_footprint(transform): + surfaces, frames = _pair(radius=2.5, center=(20., -10., 5.)) + first, second = map(_homogeneous, surfaces) + if transform == 'reverse': + first = first[::-1].copy() + elif transform == 'transpose': + first = first.swapaxes(0, 1).copy() + elif transform == 'gauge': + first = first*1e-7 + second = second*13. + budget = _budget() + result = spherical.try_spherical_overlap(first, second, ATOL, budget) + assert result is not None and len(result['overlap_regions']) == 1 + assert result._c1_resolved and result._c3_resolved + assert budget.result_fields()['complete'] + # Coverage oracle uses geometry, independent of the altered parameter chart. + from examples.ssx.ssx5_analytic_audit import distances_to_polylines + reference = _reference(np.array([20., -10., 5.]), 2.5, frames) + paths = [np.asarray(branch.curve[1]) for branch in result['branches']] + for arc in reference['arcs']: + assert distances_to_polylines(arc, paths).max() <= ATOL + for branch in result['branches']: + for q, xyz in zip(*branch.curve): + assert np.linalg.norm(eval_surface(first, *q[:2], rational=True)-xyz) <= ATOL + assert np.linalg.norm(eval_surface(second, *q[2:], rational=True)-xyz) <= ATOL + assert result['overlap_regions'][0].normal_agreement == (-1 if transform in ('reverse', 'transpose') else 1) + + +def test_small_control_and_weight_perturbations_are_judged_at_cad_accuracy(): + surfaces, frames = _pair() + first, second = map(_homogeneous, surfaces) + first[1, 1, :3] += first[1, 1, 3]*np.array([1e-5, -1e-5, .5e-5]) + second[1, 1] *= 1.+1e-5 + budget = _budget() + result = spherical.try_spherical_overlap(first, second, ATOL, budget) + assert result is not None and len(result['overlap_regions']) == 1 + assert budget.result_fields()['complete'] + + +def test_hidden_interior_bulge_cannot_pass_matching_boundary_and_witness(): + first = _homogeneous(_octant(np.eye(3))) + second = first.copy() + point = eval_surface(first, .5, .5, rational=True) + second[1, 1] += 2.*np.r_[point, 1.] + # Boundaries and the central point coincide; the interior does not. + assert np.linalg.norm(eval_surface(first, .5, .5, rational=True) + - eval_surface(second, .5, .5, rational=True)) <= ATOL + assert abs(np.linalg.norm(eval_surface(second, .25, .25, rational=True))-1.) > ATOL + budget = _budget() + assert spherical.try_spherical_overlap(first, second, ATOL, budget) is None + assert not budget.exhausted + + +def test_same_cartesian_controls_with_different_weights_need_whole_map_check(): + first = _homogeneous(_octant(np.eye(3))) + second = first.copy() + second[1, 1] *= 1.2 + np.testing.assert_allclose(first[..., :3]/first[..., 3, None], + second[..., :3]/second[..., 3, None]) + assert spherical.try_spherical_overlap(first, second, ATOL, _budget()) is None + + +def test_radius_gap_and_unsupported_shape_decline_without_absence_claim(): + first = _homogeneous(_octant(np.eye(3))) + second = _homogeneous(_octant(np.eye(3), radius=1.+4*ATOL)) + assert spherical.try_spherical_overlap(first, second, ATOL, _budget()) is None + higher_degree_shape = np.concatenate((first, first[-1:]), axis=0) + assert spherical.try_spherical_overlap(first, higher_degree_shape, ATOL, _budget()) is None + + +def test_curve_only_and_disjoint_footprints_stay_with_general_ssx(): + first = _homogeneous(_octant(np.eye(3))) + for angle in (np.pi/2, np.pi/2+.02): + second = _homogeneous(_octant(_rotation([0., 0., 1.], angle))) + assert spherical.try_spherical_overlap(first, second, ATOL, _budget()) is None + + +def test_output_denial_keeps_rims_without_dangling_region_references(): + surfaces, _ = _pair() + budget = _budget(max_output_items=2) + result = spherical.try_spherical_overlap(*map(_homogeneous, surfaces), ATOL, budget) + assert result is not None and len(result['branches']) == 2 + assert result['overlap_regions'] == [] + assert 'output_cap' in budget.result_fields()['status']['reasons'] + + +def test_work_denial_during_assembly_retains_completed_validated_rim(monkeypatch): + surfaces, _ = _pair() + budget = _budget() + original = overlap._adaptive_curved_rim + completed = [] + def stop_after_first(*args, **kwargs): + result = original(*args, **kwargs) + if result is not None and not completed: + completed.append(result) + budget.max_cells = budget.cells_processed + return result + monkeypatch.setattr(overlap, '_adaptive_curved_rim', stop_after_first) + result = spherical.try_spherical_overlap(*map(_homogeneous, surfaces), ATOL, budget) + assert budget.exhausted and result is not None + assert len(result['branches']) == 1 and result['overlap_regions'] == [] + np.testing.assert_array_equal(result['branches'][0].curve[1], completed[0]['xyz']) + + +def test_setup_denial_leaves_the_general_search_stopped_by_its_budget(): + surfaces, _ = _pair() + budget = _budget(max_cells=0) + assert spherical.try_spherical_overlap(*map(_homogeneous, surfaces), ATOL, budget) is None + assert budget.exhausted and budget.cells_processed == 0 + + +@pytest.mark.parametrize('angle', [0., np.pi/2-.02]) +def test_identical_and_thin_octants_have_the_complete_boundary(angle): + frames = (np.eye(3), _rotation([0., 0., 1.], angle)) + surfaces = tuple(_octant(frame) for frame in frames) + budget = _budget() + result = spherical.try_spherical_overlap(*map(_homogeneous, surfaces), ATOL, budget) + assert result is not None and budget.result_fields()['complete'] + assert result._c1_resolved and result._c3_resolved + _assert_region(result, surfaces, _reference(np.zeros(3), 1., frames)) + + +def test_singularity_guard_work_denial_preserves_the_validated_region(monkeypatch): + surfaces, _ = _pair() + budget = _budget() + original = spherical._resolve_spherical_c1 + def deny_guard(sources, charts, result, atol, work_budget): + work_budget.max_cells = work_budget.cells_processed + return original(sources, charts, result, atol, work_budget) + monkeypatch.setattr(spherical, '_resolve_spherical_c1', deny_guard) + result = spherical.try_spherical_overlap(*map(_homogeneous, surfaces), ATOL, budget) + assert result is not None and len(result['overlap_regions']) == 1 + assert len(result['branches']) == 6 + assert result._c3_resolved and not result._c1_resolved + assert budget.exhausted + assert 'work_budget' in budget.result_fields()['status']['reasons'] + + +def test_rank_defect_outside_the_pole_cap_keeps_c1_enabled(): + surface = _homogeneous(_octant(np.eye(3))) + budget = _budget() + chart = spherical._propose_chart(surface, budget) + result = spherical.try_spherical_overlap(surface, surface.copy(), ATOL, budget) + assert result is not None and result._c1_resolved + folded = surface.copy() + folded[-1] = folded[0] + # This actual biquadratic map has S_u=0 on u=.5, including the + # equatorial edge far from the pole. It still has a collapsed pole. + from mmcore.numeric._bezier_common import eval_surface_d1 + _, derivative_u, _ = eval_surface_d1(folded, .5, 0., rational=True) + assert np.linalg.norm(derivative_u) <= np.finfo(float).eps + assert not spherical._resolve_spherical_c1( + (folded, surface), (chart, chart), result, ATOL, _budget()) + + +def test_unrepresented_pole_cap_on_the_other_surface_keeps_c1_enabled(): + surface = _homogeneous(_octant(np.eye(3))) + chart = spherical._propose_chart(surface, _budget()) + # An unrelated corner cannot represent this pole, and the identical + # other patch provides no separating support direction. + from mmcore.numeric.intersection.ssx._ssx_substrate import SSXBranch + result = dict(overlap_regions=[object()], branches=[SSXBranch(curve=( + np.zeros((2, 4)), np.array([[1., 0., 0.], [0., 1., 0.]])))]) + assert not spherical._resolve_spherical_c1( + (surface, surface), (chart, chart), result, ATOL, _budget()) + + +def test_nonspherical_proposal_does_not_spend_the_general_solver_allowance(): + first = np.array([[[i/2, j/2, 0., 1.] for j in range(3)] for i in range(3)]) + second = np.array([[[i/2, .5, j/2-.5, 1.] for j in range(3)] for i in range(3)]) + budget = _budget(max_cells=1000) + assert spherical.try_spherical_overlap(first, second, ATOL, budget) is None + assert not budget.exhausted + from mmcore.numeric.intersection.ssx._bez_ssx5 import bez_ssx + result = bez_ssx(first, second, ATOL, rational=True, max_cells=1000) + assert result['complete'] and len(result['branches']) == 1 + xyz = np.asarray(result['branches'][0].curve[1]) + assert abs(xyz[:, 0].min()) <= ATOL and abs(xyz[:, 0].max()-1.) <= ATOL + + +def test_thin_perturbed_footprint_declines_when_chart_error_is_amplified(monkeypatch): + # Twenty-atol corner separation is resolved for the unperturbed charts + # above. A small interior perturbation can nevertheless make the common + # carrier's angular footprint uncertain at a nearly opposed pair of faces. + first = _homogeneous(_octant(np.eye(3))) + second = _homogeneous(_octant(_rotation([0., 0., 1.], np.pi/2-.02))) + first[1, 1, 2] += first[1, 1, 3]*1e-5 + def forbidden(*args, **kwargs): + raise AssertionError('uncertain footprint entered region construction') + monkeypatch.setattr(spherical, 'assemble_overlap_regions', forbidden) + budget = _budget() + assert spherical.try_spherical_overlap(first, second, ATOL, budget) is None + assert not budget.exhausted + assert budget.cell_counts.get('spherical_bound', 0) > 0 + + +@pytest.mark.parametrize('angle', [0., np.pi/2-.02]) +def test_public_region_does_not_retype_shared_pole_parameters_as_singularities(angle): + from mmcore.numeric.intersection.ssx import nurbs_ssx + from test_ssx_spherical_overlap import _assert_region, _reference + + frames = [np.eye(3), _rotation([0., 0., 1.], angle)] + surfaces = tuple(_octant(frame) for frame in frames) + result = nurbs_ssx(*surfaces, atol=ATOL) + _assert_region(result, surfaces, _reference([0., 0., 0.], 1., frames)) + assert result['complete'] and result['status']['reasons'] == [] + assert result['points'] == [] and result['singularities'] == [] From 5c5149b70fc4dbbb55e315f73baec9ce9800007a Mon Sep 17 00:00:00 2001 From: sth-v Date: Fri, 25 Sep 2026 16:43:59 +0300 Subject: [PATCH 2/3] Preserve separate contacts and original apex identity during overlap search --- .../2026-09-25-case17-curved-overlap.md | 16 ++ mmcore/numeric/intersection/csx/_bez_csx4.py | 58 +++++++- mmcore/numeric/intersection/ssx/_bez_ssx5.py | 72 ++++++--- tests/test_csx_cad_span_contacts.py | 50 +++++++ ...test_ssx_collapsed_fiber_classification.py | 139 ++++++++++++++++++ 5 files changed, 306 insertions(+), 29 deletions(-) create mode 100644 tests/test_ssx_collapsed_fiber_classification.py diff --git a/docs/superpowers/reviews/2026-09-25-case17-curved-overlap.md b/docs/superpowers/reviews/2026-09-25-case17-curved-overlap.md index c305ee85..b4c868b8 100644 --- a/docs/superpowers/reviews/2026-09-25-case17-curved-overlap.md +++ b/docs/superpowers/reviews/2026-09-25-case17-curved-overlap.md @@ -34,6 +34,22 @@ no overlap region. membership and domain-end checks. The captured cut changed from 22.54 seconds / 67,104 work units to 0.233 seconds / 1,126 units. +The first full-suite run caught two integration regressions before merge. +The candidate-interval path could mistake following a surface edge for +leaving its domain, joining two lines eight tolerances apart. Domain +pinning must establish an actual exit; touching an edge and returning to +the interior does not qualify either. Both independent negative controls +retain the separate roots. + +Changed subdivision also exposed an existing coordinate-dependent apex +classification: a point at an original collapsed cone edge could sit far +from that edge in a very thin child's rescaled parameters, acquiring a +false isolated tangent-point label. Cells now retain their original +surface context. The type check uses the fixed parameter axis and distance +to the collapsed edge's image, while retaining local checks for an interior +collapsed isoline exposed by a split. The cone generator and genuine +isolated touches remain covered by their existing geometry assertions. + Correcting the boundary still left the ordinary one-dimensional search subdividing a two-dimensional coincident area. A separate geometric reduction now handles supported rational quadratic spherical octants. diff --git a/mmcore/numeric/intersection/csx/_bez_csx4.py b/mmcore/numeric/intersection/csx/_bez_csx4.py index fa850b9f..5a868031 100644 --- a/mmcore/numeric/intersection/csx/_bez_csx4.py +++ b/mmcore/numeric/intersection/csx/_bez_csx4.py @@ -858,19 +858,63 @@ def _pinned(t_end, outward): if t_end <= edge_pad or t_end >= 1.0 - edge_pad: return True _p, u, v, _d = _member(t_end, (0.5, 0.5)) - if min(u, 1.0 - u) <= 1e-6 or min(v, 1.0 - v) <= 1e-6: - return True # The tolerance boundary and the uv-domain exit can COINCIDE # (measured on user data: d crosses atol at almost the same t # where the projected path leaves through u=1, so the # projection AT the refined boundary is still interior). - # Probe one grid step OUTWARD: a domain-clipped span clamps to - # an edge there; a genuine interior fade-out (the offset-twin - # signature, which must stay refused) does not. + # A domain pin must actually ENTER through that edge. Merely + # lying on an edge does not pin t: a curve can follow v=0 past + # two distinct isolated crossings, and an endpoint-seeded + # window between those roots must not become an overlap. probe_step = (parameter_high - parameter_low) / (len(ts) - 1.) t_probe = min(1.0, max(0.0, t_end + outward * probe_step)) - _p, u, v, _d = _member(t_probe, (u, v)) - return (min(u, 1.0 - u) <= 1e-6 or min(v, 1.0 - v) <= 1e-6) + probe_point, u_out, v_out, projected_distance = _member(t_probe, (u, v)) + exit_roundoff = max(tiny, 512.*np.finfo(float).eps*max( + 1., float(np.max(abs(probe_point))), + float(np.max(abs(points))))) + # An accepted tolerance fringe may already clamp several + # samples before its endpoint, so use the whole run for the + # interior witness rather than only the adjacent grid sample. + for coordinate, side in ((0, 0.), (0, 1.), (1, 0.), (1, 1.)): + inner = uvs[k0:k1+1, coordinate] + if side == 0. and float(inner.max()) > 0.: + sense = 1. + elif side == 1. and float(inner.min()) < 1.: + sense = -1. + else: + continue + # Test an actual boundary-constrained candidate. The + # coupled projector can stall just inside the edge + # when clipping its step worsens the free-coordinate fit; + # its final UV bits are not an active-set certificate. + boundary_uv = [u_out, v_out] + boundary_uv[coordinate] = side + surface_point, du, dv = eval_surface_d1( + S64, *boundary_uv, rational=rational) + boundary_distance = float(np.linalg.norm(probe_point-surface_point)) + if not (np.isfinite(boundary_distance) + and np.isfinite(projected_distance) + and boundary_distance <= projected_distance+exit_roundoff): + continue + # Clamping must be active in geometry too. A curve that + # only grazes an edge can remain inside on both sides; + # small positive UV values are not a domain exit. Remove + # the along-edge tangent before measuring outward motion, + # so a normal offset or projection lag cannot supply it. + across, along = (du, dv) if coordinate == 0 else (dv, du) + along_sq = float(along@along) + if along_sq <= 0.: + continue + inward = sense*(across-along*float(across@along)/along_sq) + inward_length = float(np.linalg.norm(inward)) + derivative_roundoff = 512.*np.finfo(float).eps*max( + float(np.linalg.norm(du)), float(np.linalg.norm(dv))) + if not np.isfinite(inward_length) or inward_length <= derivative_roundoff: + continue + exit_distance = float((probe_point-surface_point)@inward/inward_length) + if exit_distance < -exit_roundoff: + return True + return False if not (_pinned(t_lo, -1.0) and _pinned(t_hi, +1.0)): continue # A within-atol span remains CAD-coincident when rounding puts diff --git a/mmcore/numeric/intersection/ssx/_bez_ssx5.py b/mmcore/numeric/intersection/ssx/_bez_ssx5.py index 0182a04f..4398417f 100644 --- a/mmcore/numeric/intersection/ssx/_bez_ssx5.py +++ b/mmcore/numeric/intersection/ssx/_bez_ssx5.py @@ -562,25 +562,57 @@ def _weight_net_uniform(S) -> bool: def _on_collapsed_boundary_fiber( - S, u, v, rational=True, *, param_tol: float = 1e-10) -> bool: + S, u, v, rational=True, *, param_tol=1e-10, + distance_tol=None) -> bool: """True only for a parameter point on an identically collapsed edge. This is deliberately narrower than ``Sigma=0``: a C1 point can still be regular on the 4D Psi curve and may be a required branch seed. We suppress only a certified positive-dimensional boundary preimage. """ - eps = max(0.0, float(param_tol)) - if abs(u) <= eps and _curve_geometry_collapsed(S[0, :, :], rational): - return True - if abs(u - 1.0) <= eps and _curve_geometry_collapsed(S[-1, :, :], rational): - return True - if abs(v) <= eps and _curve_geometry_collapsed(S[:, 0, :], rational): - return True - if abs(v - 1.0) <= eps and _curve_geometry_collapsed(S[:, -1, :], rational): - return True + eps = np.maximum(0., np.broadcast_to(param_tol, (2,))) + for coordinate, tolerance, low, high in ( + (u, eps[0], S[0, :, :], S[-1, :, :]), + (v, eps[1], S[:, 0, :], S[:, -1, :])): + for side, edge in ((0., low), (1., high)): + if (not abs(coordinate - side) <= tolerance + or not _curve_geometry_collapsed(edge, rational)): + continue + if distance_tol is None: + return True + point = eval_surface(S, u, v, rational=rational) + pole = eval_curve(edge, .5, rational=rational) + distance = float(np.linalg.norm(point - pole)) + if np.isfinite(distance) and distance <= distance_tol: + return True return False +def _cell_root_on_collapsed_fiber(cell, root, atol, global_ptol): + """Recognize original apex fibers after a thin child rescales UV. + + A near-apex root can be far from the child's local parameter edge even + though its model-space image is at the original collapsed endpoint. + Check each fixed axis in the original chart and require spatial + proximity too; nearby parameters alone do not identify a CAD point. + """ + def on_fiber(surfaces, parameters, tolerances): + return any(_on_collapsed_boundary_fiber( + surface, *parameters[2*i:2*i+2], rational=True, + param_tol=tolerances[2*i:2*i+2], distance_tol=atol) + for i, surface in enumerate(surfaces)) + + # A split can also expose an interior collapsed isoline. Retain that + # local evidence even when it is not an edge of either original chart. + if on_fiber((cell.g1.surface, cell.g2.surface), np.asarray(root), + _cell_ptol4(cell, atol)): + return True + sources = getattr(cell, 'source_surfaces', None) + return sources is not None and on_fiber( + sources, _local_to_global(np.asarray(root), cell.box), + np.asarray(global_ptol)) + + def _canonicalize_collapsed_fiber_params( S, uv, anchor_uv, *, rational=True, param_tol=1e-10): """Choose the limiting branch representative on collapsed edge fibers. @@ -1295,12 +1327,7 @@ def _emit_tangent_roots(cell, atol, unify_tol, all_singularities, stuv_g = _local_to_global(np.asarray(xw), cell.box) if _stuv_in_overlap_boxes(stuv_g, overlap_boxes): continue - if (_on_collapsed_boundary_fiber( - cell.g1.surface, xw[0], xw[1], rational=True, - param_tol=float(max(local_ptol4[0], local_ptol4[1]))) - or _on_collapsed_boundary_fiber( - cell.g2.surface, xw[2], xw[3], rational=True, - param_tol=float(max(local_ptol4[2], local_ptol4[3])))): + if _cell_root_on_collapsed_fiber(cell, xw, atol, unify_tol): # Every value of the collapsed edge's free parameter maps to # this endpoint. It belongs to a C1/parameter-fiber set, never # an isolated C2 tangent point; relative normal tests are @@ -1876,12 +1903,7 @@ def _xyz(x): stuv_g = _local_to_global(np.asarray(xw), cell.box) if _stuv_in_overlap_boxes(stuv_g, overlap_boxes): continue - if (_on_collapsed_boundary_fiber( - cell.g1.surface, xw[0], xw[1], rational=True, - param_tol=float(max(ptol4[0], ptol4[1]))) - or _on_collapsed_boundary_fiber( - cell.g2.surface, xw[2], xw[3], rational=True, - param_tol=float(max(ptol4[2], ptol4[3])))): + if _cell_root_on_collapsed_fiber(cell, xw, atol, unify_tol): continue if _normals_degenerate_at(cell.g1.surface, cell.g2.surface, xw): continue # L15: Sigma=0 root — C1 candidate, not a C2 touch @@ -7086,6 +7108,9 @@ class _Cell: # remaining domain still needs searching; this flag only enables # coverage checks against geometry that has actually been retained. tangency_remainder: bool = False + # Original charts retain collapsed-edge identity when a tiny descendant + # rescales a near-apex point far from its local UV boundary. + source_surfaces: Optional[tuple] = None def _probe_children(cell): @@ -7139,6 +7164,7 @@ def _split2(T): probe_only=True, work_budget=cell.work_budget, csx_fn=cell.csx_fn, tangency_remainder=cell.tangency_remainder, + source_surfaces=cell.source_surfaces, )) return out @@ -7612,6 +7638,7 @@ def _run_top_boundary_csx(curve, surface, **kwargs): new_crossings=list(crossings), F_sq=F_sq_top, w_scale=w_scale_top, work_budget=budget, csx_fn=_run_csx, + source_surfaces=(S1_h_top, S2_h_top), ) top_cell.partitions = _build_outer_partitions(top_cell) @@ -8421,6 +8448,7 @@ def _known_tangent_crossings(cell): w_scale=cell.w_scale, work_budget=budget, csx_fn=_run_csx, tangency_remainder=cell.tangency_remainder, + source_surfaces=cell.source_surfaces, ) scell.partitions = _build_cell_partitions(scell) for c in sub_cx: diff --git a/tests/test_csx_cad_span_contacts.py b/tests/test_csx_cad_span_contacts.py index a4ce319d..7be28258 100644 --- a/tests/test_csx_cad_span_contacts.py +++ b/tests/test_csx_cad_span_contacts.py @@ -167,6 +167,56 @@ def test_two_boundary_contacts_do_not_promote_an_out_of_tolerance_bulge(): _validate_contacts(curve, surface, result['isolated']) +@pytest.mark.parametrize('side', [0., 1.]) +def test_following_a_surface_edge_does_not_pin_an_interior_tolerance_window(side): + # These two roots are close in curve parameters but eight CAD atols + # apart in model space. The entire owner curve follows a target UV + # edge, which does not delimit the interval between the two roots. + a, b = .49, .4905 + z = 32.*np.array([a*b, a*b-(a+b)/2., (1.-a)*(1.-b)]) + curve = np.column_stack((np.array([0., 8., 16.]), np.full(3, side), z)) + surface = np.array([[[16.*u, v, 0.] for v in (0., 1.)] + for u in (0., 1.)]) + result = csx.bez_csx(curve, surface, atol=1e-3, rational=False) + assert not result['budget_exhausted'] + assert result['overlaps'] == [] + _validate_contacts(curve, surface, result['isolated']) + for x in (7.84, 7.848): + assert any(np.linalg.norm(point['point']-[x, side, 0.]) <= 1e-3 + for point in result['isolated']) + + +@pytest.mark.parametrize('side', [0., 1.]) +def test_touching_a_domain_edge_without_exiting_does_not_pin_a_span(side): + # The UV path touches the same edge at .4 and .6 but stays inside on + # BOTH sides. Its outward probes are only 4e-7 from that edge; treating + # "UV near zero" as an active clamp incorrectly promoted this window. + degree = 4 + product = np.polynomial.polynomial.polyfromroots([.4, .6]) + y = np.polynomial.polynomial.polymul(product, product) + if side: + y = -y + y[0] += 1. + def coefficients(power): + power = np.pad(power, (0, degree+1-len(power))) + return [sum(power[k]*comb(i, k)/comb(degree, k) + for k in range(i+1)) for i in range(degree+1)] + curve = np.column_stack((coefficients([0., 16.]), coefficients(y), + coefficients(.05*product))) + surface = np.array([[[16.*u, v, 0.] for v in (0., 1.)] + for u in (0., 1.)]) + assert csx._tolerance_csx_overlap_certificate( + curve, surface, 1e-3, False, 1e-3/16., None, + parameter_range=(.4, .6)) is None + result = csx.bez_csx(curve, surface, atol=1e-3, rational=False) + assert result['overlaps'] == [] + assert not result['budget_exhausted'] + _validate_contacts(curve, surface, result['isolated']) + for x in (6.4, 9.6): + assert any(np.linalg.norm(point['point']-[x, side, 0.]) <= 1e-3 + for point in result['isolated']) + + @pytest.mark.parametrize('owner,axis,side', [ (owner, axis, side) for owner in (0, 1) for axis in (0, 1) for side in (0, 1)]) def test_rounded_spherical_boundary_spans_cover_the_independent_common_rim(owner, axis, side): diff --git a/tests/test_ssx_collapsed_fiber_classification.py b/tests/test_ssx_collapsed_fiber_classification.py new file mode 100644 index 00000000..f040e3aa --- /dev/null +++ b/tests/test_ssx_collapsed_fiber_classification.py @@ -0,0 +1,139 @@ +"""CAD-scale apex classification retains original chart and sheet identity.""" +from types import SimpleNamespace + +import numpy as np +import pytest + +from mmcore.numeric._bezier_common import eval_surface, eval_surface_d1, restrict_net_axis_v +from mmcore.numeric.intersection.ssx import _bez_ssx5 as ssx + +ATOL = 1e-3 +UNIT_BOX = ((0., 1.),)*4 + + +def _plane(): + return np.array([[[u, v, 0., 1.] for v in (0., 1.)] for u in (0., 1.)]) + + +def _fan(scale=1.): + # S(u,v)=(u*v,v,0), with every point of v=0 mapping to the apex. + return np.array([[[scale*u*v, scale*v, 0., 1.] + for v in (0., 1.)] for u in (0., 1.)]) + + +def _cell(sources, box=UNIT_BOX, original_context=True): + children = [] + for owner, surface in enumerate(sources): + child = surface.copy() + for axis in (0, 1): + child = restrict_net_axis_v(child, axis, *box[2*owner+axis], 0., 1.) + children.append(SimpleNamespace(surface=child)) + cell = SimpleNamespace(g1=children[0], g2=children[1], box=box) + if original_context: + cell.source_surfaces = tuple(sources) + return cell + + +@pytest.mark.parametrize('swap', [False, True]) +@pytest.mark.parametrize('transpose', [False, True]) +@pytest.mark.parametrize('apex_side', [0, 1]) +def test_thin_child_uses_original_apex_chart(swap, transpose, apex_side): + first, second = _fan(), _plane() + near, far = 1e-12, 1e-8 + parameters = np.array([.35, .5*(near+far)]) + first_box = [(0., 1.), (near, far)] + if apex_side == 1: + first = first[:, ::-1].copy() + parameters[1] = 1.-parameters[1] + first_box[1] = (1.-far, 1.-near) + if transpose: + first = first.swapaxes(0, 1).copy() + parameters = parameters[::-1].copy() + first_box.reverse() + point = eval_surface(first, *parameters, rational=True) + assert np.linalg.norm(point) <= ATOL + global_root = np.r_[parameters, point[:2]] + sources = (first, second) + box = tuple(first_box)+((0., 1.), (0., 1.)) + if swap: + sources = sources[::-1] + global_root = global_root[[2, 3, 0, 1]] + box = box[2:]+box[:2] + cell = _cell(sources, box) + local_root = ssx._global_to_local(global_root, box) + owner = 1 if swap else 0 + local_surface = (cell.g1.surface, cell.g2.surface)[owner] + local_parameters = local_root[2*owner:2*owner+2] + local_tolerance = ssx._cell_ptol4(cell, ATOL)[2*owner:2*owner+2] + # Neither child boundary is collapsed: the original apex lies just + # outside this child's interval. Its CAD neighborhood still belongs + # to that original fiber, not to a new isolated C2 touch. + assert not ssx._on_collapsed_boundary_fiber( + local_surface, *local_parameters, param_tol=local_tolerance, + distance_tol=ATOL) + assert ssx._cell_root_on_collapsed_fiber(cell, local_root, ATOL, + np.full(4, ATOL)) + + +def test_parameter_nearness_does_not_hide_a_spatially_distant_point(): + sources = (_fan(scale=1000.), _plane()) + box = ((0., 1.), (0., 1e-3), (0., 1.), (0., 1.)) + global_root = np.array([.5, .5e-3, .25, .5]) + point = eval_surface(sources[0], *global_root[:2], rational=True) + assert np.linalg.norm(point) > ATOL + assert ssx._on_collapsed_boundary_fiber(sources[0], *global_root[:2], + param_tol=np.full(2, ATOL)) + assert not ssx._cell_root_on_collapsed_fiber( + _cell(sources, box), ssx._global_to_local(global_root, box), + ATOL, np.full(4, ATOL)) + + +def test_fixed_axis_keeps_its_own_parameter_tolerance(): + # f(v)=1000*v^3 has a small spatial image near zero but large global + # speed. The u direction is independently scaled down by 1e-5. + f = np.array([0., 0., 0., 1000.]) + first = np.array([[[1e-5*u*y, y, 0., 1.] for y in f] for u in (0., 1.)]) + sources = (first, _plane()) + box = ((0., 1.), (.004, .006), (0., 1.), (0., 1.)) + q = np.array([.5, .005]) + point = eval_surface(first, *q, rational=True) + assert np.linalg.norm(point) < ATOL + global_root = np.r_[q, point[:2]] + tolerances = np.array([.1, 1e-6, ATOL, ATOL]) + assert not ssx._cell_root_on_collapsed_fiber( + _cell(sources, box), ssx._global_to_local(global_root, box), ATOL, tolerances) + + +def test_regular_other_preimage_at_the_same_xyz_is_not_the_apex_fiber(): + # S(u,v)=(v*(u-.5), v*(v-.5), 0). The boundary v=0 is collapsed, + # while (.5,.5) maps to the same XYZ with two independent derivatives. + first = np.array([[[v*(u-.5), y, 0., 1.] + for v, y in zip((0., .5, 1.), (0., -.25, .5))] + for u in (0., 1.)]) + point, du, dv = eval_surface_d1(first, .5, .5, rational=True) + assert np.linalg.norm(point) <= ATOL + assert np.linalg.norm(np.cross(du, dv)) > .1 + box = ((.4, .6), (.4, .6), (0., 1.), (0., 1.)) + global_root = np.array([.5, .5, 0., 0.]) + assert not ssx._cell_root_on_collapsed_fiber( + _cell((first, _plane()), box), ssx._global_to_local(global_root, box), + ATOL, np.full(4, ATOL)) + + +def test_child_boundary_can_expose_an_original_interior_collapsed_isoline(): + first = np.array([[[u*(v-.5), v-.5, 0., 1.] + for v in (0., 1.)] for u in (0., 1.)]) + box = ((0., 1.), (.5, 1.), (0., 1.), (0., 1.)) + global_root = np.array([.3, .5, 0., 0.]) + assert not ssx._on_collapsed_boundary_fiber(first, *global_root[:2], + param_tol=np.full(2, ATOL), + distance_tol=ATOL) + assert ssx._cell_root_on_collapsed_fiber( + _cell((first, _plane()), box), ssx._global_to_local(global_root, box), + ATOL, np.full(4, ATOL)) + + +def test_standalone_cell_keeps_local_collapsed_edge_evidence(): + cell = _cell((_fan(), _plane()), original_context=False) + assert ssx._cell_root_on_collapsed_fiber(cell, np.array([.3, 0., 0., 0.]), + ATOL, np.full(4, ATOL)) From 580b016bdc800c7649ad0e6adac910f2aa776bb0 Mon Sep 17 00:00:00 2001 From: sth-v Date: Fri, 25 Sep 2026 17:07:38 +0300 Subject: [PATCH 3/3] Distinguish isolated endpoint contacts from whole CAD overlap spans --- .../2026-09-25-case17-curved-overlap.md | 9 ++++ tests/test_ssx_cut_overlap_contacts.py | 54 ++++++++++++++----- 2 files changed, 51 insertions(+), 12 deletions(-) diff --git a/docs/superpowers/reviews/2026-09-25-case17-curved-overlap.md b/docs/superpowers/reviews/2026-09-25-case17-curved-overlap.md index b4c868b8..1155e3b7 100644 --- a/docs/superpowers/reviews/2026-09-25-case17-curved-overlap.md +++ b/docs/superpowers/reviews/2026-09-25-case17-curved-overlap.md @@ -102,6 +102,15 @@ exact-mode roots remain checked. The new narrow-strip test has an independent nonlinear parameter correspondence; its interval is only `0.0005` in curve parameters but spans fifty model tolerances. +The full suites also exposed an inappropriate storage expectation in the +cut-contact test: `C(t)=(t,0,(1-t)^2)` only meets the plane at its right +endpoint, but the test required a tolerance-fringe overlap with no isolated +point. That fringe has no domain entrance. The corrected test preserves +the isolated endpoint, and an additional curve with height scaled to +`atol/2` checks a genuine whole-domain CAD overlap retaining its paired +endpoint contact. Both controls validate geometry at the requested model +tolerance instead of imposing `1e-12` source-image equality. + Contact proposals stored inside overlap records share the result allowance with isolated roots. A denied CAD-classification step reports the work stop and retains already found geometry. Completed paired rims also survive diff --git a/tests/test_ssx_cut_overlap_contacts.py b/tests/test_ssx_cut_overlap_contacts.py index 21aaa6be..4fa5b5b6 100644 --- a/tests/test_ssx_cut_overlap_contacts.py +++ b/tests/test_ssx_cut_overlap_contacts.py @@ -10,28 +10,58 @@ @pytest.mark.parametrize("matrix", [np.eye(3), np.array([[0., 1., 0.], [0., 0., 1.], [1., 0., 0.]])]) -def test_tolerance_span_keeps_its_polished_endpoint_contact(matrix): +def test_endpoint_touch_keeps_its_cut_contact(matrix): # C(t)=(t,0,(1-t)^2) meets S(u,v)=(u,v,0) only at t=u=1,v=0. - # Its CAD tolerance neighborhood can be a valid CSX overlap span, - # but that representation must not erase the SSX cut's point seed. + # The near-contact interval starts at a normal-gap threshold, not a + # domain entrance. Following v=0 does not turn it into an overlap. curve = np.array([[0., 0., 1.], [.5, 0., 0.], [1., 0., 0.]]) @ matrix.T surface = np.array([[[0., 0., 0.], [0., 1., 0.]], [[1., 0., 0.], [1., 1., 0.]]]) @ matrix.T result = bez_csx(curve, surface, 1e-3, rational=False, max_results=8) + assert not result['budget_exhausted'] + assert len(result['isolated']) == 1 + assert result['overlaps'] == [] + contacts = list(_cut_face_contacts(result)) + assert len(contacts) == 1 + _assert_endpoint_and_valid_contacts(curve, surface, contacts, matrix, 1e-3) + + +def _assert_endpoint_and_valid_contacts(curve, surface, contacts, matrix, atol): + target = np.array([1., 0., 0.]) @ matrix.T + assert any(np.linalg.norm(point['point']-target) <= atol for point in contacts) + for point in contacts: + parameters = np.array([point[k] for k in ('t', 'u', 'v')]) + assert np.isfinite(parameters).all() + assert np.all((0. <= parameters) & (parameters <= 1.)) + assert np.linalg.norm(eval_curve(curve, point['t'], rational=False) + - point['point']) <= atol + assert np.linalg.norm(eval_surface(surface, point['u'], point['v'], + rational=False)-point['point']) <= atol + + +@pytest.mark.parametrize("matrix", [np.eye(3), + np.array([[0., 1., 0.], [0., 0., 1.], [1., 0., 0.]])]) +def test_whole_cad_overlap_retains_its_endpoint_cut_contact(matrix): + # Every point has normal gap <=atol/2. Both ends are actual curve-domain + # endpoints, so the whole span is a CAD overlap. Absorbing its isolated + # roots must retain the exact right endpoint as a paired SSX seed. + atol = 1e-3 + curve = np.array([[0., 0., .5*atol], [.5, 0., 0.], [1., 0., 0.]]) @ matrix.T + surface = np.array([[[0., 0., 0.], [0., 1., 0.]], + [[1., 0., 0.], [1., 1., 0.]]]) @ matrix.T + result = bez_csx(curve, surface, atol, rational=False) + assert not result['budget_exhausted'] assert result['isolated'] == [] assert len(result['overlaps']) == 1 + overlap = result['overlaps'][0] + for actual, expected in zip(overlap['t_range'], (0., 1.)): + assert np.linalg.norm(eval_curve(curve, actual, rational=False) + -eval_curve(curve, expected, rational=False)) <= atol contacts = list(_cut_face_contacts(result)) - assert len(contacts) == 1 - point = contacts[0] - np.testing.assert_allclose([point[k] for k in ('t', 'u', 'v')], - [1., 1., 0.], atol=1e-6) - np.testing.assert_allclose(eval_curve(curve, point['t'], rational=False), - point['point'], atol=1e-12) - np.testing.assert_allclose( - eval_surface(surface, point['u'], point['v'], rational=False), - point['point'], atol=1e-12) + assert contacts and overlap['boundary_contacts'] + _assert_endpoint_and_valid_contacts(curve, surface, contacts, matrix, atol) @pytest.mark.parametrize("max_results", [0, 1])