brep_kernel/brep/classification/

face_containment.rs

use super::*;

fn point_segment_distance(point: Vec2, start: Vec2, end: Vec2) -> f64 {
    let segment = end.sub(start);
    let length_squared = segment.dot(segment);
    if length_squared <= 1e-30 {
        return point.sub(start).length();
    }
    let parameter = (point.sub(start).dot(segment) / length_squared).clamp(0.0, 1.0);
    point.sub(start.add(segment.scale(parameter))).length()
}

#[derive(Clone, Copy, Debug, PartialEq)]
pub enum PolygonClass {
    Inside,
    Outside,
    Boundary,
}

fn point_in_polygon(point: Vec2, polygon: &[Vec2], tolerance: f64) -> PolygonClass {
    for index in 0..polygon.len() {
        if point_segment_distance(point, polygon[index], polygon[(index + 1) % polygon.len()])
            <= tolerance
        {
            return PolygonClass::Boundary;
        }
    }
    let mut inside = false;
    for index in 0..polygon.len() {
        let a = polygon[index];
        let b = polygon[(index + 1) % polygon.len()];
        if (a.y > point.y) != (b.y > point.y) {
            let crossing = a.x + (point.y - a.y) / (b.y - a.y) * (b.x - a.x);
            if crossing > point.x {
                inside = !inside;
            }
        }
    }
    if inside {
        PolygonClass::Inside
    } else {
        PolygonClass::Outside
    }
}

struct SegmentReference<'a> {
    start: Vec2,
    end: Vec2,
    curve: &'a NurbsCurve,
    parameter_start: f64,
    parameter_end: f64,
}

fn interior_knot_count(curve: &NurbsCurve) -> usize {
    let start = curve.knots[curve.degree];
    let end = curve.knots[curve.knots.len() - 1 - curve.degree];
    let mut previous = None;
    let mut count = 0;
    for &knot in &curve.knots {
        if knot <= start + 1e-12 || knot >= end - 1e-12 {
            continue;
        }
        if previous.is_none_or(|value: f64| (value - knot).abs() > 1e-12) {
            previous = Some(knot);
            count += 1;
        }
    }
    count
}

