brepkit-operations 4.6.1

//! Sectioning (slicing) solids with planes.
//!
//! Computes the cross-section of a solid at a given cutting plane,
//! producing face(s) representing the intersection.

#![allow(clippy::too_many_lines, clippy::doc_markdown)]

use brepkit_math::tolerance::Tolerance;
use brepkit_math::vec::{Point3, Vec3};
use brepkit_topology::Topology;
use brepkit_topology::edge::{Edge, EdgeCurve};
use brepkit_topology::face::{Face, FaceId, FaceSurface};
use brepkit_topology::solid::SolidId;
use brepkit_topology::vertex::Vertex;
use brepkit_topology::wire::{OrientedEdge, Wire, WireId};

use brepkit_math::nurbs::intersection::{
    IntersectionCurve, IntersectionPoint, intersect_plane_nurbs,
};
use brepkit_math::nurbs::surface::NurbsSurface;

use crate::boolean::face_polygon;
use crate::dot_normal_point;

/// How far, in spans of the domain, `mid` lies from the nearest of the
/// section's `curves` other than the one at `index`, measured round a closed
/// direction's seam the shorter way. A contact along the seam itself comes
/// back as two curves, one at each of its sides: the curve at the other
/// side is the same one, not a neighbour.
fn gap_to_others(
    curves: &[IntersectionCurve],
    index: usize,
    mid: (f64, f64),
    ((u_lo, u_hi), (v_lo, v_hi)): ((f64, f64), (f64, f64)),
    (u_closed, v_closed): (bool, bool),
) -> f64 {
    let (wide, tall) = (u_hi - u_lo, v_hi - v_lo);
    let short = |x: f64, closed: bool| if closed { x - x.round() } else { x };
    let from = |q: &IntersectionPoint| {
        (
            short((q.param1.0 - mid.0) / wide, u_closed),
            short((q.param1.1 - mid.1) / tall, v_closed),
        )
    };
    // Which side of a closed direction's seam all of a curve lies on.
    let side =
        |c: &[IntersectionPoint], (lo, hi): (f64, f64), at: fn(&IntersectionPoint) -> f64| {
            let on = |edge: f64| c.iter().all(|q| (at(q) - edge).abs() <= 1e-12 * (hi - lo));
            if on(lo) {
                Some(false)
            } else if on(hi) {
                Some(true)
            } else {
                None
            }
        };
    let ours = &curves[index].points;
    let copy = |c: &[IntersectionPoint]| {
        let across = |closed: bool, range: (f64, f64), at: fn(&IntersectionPoint) -> f64| {
            closed
                && matches!((side(ours, range, at), side(c, range, at)), (Some(a), Some(b)) if a != b)
        };
        across(u_closed, (u_lo, u_hi), |q| q.param1.0)
            || across(v_closed, (v_lo, v_hi), |q| q.param1.1)
    };
    curves
        .iter()
        .enumerate()
        .filter(|&(j, c)| j != index && !copy(&c.points))
        .map(|(_, c)| {
            let single = (c.points.len() == 1).then(|| {
                let (x, y) = from(&c.points[0]);
                x.hypot(y)
            });
            let spans = c.points.windows(2).map(|pair| {
                // The segment's far end on the near end's side of a seam.
                let (ax, ay) = from(&pair[0]);
                let (bx, by) = (
                    ax + short((pair[1].param1.0 - pair[0].param1.0) / wide, u_closed),
                    ay + short((pair[1].param1.1 - pair[0].param1.1) / tall, v_closed),
                );
                let (dx, dy) = (bx - ax, by - ay);
                let length = dx.mul_add(dx, dy * dy);
                let t = if length > 0.0 {
                    (-ax.mul_add(dx, ay * dy) / length).clamp(0.0, 1.0)
                } else {
                    0.0
                };
                t.mul_add(dx, ax).hypot(t.mul_add(dy, ay))
            });
            single
                .into_iter()
                .chain(spans)
                .fold(f64::INFINITY, f64::min)
        })
        .fold(f64::INFINITY, f64::min)
}

/// Whether the plane only touches `nurbs` along the section through
/// `points`, without crossing it.
///
/// The surface's normal lies along the plane's at every point, and across
/// the curve at its middle point the surface stays on one side of the
/// plane. Where the surface crosses with its normal along the plane's, as
/// `z = x³` crosses `z = 0`, the curve parts the regions either side. Each
/// side's probe starts a billionth of the domain across and widens tenfold
/// while the surface stays within rounding of the plane, as across
/// `z = x⁷`, but stops short of halfway to the section's nearest other
/// curve, so that it tells the side before another branch comes between;
/// a curve whose sides no probe tells stays. On a closed surface the probes
/// wrap round its seam.
fn touches_along(
    nurbs: &NurbsSurface,
    normal: Vec3,
    d: f64,
    curves: &[IntersectionCurve],
    index: usize,
) -> bool {
    let points = &curves[index].points;
    let tangent = |p: &IntersectionPoint| {
        nurbs
            .normal(p.param1.0, p.param1.1)
            .is_ok_and(|n| n.cross(normal).length() <= 1e-9 * normal.length())
    };
    if points.len() < 2 || !points.iter().all(tangent) {
        return false;
    }
    let (u_lo, u_hi) = nurbs.domain_u();
    let (v_lo, v_hi) = nurbs.domain_v();
    let (wide, tall) = (u_hi - u_lo, v_hi - v_lo);
    let k = points.len() / 2;
    let (a, b) = (
        points[k - 1].param1,
        points[(k + 1).min(points.len() - 1)].param1,
    );
    // Across the curve in parameters scaled to the domain.
    let (tu, tv) = ((b.0 - a.0) / wide, (b.1 - a.1) / tall);
    let len = tu.hypot(tv);
    if len == 0.0 {
        return false;
    }
    let mid = points[k].param1;
    let rounding = 4.0
        * f64::EPSILON
        * nurbs
            .control_points()
            .iter()
            .flatten()
            .map(|p| {
                (normal.x() * p.x()).abs()
                    + (normal.y() * p.y()).abs()
                    + (normal.z() * p.z()).abs()
                    + d.abs()
            })
            .fold(0.0, f64::max);
    // Round a closed direction's seam, onto the domain's far side.
    let wrap = |x: f64, lo: f64, hi: f64, closed: bool| {
        if closed {
            (x - lo).rem_euclid(hi - lo) + lo
        } else {
            x.clamp(lo, hi)
        }
    };
    let (u_closed, v_closed) = (nurbs.is_periodic_u(), nurbs.is_periodic_v());
    let side = |reach: f64| {
        let u = wrap(
            (reach * wide).mul_add(-tv / len, mid.0),
            u_lo,
            u_hi,
            u_closed,
        );
        let v = wrap(
            (reach * tall).mul_add(tu / len, mid.1),
            v_lo,
            v_hi,
            v_closed,
        );
        let p = nurbs.evaluate(u, v);
        normal.dot(Vec3::new(p.x(), p.y(), p.z())) - d
    };
    // No probe reaches halfway to another curve of the section, so that
    // none passes beyond it.
    let gap = gap_to_others(
        curves,
        index,
        mid,
        ((u_lo, u_hi), (v_lo, v_hi)),
        (u_closed, v_closed),
    );
    let nearest = |sign: f64| {
        [1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1]
            .into_iter()
            .take_while(|&reach| reach < 0.5 * gap)
            .map(|reach| side(sign * reach))
            .find(|f| f.abs() > rounding)
    };
    match (nearest(1.0), nearest(-1.0)) {
        (Some(left), Some(right)) => left * right > 0.0,
        _ => false,
    }
}

