brep_kernel/props/mass_properties/

polygons.rs

use super::*;

pub fn trim_polygons(face: &FaceRecord) -> Result<Vec<Vec<[f64; 2]>>, String> {
    let mut polygons = Vec::new();
    for loop_record in &face.loops {
        let mut polygon = Vec::new();
        for coedge in &loop_record.coedges {
            sample_coedge(coedge, [0.0, 0.0], &mut polygon)?;
        }
        polygons.push(polygon);
    }
    Ok(polygons)
}

/// Rebuild a sphere cap bounded by a varying full-period contact rim and the
/// collapsed opposite-winding rim at one parameter pole.  In the flat UV
/// plane both individual loops have zero enclosed area; joined across one cut
/// of the periodic cover they form the actual cap polygon.
pub(super) fn singly_periodic_sphere_cap_polygons(
    face: &FaceRecord,
    closed_u: bool,
    closed_v: bool,
    domain: [f64; 4],
) -> Result<Option<Vec<Vec<[f64; 2]>>>, String> {
    if !matches!(
        face.surface.analytic(),
        Some(crate::AnalyticSurface::Sphere { .. })
    ) || (closed_u, closed_v) != (true, false)
        || face.loops.len() != 2
    {
        return Ok(None);
    }
    let [u0, u1, v0, v1] = domain;
    let period = u1 - u0;
    let v_span = v1 - v0;
    struct Rim {
        points: Vec<[f64; 2]>,
        winding: f64,
        vmin: f64,
        vmax: f64,
    }
    let mut rims = Vec::with_capacity(2);
    for loop_record in &face.loops {
        let mut points = Vec::new();
        for coedge in &loop_record.coedges {
            let [start, end] = coedge.pcurve.domain()?;
            for index in 0..=64 {
                let parameter = start + (end - start) * index as f64 / 64.0;
                let point = coedge.pcurve.evaluate(parameter)?;
                points.push([point.x, point.y]);
            }
        }
        if points.len() < 2 {
            return Ok(None);
        }
        let winding = points.last().unwrap()[0] - points[0][0];
        let (mut vmin, mut vmax) = (f64::INFINITY, f64::NEG_INFINITY);
        for point in &points {
            vmin = vmin.min(point[1]);
            vmax = vmax.max(point[1]);
        }
        rims.push(Rim {
            points,
            winding,
            vmin,
            vmax,
        });
    }
    if rims
        .iter()
        .any(|rim| (rim.winding.abs() - period).abs() > 0.05 * period)
        || rims[0].winding * rims[1].winding >= 0.0
    {
        return Ok(None);
    }
    let flat = |rim: &Rim| rim.vmax - rim.vmin <= 1e-6 * v_span;
    let pole = match (flat(&rims[0]), flat(&rims[1])) {
        (true, false) => &rims[0],
        (false, true) => &rims[1],
        _ => return Ok(None),
    };
    let pole_v = 0.5 * (pole.vmin + pole.vmax);
    if (pole_v - v0).abs() > 1e-6 * v_span && (pole_v - v1).abs() > 1e-6 * v_span {
        return Ok(None);
    }

    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);
    };
    let mut lower_points = lower.points.clone();
    if lower.winding < 0.0 {
        lower_points.reverse();
    }
    let mut upper_points = upper.points.clone();
    if upper.winding > 0.0 {
        upper_points.reverse();
    }
    let shift = ((lower_points[0][0] - upper_points.last().unwrap()[0]) / period).round() * period;
    for point in &mut upper_points {
        point[0] += shift;
    }
    lower_points.extend(upper_points);
    Ok(Some(vec![lower_points]))
}

/// Trim polygons with implicit torus seams unwrapped onto the covering plane.
/// Returns `(polygons, unwrapped)`; `unwrapped` is true iff any loop crossed an
/// implicit seam (and therefore carries out-of-domain parameters). Only attempts
/// the unwrap on a doubly-periodic surface — single-periodic seams the raw
/// integrator already handles are left in-domain.
pub(super) fn trim_polygons_unwrapped(
    face: &FaceRecord,
    closed_u: bool,
    closed_v: bool,
    u_period: f64,
    v_period: f64,
) -> Result<(Vec<Vec<[f64; 2]>>, bool), String> {
    let mut polygons = Vec::new();
    let mut unwrapped = false;
    for loop_record in &face.loops {
        let offsets = crate::topology::loop_seam_offsets(
            &loop_record.coedges,
            closed_u,
            closed_v,
            u_period,
            v_period,
        )?;
        let mut polygon = Vec::new();
        for (index, coedge) in loop_record.coedges.iter().enumerate() {
            let offset = offsets[index];
            if offset[0] != 0.0 || offset[1] != 0.0 {
                unwrapped = true;
            }
            sample_coedge(coedge, offset, &mut polygon)?;
        }
        polygons.push(polygon);
    }
    Ok((polygons, unwrapped))
}