/// Point-in-face for a DOUBLY-periodic (torus) seam-band face, whose material
/// region the raw even-odd/tangent test below gets wrong: the constant-level
/// seam rim is a zero-area iso line and the wavy cut's own polygon encloses only
/// the thin strip it bounds, so a point in the band reads Outside (or flips on
/// the fragile nearest-tangent refinement). Reconstruct the band as one simple
/// (u,v) polygon — identical topology to the mass integrator and tessellator —
/// and test the point against it, so all three subsystems agree on the material
/// side. Returns None (keep the general path) for every other face.
fn seam_band_point_in_face(
    face: &FaceRecord,
    point: Vec2,
    tolerance: f64,
) -> Result<Option<PolygonClass>, String> {
    if face.surface.closed_directions()? != (true, true) {
        return Ok(None);
    }
    let [u0, u1] = face.surface.domain_u()?;
    let [v0, v1] = face.surface.domain_v()?;
    if crate::topology::doubly_periodic_has_only_collapsed_loops(face)? {
        return Ok(Some(PolygonClass::Inside));
    }
    let u_span = (u1 - u0).abs().max(1e-30);
    let v_span = (v1 - v0).abs().max(1e-30);
    let mut loops_uv: Vec<Vec<[f64; 2]>> = Vec::with_capacity(face.loops.len());
    for loop_record in &face.loops {
        let mut points: Vec<[f64; 2]> = Vec::new();
        for coedge in &loop_record.coedges {
            let [d0, d1] = coedge.pcurve.domain()?;
            let samples = 24;
            for k in 0..=samples {
                let t = d0 + (d1 - d0) * k as f64 / samples as f64;
                let p = coedge.pcurve.evaluate(t)?;
                points.push([p.x, p.y]);
            }
        }
        let (mut umin, mut umax, mut vmin, mut vmax) = (
            f64::INFINITY,
            f64::NEG_INFINITY,
            f64::INFINITY,
            f64::NEG_INFINITY,
        );
        for p in &points {
            umin = umin.min(p[0]);
            umax = umax.max(p[0]);
            vmin = vmin.min(p[1]);
            vmax = vmax.max(p[1]);
        }
        // Mirror `mass_properties::biperiodic_band_range`: a torus has no
        // geometric pole, so a VERTEX_LOOP whose whole pcurve trace collapses
        // to one parameter point is a zero-area puncture, not a band boundary.
        // Restrict the new behavior to faces that actually contain such an
        // extra loop; ordinary two-loop torus classification stays unchanged.
        if face.loops.len() > 2 && (umax - umin) <= 1e-3 * u_span && (vmax - vmin) <= 1e-3 * v_span
        {
            continue;
        }
        loops_uv.push(points);
    }
    // MERGED SEAM-CARRYING BAND (single loop): `insert_periodic_band_seam_edges`
    // fuses a wrapped band's two rims into ONE loop joined by seam columns, so
    // the two-loop paths below never see it, and the plain even-odd polygon
    // carries a period-magnitude hop at each rim→column junction (a phantom
    // diagonal across the domain) that misclassifies the whole middle zone.
    // Unwrap the loop onto the covering plane (`loop_seam_offsets`, the same
    // fold mass integration and tessellation use), where it IS a simple
    // polygon, and even-odd the query's period images against it. Gated on a
    // genuinely seam-crossing single loop on a doubly-periodic surface; every
    // other face falls through unchanged. Hatch: BREP_SEAM_BAND_MERGED=0.
    if face.loops.len() == 1
        && loops_uv.len() == 1
        && std::env::var("BREP_SEAM_BAND_MERGED").as_deref() != Ok("0")
    {
        let coedges = &face.loops[0].coedges;
        let offsets = crate::topology::loop_seam_offsets(coedges, true, true, u_span, v_span)?;
        if offsets.iter().any(|o| o[0] != 0.0 || o[1] != 0.0) {
            let mut polygon: Vec<Vec2> = Vec::new();
            for (coedge_index, coedge) in coedges.iter().enumerate() {
                let [d0, d1] = coedge.pcurve.domain()?;
                let samples = 24;
                for k in 0..samples {
                    let t = d0 + (d1 - d0) * k as f64 / samples as f64;
                    let p = coedge.pcurve.evaluate(t)?;
                    polygon.push(Vec2 {
                        x: p.x + offsets[coedge_index][0],
                        y: p.y + offsets[coedge_index][1],
                    });
                }
            }
            let mut best = PolygonClass::Outside;
            'images: for du in [-1.0, 0.0, 1.0] {
                for dv in [-1.0, 0.0, 1.0] {
                    let image = Vec2 {
                        x: point.x + du * u_span,
                        y: point.y + dv * v_span,
                    };
                    match point_in_polygon(image, &polygon, tolerance) {
                        PolygonClass::Boundary => {
                            best = PolygonClass::Boundary;
                            break 'images;
                        }
                        PolygonClass::Inside => best = PolygonClass::Inside,
                        PolygonClass::Outside => {}
                    }
                }
            }
            return Ok(Some(best));
        }
    }
    if loops_uv.len() != 2 {
        return Ok(None);
    }
    if let Some(band) = crate::topology::analyze_doubly_periodic_seam_band(
        &loops_uv,
        [u0, u1, v0, v1],
        face.same_sense,
    ) {
        let polygon: Vec<Vec2> =
            crate::topology::seam_band_uv_polygon(&loops_uv, [u0, u1, v0, v1], &band)
                .into_iter()
                .map(|p| Vec2 { x: p[0], y: p[1] })
                .collect();
        return Ok(Some(point_in_polygon(point, &polygon, tolerance)));
    }
    // TWO clean full-wrap rims are the companion `biperiodic_band_range` case,
    // with or without an extra collapsed puncture. Classify only the cross
    // parameter: the material is either the strip between the rims or its
    // periodic complement. The exact two-loop/full-wrap/constant-cross-level
    // gates below mirror mass integration and tessellation; every other torus
    // face retains the general classifier.
    for p_is_u in [true, false] {
        let (period, q_extent) = if p_is_u {
            (u1 - u0, v_span)
        } else {
            (v1 - v0, u_span)
        };
        if !(period > 0.0) {
            continue;
        }
        let coord = |p: &[f64; 2]| if p_is_u { (p[0], p[1]) } else { (p[1], p[0]) };
        let mut rings = Vec::with_capacity(2);
        for points in &loops_uv {
            let (mut pmin, mut pmax, mut qmin, mut qmax) = (
                f64::INFINITY,
                f64::NEG_INFINITY,
                f64::INFINITY,
                f64::NEG_INFINITY,
            );
            for p in points {
                let (periodic, cross) = coord(p);
                pmin = pmin.min(periodic);
                pmax = pmax.max(periodic);
                qmin = qmin.min(cross);
                qmax = qmax.max(cross);
            }
            if (pmax - pmin) < 0.6 * period || (qmax - qmin) > 0.05 * q_extent {
                rings.clear();
                break;
            }
            let mut net = 0.0;
            for pair in points.windows(2) {
                let mut delta = coord(&pair[1]).0 - coord(&pair[0]).0;
                if delta > 0.5 * period {
                    delta -= period;
                } else if delta < -0.5 * period {
                    delta += period;
                }
                net += delta;
            }
            let direction = if net > 0.25 * period {
                1
            } else if net < -0.25 * period {
                -1
            } else {
                0
            };
            rings.push((0.5 * (qmin + qmax), direction));
        }
        if rings.len() != 2 {
            continue;
        }
        rings.sort_by(|a, b| a.0.total_cmp(&b.0));
        let [(q_lo, lower_direction), (q_hi, upper_direction)] = rings.as_slice() else {
            unreachable!()
        };
        let inconclusive =
            *lower_direction == 0 || *upper_direction == 0 || lower_direction == upper_direction;
        let between_is_ccw_uv = if p_is_u {
            *lower_direction > 0
        } else {
            *lower_direction < 0
        };
        let complement = !inconclusive && between_is_ccw_uv != face.same_sense;
        let q = if p_is_u { point.y } else { point.x };
        if (q - q_lo).abs() <= tolerance || (q - q_hi).abs() <= tolerance {
            return Ok(Some(PolygonClass::Boundary));
        }
        let between = q > *q_lo && q < *q_hi;
        return Ok(Some(if between != complement {
            PolygonClass::Inside
        } else {
            PolygonClass::Outside
        }));
    }
    Ok(None)
}