/// Chain intersection curve points into consecutive segment pairs.
///
/// Instead of connecting only the first and last point (chord approximation),
/// this chains all intermediate points as individual segments, faithfully
/// tracing the actual intersection curve.
fn chain_curve_points(points: &[IntersectionPoint], segments: &mut Vec<(Point3, Point3)>) {
    if points.len() < 2 {
        return;
    }
    for pair in points.windows(2) {
        segments.push((pair[0].point, pair[1].point));
    }
}

/// A cross-section result: one or more planar faces on the cutting plane.
#[derive(Debug)]
pub struct Section {
    /// The face IDs of the cross-section faces in the topology.
    pub faces: Vec<FaceId>,
}

/// Compute the cross-section of a solid with a plane.
///
/// The cutting plane is defined by a point on the plane and its normal.
/// Returns the cross-section as one or more planar faces lying on the
/// cutting plane. For a simple convex solid this is typically one face;
/// for solids with holes or disconnected volumes there may be multiple.
/// A plane that misses the solid entirely yields an empty face list
/// (a successful empty result, not an error).
///
/// # Algorithm
///
/// 1. For each face of the solid, compute the intersection with the
///    cutting plane (line segments for planar faces).
/// 2. Collect all intersection segments.
/// 3. Assemble segments into closed wires.
/// 4. Create planar faces from the wires.
///
/// # Errors
///
/// Returns an error if NURBS intersection computation fails, or if
/// intersection segments exist but cannot be assembled into a closed
/// cross-section wire.
pub fn section(
    topo: &mut Topology,
    solid: SolidId,
    plane_point: Point3,
    plane_normal: Vec3,
) -> Result<Section, crate::OperationsError> {
    let tol = Tolerance::new();

    let normal = plane_normal.normalize()?;
    let d = dot_normal_point(normal, plane_point);

    let solid_data = topo.solid(solid)?;
    let all_shell_ids: Vec<_> = std::iter::once(solid_data.outer_shell())
        .chain(solid_data.inner_shells().iter().copied())
        .collect();

    let mut face_ids: Vec<FaceId> = Vec::new();
    for shell_id in &all_shell_ids {
        let shell = topo.shell(*shell_id)?;
        face_ids.extend_from_slice(shell.faces());
    }

    let mut segments: Vec<(Point3, Point3)> = Vec::new();

    for &fid in &face_ids {
        let face = topo.face(fid)?;
        match face.surface() {
            FaceSurface::Plane {
                normal: face_normal,
                d: face_d,
            } => {
                let face_normal = *face_normal;
                let face_d = *face_d;

                let verts = face_polygon(topo, fid)?;
                if let Some(seg) =
                    intersect_planar_face_with_plane(&verts, face_normal, face_d, normal, d, tol)
                {
                    segments.push(seg);
                }
            }
            FaceSurface::Nurbs(nurbs) => {
                for curve in &nurbs_section_curves(nurbs, normal, d)? {
                    chain_curve_points(&curve.points, &mut segments);
                }
            }
            FaceSurface::Cylinder(cyl) => {
                let curves =
                    brepkit_math::analytic_intersection::intersect_plane_cylinder(cyl, normal, d)?;
                for curve in &curves {
                    chain_curve_points(&curve.points, &mut segments);
                }
            }
            FaceSurface::Cone(cone) => {
                let curves =
                    brepkit_math::analytic_intersection::intersect_plane_cone(cone, normal, d)?;
                for curve in &curves {
                    chain_curve_points(&curve.points, &mut segments);
                }
            }
            FaceSurface::Sphere(sphere) => {
                let curves =
                    brepkit_math::analytic_intersection::intersect_plane_sphere(sphere, normal, d)?;
                for curve in &curves {
                    chain_curve_points(&curve.points, &mut segments);
                }
            }
            FaceSurface::Torus(torus) => {
                let curves =
                    brepkit_math::analytic_intersection::intersect_plane_torus(torus, normal, d)?;
                for curve in &curves {
                    chain_curve_points(&curve.points, &mut segments);
                }
            }
        }
    }

    // Analytic faces can share a section curve across two oppositely-oriented
    // faces — e.g. a sphere's two hemispheres each yield the full equatorial
    // circle — emitting every segment twice in opposing directions. Left in
    // place, the wire assembler chains both copies into a single zero-area
    // double loop, so collapse coincident duplicates first.
    dedup_coincident_segments(&mut segments, tol);

    // If no crossing segments were found, the cutting plane may be exactly
    // coplanar with one or more faces. In that case, extract the boundary
    // edges of those coplanar faces as the cross-section.
    if segments.is_empty() {
        let coplanar_segs = extract_coplanar_boundary(topo, &face_ids, normal, d, tol)?;
        segments = coplanar_segs;
    }

    // No crossing segments and no coplanar boundary across outer + inner
    // shells means the cutting plane lies wholly on one side of the solid:
    // a genuine miss. The empty set is a valid section, not an error.
    if segments.is_empty() {
        return Ok(Section { faces: Vec::new() });
    }

    let mut wires = assemble_wires(topo, &segments, normal, d, tol)?;

    // Fallback: if crossing-based segments didn't form closed wires,
    // try using coplanar face boundaries instead (handles the case where
    // the cutting plane exactly coincides with a face of the solid).
    if wires.is_empty() {
        let coplanar_segs = extract_coplanar_boundary(topo, &face_ids, normal, d, tol)?;
        if !coplanar_segs.is_empty() {
            wires = assemble_wires(topo, &coplanar_segs, normal, d, tol)?;
        }
    }

    if wires.is_empty() {
        return Err(crate::OperationsError::InvalidInput {
            reason: "no closed cross-section could be assembled".into(),
        });
    }

    // Group wires by containment: outer wires get their inner wires as holes.
    // Without this, a hollow shape's section would produce two separate faces
    // (outer rectangle + inner hole rectangle) instead of one face-with-hole.
    let groups = group_wires_by_containment(topo, &wires, normal);

    let mut result_faces = Vec::with_capacity(groups.len());
    for (outer, inners) in groups {
        let face = topo.add_face(Face::new(outer, inners, FaceSurface::Plane { normal, d }));
        result_faces.push(face);
    }

    Ok(Section {
        faces: result_faces,
    })
}