/// A singly-periodic (cylinder/cone) wall whose seam is IMPLICIT: the closed
/// wall is bounded by two full-wrap rim loops (iso-parametric lines at constant
/// cross-level) with NO seam edge connecting them — the two-rim-circles form
/// OCC/STEP emit for a full lateral face. This is the single-periodic analogue
/// of [`biperiodic_band_range`]: each rim is a degenerate iso line in parameter
/// space, so the raw trim integrates the band between them to ~zero (only the
/// in-band holes survive — a full cylinder wall collapses to a ~1000x-too-small
/// sliver, ABC 00000056 solid 2). We restrict to the unambiguous case where the
/// two rims sit at the cross-domain EXTREMES: the material region is then the
/// WHOLE surface (full periodic span × full cross span) minus any hole/notch
/// loops, with no between-rims-vs-complement choice to get wrong. Replacing the
/// two rim lines with a single full-domain band rectangle lets the existing
/// winding integrator fill the wall while the remaining loops subtract as holes.
///
/// Returns the rebuilt polygon set (band rectangle first, then the unwrapped
/// non-rim loops), or None when `face` is not a clean two-extreme-rim band —
/// every other periodic face keeps its existing integration path unchanged.
pub(super) fn singly_periodic_wall_polygons(
    face: &FaceRecord,
    closed_u: bool,
    closed_v: bool,
    domain: [f64; 4],
) -> Result<Option<Vec<Vec<[f64; 2]>>>, String> {
    // Exactly one periodic direction (the doubly-periodic case is handled by
    // `biperiodic_band_integral` upstream).
    if closed_u == closed_v {
        return Ok(None);
    }
    let [u0, u1, v0, v1] = domain;
    let p_is_u = closed_u;
    let period = if p_is_u { u1 - u0 } else { v1 - v0 };
    if !(period > 0.0) {
        return Ok(None);
    }
    let (q_dom_lo, q_dom_hi) = if p_is_u { (v0, v1) } else { (u0, u1) };
    let q_extent = (q_dom_hi - q_dom_lo).abs();
    if !(q_extent > 0.0) {
        return Ok(None);
    }
    // (periodic, cross) coordinate split.
    let coord = |p: [f64; 2]| if p_is_u { (p[0], p[1]) } else { (p[1], p[0]) };

    let mut rim_indices: Vec<usize> = Vec::new();
    let mut rim_levels: Vec<f64> = Vec::new();
    for (index, loop_record) in face.loops.iter().enumerate() {
        // Sample the loop's pcurves in traversal order.
        let mut points: Vec<[f64; 2]> = Vec::new();
        for coedge in &loop_record.coedges {
            let [d0, d1] = coedge.pcurve.domain()?;
            let samples = 16;
            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]);
            }
        }
        if points.len() < 2 {
            return Ok(None);
        }
        let (mut qmin, mut qmax) = (f64::INFINITY, f64::NEG_INFINITY);
        let mut net = 0.0;
        for pair in points.windows(2) {
            let (p0, q0) = coord(pair[0]);
            let (p1, _q1) = coord(pair[1]);
            qmin = qmin.min(q0);
            qmax = qmax.max(q0);
            let mut delta = p1 - p0;
            if delta > 0.5 * period {
                delta -= period;
            } else if delta < -0.5 * period {
                delta += period;
            }
            net += delta;
        }
        if let Some(last) = points.last() {
            let (_p, q) = coord(*last);
            qmin = qmin.min(q);
            qmax = qmax.max(q);
        }
        // A rim: nets a full period in the periodic direction and stays on a
        // near-constant cross-level (a true iso line).
        let is_rim = net.abs() > 0.75 * period && (qmax - qmin) <= 0.02 * q_extent;
        if is_rim {
            rim_indices.push(index);
            rim_levels.push(0.5 * (qmin + qmax));
        }
    }
    if rim_levels.len() != 2 {
        return Ok(None);
    }
    let (q_lo, q_hi) = if rim_levels[0] <= rim_levels[1] {
        (rim_levels[0], rim_levels[1])
    } else {
        (rim_levels[1], rim_levels[0])
    };
    // Require the rims at the cross-domain extremes: material is then the whole
    // wall (no between-rims/complement ambiguity). A sub-band (rims interior to
    // the cross domain) is declined here — the existing path keeps handling it.
    let tol = 0.02 * q_extent;
    if (q_lo - q_dom_lo).abs() > tol || (q_hi - q_dom_hi).abs() > tol {
        return Ok(None);
    }
    if (q_hi - q_lo) <= 0.5 * q_extent {
        return Ok(None);
    }

    // Band rectangle over the full periodic span × the rim cross-levels.
    let (r_ulo, r_uhi, r_vlo, r_vhi) = if p_is_u {
        (u0, u1, q_lo, q_hi)
    } else {
        (q_lo, q_hi, v0, v1)
    };
    let rectangle = vec![
        [r_ulo, r_vlo],
        [r_uhi, r_vlo],
        [r_uhi, r_vhi],
        [r_ulo, r_vhi],
    ];
    let mut polygons: Vec<Vec<[f64; 2]>> = vec![rectangle];
    // Keep the non-rim loops (holes / notches), unwrapping any implicit-seam hop
    // so a seam-straddling hole subtracts as one region.
    for (index, loop_record) in face.loops.iter().enumerate() {
        if rim_indices.contains(&index) {
            continue;
        }
        let offsets = crate::topology::loop_seam_offsets(
            &loop_record.coedges,
            closed_u,
            closed_v,
            u1 - u0,
            v1 - v0,
        )?;
        let mut polygon = Vec::new();
        for (coedge_index, coedge) in loop_record.coedges.iter().enumerate() {
            sample_coedge(coedge, offsets[coedge_index], &mut polygon)?;
        }
        polygons.push(polygon);
    }
    Ok(Some(polygons))
}

