use bone_types::{EdgeLabel, EdgeRole, FaceLabel, Point3, Tolerance, VertexLabel, VertexRole};

use super::super::edges::{EdgeCurve3, line_between};
use super::classify;
use super::classify::local;
use super::{Body, Face, PlanarBody, PlanarFace, PlanarLoop, PlanarVertex, Tube, as_vec};

fn is_straight_simple(face: &PlanarFace) -> bool {
    face.holes.is_empty() && face.polygon.len() >= 3 && face.curves.iter().all(is_line)
}

fn is_disk(face: &PlanarFace) -> bool {
    face.holes.is_empty()
        && face.curves.len() == 1
        && matches!(face.curves.first(), Some(EdgeCurve3::Circle(_)))
}

pub(crate) fn arrange_body(
    body: &Body,
    other: &Body,
    tolerance: Tolerance,
) -> Result<Body, super::super::BrepError> {
    let other_planar: PlanarBody = other
        .iter()
        .filter_map(|face| match face {
            Face::Planar(planar) => Some((**planar).clone()),
            Face::Tube(_) => None,
        })
        .collect();
    let other_tubes: Vec<Tube> = other
        .iter()
        .filter_map(|face| match face {
            Face::Tube(tube) => Some((**tube).clone()),
            Face::Planar(_) => None,
        })
        .collect();
    body.iter()
        .map(|face| match face {
            Face::Planar(planar) => {
                let lines = planar_lines(planar, &other_planar, &other_tubes, tolerance);
                let bores = tube_holes(planar, &other_tubes, tolerance)?;
                let (outer_bores, extras) = split_washers(planar, &bores, tolerance);
                let mut faces: Vec<Face> = arrange_planar(planar, &lines, &outer_bores, tolerance)?
                    .into_iter()
                    .map(Face::planar)
                    .collect();
                faces.extend(extras.into_iter().map(Face::planar));
                Ok(faces)
            }
            Face::Tube(tube) => arrange_tube(tube, other, tolerance),
        })
        .collect::<Result<Vec<Vec<Face>>, _>>()
        .map(|nested| nested.into_iter().flatten().collect())
}

fn arrange_planar(
    face: &PlanarFace,
    lines: &[ChordLine],
    bores: &[PlanarLoop],
    tolerance: Tolerance,
) -> Result<Vec<PlanarFace>, super::super::BrepError> {
    if is_disk(face) {
        return arrange_disk_bored(face, lines, bores, tolerance);
    }
    if is_holed_polygon(face) || (is_straight_simple(face) && !bores.is_empty()) {
        return arrange_holed(face, lines, bores, tolerance);
    }
    if is_straight_simple(face) {
        return Ok(split_with_lines(face, lines, tolerance));
    }
    Ok(vec![face.clone()])
}

fn arrange_disk_bored(
    face: &PlanarFace,
    lines: &[ChordLine],
    bores: &[PlanarLoop],
    tolerance: Tolerance,
) -> Result<Vec<PlanarFace>, super::super::BrepError> {
    if bores.is_empty() {
        return arrange_disk(face, lines, tolerance);
    }
    if !lines.is_empty() {
        return Err(super::super::BrepError::BooleanUnsupported {
            detail: super::super::BooleanGap::PartialArcProfile,
        });
    }
    let mut annulus = face.clone();
    annulus.holes = bores.to_vec();
    Ok(vec![annulus])
}

enum HoleHit {
    Bore(PlanarLoop),
    Partial,
    Clear,
}

fn tube_holes(
    face: &PlanarFace,
    tubes: &[Tube],
    tolerance: Tolerance,
) -> Result<Vec<PlanarLoop>, super::super::BrepError> {
    use super::ssi::{IntersectionCurve, SurfaceIntersection, intersect_surfaces};
    let eps = tolerance.value();
    let hits: Vec<HoleHit> = tubes
        .iter()
        .flat_map(|tube| {
            match intersect_surfaces(face.analytic(), tube.cyl.analytic(), tolerance) {
                SurfaceIntersection::Curves(curves) => curves
                    .into_iter()
                    .map(|curve| match curve {
                        IntersectionCurve::Circle(circle) => {
                            let z = super::model::project_onto_axis(
                                tube.cyl.axis_origin,
                                tube.cyl.axis,
                                circle.center(),
                            );
                            if z < tube.z0 - eps || z > tube.z1 + eps {
                                return HoleHit::Clear;
                            }
                            match circle_coverage(face, &circle, tolerance) {
                                Coverage::Full => HoleHit::Bore(circle_hole(face, tube, &circle)),
                                Coverage::Partial => HoleHit::Partial,
                                Coverage::None => HoleHit::Clear,
                            }
                        }
                        IntersectionCurve::Line(_) => HoleHit::Clear,
                    })
                    .collect::<Vec<_>>(),
                _ => Vec::new(),
            }
        })
        .collect();
    if hits.iter().any(|hit| matches!(hit, HoleHit::Partial)) {
        return Err(super::super::BrepError::BooleanUnsupported {
            detail: super::super::BooleanGap::PartialCircleImprint,
        });
    }
    Ok(hits
        .into_iter()
        .filter_map(|hit| match hit {
            HoleHit::Bore(loop_) => Some(loop_),
            _ => None,
        })
        .collect())
}