/// The plane's sections of a NURBS face that bound regions of the solid's
/// section. A curve along which the face touches the plane from one side
/// bounds none, and a curve along one of the face's own edges is that
/// edge, which its neighbour's section gives. A seam, where the face closes
/// on itself, has no neighbour: its curve stays.
fn nurbs_section_curves(
    nurbs: &NurbsSurface,
    normal: Vec3,
    d: f64,
) -> Result<Vec<IntersectionCurve>, crate::OperationsError> {
    let (u_lo, u_hi) = nurbs.domain_u();
    let (v_lo, v_hi) = nurbs.domain_v();
    let (u_seam, v_seam) = (nurbs.is_periodic_u(), nurbs.is_periodic_v());
    let on_edge = |points: &[IntersectionPoint]| {
        [(0, u_lo), (0, u_hi), (1, v_lo), (1, v_hi)]
            .iter()
            .filter(|&&(k, _)| if k == 0 { !u_seam } else { !v_seam })
            .any(|&(k, edge)| {
                let slack = 1e-9 * if k == 0 { u_hi - u_lo } else { v_hi - v_lo };
                points.iter().all(|p| {
                    let at = if k == 0 { p.param1.0 } else { p.param1.1 };
                    (at - edge).abs() <= slack
                })
            })
    };
    let curves = intersect_plane_nurbs(nurbs, normal, d, 50)?;
    let touching: Vec<bool> = (0..curves.len())
        .map(|k| touches_along(nurbs, normal, d, &curves, k))
        .collect();
    Ok(curves
        .into_iter()
        .zip(touching)
        .filter(|(c, touches)| !touches && !on_edge(&c.points))
        .map(|(c, _)| c)
        .collect())
}

/// Group section wires by containment: each outer wire collects all wires that
/// are spatially inside it as inner (hole) wires. Returns `(outer, [inner])` tuples.
///
/// A wire is "inner" when it sits inside another wire (in the section plane) and
/// is not itself inside any wire that's already inside another. Two-level nesting
/// (holes within holes within holes) is collapsed into the outermost containment;
/// CAD sections don't typically produce deeper nesting.
fn group_wires_by_containment(
    topo: &Topology,
    wires: &[WireId],
    normal: Vec3,
) -> Vec<(WireId, Vec<WireId>)> {
    if wires.len() <= 1 {
        return wires.iter().map(|&w| (w, vec![])).collect();
    }

    // Collect each wire's vertex polygon (used for both bounding-box and
    // point-in-polygon checks below).
    let polygons: Vec<Vec<Point3>> = wires.iter().map(|&wid| wire_vertices(topo, wid)).collect();

    // For each wire, find the smallest wire that strictly contains it.
    // A wire with no strict container is itself an outer wire.
    let mut parent: Vec<Option<usize>> = vec![None; wires.len()];
    for (i, poly_i) in polygons.iter().enumerate() {
        if poly_i.is_empty() {
            continue;
        }
        let sample = poly_i[0];
        let mut best_parent: Option<usize> = None;
        let mut best_area = f64::INFINITY;
        for (j, poly_j) in polygons.iter().enumerate() {
            if i == j || poly_j.len() < 3 {
                continue;
            }
            if crate::distance::point_in_polygon_3d(&sample, poly_j, &normal) {
                let area = polygon_area_3d(poly_j, normal);
                if area < best_area {
                    best_area = area;
                    best_parent = Some(j);
                }
            }
        }
        parent[i] = best_parent;
    }

    // Build groups: each outer wire (no parent) collects its direct children.
    let mut groups: Vec<(WireId, Vec<WireId>)> = Vec::new();
    let mut group_of: std::collections::HashMap<usize, usize> = std::collections::HashMap::new();
    for (i, &wid) in wires.iter().enumerate() {
        if parent[i].is_none() {
            group_of.insert(i, groups.len());
            groups.push((wid, vec![]));
        }
    }
    for (i, &wid) in wires.iter().enumerate() {
        if let Some(p) = parent[i] {
            // Walk up to the outermost ancestor (in case of nested holes).
            let mut top = p;
            while let Some(next) = parent[top] {
                top = next;
            }
            if let Some(&group_idx) = group_of.get(&top) {
                groups[group_idx].1.push(wid);
            }
        }
    }

    groups
}

fn wire_vertices(topo: &Topology, wire_id: WireId) -> Vec<Point3> {
    let Ok(wire) = topo.wire(wire_id) else {
        return vec![];
    };
    let mut pts = Vec::with_capacity(wire.edges().len());
    for oe in wire.edges() {
        let Ok(edge) = topo.edge(oe.edge()) else {
            continue;
        };
        let vid = if oe.is_forward() {
            edge.start()
        } else {
            edge.end()
        };
        if let Ok(v) = topo.vertex(vid) {
            pts.push(v.point());
        }
    }
    pts
}

