use std::collections::HashMap;

use bone_types::{
    EdgeLabel, FaceLabel, Length, Plane3, Point3, Tolerance, UnitVec3, Vec3, VertexLabel,
};
use truck_modeling::{Curve, EdgeID, Surface};
use uom::si::length::millimeter;

use super::build;
use super::convert::point_from_truck;
use super::edges::EdgeCurve3;
use super::{BrepError, BrepSolid, LabelKind};
use classify::FaceClass;

pub mod classify;
pub mod cyl;
pub mod decompose;
pub mod imprint;
pub mod model;
pub mod ssi;
pub mod stitch;

pub(crate) use model::{Body, Cyl, Face, PlanarBody, PlanarFace, PlanarLoop, PlanarVertex, Tube};

#[derive(Copy, Clone, Debug, PartialEq, Eq)]
pub enum BooleanOp {
    Union,
    Difference,
    Intersection,
}

pub fn boolean(
    a: &BrepSolid,
    b: &BrepSolid,
    op: BooleanOp,
    tolerance: Tolerance,
) -> Result<BrepSolid, BrepError> {
    if op == BooleanOp::Difference
        && let Some(tool) = cyl::recognize_cylinder(b, tolerance)
    {
        return match cyl::cylinder_cut(a, &tool, tolerance) {
            Ok(result) => Ok(result),
            Err(fast_error) => fall_back(a, b, op, tolerance, fast_error),
        };
    }
    general_boolean(a, b, op, tolerance)
}

fn fall_back(
    a: &BrepSolid,
    b: &BrepSolid,
    op: BooleanOp,
    tolerance: Tolerance,
    fast_error: BrepError,
) -> Result<BrepSolid, BrepError> {
    match general_boolean(a, b, op, tolerance) {
        Ok(result) if result.validate(tolerance).is_ok() => Ok(result),
        _ => Err(fast_error),
    }
}

fn general_boolean(
    a: &BrepSolid,
    b: &BrepSolid,
    op: BooleanOp,
    tolerance: Tolerance,
) -> Result<BrepSolid, BrepError> {
    let soup_a = classify::SolidSoup::new(a, tolerance)?;
    let soup_b = classify::SolidSoup::new(b, tolerance)?;
    let body_a = extract_body(a, &soup_a, tolerance)?;
    let body_b = extract_body(b, &soup_b, tolerance)?;
    let arranged_a = imprint::arrange_body(&body_a, &body_b, tolerance)?;
    let arranged_b = imprint::arrange_body(&body_b, &body_a, tolerance)?;
    let kept = select(op, &arranged_a, &arranged_b, &soup_a, &soup_b, tolerance)?;
    let kept = imprint::conform_body(&kept, tolerance);
    let result = stitch::assemble_result(&kept, tolerance)?;
    if volume_plausible(&result, a, b, op) {
        Ok(result)
    } else {
        Err(BrepError::BooleanIndeterminate)
    }
}

fn is_curved(solid: &BrepSolid) -> bool {
    solid
        .arena
        .solid()
        .face_iter()
        .any(|face| !matches!(face.surface(), Surface::Plane(_)))
}

fn mesh_volume(solid: &BrepSolid) -> Option<f64> {
    let mesh = solid
        .tessellate(
            bone_types::ChordHeightTolerance::from_mm(0.01),
            bone_types::AngleTolerance::from_radians(0.05),
        )
        .ok()?;
    Some(
        mesh.faces()
            .iter()
            .flat_map(|slab| {
                let positions = slab.positions().to_vec();
                slab.triangles().iter().map(move |tri| {
                    let corner = |index: u32| positions[index as usize].coords_mm();
                    let (ax, ay, az) = corner(tri[0]);
                    let (bx, by, bz) = corner(tri[1]);
                    let (cx, cy, cz) = corner(tri[2]);
                    let cross = (by * cz - bz * cy, bz * cx - bx * cz, bx * cy - by * cx);
                    (ax * cross.0 + ay * cross.1 + az * cross.2) / 6.0
                })
            })
            .sum(),
    )
}

fn volume_plausible(result: &BrepSolid, a: &BrepSolid, b: &BrepSolid, op: BooleanOp) -> bool {
    if !is_curved(a) && !is_curved(b) {
        return true;
    }
    let (Some(vr), Some(va), Some(vb)) = (mesh_volume(result), mesh_volume(a), mesh_volume(b))
    else {
        return true;
    };
    let slack = 0.02 * va.abs().max(vb.abs()) + 1.0e-3;
    match op {
        BooleanOp::Difference => vr <= va + slack,
        BooleanOp::Intersection => vr <= va.min(vb) + slack,
        BooleanOp::Union => vr <= va + vb + slack,
    }
}

pub(crate) fn as_vec(u: UnitVec3) -> Vec3 {
    let (x, y, z) = u.components();
    Vec3::from_mm(x, y, z)
}

pub(crate) fn unit(v: Vec3, tolerance: Tolerance) -> Option<UnitVec3> {
    let (x, y, z) = v.coords_mm();
    UnitVec3::try_from_components(x, y, z, tolerance).ok()
}