fn circle_coverage(
    face: &PlanarFace,
    circle: &crate::circle3::Circle3,
    tolerance: Tolerance,
) -> Coverage {
    let Some((rx, ry)) = super::model::canonical_basis(circle.normal(), tolerance) else {
        return Coverage::None;
    };
    let center = circle.center();
    let radius = circle.radius_mm();
    let inside = (0..COVERAGE_SAMPLES)
        .filter(|index| {
            let theta = core::f64::consts::TAU * f64::from(*index) / f64::from(COVERAGE_SAMPLES);
            let point =
                center + as_vec(rx) * (radius * theta.cos()) + as_vec(ry) * (radius * theta.sin());
            classify::region_interior(face, point, tolerance)
        })
        .count();
    coverage_of(inside, usize::from(COVERAGE_SAMPLES))
}

fn circle_hole(face: &PlanarFace, tube: &Tube, circle: &crate::circle3::Circle3) -> PlanarLoop {
    PlanarLoop {
        verts: vec![PlanarVertex {
            position: circle.center(),
            label: bone_types::VertexLabel {
                feature: face.label.feature,
                role: bone_types::VertexRole::Cut {
                    key: super::stitch::cut_key(circle.center()),
                },
            },
        }],
        edges: vec![super::model::elabel_intersection(face.label, tube.label)],
        curves: vec![EdgeCurve3::Circle(*circle)],
    }
}

fn bore_circle(loop_: &PlanarLoop) -> Option<&crate::circle3::Circle3> {
    match loop_.curves.first() {
        Some(EdgeCurve3::Circle(circle)) => Some(circle),
        _ => None,
    }
}

fn coaxial(a: &crate::circle3::Circle3, b: &crate::circle3::Circle3, eps: f64) -> bool {
    (a.center() - b.center()).norm_mm() <= eps && a.normal().dot(b.normal()).abs() >= 1.0 - eps
}

fn split_washers(
    face: &PlanarFace,
    bores: &[PlanarLoop],
    tolerance: Tolerance,
) -> (Vec<PlanarLoop>, Vec<PlanarFace>) {
    let eps = tolerance.value();
    let circles: Vec<Option<crate::circle3::Circle3>> = bores
        .iter()
        .map(|bore| bore_circle(bore).copied())
        .collect();
    let groups: Vec<Vec<usize>> = (0..bores.len()).fold(Vec::new(), |mut groups, index| {
        let slot = circles[index].and_then(|here| {
            groups.iter_mut().find(|group| {
                group
                    .iter()
                    .any(|&other| circles[other].is_some_and(|seen| coaxial(&here, &seen, eps)))
            })
        });
        match slot {
            Some(group) => group.push(index),
            None => groups.push(vec![index]),
        }
        groups
    });
    groups.iter().fold(
        (Vec::new(), Vec::new()),
        |(mut holes, mut extras), group| {
            let mut ordered = group.clone();
            ordered.sort_by(|&a, &b| {
                let ra = circles[a].map_or(0.0, crate::circle3::Circle3::radius_mm);
                let rb = circles[b].map_or(0.0, crate::circle3::Circle3::radius_mm);
                ra.total_cmp(&rb)
            });
            match ordered.as_slice() {
                [inner, outer] if circles[*inner].is_some() && circles[*outer].is_some() => {
                    holes.push(bores[*outer].clone());
                    let island = &bores[*inner];
                    extras.push(PlanarFace {
                        label: face.label,
                        plane: face.plane,
                        polygon: island.verts.clone(),
                        edges: island.edges.clone(),
                        curves: island.curves.clone(),
                        holes: Vec::new(),
                    });
                    extras.push(PlanarFace {
                        label: face.label,
                        plane: face.plane,
                        polygon: bores[*outer].verts.clone(),
                        edges: bores[*outer].edges.clone(),
                        curves: bores[*outer].curves.clone(),
                        holes: vec![bores[*inner].clone()],
                    });
                }
                _ => ordered
                    .iter()
                    .for_each(|&index| holes.push(bores[index].clone())),
            }
            (holes, extras)
        },
    )
}

fn planar_lines(
    face: &PlanarFace,
    other_planar: &PlanarBody,
    other_tubes: &[Tube],
    tolerance: Tolerance,
) -> Vec<ChordLine> {
    let mut lines = chord_lines(face, other_planar, tolerance);
    lines.extend(tube_chord_lines(face, other_tubes, tolerance));
    dedup_lines(lines, tolerance)
}