/// Point-in-face for a SINGLY-PERIODIC face with a SEAM-STRADDLING loop: a
/// trim loop crossing the u seam is stored with its samples folded into the
/// domain, so its flat polygon is TORN at the seam and plain even-odd
/// misclassifies near the tear (STEP 00000585-family: face 519's big inner
/// loop straddles the seam and a REAL section sub-segment between two accepted
/// pieces read Outside in `process_curve`'s midpoint gate — trial 174's
/// missing middle segment). Mirror fragment.rs's proven seam-straddling
/// fallback: UNWRAP every loop by accumulating period-folded deltas (a
/// zero-winding loop closes in the covering plane), then even-odd the query at
/// its −P/0/+P period images across all unwrapped loops. Fires only when the
/// face is closed in exactly u, has >2 loops, and at least one loop straddles
/// the seam (≥1 fold jump); every other face returns None and keeps the plain
/// classifier bit-identically. Loops that WIND the period (net winding ≠ 0)
/// do not close under unwrapping — bail to the plain path rather than guess.
/// Escape hatch `BREP_HORIZON_CONTAINMENT=0`.
fn wrapped_horizon_point_in_face(
    face: &FaceRecord,
    point: Vec2,
    tolerance: f64,
) -> Result<Option<PolygonClass>, String> {
    if face.surface.closed_directions()? != (true, false) || face.loops.is_empty() {
        return Ok(None);
    }
    if std::env::var("BREP_HORIZON_CONTAINMENT").as_deref() == Ok("0") {
        return Ok(None);
    }
    // SINGLE/DOUBLE-LOOP straddlers (t91's face 675: ONE loop drawn in-domain
    // that hops the u-seam twice, trimming a sliver strip AT the seam) used
    // to be excluded by a `loops.len() <= 2` gate and fell to the plain
    // even-odd, whose verdict on the hopped polygon is INVERTED (the strip
    // interior read Outside, far azimuths read Inside) — minting bogus
    // far-azimuth section pieces and discarding the real ones. The unwrap +
    // period-image even-odd below is loop-count-agnostic, so the gate now
    // only requires a non-empty loop set; non-straddling loops still bail to
    // the plain path unchanged. Escape hatch: BREP_HORIZON_SINGLE_LOOP=0
    // restores the old minimum.
    if face.loops.len() <= 2
        && std::env::var("BREP_HORIZON_SINGLE_LOOP").as_deref() == Ok("0")
    {
        return Ok(None);
    }
    let [u0, u1] = face.surface.domain_u()?;
    let u_period = u1 - u0;
    if u_period <= 0.0 {
        return Ok(None);
    }
    let debug = std::env::var("BREP_DEBUG_HORIZON").is_ok();
    let mut loops: Vec<Vec<Vec2>> = Vec::new();
    let mut straddling = false;
    for loop_record in &face.loops {
        let mut points: Vec<Vec2> = Vec::new();
        for coedge in &loop_record.coedges {
            let curve = &coedge.pcurve;
            let [start, end] = curve.domain()?;
            let sample_count =
                2usize.max((interior_knot_count(curve) + 1) * (curve.degree + 1) * 4);
            for index in 0..sample_count {
                let parameter = start + (end - start) * index as f64 / sample_count as f64;
                let evaluated = curve.evaluate(parameter)?;
                points.push(Vec2 {
                    x: evaluated.x,
                    y: evaluated.y,
                });
            }
        }
        if points.len() < 3 {
            continue;
        }
        // Unwrap into the covering plane: place each sample at the previous
        // one plus the period-folded delta. A seam-straddling loop becomes a
        // continuous closed polygon; a WINDING loop does not close — bail.
        let mut unwrapped = Vec::with_capacity(points.len());
        let mut jumps = 0usize;
        let mut cursor = points[0];
        unwrapped.push(cursor);
        for pair in points.windows(2) {
            let mut du = pair[1].x - pair[0].x;
            let folded = du - u_period * (du / u_period).round();
            if (du - folded).abs() > 0.25 * u_period {
                jumps += 1;
            }
            du = folded;
            cursor = Vec2 {
                x: cursor.x + du,
                y: pair[1].y,
            };
            unwrapped.push(cursor);
        }
        let closure = (unwrapped[0].x - unwrapped[unwrapped.len() - 1].x).abs();
        if closure > 0.25 * u_period {
            if debug {
                eprintln!(
                    "horizon: face {} loop winds the period (closure {closure:.3}) — bail",
                    face.id
                );
            }
            return Ok(None); // winding loop: unwrapping cannot close it
        }
        if jumps > 0 {
            straddling = true;
        }
        loops.push(unwrapped);
    }
    if !straddling || loops.is_empty() {
        return Ok(None);
    }
    // Even-odd on the QUOTIENT cylinder: each unwrapped loop is tested at
    // every period image of the query and contributes its own parity. The old
    // combining ("any single image with an odd crossing count ⇒ Inside") is
    // wrong whenever two loops live in DIFFERENT period frames after
    // unwrapping — t181/00000312#p7: face 519's seam-straddling HOLE unwraps
    // to u∈[−0.039, 0.039] while the outer rectangle spans [0, 1], so a query
    // inside the hole at u≈0.986 is contained by the outer loop at shift 0
    // (odd ⇒ Inside) and by the hole at shift −P (odd ⇒ Inside again); the
    // true even count (outer + hole = 2 ⇒ Outside) is never assembled at any
    // single image. Fragment selection then kept a PHANTOM half-hole region
    // whose boundary (half the hole loop + a derived chord) minted same-sense
    // coedges at assembly. Summing per-loop parities across images restores
    // the covering-space even-odd. Escape hatch: BREP_HORIZON_CROSS_FRAME=0
    // restores the old per-image combining.
    if std::env::var("BREP_HORIZON_CROSS_FRAME").as_deref() == Ok("0") {
        let mut best: Option<PolygonClass> = None;
        for shift in [-u_period, 0.0, u_period] {
            let image = Vec2 {
                x: point.x + shift,
                y: point.y,
            };
            let mut crossings = 0usize;
            for polygon in &loops {
                match point_in_polygon(image, polygon, tolerance) {
                    PolygonClass::Boundary => return Ok(Some(PolygonClass::Boundary)),
                    PolygonClass::Inside => crossings += 1,
                    PolygonClass::Outside => {}
                }
            }
            if crossings % 2 == 1 {
                best = Some(PolygonClass::Inside);
            } else if best.is_none() {
                best = Some(PolygonClass::Outside);
            }
        }
        if debug {
            eprintln!(
                "horizon: face {} loops={} straddling (legacy per-image) -> {:?}",
                face.id,
                loops.len(),
                best
            );
        }
        return Ok(best);
    }
    let mut crossings = 0usize;
    for polygon in &loops {
        let mut image_hits = 0usize;
        for shift in [-u_period, 0.0, u_period] {
            let image = Vec2 {
                x: point.x + shift,
                y: point.y,
            };
            match point_in_polygon(image, polygon, tolerance) {
                PolygonClass::Boundary => return Ok(Some(PolygonClass::Boundary)),
                PolygonClass::Inside => image_hits += 1,
                PolygonClass::Outside => {}
            }
        }
        crossings += image_hits % 2;
    }
    let class = if crossings % 2 == 1 {
        PolygonClass::Inside
    } else {
        PolygonClass::Outside
    };
    if debug {
        eprintln!(
            "horizon: face {} loops={} straddling -> {:?}",
            face.id,
            loops.len(),
            class
        );
    }
    Ok(Some(class))
}