/// Unsigned area of a 3D planar polygon (projected to dominant axis plane).
fn polygon_area_3d(polygon: &[Point3], normal: Vec3) -> f64 {
    if polygon.len() < 3 {
        return 0.0;
    }
    let ax = normal.x().abs();
    let ay = normal.y().abs();
    let az = normal.z().abs();
    let to_2d = |p: Point3| -> (f64, f64) {
        if az >= ax && az >= ay {
            (p.x(), p.y())
        } else if ay >= ax {
            (p.x(), p.z())
        } else {
            (p.y(), p.z())
        }
    };
    let mut sum = 0.0;
    let n = polygon.len();
    for i in 0..n {
        let (x0, y0) = to_2d(polygon[i]);
        let (x1, y1) = to_2d(polygon[(i + 1) % n]);
        sum += x0 * y1 - x1 * y0;
    }
    (sum * 0.5).abs()
}

type EdgeKey = ((i64, i64, i64), (i64, i64, i64));

/// Quantize a point onto an integer lattice for tolerant endpoint matching.
fn quantize_point(p: Point3, tol: Tolerance) -> (i64, i64, i64) {
    let scale = 1.0 / (tol.linear * 10.0);
    (
        (p.x() * scale).round() as i64,
        (p.y() * scale).round() as i64,
        (p.z() * scale).round() as i64,
    )
}

/// Build an orientation-insensitive key for a segment from its quantized
/// endpoints, so `(a, b)` and `(b, a)` collide.
fn make_edge_key(a: Point3, b: Point3, tol: Tolerance) -> EdgeKey {
    let qa = quantize_point(a, tol);
    let qb = quantize_point(b, tol);
    if qa <= qb { (qa, qb) } else { (qb, qa) }
}

/// Collapse geometrically-coincident segments, ignoring orientation.
///
/// A section curve shared by two oppositely-oriented faces is emitted once in
/// each direction (see the sphere-hemisphere case at the call site). Keying
/// each segment by its unordered quantized endpoints keeps a single copy.
/// Distinct polygon edges share at most one endpoint, so real edges of the
/// section outline are never merged.
fn dedup_coincident_segments(segments: &mut Vec<(Point3, Point3)>, tol: Tolerance) {
    use std::collections::HashSet;

    let mut seen: HashSet<EdgeKey> = HashSet::new();
    segments.retain(|&(a, b)| seen.insert(make_edge_key(a, b, tol)));
}

/// Extract boundary edges of faces coplanar with the cutting plane.
///
/// For each coplanar face, its edges are boundary edges if they are not shared
/// with another coplanar face. These boundary edges form the cross-section
/// outline where the solid meets the cutting plane.
fn extract_coplanar_boundary(
    topo: &Topology,
    face_ids: &[FaceId],
    cut_normal: Vec3,
    cut_d: f64,
    tol: Tolerance,
) -> Result<Vec<(Point3, Point3)>, crate::OperationsError> {
    use std::collections::HashMap;

    // Use a relaxed tolerance for coplanar detection to handle
    // floating-point precision differences across platforms (e.g. WASM).
    let coplanar_tol = tol.linear * 100.0;

    let mut coplanar_faces = Vec::new();
    for &fid in face_ids {
        let verts = face_polygon(topo, fid)?;
        if verts.len() >= 3
            && verts
                .iter()
                .all(|v| (dot_normal_point(cut_normal, *v) - cut_d).abs() < coplanar_tol)
        {
            coplanar_faces.push(fid);
        }
    }

    if coplanar_faces.is_empty() {
        return Ok(Vec::new());
    }

    // Collect all edges of coplanar faces. Each edge is represented by a pair
    // of quantized endpoint coordinates (to handle floating-point matching).
    // An edge shared by two coplanar faces appears twice and is internal.
    // An edge appearing once is a boundary edge.
    let mut edge_counts: HashMap<EdgeKey, (Point3, Point3, usize)> = HashMap::new();

    for &fid in &coplanar_faces {
        let verts = face_polygon(topo, fid)?;
        let n = verts.len();
        for i in 0..n {
            let a = verts[i];
            let b = verts[(i + 1) % n];
            let key = make_edge_key(a, b, tol);
            edge_counts
                .entry(key)
                .and_modify(|e| e.2 += 1)
                .or_insert((a, b, 1));
        }
    }

    // Boundary edges appear exactly once.
    let boundary: Vec<(Point3, Point3)> = edge_counts
        .into_values()
        .filter(|(_, _, count)| *count == 1)
        .map(|(a, b, _)| (a, b))
        .collect();

    Ok(boundary)
}

/// Intersect a planar face polygon with a cutting plane.
///
/// Returns the line segment (if any) where the cutting plane crosses
/// the face polygon. Uses the Sutherland-Hodgman approach: classify
/// each vertex as above/below the plane, find edge crossings.
fn intersect_planar_face_with_plane(
    verts: &[Point3],
    _face_normal: Vec3,
    _face_d: f64,
    cut_normal: Vec3,
    cut_d: f64,
    tol: Tolerance,
) -> Option<(Point3, Point3)> {
    let n = verts.len();
    if n < 3 {
        return None;
    }

    let dists: Vec<f64> = verts
        .iter()
        .map(|v| dot_normal_point(cut_normal, *v) - cut_d)
        .collect();

    // If all vertices lie on the cutting plane the face is coplanar —
    // the cross-section boundary comes from adjacent faces, not this one.
    // Use a relaxed tolerance (100× linear) for cross-platform consistency.
    let coplanar_tol = tol.linear * 100.0;
    if dists.iter().all(|d| d.abs() < coplanar_tol) {
        return None;
    }

    let mut crossings = Vec::new();

    for i in 0..n {
        let j = (i + 1) % n;
        let di = dists[i];
        let dj = dists[j];

        if di.abs() < tol.linear {
            crossings.push(verts[i]);
            continue;
        }

        // Check for sign change (edge crosses plane).
        if (di > tol.linear && dj < -tol.linear) || (di < -tol.linear && dj > tol.linear) {
            let t = di / (di - dj);
            let pi = verts[i];
            let pj = verts[j];
            let ix = Point3::new(
                (pj.x() - pi.x()).mul_add(t, pi.x()),
                (pj.y() - pi.y()).mul_add(t, pi.y()),
                (pj.z() - pi.z()).mul_add(t, pi.z()),
            );
            crossings.push(ix);
        }
    }

    let mut unique = Vec::new();
    for p in &crossings {
        if !unique
            .iter()
            .any(|q: &Point3| (*p - *q).length_squared() < tol.linear * tol.linear)
        {
            unique.push(*p);
        }
    }

    if unique.len() >= 2 {
        Some((unique[0], unique[1]))
    } else {
        None
    }
}