fn split_with_lines(
    face: &PlanarFace,
    lines: &[ChordLine],
    tolerance: Tolerance,
) -> Vec<PlanarFace> {
    if lines.is_empty() {
        vec![face.clone()]
    } else {
        split_face(face, lines, tolerance)
    }
}

fn tube_chord_lines(face: &PlanarFace, tubes: &[Tube], tolerance: Tolerance) -> Vec<ChordLine> {
    use super::ssi::{IntersectionCurve, SurfaceIntersection, intersect_surfaces};
    let eps = tolerance.value();
    let x_axis = as_vec(face.plane.x_axis());
    let y_axis = as_vec(face.plane.y_axis());
    tubes
        .iter()
        .flat_map(|tube| {
            match intersect_surfaces(face.analytic(), tube.cyl.analytic(), tolerance) {
                SurfaceIntersection::Curves(curves) => curves
                    .iter()
                    .filter_map(|curve| match curve {
                        IntersectionCurve::Line(line) => {
                            let origin = local(face, line.point());
                            let dir3 = as_vec(line.direction());
                            let dir = (dir3.dot_mm2(x_axis), dir3.dot_mm2(y_axis));
                            let normal = (-dir.1, dir.0);
                            let scale = normal.0.hypot(normal.1);
                            if scale <= eps {
                                return None;
                            }
                            let (a, b) = (normal.0 / scale, normal.1 / scale);
                            let c = -(a * origin.0 + b * origin.1);
                            let (a, b, c) = if a < -eps || (a.abs() <= eps && b < 0.0) {
                                (-a, -b, -c)
                            } else {
                                (a, b, c)
                            };
                            Some(ChordLine {
                                a,
                                b,
                                c,
                                cut: tube.label,
                            })
                        }
                        IntersectionCurve::Circle(_) => None,
                    })
                    .collect::<Vec<_>>(),
                _ => Vec::new(),
            }
        })
        .collect()
}

fn split_range(lo: f64, hi: f64, cuts: &[f64], eps: f64) -> Vec<(f64, f64)> {
    let mut points = vec![lo];
    cuts.iter()
        .filter(|c| **c > lo + eps && **c < hi - eps)
        .for_each(|c| points.push(*c));
    points.push(hi);
    points.sort_by(f64::total_cmp);
    let unique = points
        .into_iter()
        .fold(Vec::<f64>::new(), |mut acc, value| {
            if acc.last().is_none_or(|last| (value - last).abs() > eps) {
                acc.push(value);
            }
            acc
        });
    unique.windows(2).map(|pair| (pair[0], pair[1])).collect()
}

fn split_circle(t0: f64, t1: f64, cuts: &[f64], eps: f64) -> Vec<(f64, f64)> {
    let period = t1 - t0;
    let mut sorted: Vec<f64> = cuts
        .iter()
        .map(|c| t0 + (c - t0).rem_euclid(period))
        .collect();
    sorted.sort_by(f64::total_cmp);
    let unique = sorted
        .into_iter()
        .fold(Vec::<f64>::new(), |mut acc, value| {
            if acc.last().is_none_or(|last| (value - last).abs() > eps) {
                acc.push(value);
            }
            acc
        });
    if unique.len() < 2 {
        return vec![(t0, t1)];
    }
    let count = unique.len();
    (0..count)
        .map(|i| {
            let a = unique[i];
            let b = if i + 1 < count {
                unique[i + 1]
            } else {
                unique[0] + period
            };
            (a, b)
        })
        .collect()
}

enum Coverage {
    Full,
    Partial,
    None,
}

enum TubeCut {
    Z(f64),
    Theta(f64),
    Partial,
}

const COVERAGE_SAMPLES: u16 = 16;

fn coverage_of(inside: usize, total: usize) -> Coverage {
    if inside == 0 {
        Coverage::None
    } else if inside == total {
        Coverage::Full
    } else {
        Coverage::Partial
    }
}

fn ring_coverage(tube: &Tube, planar: &PlanarFace, z: f64, tolerance: Tolerance) -> Coverage {
    let span = tube.theta1 - tube.theta0;
    let inside = (0..COVERAGE_SAMPLES)
        .filter(|index| {
            let theta = tube.theta0 + span * f64::from(*index) / f64::from(COVERAGE_SAMPLES);
            classify::region_contains(planar, tube.cyl.point(theta, z), tolerance)
        })
        .count();
    coverage_of(inside, usize::from(COVERAGE_SAMPLES))
}