/// A DOUBLY-periodic (torus) face whose trim is a full-wrap BAND bounded by a
/// constant-cross-level seam RIM (at a cross-domain extreme) and a single-valued
/// WAVY CUT. The raw winding integrator collapses it to a sliver (the rim is a
/// zero-area iso line and the cut's own polygon encloses only the thin strip it
/// bounds), so — mirroring [`singly_periodic_wall_polygons`] for the singly-
/// periodic wall — rebuild the material region as ONE simple (u,v) polygon:
/// the wavy cut as one boundary, the flat rim iso-line (at the OPPOSITE cross
/// extreme, which is the same physical seam circle) as the other, closed by the
/// two private wrap-seam edges. The winding integrator then fills the band.
///
/// Returns None for anything that is not this exact signature (see
/// [`crate::topology::analyze_doubly_periodic_seam_band`]) — every other face
/// keeps its existing integration path bit-for-bit.
pub(super) fn doubly_periodic_seam_band_polygons(
    face: &FaceRecord,
    domain: [f64; 4],
) -> Result<Option<Vec<Vec<[f64; 2]>>>, String> {
    let (closed_u, closed_v) = face.surface.closed_directions()?;
    if !(closed_u && closed_v) {
        return Ok(None);
    }
    // Sample each loop's pcurves in traversal order.
    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]);
            }
        }
        loops_uv.push(points);
    }
    let Some(band) =
        crate::topology::analyze_doubly_periodic_seam_band(&loops_uv, domain, face.same_sense)
    else {
        return Ok(None);
    };
    let polygon = crate::topology::seam_band_uv_polygon(&loops_uv, domain, &band);
    Ok(Some(vec![polygon]))
}

pub(super) fn is_untrimmed(face: &FaceRecord) -> Result<bool, String> {
    if face.loops.len() != 1 {
        return Ok(false);
    }
    let ku = crate::KnotVector::new(face.surface.knots_u.clone(), face.surface.degree_u)?;
    let kv = crate::KnotVector::new(face.surface.knots_v.clone(), face.surface.degree_v)?;
    let [u0, u1] = ku.domain();
    let [v0, v1] = kv.domain();
    let domain_area = (u1 - u0) * (v1 - v0);
    Ok((parameter_space_area(face)?.abs() - domain_area).abs() <= 1e-6 * domain_area)
}