/// Assemble intersection segments into closed wires.
///
/// Chains segments endpoint-to-endpoint using spatial proximity,
/// producing one or more closed wires.
fn assemble_wires(
    topo: &mut Topology,
    segments: &[(Point3, Point3)],
    _normal: Vec3,
    _d: f64,
    tol: Tolerance,
) -> Result<Vec<WireId>, crate::OperationsError> {
    if segments.is_empty() {
        return Ok(vec![]);
    }

    let mut remaining: Vec<(Point3, Point3)> = segments.to_vec();
    let mut wires = Vec::new();

    // Chaining tolerance must be tight enough to keep separate loops apart
    // (e.g., the outer perimeter and the inner hole of a hollow shape) while
    // still tolerating the small endpoint drift produced by tessellated NURBS
    // intersections. `tol.linear * 1000` (≈ 1e-4) balances both.
    //
    // Earlier this defaulted to 50% of the average segment length, which for
    // a hollow 20×20 box section was ≈ 15 — wildly generous, bridging the
    // outer and inner wires into one self-intersecting polygon.
    let chain_tol = tol.linear * 1000.0;

    while !remaining.is_empty() {
        let first = remaining.remove(0);
        let mut chain: Vec<Point3> = vec![first.0, first.1];

        let mut changed = true;
        while changed {
            changed = false;
            let chain_end = chain[chain.len() - 1];
            let threshold_sq = chain_tol * chain_tol;

            let mut best_idx = None;
            let mut best_dist = threshold_sq;
            let mut best_forward = true;

            for i in 0..remaining.len() {
                let (a, b) = remaining[i];
                let dist_a = (a - chain_end).length_squared();
                let dist_b = (b - chain_end).length_squared();

                if dist_a < best_dist {
                    best_dist = dist_a;
                    best_idx = Some(i);
                    best_forward = true;
                }
                if dist_b < best_dist {
                    best_dist = dist_b;
                    best_idx = Some(i);
                    best_forward = false;
                }
            }

            if let Some(idx) = best_idx {
                let (a, b) = remaining.remove(idx);
                chain.push(if best_forward { b } else { a });
                changed = true;
            }
        }

        if chain.len() < 3 {
            continue;
        }

        let start = chain[0];
        let end = chain[chain.len() - 1];
        let closed = (start - end).length_squared() < chain_tol * chain_tol;

        if !closed {
            continue;
        }

        if chain.len() > 3 {
            chain.pop();
        }

        let n = chain.len();
        let vert_ids: Vec<_> = chain
            .iter()
            .map(|&p| topo.add_vertex(Vertex::new(p, tol.linear)))
            .collect();

        let mut oriented_edges = Vec::with_capacity(n);
        for i in 0..n {
            let j = (i + 1) % n;
            let edge = topo.add_edge(Edge::new(vert_ids[i], vert_ids[j], EdgeCurve::Line));
            oriented_edges.push(OrientedEdge::new(edge, true));
        }

        let wire = Wire::new(oriented_edges, true).map_err(crate::OperationsError::Topology)?;
        wires.push(topo.add_wire(wire));
    }

    Ok(wires)
}

#[cfg(test)]
mod tests {
    #![allow(clippy::unwrap_used)]

    use brepkit_math::vec::{Point3, Vec3};
    use brepkit_topology::Topology;
    use brepkit_topology::test_utils::make_unit_cube_manifold;

    use super::*;

    #[test]
    fn section_cube_at_half_height() {
        let mut topo = Topology::new();
        let cube = make_unit_cube_manifold(&mut topo);

        // Cut with a horizontal plane at z=0.5.
        let result = section(
            &mut topo,
            cube,
            Point3::new(0.0, 0.0, 0.5),
            Vec3::new(0.0, 0.0, 1.0),
        )
        .unwrap();

        assert_eq!(
            result.faces.len(),
            1,
            "should produce one cross-section face"
        );

        // The cross section of a unit cube at z=0.5 should be a unit square.
        let area = crate::measure::face_area(&topo, result.faces[0], 0.1).unwrap();
        assert!(
            (area - 1.0).abs() < 1e-6,
            "cross-section area should be ~1.0, got {area}"
        );
    }

    #[test]
    fn section_cube_at_quarter_height() {
        let mut topo = Topology::new();
        let cube = make_unit_cube_manifold(&mut topo);

        let result = section(
            &mut topo,
            cube,
            Point3::new(0.0, 0.0, 0.25),
            Vec3::new(0.0, 0.0, 1.0),
        )
        .unwrap();

        assert_eq!(result.faces.len(), 1);

        // Still a 1×1 square at any height for a cube.
        let area = crate::measure::face_area(&topo, result.faces[0], 0.1).unwrap();
        assert!(
            (area - 1.0).abs() < 1e-6,
            "cross-section area should be ~1.0, got {area}"
        );
    }

    #[test]
    fn section_cube_along_x() {
        let mut topo = Topology::new();
        let cube = make_unit_cube_manifold(&mut topo);

        // Cut with a vertical plane at x=0.5.
        let result = section(
            &mut topo,
            cube,
            Point3::new(0.5, 0.0, 0.0),
            Vec3::new(1.0, 0.0, 0.0),
        )
        .unwrap();

        assert_eq!(result.faces.len(), 1);

        // Cross-section of unit cube at x=0.5 is a 1×1 square in YZ.
        let area = crate::measure::face_area(&topo, result.faces[0], 0.1).unwrap();
        assert!(
            (area - 1.0).abs() < 1e-6,
            "cross-section area should be ~1.0, got {area}"
        );
    }