fn generator_coverage(
    tube: &Tube,
    planar: &PlanarFace,
    theta: f64,
    tolerance: Tolerance,
) -> Coverage {
    let span = tube.z1 - tube.z0;
    let inside = (0..=COVERAGE_SAMPLES)
        .filter(|index| {
            let z = tube.z0 + span * f64::from(*index) / f64::from(COVERAGE_SAMPLES);
            classify::region_contains(planar, tube.cyl.point(theta, z), tolerance)
        })
        .count();
    coverage_of(inside, usize::from(COVERAGE_SAMPLES) + 1)
}

fn tube_cut(
    tube: &Tube,
    planar: &PlanarFace,
    curve: &super::ssi::IntersectionCurve,
    tolerance: Tolerance,
) -> Option<TubeCut> {
    use super::ssi::IntersectionCurve;
    let eps = tolerance.value();
    match curve {
        IntersectionCurve::Circle(circle) => {
            let z = super::model::project_onto_axis(
                tube.cyl.axis_origin,
                tube.cyl.axis,
                circle.center(),
            );
            if z <= tube.z0 + eps || z >= tube.z1 - eps {
                return None;
            }
            match ring_coverage(tube, planar, z, tolerance) {
                Coverage::Full => Some(TubeCut::Z(z)),
                Coverage::Partial => Some(TubeCut::Partial),
                Coverage::None => None,
            }
        }
        IntersectionCurve::Line(line) => {
            let (theta, _) = tube.cyl.theta_z(line.point());
            match generator_coverage(tube, planar, theta, tolerance) {
                Coverage::Full => Some(TubeCut::Theta(theta)),
                Coverage::Partial => Some(TubeCut::Partial),
                Coverage::None => None,
            }
        }
    }
}

fn arrange_tube(
    tube: &Tube,
    other: &Body,
    tolerance: Tolerance,
) -> Result<Vec<Face>, super::super::BrepError> {
    use super::ssi::{SurfaceIntersection, intersect_surfaces};
    let eps = tolerance.value();
    let cuts: Vec<TubeCut> = other
        .iter()
        .filter_map(|face| match face {
            Face::Planar(planar) => Some(planar),
            Face::Tube(_) => None,
        })
        .flat_map(|planar| {
            match intersect_surfaces(tube.cyl.analytic(), planar.analytic(), tolerance) {
                SurfaceIntersection::Curves(curves) => curves
                    .iter()
                    .filter_map(|curve| tube_cut(tube, planar, curve, tolerance))
                    .collect::<Vec<_>>(),
                _ => Vec::new(),
            }
        })
        .collect();
    if cuts.iter().any(|cut| matches!(cut, TubeCut::Partial)) {
        return Err(super::super::BrepError::BooleanUnsupported {
            detail: super::super::BooleanGap::PartialArcProfile,
        });
    }
    let z_cuts: Vec<f64> = cuts
        .iter()
        .filter_map(|cut| match cut {
            TubeCut::Z(z) => Some(*z),
            _ => None,
        })
        .collect();
    let theta_cuts: Vec<f64> = cuts
        .iter()
        .filter_map(|cut| match cut {
            TubeCut::Theta(theta) => Some(*theta),
            _ => None,
        })
        .collect();
    let z_bands = split_range(tube.z0, tube.z1, &z_cuts, eps);
    let theta_bands = split_circle(tube.theta0, tube.theta1, &theta_cuts, eps);
    Ok(z_bands
        .iter()
        .flat_map(|&(za, zb)| {
            theta_bands.iter().map(move |&(ta, tb)| {
                let mut piece = tube.clone();
                piece.z0 = za;
                piece.z1 = zb;
                piece.theta0 = ta;
                piece.theta1 = tb;
                piece.full = (tb - ta - core::f64::consts::TAU).abs() <= eps;
                if (za - tube.z0).abs() > eps {
                    piece.rim_bottom = super::model::elabel_intersection(tube.label, tube.label);
                }
                if (zb - tube.z1).abs() > eps {
                    piece.rim_top = super::model::elabel_intersection(tube.label, tube.label);
                }
                Face::tube(piece)
            })
        })
        .collect())
}

fn is_holed_polygon(face: &PlanarFace) -> bool {
    !face.holes.is_empty() && face.polygon.len() >= 3 && face.curves.iter().all(is_line)
}

fn outer_only(face: &PlanarFace) -> PlanarFace {
    PlanarFace::straight(
        face.label,
        face.plane,
        face.polygon.clone(),
        face.edges.clone(),
    )
}

fn point_in_polygon(locals: &[(f64, f64)], p: (f64, f64)) -> bool {
    let count = locals.len();
    (0..count).fold(false, |inside, index| {
        let (xi, yi) = locals[index];
        let (xj, yj) = locals[(index + count - 1) % count];
        let crosses = (yi > p.1) != (yj > p.1) && p.0 < (xj - xi) * (p.1 - yi) / (yj - yi) + xi;
        inside ^ crosses
    })
}