/// Point-in-face for a SINGLY-PERIODIC (closed-u) FULL-BAND STRIP whose two
/// full-wrap rim loops are drawn as COVERING-PLANE pcurves running OUT OF the
/// u-domain by up to a whole period (case 08:10-genus face 1513, the 4th seam
/// representation after helmet-unwrapped / band-complement / 675-hops): the
/// raw even-odd sees u-winding loop polygons whose parity is meaningless, so
/// nearly the whole material strip reads Outside, the marched clip crumbles
/// the section curves to sub-mm bits, and the neighbours strand one-use.
///
/// The material region is the strip BETWEEN the two rims: on an open-v
/// surface a 2-loop face whose loops both wrap the full period has no other
/// representable region. The verdict is a sampled per-u hull — each rim is a
/// single periodic polyline v = rim(u) in the covering plane, and the query
/// point is between the rims iff an upward (+v) vertical ray at the query's
/// u (folded modulo the period into each rim's own one-period window) crosses
/// the two rim polylines an odd number of times in total. NO flat-level
/// approximation is used (the rims are scalloped: 1513's upper rim is flat at
/// v=1 on parts of u but dips to v=0.745 elsewhere — a v-between-extremes
/// rule is provably wrong on it).
///
/// STRUCTURAL gates (all must hold; anything else returns None and keeps the
/// previous classification bit-identically):
///   - surface closed in exactly u; exactly 2 loops;
///   - at least one loop's samples run OUT of the u-domain by >1e-3 period
///     (the covering-plane discriminator — in-domain 2-rim bands keep their
///     current path);
///   - each loop's unwrapped net u-travel is exactly ±one period and its v
///     closes (a true full-wrap rim, not a partial arc);
///   - opposite travel directions, disjoint v-hulls, and the material-left
///     winding rule confirming the BETWEEN strip (a complement verdict cannot
///     be represented on an open-v surface — decline, never guess).
/// Escape hatch: BREP_COVERING_RIM_STRIP=0.
fn covering_rim_strip_point_in_face(
    face: &FaceRecord,
    point: Vec2,
    tolerance: f64,
) -> Result<Option<PolygonClass>, String> {
    if face.surface.closed_directions()? != (true, false) || face.loops.len() != 2 {
        return Ok(None);
    }
    if std::env::var("BREP_COVERING_RIM_STRIP").as_deref() == Ok("0") {
        return Ok(None);
    }
    let [u0, u1] = face.surface.domain_u()?;
    let [v0, v1] = face.surface.domain_v()?;
    let period = u1 - u0;
    if !(period > 0.0) {
        return Ok(None);
    }
    let v_span = (v1 - v0).abs().max(1e-30);
    struct Rim {
        points: Vec<Vec2>,
        vmin: f64,
        vmax: f64,
        ascending: bool,
    }
    let mut rims: Vec<Rim> = Vec::with_capacity(2);
    let mut out_of_domain = false;
    for loop_record in &face.loops {
        let mut points: Vec<Vec2> = Vec::new();
        for coedge in &loop_record.coedges {
            let curve = &coedge.pcurve;
            let [start, end] = curve.domain()?;
            let sample_count =
                2usize.max((interior_knot_count(curve) + 1) * (curve.degree + 1) * 4);
            for index in 0..=sample_count {
                let parameter = start + (end - start) * index as f64 / sample_count as f64;
                let evaluated = curve.evaluate(parameter)?;
                points.push(Vec2 {
                    x: evaluated.x,
                    y: evaluated.y,
                });
            }
        }
        if points.len() < 3 {
            return Ok(None);
        }
        // Covering-plane discriminator: this lane exists for loops drawn PAST
        // the domain; everything in-domain keeps the existing classifiers.
        for p in &points {
            if p.x < u0 - 1e-3 * period || p.x > u1 + 1e-3 * period {
                out_of_domain = true;
            }
        }
        // Unwrap into the covering plane (period-folded deltas — a no-op for
        // an already-continuous covering-plane representation, and the same
        // fold the horizon lane uses for in-domain seam hops).
        let mut unwrapped: Vec<Vec2> = Vec::with_capacity(points.len());
        let mut cursor = points[0];
        unwrapped.push(cursor);
        for pair in points.windows(2) {
            let du = pair[1].x - pair[0].x;
            let folded = du - period * (du / period).round();
            cursor = Vec2 {
                x: cursor.x + folded,
                y: pair[1].y,
            };
            unwrapped.push(cursor);
        }
        let net = unwrapped[unwrapped.len() - 1].x - unwrapped[0].x;
        if (net.abs() - period).abs() > 2e-2 * period {
            return Ok(None); // not a clean single full wrap
        }
        if (unwrapped[unwrapped.len() - 1].y - unwrapped[0].y).abs() > 1e-3 * v_span {
            return Ok(None); // rim does not close in v
        }
        let (mut vmin, mut vmax) = (f64::INFINITY, f64::NEG_INFINITY);
        for p in &unwrapped {
            vmin = vmin.min(p.y);
            vmax = vmax.max(p.y);
        }
        rims.push(Rim {
            points: unwrapped,
            vmin,
            vmax,
            ascending: net > 0.0,
        });
    }
    if !out_of_domain {
        return Ok(None);
    }
    if rims[0].ascending == rims[1].ascending {
        return Ok(None); // rims must traverse opposite u directions
    }
    let (lower, upper) = if rims[0].vmax <= rims[1].vmin {
        (&rims[0], &rims[1])
    } else if rims[1].vmax <= rims[0].vmin {
        (&rims[1], &rims[0])
    } else {
        return Ok(None); // interleaved v-hulls: not a clean strip
    };
    if upper.vmin - lower.vmax <= 1e-6 * v_span {
        return Ok(None);
    }
    // Material-left rule (same convention as the doubly-periodic band lanes):
    // the between strip is CCW in uv iff the LOWER rim travels +u; that must
    // match the face sense, else the material would be the complement — which
    // an open-v surface cannot represent. Decline rather than guess.
    if (lower.ascending) != face.same_sense {
        return Ok(None);
    }
    // Boundary: proximity to either rim polyline at the query's period images.
    for rim in [lower, upper] {
        for shift in [-period, 0.0, period] {
            let image = Vec2 {
                x: point.x + shift,
                y: point.y,
            };
            for pair in rim.points.windows(2) {
                if point_segment_distance(image, pair[0], pair[1]) <= tolerance {
                    return Ok(Some(PolygonClass::Boundary));
                }
            }
        }
    }
    // Sampled per-u hull: fold the query into each rim's own one-period
    // window and count upward-ray crossings; odd total = between the rims.
    let mut crossings = 0usize;
    for rim in [lower, upper] {
        let window_base = rim.points[0].x.min(rim.points[rim.points.len() - 1].x);
        let x = window_base + (point.x - window_base).rem_euclid(period);
        for pair in rim.points.windows(2) {
            let (a, b) = (pair[0], pair[1]);
            if (a.x > x) != (b.x > x) {
                let v_cross = a.y + (x - a.x) / (b.x - a.x) * (b.y - a.y);
                if v_cross > point.y {
                    crossings += 1;
                }
            }
        }
    }
    Ok(Some(if crossings % 2 == 1 {
        PolygonClass::Inside
    } else {
        PolygonClass::Outside
    }))
}