    /// Section a box-minus-sphere at z=0.
    ///
    /// Box is (0,0,0)-(20,20,20), sphere at origin r=12.
    /// At z=0, the sphere intersects the box creating a circular cutout.
    /// Cross-section = 20×20 rectangle minus circle of radius 12
    ///   (but sphere center is at origin and box starts at 0, so the
    ///   sphere only clips a quarter-circle at the z=0 face corner).
    ///
    /// The exact area depends on how much of the sphere lies within the
    /// box at this plane. At minimum, we verify the section succeeds and
    /// produces a face with positive area less than the full 20×20 = 400.
    #[test]
    fn section_after_boolean_cut() {
        let mut topo = Topology::new();
        let b = crate::primitives::make_box(&mut topo, 20.0, 20.0, 20.0).unwrap();
        let s = crate::primitives::make_sphere(&mut topo, 12.0, 16).unwrap();

        let solid =
            crate::boolean::boolean(&mut topo, crate::boolean::BooleanOp::Cut, b, s).unwrap();

        let sec = section(
            &mut topo,
            solid,
            Point3::new(0.0, 0.0, 0.0),
            Vec3::new(0.0, 0.0, 1.0),
        )
        .unwrap();

        assert!(!sec.faces.is_empty(), "should produce at least one face");

        // The cross-section area should be positive and less than the full
        // box face (400). The sphere removes a quarter-disc of radius 12
        // from the corner, so area ≈ 400 - π(144)/4 ≈ 400 - 113.1 ≈ 286.9.
        let total_area: f64 = sec
            .faces
            .iter()
            .map(|&fid| crate::measure::face_area(&topo, fid, 0.1).unwrap())
            .sum();
        assert!(
            total_area > 200.0,
            "section area should be > 200 (box face minus sphere), got {total_area:.2}"
        );
        assert!(
            total_area < 400.0,
            "section area should be < 400 (full box face), got {total_area:.2}"
        );
    }

    #[test]
    fn section_plane_misses_solid() {
        let mut topo = Topology::new();
        let cube = make_unit_cube_manifold(&mut topo);

        // Plane above the cube — misses entirely.
        let result = section(
            &mut topo,
            cube,
            Point3::new(0.0, 0.0, 5.0),
            Vec3::new(0.0, 0.0, 1.0),
        )
        .unwrap();
        assert!(
            result.faces.is_empty(),
            "plane above cube should produce an empty section"
        );
    }

    #[test]
    fn section_plane_flush_with_top_face() {
        let mut topo = Topology::new();
        let cube = make_unit_cube_manifold(&mut topo);

        // Plane coincident with the cube's top face at z=1 — the coplanar
        // boundary fallback must yield the face outline, not a miss.
        let result = section(
            &mut topo,
            cube,
            Point3::new(0.0, 0.0, 1.0),
            Vec3::new(0.0, 0.0, 1.0),
        )
        .unwrap();
        assert_eq!(
            result.faces.len(),
            1,
            "flush plane should yield the face outline"
        );
        let area = crate::measure::face_area(&topo, result.faces[0], 0.1).unwrap();
        assert!(
            (area - 1.0).abs() < 1e-6,
            "flush-face section area should be ~1.0, got {area}"
        );
    }

    #[test]
    fn section_extruded_box() {
        let mut topo = Topology::new();
        let solid = crate::primitives::make_box(&mut topo, 2.0, 3.0, 4.0).unwrap();

        // Box extends from (0,0,0) to (2,3,4). Cut at z=2 (middle of the box).
        let result = section(
            &mut topo,
            solid,
            Point3::new(0.0, 0.0, 2.0),
            Vec3::new(0.0, 0.0, 1.0),
        )
        .unwrap();

        assert_eq!(result.faces.len(), 1);

        let area = crate::measure::face_area(&topo, result.faces[0], 0.1).unwrap();
        // 2×3 = 6
        assert!(
            (area - 6.0).abs() < 1e-6,
            "cross-section area should be ~6.0, got {area}"
        );
    }

    /// Diagonal section: plane x + y = 0.8, normal = (1,1,0).
    ///
    /// This plane intersects the unit cube at:
    ///   front (y=0): x=0.8 → line (0.8,0,0)-(0.8,0,1)
    ///   left (x=0): y=0.8 → line (0,0.8,0)-(0,0.8,1)
    ///   top (z=1): x+y=0.8 → line (0.8,0,1)-(0,0.8,1)
    ///   bottom (z=0): x+y=0.8 → line (0.8,0,0)-(0,0.8,0)
    ///
    /// The cross-section is a rectangle with:
    ///   width = distance between the two vertical lines
    ///         = |(0.8,0)-(0,0.8)| = √(0.64+0.64) = 0.8√2
    ///   height = 1.0 (z extent)
    ///   area = 0.8√2 ≈ 1.1314
    #[test]
    fn section_diagonal_plane() {
        let mut topo = Topology::new();
        let cube = make_unit_cube_manifold(&mut topo);

        let result = section(
            &mut topo,
            cube,
            Point3::new(0.4, 0.4, 0.5),
            Vec3::new(1.0, 1.0, 0.0),
        );

        assert!(
            result.is_ok(),
            "diagonal plane should intersect cube: {:?}",
            result.err()
        );
        let sec = result.unwrap();
        assert_eq!(sec.faces.len(), 1);

        // Area = width × height = 0.8√2 × 1.0 ≈ 1.1314
        let area = crate::measure::face_area(&topo, sec.faces[0], 0.01).unwrap();
        let expected = 0.8 * std::f64::consts::SQRT_2;
        let rel_err = (area - expected).abs() / expected;
        assert!(
            rel_err < 1e-4,
            "diagonal section area should be 0.8√2 ≈ {expected:.4}, got {area:.4} \
             (rel_err={rel_err:.2e})"
        );
    }

    /// Section of a cylinder at mid-height → circular cross-section.
    /// Cylinder r=5, h=10, section at z=5 → circle area = πr² = 25π ≈ 78.54.
    #[test]
    fn section_cylinder_at_midheight() {
        let mut topo = Topology::new();
        let solid = crate::primitives::make_cylinder(&mut topo, 5.0, 10.0).unwrap();

        let result = section(
            &mut topo,
            solid,
            Point3::new(0.0, 0.0, 5.0),
            Vec3::new(0.0, 0.0, 1.0),
        );

        assert!(
            result.is_ok(),
            "section of cylinder should succeed: {:?}",
            result.err()
        );
        let sec = result.unwrap();
        assert!(!sec.faces.is_empty(), "should produce at least one face");

        let total_area: f64 = sec
            .faces
            .iter()
            .map(|&fid| crate::measure::face_area(&topo, fid, 0.01).unwrap())
            .sum();
        // Circle area = πr² = 25π ≈ 78.54
        let expected = std::f64::consts::PI * 25.0;
        let rel_err = (total_area - expected).abs() / expected;
        assert!(
            rel_err < 0.05,
            "cylinder section area should be πr² = {expected:.2}, got {total_area:.2} \
             (rel_err={rel_err:.2e})"
        );
    }