/// Symmetric degree-5 triangle cubature (barycentric points and weights,
/// normalized to unit area) — Golovanov §8.10 Table 8.10.2.
pub(super) const TRIANGLE_CUBATURE: [([f64; 3], f64); 7] = [
    ([1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0], 0.225),
    (
        [0.101286507323456, 0.101286507323456, 0.797426985353087],
        0.125939180544827,
    ),
    (
        [0.101286507323456, 0.797426985353087, 0.101286507323456],
        0.125939180544827,
    ),
    (
        [0.797426985353087, 0.101286507323456, 0.101286507323456],
        0.125939180544827,
    ),
    (
        [0.470142064105115, 0.470142064105115, 0.059715871789770],
        0.132394152788506,
    ),
    (
        [0.470142064105115, 0.059715871789770, 0.470142064105115],
        0.132394152788506,
    ),
    (
        [0.059715871789770, 0.470142064105115, 0.470142064105115],
        0.132394152788506,
    ),
];

pub(super) fn polygon_signed_area(polygon: &[[f64; 2]]) -> f64 {
    let mut area = 0.0;
    for index in 0..polygon.len() {
        let a = polygon[index];
        let b = polygon[(index + 1) % polygon.len()];
        area += a[0] * b[1] - b[0] * a[1];
    }
    area * 0.5
}

/// Sutherland–Hodgman clip of one loop against an axis-aligned cell.  A
/// loop that fully encloses the cell clips to the whole cell; a disjoint
/// loop clips to nothing — so clipping every loop handles containment,
/// holes, and partial coverage uniformly, with the loop's winding as the
/// contribution sign.
pub(super) fn clip_polygon_to_cell(polygon: &[[f64; 2]], cell: [f64; 4]) -> Vec<[f64; 2]> {
    let [u_low, u_high, v_low, v_high] = cell;
    let mut current = polygon.to_vec();
    // (axis, bound, keep_below)
    for (axis, bound, keep_below) in [
        (0, u_low, false),
        (0, u_high, true),
        (1, v_low, false),
        (1, v_high, true),
    ] {
        if current.len() < 3 {
            return Vec::new();
        }
        let inside = |point: &[f64; 2]| {
            if keep_below {
                point[axis] <= bound
            } else {
                point[axis] >= bound
            }
        };
        let mut next = Vec::with_capacity(current.len() + 4);
        for index in 0..current.len() {
            let a = current[index];
            let b = current[(index + 1) % current.len()];
            let a_in = inside(&a);
            let b_in = inside(&b);
            if a_in {
                next.push(a);
            }
            if a_in != b_in {
                let t = (bound - a[axis]) / (b[axis] - a[axis]);
                let mut crossing = [0.0; 2];
                crossing[axis] = bound;
                crossing[1 - axis] = a[1 - axis] + t * (b[1 - axis] - a[1 - axis]);
                next.push(crossing);
            }
        }
        current = next;
    }
    current
}

pub(super) fn cell_breaks(knots: &[f64], degree: usize, low: f64, high: f64) -> Vec<f64> {
    let mut breaks = vec![low, high];
    breaks.extend(interior_knots(knots, degree));
    breaks.sort_by(f64::total_cmp);
    breaks.dedup_by(|a, b| (*a - *b).abs() <= (high - low) * 1e-12);
    // Subdivide spans so the base grid has at least MINIMUM_CELLS cells
    // across the domain; boundary cells then refine adaptively.
    const MINIMUM_CELLS: usize = 8;
    let target = (high - low) / MINIMUM_CELLS as f64;
    let mut refined = Vec::with_capacity(breaks.len() * 2);
    for pair in breaks.windows(2) {
        refined.push(pair[0]);
        let span = pair[1] - pair[0];
        let pieces = (span / target).ceil().max(1.0) as usize;
        for piece in 1..pieces {
            refined.push(pair[0] + span * piece as f64 / pieces as f64);
        }
    }
    refined.push(high);
    refined
}

/// Tile one period's `cell_breaks` across `[lo, hi]` for a periodic direction
/// whose unwrapped trim ran past the domain end. Shifted copies of the base
/// breaks keep the knot-conforming grid on every covered period; a direction
/// that stayed in-domain yields exactly the base breaks.
pub(super) fn tiled_breaks(base: &[f64], lo: f64, hi: f64, period: f64) -> Vec<f64> {
    if period <= 0.0 || base.len() < 2 {
        return base.to_vec();
    }
    let origin = base[0];
    let first = ((lo - origin) / period).floor() as i64;
    let last = ((hi - origin) / period).ceil() as i64;
    let mut breaks = Vec::new();
    for tile in first..=last {
        let shift = tile as f64 * period;
        for &value in base {
            let shifted = value + shift;
            if breaks.last().map_or(true, |&previous: &f64| {
                (shifted - previous).abs() > 1e-12 * period
            }) {
                breaks.push(shifted);
            }
        }
    }
    breaks
}