pub(crate) fn basis(normal: UnitVec3, tolerance: Tolerance) -> Option<(UnitVec3, UnitVec3)> {
    let seed = if normal.dot(UnitVec3::x_axis()).abs() <= 0.9 {
        UnitVec3::x_axis()
    } else {
        UnitVec3::y_axis()
    };
    let x = normal.cross(seed, tolerance).ok()?;
    let y = normal.cross(x, tolerance).ok()?;
    Some((x, y))
}

fn newell(points: &[Point3]) -> Vec3 {
    let count = points.len();
    (0..count).fold(Vec3::from_mm(0.0, 0.0, 0.0), |acc, i| {
        let (x0, y0, z0) = points[i].coords_mm();
        let (x1, y1, z1) = points[(i + 1) % count].coords_mm();
        acc + Vec3::from_mm(
            (y0 - y1) * (z0 + z1),
            (z0 - z1) * (x0 + x1),
            (x0 - x1) * (y0 + y1),
        )
    })
}

pub(crate) fn extract_body(
    solid: &BrepSolid,
    soup: &classify::SolidSoup,
    tolerance: Tolerance,
) -> Result<Body, BrepError> {
    let vertex_labels: Vec<(Point3, VertexLabel)> = solid
        .iter_vertices()
        .map(|vertex| (vertex.position(), vertex.label()))
        .collect();
    let edges = edge_index(solid);
    let mut body: Vec<Face> = solid
        .arena
        .solid()
        .face_iter()
        .filter(|face| matches!(face.surface(), Surface::Plane(_)))
        .map(|face| extract_planar_face(solid, face, &vertex_labels, &edges, tolerance))
        .collect::<Result<Vec<Vec<Face>>, _>>()?
        .into_iter()
        .flatten()
        .collect();
    let cap_specs: std::collections::BTreeMap<FaceLabel, Annulus> = solid
        .arena
        .solid()
        .face_iter()
        .filter_map(|face| {
            let annulus = rev_annulus(&face.surface(), tolerance)?;
            Some((bone_label(solid, face)?, annulus))
        })
        .collect();
    let caps = cap_specs
        .into_iter()
        .map(|(label, annulus)| synth_annular_face(solid, label, &annulus, soup, tolerance))
        .collect::<Result<Vec<_>, _>>()?;
    body.extend(caps);
    let tube_labels: std::collections::BTreeSet<FaceLabel> = solid
        .arena
        .solid()
        .face_iter()
        .filter(|face| !matches!(face.surface(), Surface::Plane(_)))
        .filter(|face| rev_annulus(&face.surface(), tolerance).is_none())
        .filter_map(|face| bone_label(solid, face))
        .collect();
    let tubes = tube_labels
        .into_iter()
        .map(|label| extract_tube_for_label(solid, label, soup, tolerance).map(Face::tube))
        .collect::<Result<Vec<_>, _>>()?;
    body.extend(tubes);
    Ok(body)
}

fn extract_tube_for_label(
    solid: &BrepSolid,
    label: FaceLabel,
    soup: &classify::SolidSoup,
    tolerance: Tolerance,
) -> Result<Tube, BrepError> {
    extract_tube_from_bone(solid, label, soup, tolerance).or_else(|_| {
        let face = solid
            .arena
            .solid()
            .face_iter()
            .find(|face| bone_label(solid, face) == Some(label))
            .ok_or(BrepError::MissingLabel {
                kind: LabelKind::Face,
            })?;
        tube_from_surface(solid, face, label, soup, tolerance)
    })
}

fn bone_label(solid: &BrepSolid, face: &truck_modeling::Face) -> Option<FaceLabel> {
    solid
        .arena
        .face_index()
        .get(&face.id())
        .and_then(|brep_id| solid.faces.get(*brep_id))
        .map(super::BrepFace::label)
}

pub(crate) fn extract_planar_body(
    solid: &BrepSolid,
    tolerance: Tolerance,
) -> Result<PlanarBody, BrepError> {
    let soup = classify::SolidSoup::new(solid, tolerance)?;
    extract_body(solid, &soup, tolerance)?
        .into_iter()
        .map(|face| match face {
            Face::Planar(planar) => Ok(*planar),
            Face::Tube(_) => Err(unsupported(super::BooleanGap::CurvedFace)),
        })
        .collect()
}

fn edge_index(solid: &BrepSolid) -> HashMap<EdgeID, (EdgeLabel, EdgeCurve3)> {
    solid
        .iter_edges()
        .flat_map(|edge| {
            let label = edge.label();
            let curve = edge.curve().clone();
            edge.handles
                .iter()
                .map(move |handle| (solid.arena.edge(*handle).id(), (label, curve.clone())))
        })
        .collect()
}

pub(crate) struct RawLoop {
    pub verts: Vec<PlanarVertex>,
    pub edges: Vec<EdgeLabel>,
    pub curves: Vec<EdgeCurve3>,
}

type EdgeIndex = HashMap<EdgeID, (EdgeLabel, EdgeCurve3)>;

fn extract_planar_face(
    solid: &BrepSolid,
    face: &truck_modeling::Face,
    vertex_labels: &[(Point3, VertexLabel)],
    edges: &EdgeIndex,
    tolerance: Tolerance,
) -> Result<Vec<Face>, BrepError> {
    let label = bone_label(solid, face).ok_or(BrepError::MissingLabel {
        kind: LabelKind::Face,
    })?;
    extract_planar(label, face, vertex_labels, edges, tolerance)
}