    /// Section of a sphere through its center → a disk bounded by the great
    /// circle. A bare sphere is two hemisphere faces that each yield the full
    /// equatorial circle, so the section must collapse the duplicate into one
    /// disk rather than a zero-area double loop (#860). r=12 at z=0 → πr² =
    /// 144π ≈ 452.39.
    #[test]
    fn section_sphere_through_center() {
        let mut topo = Topology::new();
        let solid = crate::primitives::make_sphere(&mut topo, 12.0, 16).unwrap();

        let sec = section(
            &mut topo,
            solid,
            Point3::new(0.0, 0.0, 0.0),
            Vec3::new(0.0, 0.0, 1.0),
        )
        .unwrap();

        assert_eq!(sec.faces.len(), 1, "sphere section should be a single disk");

        let total_area: f64 = sec
            .faces
            .iter()
            .map(|&fid| crate::measure::face_area(&topo, fid, 0.01).unwrap())
            .sum();
        let expected = std::f64::consts::PI * 144.0;
        let rel_err = (total_area - expected).abs() / expected;
        assert!(
            rel_err < 0.05,
            "sphere great-circle section area should be πr² = {expected:.2}, got \
             {total_area:.2} (rel_err={rel_err:.2e})"
        );
    }

    /// Section of a sphere off-center → a smaller disk. r=12 at z=5 →
    /// circle radius √(144−25)=√119, area = 119π ≈ 373.85.
    #[test]
    fn section_sphere_off_center() {
        let mut topo = Topology::new();
        let solid = crate::primitives::make_sphere(&mut topo, 12.0, 16).unwrap();

        let sec = section(
            &mut topo,
            solid,
            Point3::new(0.0, 0.0, 5.0),
            Vec3::new(0.0, 0.0, 1.0),
        )
        .unwrap();

        assert_eq!(sec.faces.len(), 1, "sphere section should be a single disk");

        let total_area: f64 = sec
            .faces
            .iter()
            .map(|&fid| crate::measure::face_area(&topo, fid, 0.01).unwrap())
            .sum();
        let expected = std::f64::consts::PI * 119.0;
        let rel_err = (total_area - expected).abs() / expected;
        assert!(
            rel_err < 0.05,
            "sphere off-center section area should be {expected:.2}, got {total_area:.2} \
             (rel_err={rel_err:.2e})"
        );
    }

    /// The Bernstein coefficients of the polynomial whose power coefficients
    /// are `power`, lowest first.
    fn bernstein(power: &[f64]) -> Vec<f64> {
        let n = power.len() - 1;
        let choose =
            |n: usize, k: usize| (0..k).fold(1.0, |c, m| c * (n - m) as f64 / (m + 1) as f64);
        (0..=n)
            .map(|i| {
                (0..=i)
                    .map(|k| choose(i, k) / choose(n, k) * power[k])
                    .sum()
            })
            .collect()
    }

    /// A sheet `z = (x − a)³(x − b)` over `x` from 0 to 1, crossing the plane
    /// `z = 0` along `x = a` with its normal along the plane's and along
    /// `x = b` 0.0008 away, bounds the section along both: the side test's
    /// probes start close enough not to reach past `x = b`.
    #[test]
    fn a_nurbs_face_crossing_the_plane_flat_beside_another_crossing_bounds_both() {
        // Round a node of the 50-sample grid, so that the grid sees both.
        let (a, b) = (24.0 / 49.0 - 0.0004, 24.0 / 49.0 + 0.0004);
        let power = [
            a * a * a * b,
            -a * a * (a + 3.0 * b),
            3.0 * a * (a + b),
            -(3.0 * a + b),
            1.0,
        ];
        let rows = bernstein(&power)
            .into_iter()
            .enumerate()
            .map(|(i, z)| {
                let x = i as f64 / 4.0;
                vec![Point3::new(x, -1.0, z), Point3::new(x, 1.0, z)]
            })
            .collect();
        let sheet = NurbsSurface::new(
            4,
            1,
            [vec![0.0; 5], vec![1.0; 5]].concat(),
            vec![0.0, 0.0, 1.0, 1.0],
            rows,
            vec![vec![1.0; 2]; 5],
        )
        .unwrap();
        let curves = nurbs_section_curves(&sheet, Vec3::new(0.0, 0.0, 1.0), 0.0).unwrap();
        assert_eq!(curves.len(), 2);
    }

    /// A sheet stretched a thousandfold along `x`, `z = (u − a)(u − b)` with
    /// `x = 1000u`, crosses `z = 0` along two lines 8e-7 apart in `u`, where
    /// it lies within 1e-9 of flat: no probe of the side test reaches past
    /// the other line, and the section keeps both.
    #[test]
    fn a_nurbs_face_crossing_the_plane_twice_close_by_bounds_both() {
        // Round a node of the 50-sample grid, so that the grid sees both.
        let (a, b) = (24.0 / 49.0 - 4e-7, 24.0 / 49.0 + 4e-7);
        let heights = [
            a * b,
            0.5f64.mul_add(-(a + b), a * b),
            1.0 - (a + b) + a * b,
        ];
        let rows = heights
            .iter()
            .enumerate()
            .map(|(i, &z)| {
                let x = 500.0 * i as f64;
                vec![Point3::new(x, -1.0, z), Point3::new(x, 1.0, z)]
            })
            .collect();
        let sheet = NurbsSurface::new(
            2,
            1,
            vec![0.0, 0.0, 0.0, 1.0, 1.0, 1.0],
            vec![0.0, 0.0, 1.0, 1.0],
            rows,
            vec![vec![1.0; 2]; 3],
        )
        .unwrap();
        let curves = nurbs_section_curves(&sheet, Vec3::new(0.0, 0.0, 1.0), 0.0).unwrap();
        assert_eq!(curves.len(), 2);
    }