fn arrange_holed(
    face: &PlanarFace,
    lines: &[ChordLine],
    bores: &[PlanarLoop],
    tolerance: Tolerance,
) -> Result<Vec<PlanarFace>, super::super::BrepError> {
    let eps = tolerance.value();
    let all_holes: Vec<PlanarLoop> = face.holes.iter().chain(bores).cloned().collect();
    let base = outer_only(face);
    let crosses_hole = all_holes.iter().any(|hole| match hole.curves.first() {
        Some(EdgeCurve3::Circle(circle)) => {
            let center = local(face, circle.center());
            lines
                .iter()
                .any(|line| line.side(center).abs() < circle.radius_mm() - eps)
        }
        _ => false,
    });
    if crosses_hole {
        return Err(super::super::BrepError::BooleanUnsupported {
            detail: super::super::BooleanGap::PartialCircleImprint,
        });
    }
    let cells = split_with_lines(&base, lines, tolerance);
    Ok(cells
        .into_iter()
        .map(|mut cell| {
            let locals: Vec<(f64, f64)> = cell
                .polygon
                .iter()
                .map(|v| local(&cell, v.position))
                .collect();
            cell.holes = all_holes
                .iter()
                .filter(|hole| match hole.curves.first() {
                    Some(EdgeCurve3::Circle(circle)) => {
                        point_in_polygon(&locals, local(&cell, circle.center()))
                    }
                    _ => hole
                        .verts
                        .first()
                        .is_some_and(|v| point_in_polygon(&locals, local(&cell, v.position))),
                })
                .cloned()
                .collect();
            cell
        })
        .collect())
}

fn arrange_disk(
    face: &PlanarFace,
    lines: &[ChordLine],
    tolerance: Tolerance,
) -> Result<Vec<PlanarFace>, super::super::BrepError> {
    let Some(EdgeCurve3::Circle(circle)) = face.curves.first() else {
        return Ok(vec![face.clone()]);
    };
    let circle = *circle;
    let eps = tolerance.value();
    let radius = circle.radius_mm();
    let center_local = local(face, circle.center());
    let crossing: Vec<&ChordLine> = lines
        .iter()
        .filter(|line| line.side(center_local).abs() < radius - eps)
        .collect();
    match crossing.as_slice() {
        [] => Ok(vec![face.clone()]),
        [line] => Ok(split_disk(face, &circle, line, center_local, tolerance)),
        _ => Err(super::super::BrepError::BooleanUnsupported {
            detail: super::super::BooleanGap::PartialArcProfile,
        }),
    }
}

fn split_disk(
    face: &PlanarFace,
    circle: &crate::circle3::Circle3,
    line: &ChordLine,
    center_local: (f64, f64),
    tolerance: Tolerance,
) -> Vec<PlanarFace> {
    use bone_types::Angle;
    use uom::si::angle::radian;
    let radius = circle.radius_mm();
    let side = line.side(center_local);
    let foot = (
        center_local.0 - side * line.a,
        center_local.1 - side * line.b,
    );
    let half = (radius * radius - side * side).max(0.0).sqrt();
    let dir = (-line.b, line.a);
    let p1l = (foot.0 + half * dir.0, foot.1 + half * dir.1);
    let p2l = (foot.0 - half * dir.0, foot.1 - half * dir.1);
    let p1 = face.plane.point_at_local(p1l.0, p1l.1, 0.0);
    let p2 = face.plane.point_at_local(p2l.0, p2l.1, 0.0);
    let Some((rx, ry)) = super::model::canonical_basis(circle.normal(), tolerance) else {
        return vec![face.clone()];
    };
    let center = circle.center();
    let angle_of = |p: Point3| {
        let d = p - center;
        d.dot_mm2(as_vec(ry)).atan2(d.dot_mm2(as_vec(rx)))
    };
    let theta1 = angle_of(p1);
    let theta2 = angle_of(p2);
    let sweep_a = (theta2 - theta1).rem_euclid(core::f64::consts::TAU);
    let plane = bone_types::Plane3::new_unchecked(center, rx, ry);
    let arc = |start: f64, sweep: f64| {
        crate::arc3::Arc3::new(
            plane,
            circle.radius(),
            Angle::new::<radian>(start),
            Angle::new::<radian>(sweep),
            tolerance,
        )
        .ok()
        .map(EdgeCurve3::Arc)
    };
    let (Some(arc_a), Some(arc_b)) = (
        arc(theta1, sweep_a),
        arc(theta2, core::f64::consts::TAU - sweep_a),
    ) else {
        return vec![face.clone()];
    };
    let vlabel = super::model::vlabel_intersection([face.label, line.cut, line.cut]);
    let rim_label = face.edges[0];
    let chord_label = super::model::elabel_intersection(face.label, line.cut);
    let pv = |position| PlanarVertex {
        position,
        label: vlabel,
    };
    let segment = |verts: Vec<PlanarVertex>, curves: Vec<EdgeCurve3>| PlanarFace {
        label: face.label,
        plane: face.plane,
        polygon: verts,
        edges: vec![rim_label, chord_label],
        curves,
        holes: Vec::new(),
    };
    vec![
        segment(vec![pv(p1), pv(p2)], vec![arc_a, line_between(p2, p1)]),
        segment(vec![pv(p2), pv(p1)], vec![arc_b, line_between(p1, p2)]),
    ]
}