fn is_line(curve: &EdgeCurve3) -> bool {
    matches!(curve, EdgeCurve3::Line(_))
}

fn extract_planar(
    label: FaceLabel,
    face: &truck_modeling::Face,
    vertex_labels: &[(Point3, VertexLabel)],
    edges: &EdgeIndex,
    tolerance: Tolerance,
) -> Result<Vec<Face>, BrepError> {
    let loops = raw_loops(face, vertex_labels, edges, tolerance)?;
    let outer_index = (0..loops.len())
        .max_by(|&a, &b| loop_area(&loops[a]).total_cmp(&loop_area(&loops[b])))
        .ok_or(BrepError::BooleanIndeterminate)?;
    let outer = &loops[outer_index];
    let outer_positions: Vec<Point3> = outer.verts.iter().map(|vertex| vertex.position).collect();
    let plane = plane_from_loop(face, &outer_positions, tolerance)?;

    let curved = loops
        .iter()
        .any(|item| item.curves.iter().any(|curve| !is_line(curve)));
    if curved {
        let holes = loops
            .iter()
            .enumerate()
            .filter(|(index, _)| *index != outer_index)
            .map(|(_, item)| PlanarLoop {
                verts: item.verts.clone(),
                edges: item.edges.clone(),
                curves: item.curves.clone(),
            })
            .collect();
        return Ok(vec![Face::planar(PlanarFace {
            label,
            plane,
            polygon: outer.verts.clone(),
            edges: outer.edges.clone(),
            curves: outer.curves.clone(),
            holes,
        })]);
    }

    if loops.len() == 1 && decompose::loop_is_convex(&outer.verts, &plane, tolerance) {
        return Ok(vec![Face::planar(PlanarFace::straight(
            label,
            plane,
            outer.verts.clone(),
            outer.edges.clone(),
        ))]);
    }
    let face_loops: Vec<decompose::FaceLoop> = loops
        .iter()
        .map(|item| decompose::FaceLoop {
            verts: item.verts.clone(),
            edges: item.edges.clone(),
        })
        .collect();
    Ok(decompose::decompose(label, &plane, &face_loops, tolerance)
        .into_iter()
        .map(Face::planar)
        .collect())
}

fn revolution_line(
    surface: &Surface,
    tolerance: Tolerance,
) -> Option<(Point3, Point3, Point3, UnitVec3)> {
    let Surface::RevolutedCurve(rev) = surface else {
        return None;
    };
    let Curve::Line(line) = rev.entity_curve() else {
        return None;
    };
    let origin = point_from_truck(rev.origin());
    let ax = rev.axis();
    let axis = model::canonical_axis(unit(Vec3::from_mm(ax.x, ax.y, ax.z), tolerance)?);
    Some((
        point_from_truck(line.0),
        point_from_truck(line.1),
        origin,
        axis,
    ))
}

fn cyl_from_revolution(surface: &Surface, tolerance: Tolerance) -> Option<(Cyl, f64, f64)> {
    let (p0, p1, origin, axis) = revolution_line(surface, tolerance)?;
    let dir = unit(p1 - p0, tolerance)?;
    if model::dot(dir, axis).abs() < 1.0 - tolerance.value() {
        return None;
    }
    let radius_mm = model::radial_distance(origin, axis, p0);
    if radius_mm <= tolerance.value() {
        return None;
    }
    let z0 = model::project_onto_axis(origin, axis, p0);
    let z1 = model::project_onto_axis(origin, axis, p1);
    let cyl = Cyl::new(
        origin,
        axis,
        Length::new::<millimeter>(radius_mm),
        tolerance,
    )?;
    Some((cyl, z0.min(z1), z0.max(z1)))
}

struct Annulus {
    center: Point3,
    axis: UnitVec3,
    inner_r: f64,
    outer_r: f64,
}

fn rev_annulus(surface: &Surface, tolerance: Tolerance) -> Option<Annulus> {
    let (p0, p1, origin, axis) = revolution_line(surface, tolerance)?;
    let z0 = model::project_onto_axis(origin, axis, p0);
    let z1 = model::project_onto_axis(origin, axis, p1);
    if (z1 - z0).abs() > tolerance.value() {
        return None;
    }
    let r0 = model::radial_distance(origin, axis, p0);
    let r1 = model::radial_distance(origin, axis, p1);
    let inner_r = r0.min(r1);
    let outer_r = r0.max(r1);
    if outer_r - inner_r <= tolerance.value() {
        return None;
    }
    Some(Annulus {
        center: origin + as_vec(axis) * z0,
        axis,
        inner_r,
        outer_r,
    })
}

const PROBE_FRACTIONS: [f64; 5] = [0.4, 0.2, 0.08, 0.03, 0.012];

fn wall_outward(soup: &classify::SolidSoup, mid: Point3, radial: UnitVec3, radius_mm: f64) -> bool {
    PROBE_FRACTIONS
        .iter()
        .find_map(|&fraction| {
            let offset = (radius_mm * fraction).max(0.05);
            let inner = soup.inside(mid - as_vec(radial) * offset)?;
            let outer = soup.inside(mid + as_vec(radial) * offset)?;
            (inner != outer).then_some(inner)
        })
        .unwrap_or_else(|| {
            soup.inside(mid - as_vec(radial) * (0.4 * radius_mm))
                .unwrap_or(true)
        })
}