    /// Across a closed direction's seam, the side test's gap to another
    /// curve of the section is measured the shorter way round, passing over
    /// a curve's own copy across the seam but not a close neighbour.
    #[test]
    fn the_gap_to_another_curve_runs_round_a_closed_seam() {
        use brepkit_math::nurbs::intersection::{IntersectionCurve as IC, IntersectionPoint as IP};
        // Two lines of the section on either side of a seam at `u = 0`, at
        // `u = 0.9996` and `u = 0.0004`.
        let line = |u: f64| {
            let points: Vec<_> = [0.2, 0.5, 0.8]
                .map(|v| IP {
                    point: Point3::new(u, v, 0.0),
                    param1: (u, v),
                    param2: (0.0, 0.0),
                })
                .to_vec();
            IC {
                curve: brepkit_math::nurbs::curve::NurbsCurve::new(
                    1,
                    vec![0.0, 0.0, 1.0, 1.0],
                    vec![points[0].point, points[2].point],
                    vec![1.0, 1.0],
                )
                .unwrap(),
                points,
            }
        };
        let curves = [line(0.9996), line(0.0004)];
        let round = gap_to_others(
            &curves,
            0,
            (0.9996, 0.5),
            ((0.0, 1.0), (0.0, 1.0)),
            (true, false),
        );
        let open = gap_to_others(
            &curves,
            0,
            (0.9996, 0.5),
            ((0.0, 1.0), (0.0, 1.0)),
            (false, false),
        );
        assert!((round - 0.0008).abs() < 1e-12, "{round}");
        assert!((open - 0.9992).abs() < 1e-12, "{open}");
        // A line along the seam comes back at `u = 1` as well as `u = 0`,
        // and that copy is no neighbour; lines 8e-10 apart are, astride the
        // seam or on an open surface.
        let copy = [line(0.0), line(1.0)];
        let gap = gap_to_others(
            &copy,
            0,
            (0.0, 0.5),
            ((0.0, 1.0), (0.0, 1.0)),
            (true, false),
        );
        assert!(gap.is_infinite(), "{gap}");
        let astride = [line(4e-10), line(1.0 - 4e-10)];
        let gap = gap_to_others(
            &astride,
            0,
            (4e-10, 0.5),
            ((0.0, 1.0), (0.0, 1.0)),
            (true, false),
        );
        assert!((gap - 8e-10).abs() < 1e-15, "{gap}");
        let close = [line(0.5), line(0.5 + 8e-10)];
        let gap = gap_to_others(
            &close,
            0,
            (0.5, 0.5),
            ((0.0, 1.0), (0.0, 1.0)),
            (false, false),
        );
        assert!((gap - 8e-10).abs() < 1e-15, "{gap}");
    }

    /// A closed rational NURBS cylinder of radius 1 about `z`, whose seam
    /// runs along `(1, 0, z)`, touches the plane `x = 1` along its seam: the
    /// side test's probes wrap round the seam, and the contact bounds no
    /// region of the section.
    #[test]
    fn a_closed_nurbs_face_touching_the_plane_along_its_seam_bounds_nothing() {
        let w = std::f64::consts::FRAC_1_SQRT_2;
        let ring = [
            (1.0, 0.0, 1.0),
            (1.0, 1.0, w),
            (0.0, 1.0, 1.0),
            (-1.0, 1.0, w),
            (-1.0, 0.0, 1.0),
            (-1.0, -1.0, w),
            (0.0, -1.0, 1.0),
            (1.0, -1.0, w),
            (1.0, 0.0, 1.0),
        ];
        let cylinder = NurbsSurface::new(
            2,
            1,
            vec![
                0.0, 0.0, 0.0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1.0, 1.0, 1.0,
            ],
            vec![0.0, 0.0, 1.0, 1.0],
            ring.iter()
                .map(|&(x, y, _)| vec![Point3::new(x, y, 0.0), Point3::new(x, y, 1.0)])
                .collect(),
            ring.iter().map(|&(_, _, w)| vec![w, w]).collect(),
        )
        .unwrap();
        assert!(cylinder.is_periodic_u());
        let curves = nurbs_section_curves(&cylinder, Vec3::new(1.0, 0.0, 0.0), 1.0).unwrap();
        assert!(curves.is_empty(), "{} curves", curves.len());
    }

    /// A NURBS face crossing the plane where its normal lies along the
    /// plane's, as `z = x³` and `z = x⁷` cross `z = 0` along `x = 0`, bounds
    /// the section there; a trough whose floor only touches the plane does
    /// not.
    #[test]
    fn a_nurbs_face_crossing_the_plane_where_it_lies_flat_bounds_the_section() {
        let sheet = |rows: Vec<Vec<Point3>>| {
            let degree = rows.len() - 1;
            let knots = [vec![0.0; degree + 1], vec![1.0; degree + 1]].concat();
            NurbsSurface::new(
                degree,
                1,
                knots,
                vec![0.0, 0.0, 1.0, 1.0],
                rows,
                vec![vec![1.0; 2]; degree + 1],
            )
            .unwrap()
        };
        let row = |x: f64, z: f64| vec![Point3::new(x, -1.0, z), Point3::new(x, 1.0, z)];
        let up = Vec3::new(0.0, 0.0, 1.0);
        for power in [3_u32, 7] {
            // `(2u − 1)ⁿ` has Bernstein coefficients `(−1)ⁿ⁻ᵏ`.
            let odd = sheet(
                (0..=power)
                    .map(|k| {
                        let x = 2.0f64.mul_add(f64::from(k) / f64::from(power), -1.0);
                        row(x, if (power - k) % 2 == 0 { 1.0 } else { -1.0 })
                    })
                    .collect(),
            );
            assert_eq!(
                nurbs_section_curves(&odd, up, 0.0).unwrap().len(),
                1,
                "x^{power}"
            );
        }
        let trough = sheet(vec![row(-1.0, 1.0), row(0.0, -1.0), row(1.0, 1.0)]);
        assert!(nurbs_section_curves(&trough, up, 0.0).unwrap().is_empty());
    }
}