pub(crate) fn conform_body(body: &Body, tolerance: Tolerance) -> Body {
    let planar: Vec<PlanarFace> = body
        .iter()
        .filter_map(|face| match face {
            Face::Planar(planar) => Some((**planar).clone()),
            Face::Tube(_) => None,
        })
        .collect();
    let conformed = conform(&planar, tolerance);
    let mut next = conformed.into_iter();
    body.iter()
        .map(|face| match face {
            Face::Planar(_) => next.next().map_or_else(|| face.clone(), Face::planar),
            Face::Tube(tube) => Face::Tube(tube.clone()),
        })
        .collect()
}

#[derive(Copy, Clone)]
struct ChordLine {
    a: f64,
    b: f64,
    c: f64,
    cut: FaceLabel,
}

impl ChordLine {
    fn side(&self, point: (f64, f64)) -> f64 {
        self.a * point.0 + self.b * point.1 + self.c
    }
}

fn chord_lines(face: &PlanarFace, other: &PlanarBody, tolerance: Tolerance) -> Vec<ChordLine> {
    let found: Vec<ChordLine> = other
        .iter()
        .filter_map(|peer| face_imprint(face, peer, tolerance))
        .collect();
    dedup_lines(found, tolerance)
}

fn dedup_lines(mut lines: Vec<ChordLine>, tolerance: Tolerance) -> Vec<ChordLine> {
    let eps = tolerance.value();
    lines.sort_by(|a, b| {
        a.a.total_cmp(&b.a)
            .then(a.b.total_cmp(&b.b))
            .then(a.c.total_cmp(&b.c))
            .then(a.cut.cmp(&b.cut))
    });
    lines.into_iter().fold(Vec::new(), |mut kept, line| {
        let duplicate = kept.iter().any(|seen: &ChordLine| {
            (seen.a - line.a).abs() <= eps
                && (seen.b - line.b).abs() <= eps
                && (seen.c - line.c).abs() <= eps
        });
        if !duplicate {
            kept.push(line);
        }
        kept
    })
}

fn region_span(
    face: &PlanarFace,
    origin: (f64, f64),
    dir: (f64, f64),
    eps: f64,
) -> Option<(f64, f64)> {
    if let Some(EdgeCurve3::Circle(circle)) = face.curves.first()
        && face.polygon.len() == 1
    {
        let center = local(face, circle.center());
        let r = circle.radius_mm();
        let fx = origin.0 - center.0;
        let fy = origin.1 - center.1;
        let qa = dir.0 * dir.0 + dir.1 * dir.1;
        let qb = 2.0 * (fx * dir.0 + fy * dir.1);
        let qc = fx * fx + fy * fy - r * r;
        let disc = qb * qb - 4.0 * qa * qc;
        if disc <= 0.0 || qa <= eps {
            return None;
        }
        let root = disc.sqrt();
        let t0 = (-qb - root) / (2.0 * qa);
        let t1 = (-qb + root) / (2.0 * qa);
        return Some((t0.min(t1), t0.max(t1)));
    }
    let locals: Vec<(f64, f64)> = face
        .polygon
        .iter()
        .map(|v| local(face, v.position))
        .collect();
    clip_to_convex(&locals, origin, dir, eps)
}

fn face_imprint(face: &PlanarFace, peer: &PlanarFace, tolerance: Tolerance) -> Option<ChordLine> {
    let eps = tolerance.value();
    let peer_normal = as_vec(peer.plane.normal());
    let x_axis = as_vec(face.plane.x_axis());
    let y_axis = as_vec(face.plane.y_axis());
    let a = x_axis.dot_mm2(peer_normal);
    let b = y_axis.dot_mm2(peer_normal);
    let c = (face.plane.origin() - peer.plane.origin()).dot_mm2(peer_normal);
    let scale = a.hypot(b);
    if scale <= eps {
        return None;
    }

    let (a, b, c) = if a < -eps || (a.abs() <= eps && b < 0.0) {
        (-a / scale, -b / scale, -c / scale)
    } else {
        (a / scale, b / scale, c / scale)
    };
    let dir = (-b, a);
    let point = if a.abs() >= b.abs() {
        (-c / a, 0.0)
    } else {
        (0.0, -c / b)
    };

    let face_span = region_span(face, point, dir, eps)?;

    let anchor = face.plane.point_at_local(point.0, point.1, 0.0);
    let heading = x_axis * dir.0 + y_axis * dir.1;
    let peer_origin = local(peer, anchor);
    let peer_dir = (
        heading.dot_mm2(as_vec(peer.plane.x_axis())),
        heading.dot_mm2(as_vec(peer.plane.y_axis())),
    );
    let peer_span = region_span(peer, peer_origin, peer_dir, eps)?;

    let lo = face_span.0.max(peer_span.0);
    let hi = face_span.1.min(peer_span.1);
    (hi - lo > eps).then_some(ChordLine {
        a,
        b,
        c,
        cut: peer.label,
    })
}