fn synth_tube_vertex(label: FaceLabel, point: Point3) -> VertexLabel {
    VertexLabel {
        feature: label.feature,
        role: bone_types::VertexRole::Cut {
            key: stitch::cut_key(point),
        },
    }
}

fn wall_label_at(
    solid: &BrepSolid,
    center: Point3,
    radius_mm: f64,
    tolerance: Tolerance,
) -> FaceLabel {
    solid
        .arena
        .solid()
        .face_iter()
        .find(|face| {
            cyl_from_revolution(&face.surface(), tolerance).is_some_and(|(cyl, lo, hi)| {
                (cyl.radius_mm() - radius_mm).abs() <= tolerance.value() && {
                    let z = model::project_onto_axis(cyl.axis_origin, cyl.axis, center);
                    z >= lo - tolerance.value() && z <= hi + tolerance.value()
                }
            })
        })
        .and_then(|face| bone_label(solid, face))
        .unwrap_or(FaceLabel {
            feature: center_feature(solid),
            role: bone_types::FaceRole::Imported {
                ordinal: bone_types::ImportOrdinal::new(0),
            },
        })
}

fn center_feature(solid: &BrepSolid) -> bone_types::FeatureId {
    solid
        .iter_faces()
        .next()
        .map_or_else(bone_types::FeatureId::default, |face| face.label().feature)
}

fn synth_circle_loop(
    label: FaceLabel,
    rim: FaceLabel,
    plane: Plane3,
    radius_mm: f64,
    tolerance: Tolerance,
) -> Option<(Vec<PlanarVertex>, Vec<EdgeLabel>, Vec<EdgeCurve3>)> {
    let circle =
        crate::circle3::Circle3::new(plane, Length::new::<millimeter>(radius_mm), tolerance)
            .ok()?;
    let point = plane.origin() + as_vec(plane.x_axis()) * radius_mm;
    let vertex = PlanarVertex {
        position: point,
        label: synth_tube_vertex(label, point),
    };
    Some((
        vec![vertex],
        vec![model::elabel_intersection(label, rim)],
        vec![EdgeCurve3::Circle(circle)],
    ))
}

fn synth_annular_face(
    solid: &BrepSolid,
    label: FaceLabel,
    annulus: &Annulus,
    soup: &classify::SolidSoup,
    tolerance: Tolerance,
) -> Result<Face, BrepError> {
    let mid_r = f64::midpoint(annulus.inner_r, annulus.outer_r);
    let (px, _) = basis(annulus.axis, tolerance).ok_or(BrepError::BooleanIndeterminate)?;
    let on_face = annulus.center + as_vec(px) * mid_r;
    let outward_axis = if soup
        .inside(on_face + as_vec(annulus.axis) * (10.0 * tolerance.value()))
        .unwrap_or(false)
    {
        annulus.axis.reversed()
    } else {
        annulus.axis
    };
    let (x, y) = basis(outward_axis, tolerance).ok_or(BrepError::BooleanIndeterminate)?;
    let plane = Plane3::new_unchecked(annulus.center, x, y);
    let outer_rim = wall_label_at(solid, annulus.center, annulus.outer_r, tolerance);
    let inner_rim = wall_label_at(solid, annulus.center, annulus.inner_r, tolerance);
    let (verts, edges, curves) =
        synth_circle_loop(label, outer_rim, plane, annulus.outer_r, tolerance)
            .ok_or(BrepError::BooleanIndeterminate)?;
    let (hverts, hedges, hcurves) =
        synth_circle_loop(label, inner_rim, plane, annulus.inner_r, tolerance)
            .ok_or(BrepError::BooleanIndeterminate)?;
    Ok(Face::planar(PlanarFace {
        label,
        plane,
        polygon: verts,
        edges,
        curves,
        holes: vec![PlanarLoop {
            verts: hverts,
            edges: hedges,
            curves: hcurves,
        }],
    }))
}

fn tube_from_surface(
    solid: &BrepSolid,
    face: &truck_modeling::Face,
    label: FaceLabel,
    soup: &classify::SolidSoup,
    tolerance: Tolerance,
) -> Result<Tube, BrepError> {
    let eps = tolerance.value();
    let (cyl, z0, z1) = cyl_from_revolution(&face.surface(), tolerance)
        .ok_or(unsupported(super::BooleanGap::CurvedFace))?;
    if (z1 - z0) <= eps {
        return Err(unsupported(super::BooleanGap::CurvedFace));
    }
    let lo_center = cyl.axis_origin + as_vec(cyl.axis) * z0;
    let hi_center = cyl.axis_origin + as_vec(cyl.axis) * z1;
    let radius_mm = cyl.radius_mm();
    let rim_bottom = model::elabel_intersection(
        label,
        cap_label_at(solid, lo_center, cyl.axis, tolerance).unwrap_or(label),
    );
    let rim_top = model::elabel_intersection(
        label,
        cap_label_at(solid, hi_center, cyl.axis, tolerance).unwrap_or(label),
    );
    let seam = scan_tube_faces(solid, label)
        .and_then(|scan| scan.seam)
        .unwrap_or_else(|| model::elabel_intersection(label, label));
    let radial = cyl
        .radial(0.0, tolerance)
        .ok_or(BrepError::BooleanIndeterminate)?;
    let mid = cyl.point(0.0, f64::midpoint(z0, z1));
    Ok(Tube {
        label,
        cyl,
        theta0: 0.0,
        theta1: core::f64::consts::TAU,
        z0,
        z1,
        full: true,
        rim_bottom,
        rim_top,
        seam,
        vbot: synth_tube_vertex(label, lo_center),
        vtop: synth_tube_vertex(label, hi_center),
        outward: wall_outward(soup, mid, radial, radius_mm),
    })
}