pub fn parameter_point_in_face(
    face: &FaceRecord,
    point: Vec2,
    tolerance: f64,
) -> Result<PolygonClass, String> {
    if let Some(class) = seam_band_point_in_face(face, point, tolerance)? {
        return Ok(class);
    }
    if let Some(class) = wrapped_horizon_point_in_face(face, point, tolerance)? {
        return Ok(class);
    }
    if let Some(class) = covering_rim_strip_point_in_face(face, point, tolerance)? {
        return Ok(class);
    }
    let mut crossings = 0;
    let mut boundary = false;
    let mut nearest: Option<(SegmentReference<'_>, f64)> = None;
    for loop_record in &face.loops {
        let mut polygon = Vec::new();
        let mut segments = Vec::new();
        for coedge in &loop_record.coedges {
            let curve = &coedge.pcurve;
            let [start, end] = curve.domain()?;
            let sample_count =
                2usize.max((interior_knot_count(curve) + 1) * (curve.degree + 1) * 4);
            for index in 0..sample_count {
                let parameter = start + (end - start) * index as f64 / sample_count as f64;
                let evaluated = curve.evaluate(parameter)?;
                polygon.push(Vec2 {
                    x: evaluated.x,
                    y: evaluated.y,
                });
                let parameter_end = if index + 1 < sample_count {
                    start + (end - start) * (index + 1) as f64 / sample_count as f64
                } else {
                    end
                };
                segments.push(SegmentReference {
                    start: Vec2 {
                        x: evaluated.x,
                        y: evaluated.y,
                    },
                    end: Vec2 { x: 0.0, y: 0.0 },
                    curve,
                    parameter_start: parameter,
                    parameter_end,
                });
            }
        }
        for index in 0..polygon.len() {
            segments[index].end = polygon[(index + 1) % polygon.len()];
        }
        match point_in_polygon(point, &polygon, tolerance) {
            PolygonClass::Boundary => boundary = true,
            PolygonClass::Inside => crossings += 1,
            PolygonClass::Outside => {}
        }
        for segment in segments {
            let distance = point_segment_distance(point, segment.start, segment.end);
            if nearest
                .as_ref()
                .is_none_or(|(_, nearest_distance)| distance < *nearest_distance)
            {
                nearest = Some((segment, distance));
            }
        }
    }
    if boundary {
        return Ok(PolygonClass::Boundary);
    }
    let parity = if crossings % 2 == 1 {
        PolygonClass::Inside
    } else {
        PolygonClass::Outside
    };
    let Some((segment, distance)) = nearest else {
        return Ok(parity);
    };
    let segment_length = segment.end.sub(segment.start).length();
    if distance > segment_length || segment_length <= 0.0 {
        return Ok(parity);
    }
    let chord_parameter = segment.parameter_start
        + (segment.parameter_end - segment.parameter_start)
            * (point.sub(segment.start).dot(segment.end.sub(segment.start))
                / (segment_length * segment_length))
                .clamp(0.0, 1.0);
    let [domain_start, domain_end] = segment.curve.domain()?;
    let mut parameter = chord_parameter;
    for _ in 0..12 {
        let derivatives = segment.curve.derivatives(parameter, 2)?;
        let on_curve = Vec2 {
            x: derivatives[0].x,
            y: derivatives[0].y,
        };
        let tangent = Vec2 {
            x: derivatives[1].x,
            y: derivatives[1].y,
        };
        let second = Vec2 {
            x: derivatives[2].x,
            y: derivatives[2].y,
        };
        let residual = on_curve.sub(point);
        let denominator = tangent.dot(tangent) + residual.dot(second);
        if denominator.abs() < 1e-30 {
            break;
        }
        let step = -residual.dot(tangent) / denominator;
        parameter = (parameter + step).clamp(domain_start, domain_end);
        if step.abs() < 1e-14 * (domain_end - domain_start + 1.0) {
            break;
        }
    }
    let margin = 1e-9 * (domain_end - domain_start);
    if parameter > domain_start + margin && parameter < domain_end - margin {
        let derivatives = segment.curve.derivatives(parameter, 1)?;
        let on_curve = Vec2 {
            x: derivatives[0].x,
            y: derivatives[0].y,
        };
        let tangent = Vec2 {
            x: derivatives[1].x,
            y: derivatives[1].y,
        };
        let offset = point.sub(on_curve);
        if offset.length() <= tolerance {
            return Ok(PolygonClass::Boundary);
        }
        let cross = tangent.x * offset.y - tangent.y * offset.x;
        if cross.abs() > 1e-30 {
            return Ok(if (cross > 0.0) == face.same_sense {
                PolygonClass::Inside
            } else {
                PolygonClass::Outside
            });
        }
    }
    Ok(parity)
}