fn clip_to_convex(
    poly: &[(f64, f64)],
    origin: (f64, f64),
    dir: (f64, f64),
    eps: f64,
) -> Option<(f64, f64)> {
    let area = signed_area(poly);
    if area.abs() <= eps * eps {
        return None;
    }
    let orient = area.signum();
    let count = poly.len();
    let bounds = (0..count).try_fold((f64::NEG_INFINITY, f64::INFINITY), |(lo, hi), index| {
        let tail = poly[index];
        let head = poly[(index + 1) % count];
        let edge = (head.0 - tail.0, head.1 - tail.1);
        let inward = (-edge.1 * orient, edge.0 * orient);
        let inward_len = inward.0.hypot(inward.1);
        let c0 = (origin.0 - tail.0) * inward.0 + (origin.1 - tail.1) * inward.1;
        let c1 = dir.0 * inward.0 + dir.1 * inward.1;
        if c1.abs() <= eps * inward_len {
            (c0 >= -eps * inward_len).then_some((lo, hi))
        } else if c1 > 0.0 {
            Some((lo.max(-c0 / c1), hi))
        } else {
            Some((lo, hi.min(-c0 / c1)))
        }
    });
    bounds.filter(|(lo, hi)| hi - lo > eps)
}

fn signed_area(poly: &[(f64, f64)]) -> f64 {
    let count = poly.len();
    (0..count)
        .map(|index| {
            let (x0, y0) = poly[index];
            let (x1, y1) = poly[(index + 1) % count];
            x0 * y1 - x1 * y0
        })
        .sum::<f64>()
        / 2.0
}

#[derive(Clone)]
struct CellVertex {
    local: (f64, f64),
    position: Point3,
    label: VertexLabel,
}

#[derive(Clone)]
struct Cell {
    verts: Vec<CellVertex>,
    edges: Vec<EdgeLabel>,
}

fn split_face(face: &PlanarFace, lines: &[ChordLine], tolerance: Tolerance) -> Vec<PlanarFace> {
    let initial = Cell {
        verts: face
            .polygon
            .iter()
            .map(|vertex| CellVertex {
                local: local(face, vertex.position),
                position: vertex.position,
                label: vertex.label,
            })
            .collect(),
        edges: face.edges.clone(),
    };
    let cells = lines.iter().fold(vec![initial], |cells, line| {
        cells
            .into_iter()
            .flat_map(|cell| split_cell(&cell, line, face, tolerance))
            .collect()
    });
    cells
        .into_iter()
        .map(|cell| {
            PlanarFace::straight(
                face.label,
                face.plane,
                cell.verts
                    .iter()
                    .map(|vertex| PlanarVertex {
                        position: vertex.position,
                        label: vertex.label,
                    })
                    .collect(),
                cell.edges,
            )
        })
        .collect()
}

struct RingEntry {
    vertex: CellVertex,
    edge: EdgeLabel,
    on_chord: bool,
}

fn split_cell(cell: &Cell, line: &ChordLine, face: &PlanarFace, tolerance: Tolerance) -> Vec<Cell> {
    let eps = tolerance.value();
    let sides: Vec<f64> = cell
        .verts
        .iter()
        .map(|vertex| line.side(vertex.local))
        .collect();
    let has_pos = sides.iter().any(|side| *side > eps);
    let has_neg = sides.iter().any(|side| *side < -eps);
    if !(has_pos && has_neg) {
        return vec![cell.clone()];
    }

    let chord = chord_label(face.label, line.cut);
    let count = cell.verts.len();
    let mut positive: Vec<RingEntry> = Vec::new();
    let mut negative: Vec<RingEntry> = Vec::new();
    (0..count).for_each(|index| {
        let vertex = cell.verts[index].clone();
        let side = sides[index];
        let edge = cell.edges[index];
        let on_chord = side.abs() <= eps;
        if side >= -eps {
            positive.push(RingEntry {
                vertex: vertex.clone(),
                edge,
                on_chord,
            });
        }
        if side <= eps {
            negative.push(RingEntry {
                vertex: vertex.clone(),
                edge,
                on_chord,
            });
        }
        let next = (index + 1) % count;
        let next_side = sides[next];
        if (side > eps && next_side < -eps) || (side < -eps && next_side > eps) {
            let crossing =
                interpolate(&cell.verts[index], &cell.verts[next], side, next_side, face);
            positive.push(RingEntry {
                vertex: crossing.clone(),
                edge,
                on_chord: true,
            });
            negative.push(RingEntry {
                vertex: crossing,
                edge,
                on_chord: true,
            });
        }
    });

    [positive, negative]
        .iter()
        .filter_map(|ring| close_ring(ring, chord, eps))
        .collect()
}