fn cap_label_at(
    solid: &BrepSolid,
    center: Point3,
    axis: UnitVec3,
    tolerance: Tolerance,
) -> Option<FaceLabel> {
    solid.iter_faces().find_map(|face| {
        let plane = Plane3::from(solid.face_plane_basis(face.id())?);
        let perpendicular = plane.normal().dot(axis).abs() >= 1.0 - tolerance.value();
        let through = (plane.origin() - center).dot_mm2(as_vec(axis)).abs() <= tolerance.value();
        (perpendicular && through).then(|| face.label())
    })
}

fn plane_from_loop(
    face: &truck_modeling::Face,
    outer_positions: &[Point3],
    tolerance: Tolerance,
) -> Result<Plane3, BrepError> {
    if outer_positions.len() >= 3
        && let Ok(plane) = plane_from_polygon(outer_positions, tolerance)
    {
        return Ok(plane);
    }
    let basis = build::face_plane_basis(face).ok_or(BrepError::BooleanIndeterminate)?;
    Ok(Plane3::from(basis))
}

struct RimEdge {
    center: Point3,
    radius_mm: f64,
    sweep: f64,
    edge: EdgeLabel,
    vertex: VertexLabel,
}

pub(crate) fn rim_geometry(curve: &EdgeCurve3) -> Option<(Point3, f64, f64)> {
    use uom::si::angle::radian;
    match curve {
        EdgeCurve3::Circle(circle) => {
            Some((circle.center(), circle.radius_mm(), core::f64::consts::TAU))
        }
        EdgeCurve3::Arc(arc) => Some((
            arc.center(),
            arc.radius_mm(),
            arc.sweep_angle().get::<radian>().abs(),
        )),
        EdgeCurve3::Line(_) => None,
    }
}

struct RimScan {
    rims: Vec<RimEdge>,
    seam: Option<EdgeLabel>,
}

fn scan_tube_faces(solid: &BrepSolid, label: FaceLabel) -> Option<RimScan> {
    let loop_map: HashMap<_, _> = solid.iter_loops().map(|item| (item.id(), item)).collect();
    let edge_map: HashMap<_, _> = solid.iter_edges().map(|edge| (edge.id(), edge)).collect();
    let vlabel: HashMap<_, _> = solid
        .iter_vertices()
        .map(|vertex| (vertex.id(), vertex.label()))
        .collect();
    let faces: Vec<_> = solid
        .iter_faces()
        .filter(|face| face.label() == label)
        .collect();
    if faces.is_empty() {
        return None;
    }
    let loop_edges = |loop_id: &_| loop_map.get(loop_id).map(|item| item.edges());
    let rims = faces
        .iter()
        .flat_map(|face| face.loops().iter())
        .filter_map(loop_edges)
        .flatten()
        .filter_map(|edge_id| {
            let edge = edge_map.get(edge_id)?;
            let (center, radius_mm, sweep) = rim_geometry(edge.curve())?;
            let vertex = vlabel.get(&edge.vertices()[0]).copied()?;
            Some(RimEdge {
                center,
                radius_mm,
                sweep,
                edge: edge.label(),
                vertex,
            })
        })
        .collect();
    let seam = faces
        .iter()
        .flat_map(|face| face.loops().iter())
        .filter_map(loop_edges)
        .flatten()
        .filter_map(|edge_id| {
            let edge = edge_map.get(edge_id)?;
            matches!(edge.curve(), EdgeCurve3::Line(_)).then(|| edge.label())
        })
        .min();
    Some(RimScan { rims, seam })
}

fn extract_tube_from_bone(
    solid: &BrepSolid,
    label: FaceLabel,
    soup: &classify::SolidSoup,
    tolerance: Tolerance,
) -> Result<Tube, BrepError> {
    let eps = tolerance.value();
    let RimScan { rims, seam } = scan_tube_faces(solid, label).ok_or(BrepError::MissingLabel {
        kind: LabelKind::Face,
    })?;
    let centers: Vec<Point3> = rims.iter().fold(Vec::new(), |mut acc, rim| {
        if !acc.iter().any(|seen| (*seen - rim.center).norm_mm() <= eps) {
            acc.push(rim.center);
        }
        acc
    });
    let [center_a, center_b] = centers.as_slice() else {
        return Err(unsupported(super::BooleanGap::CurvedFace));
    };
    let ring_of = |center: Point3| -> Vec<&RimEdge> {
        rims.iter()
            .filter(|rim| (rim.center - center).norm_mm() <= eps)
            .collect()
    };
    let ring_a = ring_of(*center_a);
    let ring_b = ring_of(*center_b);
    let radius_mm = ring_a.first().map_or(0.0, |rim| rim.radius_mm);
    let consistent = ring_a
        .iter()
        .chain(&ring_b)
        .all(|rim| (rim.radius_mm - radius_mm).abs() <= eps);
    let full = (ring_a.iter().map(|rim| rim.sweep).sum::<f64>() - core::f64::consts::TAU).abs()
        <= 1.0e-3
        && (ring_b.iter().map(|rim| rim.sweep).sum::<f64>() - core::f64::consts::TAU).abs()
            <= 1.0e-3;
    if !consistent || !full {
        return Err(unsupported(super::BooleanGap::CurvedFace));
    }
    let raw = unit(*center_b - *center_a, tolerance).ok_or(BrepError::BooleanIndeterminate)?;
    let axis = model::canonical_axis(raw);
    let (lower, upper) = if model::project_onto_axis(*center_a, axis, *center_a)
        <= model::project_onto_axis(*center_a, axis, *center_b)
    {
        (&ring_a, &ring_b)
    } else {
        (&ring_b, &ring_a)
    };
    let lower_center = lower.first().map_or(*center_a, |rim| rim.center);
    let upper_center = upper.first().map_or(*center_b, |rim| rim.center);
    let height = model::project_onto_axis(lower_center, axis, upper_center);
    if height <= eps {
        return Err(unsupported(super::BooleanGap::CurvedFace));
    }
    let radius = Length::new::<millimeter>(radius_mm);
    let cyl =
        Cyl::new(lower_center, axis, radius, tolerance).ok_or(BrepError::BooleanIndeterminate)?;
    let least_edge = |ring: &[&RimEdge]| ring.iter().map(|rim| rim.edge).min();
    let least_vertex = |ring: &[&RimEdge]| ring.iter().map(|rim| rim.vertex).min();
    let rim_bottom = least_edge(lower).ok_or(BrepError::BooleanIndeterminate)?;
    let rim_top = least_edge(upper).ok_or(BrepError::BooleanIndeterminate)?;
    let vbot = least_vertex(lower).ok_or(BrepError::BooleanIndeterminate)?;
    let vtop = least_vertex(upper).ok_or(BrepError::BooleanIndeterminate)?;
    let seam_label = seam.unwrap_or_else(|| model::elabel_intersection(label, label));
    let radial = cyl
        .radial(0.0, tolerance)
        .ok_or(BrepError::BooleanIndeterminate)?;
    let mid = cyl.point(0.0, height / 2.0);
    let outward = wall_outward(soup, mid, radial, radius_mm);
    Ok(Tube {
        label,
        cyl,
        theta0: 0.0,
        theta1: core::f64::consts::TAU,
        z0: 0.0,
        z1: height,
        full: true,
        rim_bottom,
        rim_top,
        seam: seam_label,
        vbot,
        vtop,
        outward,
    })
}

fn is_circle(curve: &EdgeCurve3) -> bool {
    matches!(curve, EdgeCurve3::Circle(_))
}

fn mergeable(
    last: &(PlanarVertex, EdgeLabel, EdgeCurve3),
    next: &(PlanarVertex, EdgeLabel, EdgeCurve3),
) -> bool {
    last.1 == next.1
        && ((is_line(&last.2) && is_line(&next.2)) || (is_circle(&last.2) && is_circle(&next.2)))
}

fn collapse_by_label(
    parts: &[(PlanarVertex, EdgeLabel, EdgeCurve3)],
) -> Vec<(PlanarVertex, EdgeLabel, EdgeCurve3)> {
    let merged = parts.iter().fold(
        Vec::<(PlanarVertex, EdgeLabel, EdgeCurve3)>::new(),
        |mut acc, part| {
            if !acc.last().is_some_and(|last| mergeable(last, part)) {
                acc.push(part.clone());
            }
            acc
        },
    );
    let wraps = merged.len() > 1
        && merged
            .first()
            .zip(merged.last())
            .is_some_and(|(first, last)| mergeable(last, first));
    if wraps {
        merged[..merged.len() - 1].to_vec()
    } else {
        merged
    }
}

fn raw_loops(
    face: &truck_modeling::Face,
    vertex_labels: &[(Point3, VertexLabel)],
    edges: &EdgeIndex,
    tolerance: Tolerance,
) -> Result<Vec<RawLoop>, BrepError> {
    face.boundaries()
        .iter()
        .map(|wire| extract_loop(wire, vertex_labels, edges, tolerance))
        .collect()
}

fn loop_area(face_loop: &RawLoop) -> f64 {
    let points: Vec<Point3> = face_loop
        .curves
        .iter()
        .flat_map(|curve| {
            use crate::curve3::Curve3;
            curve.tessellate(bone_types::ChordHeightTolerance::from_mm(0.05))
        })
        .collect();
    newell(&points).norm_mm()
}