fn close_ring(ring: &[RingEntry], chord: EdgeLabel, eps: f64) -> Option<Cell> {
    if ring.len() < 3 {
        return None;
    }
    let count = ring.len();
    let verts: Vec<CellVertex> = ring.iter().map(|entry| entry.vertex.clone()).collect();
    let edges: Vec<EdgeLabel> = (0..count)
        .map(|index| {
            let next = (index + 1) % count;
            if ring[index].on_chord && ring[next].on_chord {
                chord
            } else {
                ring[index].edge
            }
        })
        .collect();
    let local_ring: Vec<(f64, f64)> = verts.iter().map(|vertex| vertex.local).collect();
    (signed_area(&local_ring).abs() > eps * eps).then_some(Cell { verts, edges })
}

fn interpolate(
    tail: &CellVertex,
    head: &CellVertex,
    tail_side: f64,
    head_side: f64,
    face: &PlanarFace,
) -> CellVertex {
    let t = tail_side / (tail_side - head_side);
    let local = (
        tail.local.0 + (head.local.0 - tail.local.0) * t,
        tail.local.1 + (head.local.1 - tail.local.1) * t,
    );
    CellVertex {
        local,
        position: face.plane.point_at_local(local.0, local.1, 0.0),
        label: VertexLabel {
            feature: face.label.feature,
            role: VertexRole::Intersection {
                faces: [face.label, face.label, face.label],
            },
        },
    }
}

fn chord_label(face: FaceLabel, cut: FaceLabel) -> EdgeLabel {
    EdgeLabel {
        feature: face.feature.min(cut.feature),
        role: EdgeRole::intersection(face, cut),
    }
}

pub(crate) fn conform(faces: &[PlanarFace], tolerance: Tolerance) -> Vec<PlanarFace> {
    let nodes = collect_nodes(faces, tolerance);
    faces
        .iter()
        .map(|face| {
            if face.holes.is_empty() && face.curves.iter().all(is_line) {
                insert_nodes(face, &nodes, tolerance)
            } else {
                face.clone()
            }
        })
        .collect()
}

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

fn collect_nodes(faces: &[PlanarFace], tolerance: Tolerance) -> Vec<PlanarVertex> {
    let eps = tolerance.value();
    faces.iter().flat_map(|face| face.polygon.iter()).fold(
        Vec::<PlanarVertex>::new(),
        |mut nodes, vertex| {
            let seen = nodes
                .iter()
                .any(|node| (node.position - vertex.position).norm_mm() <= eps);
            if !seen {
                nodes.push(*vertex);
            }
            nodes
        },
    )
}

fn insert_nodes(face: &PlanarFace, nodes: &[PlanarVertex], tolerance: Tolerance) -> PlanarFace {
    let count = face.polygon.len();
    let (polygon, edges) = (0..count).fold(
        (Vec::<PlanarVertex>::new(), Vec::<EdgeLabel>::new()),
        |(mut polygon, mut edges), index| {
            let tail = face.polygon[index];
            let head = face.polygon[(index + 1) % count];
            let edge = face.edges[index];
            polygon.push(tail);
            edges.push(edge);
            let mut between: Vec<(f64, PlanarVertex)> = nodes
                .iter()
                .filter_map(|node| {
                    on_segment(tail.position, head.position, node.position, tolerance)
                        .map(|t| (t, *node))
                })
                .collect();
            between.sort_by(|(a, _), (b, _)| a.total_cmp(b));
            between.into_iter().for_each(|(_, node)| {
                polygon.push(node);
                edges.push(edge);
            });
            (polygon, edges)
        },
    );
    PlanarFace::straight(face.label, face.plane, polygon, edges)
}

fn on_segment(tail: Point3, head: Point3, point: Point3, tolerance: Tolerance) -> Option<f64> {
    let eps = tolerance.value();
    let span = head - tail;
    let length = span.norm_mm();
    if length <= eps {
        return None;
    }
    let offset = point - tail;
    let t = offset.dot_mm2(span) / (length * length);
    let along = t * length;
    if along <= eps || along >= length - eps {
        return None;
    }
    let projected = tail + span * t;
    ((point - projected).norm_mm() <= eps).then_some(t)
}