fn extract_loop(
    wire: &truck_modeling::Wire,
    vertex_labels: &[(Point3, VertexLabel)],
    edges: &EdgeIndex,
    tolerance: Tolerance,
) -> Result<RawLoop, BrepError> {
    let parts = wire
        .edge_iter()
        .map(|edge| {
            let front = point_from_truck(edge.front().point());
            let back = point_from_truck(edge.back().point());
            let (edge_label, curve) =
                edges
                    .get(&edge.id())
                    .cloned()
                    .ok_or(BrepError::MissingLabel {
                        kind: LabelKind::Edge,
                    })?;
            let label = if matches!(curve, EdgeCurve3::Circle(_)) {
                nearest_label(front, vertex_labels)
            } else {
                label_at(front, vertex_labels, tolerance)
            }
            .ok_or(BrepError::MissingLabel {
                kind: LabelKind::Vertex,
            })?;
            let vertex = PlanarVertex {
                position: front,
                label,
            };
            let curve = model::orient_curve(curve, front, back, tolerance);
            Ok::<_, BrepError>((vertex, edge_label, curve))
        })
        .collect::<Result<Vec<_>, _>>()?;
    let collapsed = collapse_by_label(&parts);
    let verts: Vec<PlanarVertex> = collapsed.iter().map(|(v, _, _)| *v).collect();
    let labels: Vec<EdgeLabel> = collapsed.iter().map(|(_, e, _)| *e).collect();
    let curves: Vec<EdgeCurve3> = collapsed.iter().map(|(_, _, c)| c.clone()).collect();
    let (verts, labels, curves) = normalize_circle_loop(verts, labels, curves, tolerance);
    Ok(RawLoop {
        verts,
        edges: labels,
        curves,
    })
}

fn normalize_circle_loop(
    verts: Vec<PlanarVertex>,
    labels: Vec<EdgeLabel>,
    curves: Vec<EdgeCurve3>,
    tolerance: Tolerance,
) -> (Vec<PlanarVertex>, Vec<EdgeLabel>, Vec<EdgeCurve3>) {
    use crate::circle3::Circle3;
    use uom::si::angle::radian;
    let arcs: Option<Vec<(Point3, UnitVec3, f64, f64)>> = curves
        .iter()
        .map(|curve| match curve {
            EdgeCurve3::Arc(arc) => Some((
                arc.center(),
                arc.normal(),
                arc.radius_mm(),
                arc.sweep_angle().get::<radian>().abs(),
            )),
            _ => None,
        })
        .collect();
    let Some(arcs) = arcs.filter(|arcs| arcs.len() >= 2) else {
        return (verts, labels, curves);
    };
    let eps = tolerance.value();
    let (center, normal, radius_mm, _) = arcs[0];
    let consistent = arcs.iter().all(|(c, n, r, _)| {
        (*c - center).norm_mm() <= eps
            && n.dot(normal).abs() >= 1.0 - eps
            && (r - radius_mm).abs() <= eps
    });
    let total: f64 = arcs.iter().map(|(_, _, _, sweep)| sweep).sum();
    if !consistent || (total - core::f64::consts::TAU).abs() > 1.0e-3 {
        return (verts, labels, curves);
    }
    let Some((x, y)) = basis(normal, tolerance) else {
        return (verts, labels, curves);
    };
    let plane = Plane3::new_unchecked(center, x, y);
    let Ok(circle) = Circle3::new(plane, Length::new::<millimeter>(radius_mm), tolerance) else {
        return (verts, labels, curves);
    };
    let vertex = verts.first().copied().unwrap_or(PlanarVertex {
        position: center + as_vec(x) * radius_mm,
        label: VertexLabel {
            feature: labels[0].feature,
            role: bone_types::VertexRole::Cut {
                key: stitch::cut_key(center),
            },
        },
    });
    (
        vec![vertex],
        vec![labels[0]],
        vec![EdgeCurve3::Circle(circle)],
    )
}

fn nearest_label(position: Point3, vertex_labels: &[(Point3, VertexLabel)]) -> Option<VertexLabel> {
    vertex_labels
        .iter()
        .min_by(|(a, _), (b, _)| {
            (*a - position)
                .norm_mm()
                .total_cmp(&(*b - position).norm_mm())
        })
        .map(|(_, label)| *label)
}

fn label_at(
    position: Point3,
    vertex_labels: &[(Point3, VertexLabel)],
    tolerance: Tolerance,
) -> Option<VertexLabel> {
    vertex_labels
        .iter()
        .filter(|(candidate, _)| (*candidate - position).norm_mm() <= tolerance.value())
        .min_by(|(a, _), (b, _)| {
            let da = (*a - position).norm_mm();
            let db = (*b - position).norm_mm();
            da.partial_cmp(&db).unwrap_or(core::cmp::Ordering::Equal)
        })
        .map(|(_, label)| *label)
}

pub(crate) fn plane_from_polygon(
    points: &[Point3],
    tolerance: Tolerance,
) -> Result<Plane3, BrepError> {
    let normal = unit(newell(points), tolerance).ok_or(BrepError::BooleanIndeterminate)?;
    let (x, y) = basis(normal, tolerance).ok_or(BrepError::BooleanIndeterminate)?;
    let origin = points
        .first()
        .copied()
        .ok_or(BrepError::BooleanIndeterminate)?;
    Plane3::new(origin, x, y, tolerance).map_err(|_| BrepError::BooleanIndeterminate)
}

fn unsupported(detail: super::BooleanGap) -> BrepError {
    BrepError::BooleanUnsupported { detail }
}

fn reverse_face(face: &Face) -> Face {
    match face {
        Face::Planar(planar) => Face::planar(reverse_planar(planar)),
        Face::Tube(tube) => {
            let mut reversed = tube.clone();
            reversed.outward = !reversed.outward;
            std::mem::swap(&mut reversed.rim_bottom, &mut reversed.rim_top);
            std::mem::swap(&mut reversed.vbot, &mut reversed.vtop);
            Face::Tube(reversed)
        }
    }
}

fn reverse_loop(item: &PlanarLoop) -> PlanarLoop {
    let count = item.edges.len();
    PlanarLoop {
        verts: item.verts.iter().rev().copied().collect(),
        edges: (0..count)
            .map(|i| item.edges[(count + count - 2 - i) % count])
            .collect(),
        curves: (0..count)
            .map(|i| model::reverse_curve(&item.curves[(count + count - 2 - i) % count]))
            .collect(),
    }
}

fn reverse_planar(face: &PlanarFace) -> PlanarFace {
    let flipped = Plane3::new_unchecked(
        face.plane.origin(),
        face.plane.y_axis(),
        face.plane.x_axis(),
    );
    let polygon = face.polygon.iter().rev().copied().collect();
    let count = face.edges.len();
    let edges = (0..count)
        .map(|i| face.edges[(count + count - 2 - i) % count])
        .collect();
    let curves = (0..count)
        .map(|i| model::reverse_curve(&face.curves[(count + count - 2 - i) % count]))
        .collect();
    let holes = face.holes.iter().map(reverse_loop).collect();
    PlanarFace {
        label: face.label,
        plane: flipped,
        polygon,
        edges,
        curves,
        holes,
    }
}

fn classify_all(
    body: &Body,
    other: &Body,
    soup: &classify::SolidSoup,
    tolerance: Tolerance,
) -> Result<Vec<FaceClass>, BrepError> {
    body.iter()
        .map(|face| classify::classify_face(face, other, soup, tolerance))
        .collect::<Option<_>>()
        .ok_or(BrepError::BooleanIndeterminate)
}

fn select(
    op: BooleanOp,
    body_a: &Body,
    body_b: &Body,
    soup_a: &classify::SolidSoup,
    soup_b: &classify::SolidSoup,
    tolerance: Tolerance,
) -> Result<Vec<Face>, BrepError> {
    let classes_a = classify_all(body_a, body_b, soup_b, tolerance)?;
    let classes_b = classify_all(body_b, body_a, soup_a, tolerance)?;
    let mut kept = match op {
        BooleanOp::Union => select_union(body_a, &classes_a, body_b, &classes_b)?,
        BooleanOp::Intersection => select_intersection(body_a, &classes_a, body_b, &classes_b),
        BooleanOp::Difference => select_difference(body_a, &classes_a, body_b, &classes_b)?,
    };
    kept.sort_by_key(Face::label);
    Ok(kept)
}

fn keep_where(body: &Body, classes: &[FaceClass], keep: impl Fn(FaceClass) -> bool) -> Vec<Face> {
    body.iter()
        .zip(classes)
        .filter(|(_, class)| keep(**class))
        .map(|(face, _)| face.clone())
        .collect()
}

fn select_union(
    body_a: &Body,
    classes_a: &[FaceClass],
    body_b: &Body,
    classes_b: &[FaceClass],
) -> Result<Vec<Face>, BrepError> {
    let touching = classes_a
        .iter()
        .chain(classes_b)
        .any(|class| *class != FaceClass::Out);
    if !touching {
        return Err(BrepError::BooleanDisjoint);
    }
    let from_a = keep_where(body_a, classes_a, |class| {
        matches!(class, FaceClass::Out | FaceClass::OnSame)
    });
    let from_b = keep_where(body_b, classes_b, |class| class == FaceClass::Out);
    Ok(from_a.into_iter().chain(from_b).collect())
}

fn select_intersection(
    body_a: &Body,
    classes_a: &[FaceClass],
    body_b: &Body,
    classes_b: &[FaceClass],
) -> Vec<Face> {
    let from_a = keep_where(body_a, classes_a, |class| {
        matches!(class, FaceClass::In | FaceClass::OnSame)
    });
    let from_b = keep_where(body_b, classes_b, |class| class == FaceClass::In);
    from_a.into_iter().chain(from_b).collect()
}

fn select_difference(
    body_a: &Body,
    classes_a: &[FaceClass],
    body_b: &Body,
    classes_b: &[FaceClass],
) -> Result<Vec<Face>, BrepError> {
    let touching = classes_a
        .iter()
        .chain(classes_b)
        .any(|class| *class != FaceClass::Out);
    if !touching {
        return Err(BrepError::BooleanDisjoint);
    }
    let from_a = keep_where(body_a, classes_a, |class| {
        matches!(class, FaceClass::Out | FaceClass::OnOpposite)
    });
    let from_b = body_b
        .iter()
        .zip(classes_b)
        .filter(|(_, class)| **class == FaceClass::In)
        .map(|(face, _)| reverse_face(face));
    Ok(from_a.into_iter().chain(from_b).collect())
}
