use std::collections::BTreeSet;
use std::collections::HashMap;
use std::collections::HashSet;

use bone_types::{
    Axis3, BrepEdgeId, BrepFaceId, BrepLoopId, BrepShellId, BrepVertexId, CreaseAngle, EdgeLabel,
    EdgeRole, FaceFingerprint, FaceLabel, Parameter, Plane3, Point3, SideKind, SketchEntityId,
    SketchPlaneBasis, Tolerance, UnitVec3, Vec3, VertexLabel, VertexRole,
};
use slotmap::{Key, SlotMap};
use truck_modeling::{
    BoundedCurve, Edge, EdgeID, FaceID, ParameterDivision1D, ParametricCurve, ParametricSurface,
    ParametricSurface3D, SPHint2D, SearchNearestParameter, Shell, Solid, Surface, Transform,
    VertexID,
};

use super::convert::{point_from_truck, try_unit_from_truck};
use super::edges::{EdgeCurve3, crease_from_normals, line_between};
use super::{BrepEdge, BrepError, BrepFace, BrepLoop, BrepShell, BrepSolid, BrepVertex, LabelKind};
use crate::curve3::Curve3;

#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord)]
pub(crate) struct BoundaryIndex(usize);

#[derive(Clone, Copy, Debug)]
pub(crate) struct EdgeArenaHandle(usize);

#[derive(Clone)]
pub(crate) struct Arena {
    solid: Solid,
    edges: Vec<Edge>,
    face_index: HashMap<FaceID, BrepFaceId>,
}

impl Arena {
    pub(crate) fn boundary(&self, index: BoundaryIndex) -> &Shell {
        &self.solid.boundaries()[index.0]
    }

    pub(crate) fn edge(&self, handle: EdgeArenaHandle) -> &Edge {
        &self.edges[handle.0]
    }

    pub(crate) fn solid(&self) -> &Solid {
        &self.solid
    }

    pub(crate) fn face_index(&self) -> &HashMap<FaceID, BrepFaceId> {
        &self.face_index
    }

    pub(crate) fn truck_face(&self, target: BrepFaceId) -> Option<&truck_modeling::Face> {
        self.truck_faces_of(target).next()
    }

    pub(crate) fn truck_faces_of(
        &self,
        target: BrepFaceId,
    ) -> impl Iterator<Item = &truck_modeling::Face> {
        self.solid
            .face_iter()
            .filter(move |face| self.face_index.get(&face.id()) == Some(&target))
    }
}

const FACE_PLANE_TOLERANCE: Tolerance = Tolerance::new(1.0e-9);

pub(crate) fn face_plane_basis(face: &truck_modeling::Face) -> Option<SketchPlaneBasis> {
    let Surface::Plane(plane) = face.surface() else {
        return None;
    };
    let origin = point_from_truck(plane.origin());
    let x = try_unit_from_truck(plane.u_axis(), FACE_PLANE_TOLERANCE)?;
    let normal = try_unit_from_truck(plane.normal(), FACE_PLANE_TOLERANCE)?;
    let outward = if face.orientation() {
        normal
    } else {
        normal.reversed()
    };
    let y = outward.cross(x, FACE_PLANE_TOLERANCE).ok()?;
    SketchPlaneBasis::new(origin, x, y, FACE_PLANE_TOLERANCE).ok()
}

pub(crate) fn face_fingerprint(face: &truck_modeling::Face) -> Option<FaceFingerprint> {
    let basis = face_plane_basis(face)?;
    Some(FaceFingerprint {
        plane: Plane3::from(basis),
        centroid: point_from_truck(super::persist::face_centroid(face)),
    })
}

pub(crate) fn revolution_axis(face: &truck_modeling::Face) -> Option<Axis3> {
    let Surface::RevolutedCurve(processor) = face.surface() else {
        return None;
    };
    let placement = *processor.transform();
    let revolved = processor.entity();
    let origin = point_from_truck(placement.transform_point(revolved.origin()));
    let direction = try_unit_from_truck(
        placement.transform_vector(revolved.axis()),
        FACE_PLANE_TOLERANCE,
    )?;
    Some(Axis3::new(origin, direction))
}

const LENGTH_DIVISION_TOLERANCE: Tolerance = Tolerance::new(1.0e-4);

pub(crate) fn edge_points(edge: &Edge) -> Vec<truck_modeling::Point3> {
    let curve = edge.curve();
    let (_, points) =
        curve.parameter_division(curve.range_tuple(), LENGTH_DIVISION_TOLERANCE.value());
    points
}

pub(crate) fn edge_length(edge: &Edge) -> f64 {
    edge_points(edge)
        .windows(2)
        .map(|pair| {
            let span = pair[1] - pair[0];
            span.x.hypot(span.y).hypot(span.z)
        })
        .sum()
}

#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, PartialOrd, Ord)]
pub(crate) struct FragmentId(u16);

impl FragmentId {
    pub(crate) const WHOLE: Self = Self(0);

    pub(crate) const fn new(index: u16) -> Self {
        Self(index)
    }
}

type FaceKey = (FaceLabel, FragmentId);
type EdgeKey = (EdgeLabel, FragmentId);

pub(crate) struct SolidLabeling {
    pub(crate) faces: HashMap<FaceID, FaceLabel>,
    pub(crate) edges: HashMap<EdgeID, EdgeLabel>,
    pub(crate) vertices: HashMap<VertexID, VertexLabel>,
    pub(crate) closed_curves: HashSet<SketchEntityId>,
    pub(crate) edge_curves: HashMap<EdgeID, EdgeCurve3>,
    pub(crate) face_fragments: HashMap<FaceID, FragmentId>,
    pub(crate) edge_fragments: HashMap<EdgeID, FragmentId>,
    pub(crate) vertex_fragments: HashMap<VertexID, FragmentId>,
}

fn suppressed_edge(label: EdgeLabel, closed: &HashSet<SketchEntityId>) -> bool {
    matches!(label.role, EdgeRole::SideEdge { from, side: SideKind::Seam } if !closed.contains(&from))
}

fn suppressed_vertex(label: VertexLabel, closed: &HashSet<SketchEntityId>) -> bool {
    matches!(
        label.role,
        VertexRole::StartCapVertex { from, side: SideKind::Seam }
            | VertexRole::EndCapVertex { from, side: SideKind::Seam }
            | VertexRole::Seam { from }
        if !closed.contains(&from)
    )
}

pub(crate) fn assemble(solid: Solid, labeling: &SolidLabeling) -> Result<BrepSolid, BrepError> {
    let vertices = build_vertices(&solid, labeling)?;
    let truck_faces: HashMap<FaceID, truck_modeling::Face> = solid
        .boundaries()
        .iter()
        .flat_map(Shell::face_iter)
        .map(|face| (face.id(), face.clone()))
        .collect();
    let truck_face_adj = gather_truck_face_adjacency(&solid);
    let edges = build_edges(
        &solid,
        labeling,
        &vertices.by_truck,
        &vertices.positions,
        &truck_face_adj,
        &truck_faces,
    )?;
    let faces = build_faces(&solid, labeling, &edges)?;

    let shell_order = ordered_by(&faces.shells, |shell| shell.boundary_index);
    let face_order = ordered_by(&faces.faces, BrepFace::label);
    let edge_order = ordered_by(&edges.map, BrepEdge::label);
    let vertex_order = ordered_by(&vertices.map, BrepVertex::label);
    let loop_order: Vec<BrepLoopId> = face_order
        .iter()
        .flat_map(|face_id| faces.faces[*face_id].loops().iter().copied())
        .collect();

    let reattach = super::persist::capture(&solid, labeling)?;

    Ok(BrepSolid {
        arena: Arena {
            solid,
            edges: edges.arena,
            face_index: faces.face_index,
        },
        shells: faces.shells,
        faces: faces.faces,
        loops: faces.loops,
        edges: edges.map,
        vertices: vertices.map,
        shell_order,
        face_order,
        loop_order,
        edge_order,
        vertex_order,
        edge_faces: faces.edge_faces,
        reattach,
    })
}

fn gather_truck_face_adjacency(solid: &Solid) -> HashMap<EdgeID, Vec<FaceID>> {
    solid.boundaries().iter().flat_map(Shell::face_iter).fold(
        HashMap::<EdgeID, Vec<FaceID>>::new(),
        |mut adj, face| {
            face.boundaries().iter().for_each(|wire| {
                wire.edge_iter().for_each(|edge| {
                    adj.entry(edge.id()).or_default().push(face.id());
                });
            });
            adj
        },
    )
}

fn ordered_by<K: Key, V, T: Ord>(map: &SlotMap<K, V>, key: impl Fn(&V) -> T) -> Vec<K> {
    let mut keys: Vec<K> = map.keys().collect();
    keys.sort_by_cached_key(|k| key(&map[*k]));
    keys
}

struct VertexTable {
    map: SlotMap<BrepVertexId, BrepVertex>,
    by_truck: HashMap<VertexID, BrepVertexId>,
    positions: HashMap<BrepVertexId, Point3>,
}

fn build_vertices(solid: &Solid, labeling: &SolidLabeling) -> Result<VertexTable, BrepError> {
    let (map, _by_label, by_truck, positions) = solid.vertex_iter().try_fold(
        (
            SlotMap::with_key(),
            HashMap::<(VertexLabel, FragmentId), BrepVertexId>::new(),
            HashMap::<VertexID, BrepVertexId>::new(),
            HashMap::<BrepVertexId, Point3>::new(),
        ),
        |(mut map, mut by_label, mut by_truck, mut positions), vertex| {
            if by_truck.contains_key(&vertex.id()) {
                return Ok::<_, BrepError>((map, by_label, by_truck, positions));
            }
            let label = *labeling
                .vertices
                .get(&vertex.id())
                .ok_or(BrepError::MissingLabel {
                    kind: LabelKind::Vertex,
                })?;
            if suppressed_vertex(label, &labeling.closed_curves) {
                return Ok((map, by_label, by_truck, positions));
            }
            let fragment = labeling
                .vertex_fragments
                .get(&vertex.id())
                .copied()
                .unwrap_or(FragmentId::WHOLE);
            let key = (label, fragment);
            let position = point_from_truck(vertex.point());
            let brep_id = by_label.get(&key).copied().unwrap_or_else(|| {
                let id = map.insert_with_key(|id| BrepVertex {
                    id,
                    label,
                    position,
                });
                by_label.insert(key, id);
                id
            });
            by_truck.insert(vertex.id(), brep_id);
            positions.entry(brep_id).or_insert(position);
            Ok((map, by_label, by_truck, positions))
        },
    )?;
    Ok(VertexTable {
        map,
        by_truck,
        positions,
    })
}

struct EdgeTable {
    map: SlotMap<BrepEdgeId, BrepEdge>,
    by_truck: HashMap<EdgeID, BrepEdgeId>,
    arena: Vec<Edge>,
}

struct EdgeAccum {
    handles: Vec<EdgeArenaHandle>,
    incident: Vec<BrepVertexId>,
    analytic: Option<EdgeCurve3>,
}

fn build_edges(
    solid: &Solid,
    labeling: &SolidLabeling,
    vertex_dedup: &HashMap<VertexID, BrepVertexId>,
    vertex_positions: &HashMap<BrepVertexId, Point3>,
    truck_face_adj: &HashMap<EdgeID, Vec<FaceID>>,
    truck_faces: &HashMap<FaceID, truck_modeling::Face>,
) -> Result<EdgeTable, BrepError> {
    let (accum, arena, truck_label) = solid.edge_iter().try_fold(
        (
            HashMap::<EdgeKey, EdgeAccum>::new(),
            Vec::<Edge>::new(),
            HashMap::<EdgeID, EdgeKey>::new(),
        ),
        |(mut accum, mut arena, mut truck_label), edge| {
            if truck_label.contains_key(&edge.id()) {
                return Ok::<_, BrepError>((accum, arena, truck_label));
            }
            let label = *labeling
                .edges
                .get(&edge.id())
                .ok_or(BrepError::MissingLabel {
                    kind: LabelKind::Edge,
                })?;
            if suppressed_edge(label, &labeling.closed_curves) {
                return Ok((accum, arena, truck_label));
            }
            let fragment = labeling
                .edge_fragments
                .get(&edge.id())
                .copied()
                .unwrap_or(FragmentId::WHOLE);
            let key = (label, fragment);
            let handle = EdgeArenaHandle(arena.len());
            arena.push(edge.clone());
            let accumulator = accum.entry(key).or_insert_with(|| EdgeAccum {
                handles: Vec::new(),
                incident: Vec::new(),
                analytic: None,
            });
            accumulator.handles.push(handle);
            if accumulator.analytic.is_none() {
                accumulator.analytic = labeling.edge_curves.get(&edge.id()).cloned();
            }
            [edge.front().id(), edge.back().id()]
                .into_iter()
                .filter_map(|truck_vertex| vertex_dedup.get(&truck_vertex).copied())
                .for_each(|brep_vertex| accumulator.incident.push(brep_vertex));
            truck_label.insert(edge.id(), key);
            Ok((accum, arena, truck_label))
        },
    )?;

    let mut ordered_edges: Vec<(EdgeKey, EdgeAccum)> = accum.into_iter().collect();
    ordered_edges.sort_by_key(|(key, _)| *key);
    let (map, label_to_id) = ordered_edges.into_iter().fold(
        (SlotMap::with_key(), HashMap::<EdgeKey, BrepEdgeId>::new()),
        |(mut map, mut label_to_id), (key, slot)| {
            let primary = &arena[slot.handles[0].0];
            let curve = resolve_curve(slot.analytic, primary);
            let vertices = ordered_endpoints(&slot.incident, &curve, vertex_positions);
            let crease = compute_crease(truck_face_adj.get(&primary.id()), truck_faces, primary);
            let id = map.insert_with_key(|id| BrepEdge {
                id,
                label: key.0,
                handles: slot.handles,
                vertices,
                curve,
                crease,
            });
            label_to_id.insert(key, id);
            (map, label_to_id)
        },
    );

    let by_truck = truck_label
        .into_iter()
        .map(|(truck_id, key)| (truck_id, label_to_id[&key]))
        .collect();

    Ok(EdgeTable {
        map,
        by_truck,
        arena,
    })
}

fn resolve_curve(analytic: Option<EdgeCurve3>, primary: &Edge) -> EdgeCurve3 {
    analytic.unwrap_or_else(|| {
        let front = point_from_truck(primary.absolute_front().point());
        let back = point_from_truck(primary.absolute_back().point());
        line_between(front, back)
    })
}

fn compute_crease(
    adjacent: Option<&Vec<FaceID>>,
    truck_faces: &HashMap<FaceID, truck_modeling::Face>,
    edge: &Edge,
) -> CreaseAngle {
    let Some(face_ids) = adjacent else {
        return CreaseAngle::FLAT;
    };
    let midpoint = edge_midpoint(edge);
    let normals: Vec<UnitVec3> = face_ids
        .iter()
        .filter_map(|face_id| truck_faces.get(face_id))
        .filter_map(|face| face_normal_at(face, midpoint))
        .collect();
    match normals.as_slice() {
        [n1, n2] => crease_from_normals(*n1, *n2),
        _ => CreaseAngle::FLAT,
    }
}

fn edge_midpoint(edge: &Edge) -> truck_modeling::Point3 {
    let curve = edge.curve();
    let (t0, t1) = curve.range_tuple();
    curve.subs(f64::midpoint(t0, t1))
}

fn face_normal_at(face: &truck_modeling::Face, point: truck_modeling::Point3) -> Option<UnitVec3> {
    face_normal_and_distance(face, point).map(|(normal, _)| normal)
}

const NEAREST_SEED_DIVISION: u32 = 8;
const NEAREST_NEWTON_TRIALS: usize = 100;

fn nearest_surface_parameter(
    surface: &Surface,
    point: truck_modeling::Point3,
) -> Option<(f64, f64)> {
    let hint = match surface.try_range_tuple() {
        (Some((u0, u1)), Some((v0, v1))) => {
            let query = point_from_truck(point);
            let div = NEAREST_SEED_DIVISION;
            let axis =
                |lo: f64, hi: f64, step: u32| lo + (hi - lo) * f64::from(step) / f64::from(div);
            (0..=div)
                .flat_map(|i| (0..=div).map(move |j| (i, j)))
                .map(|(i, j)| {
                    let (u, v) = (axis(u0, u1, i), axis(v0, v1, j));
                    (
                        u,
                        v,
                        (point_from_truck(surface.subs(u, v)) - query).norm_mm(),
                    )
                })
                .min_by(|a, b| a.2.total_cmp(&b.2))
                .map_or(SPHint2D::None, |(u, v, _)| SPHint2D::Parameter(u, v))
        }
        _ => SPHint2D::None,
    };
    surface.search_nearest_parameter(point, hint, NEAREST_NEWTON_TRIALS)
}

pub(crate) fn face_normal_and_distance(
    face: &truck_modeling::Face,
    point: truck_modeling::Point3,
) -> Option<(UnitVec3, f64)> {
    let surface = face.surface();
    let query = point_from_truck(point);
    let (u, v) = nearest_surface_parameter(&surface, point)?;
    let mut normal = surface.normal(u, v);
    if !face.orientation() {
        normal = -normal;
    }
    let unit = try_unit_from_truck(normal, Tolerance::new(1.0e-9))?;
    Some((
        unit,
        (point_from_truck(surface.subs(u, v)) - query).norm_mm(),
    ))
}

const FACE_BOUNDARY_CHORD_MM: f64 = 1.0e-2;

pub(crate) fn face_normal_and_bounded_distance(
    face: &truck_modeling::Face,
    point: truck_modeling::Point3,
) -> Option<(UnitVec3, f64)> {
    let surface = face.surface();
    let query = point_from_truck(point);
    let (u, v) = nearest_surface_parameter(&surface, point)?;
    let mut normal = surface.normal(u, v);
    if !face.orientation() {
        normal = -normal;
    }
    let unit = try_unit_from_truck(normal, Tolerance::new(1.0e-9))?;
    let foot = point_from_truck(surface.subs(u, v));
    let perpendicular = (foot - query).norm_mm();
    Some((
        unit,
        bounded_face_distance(face, query, foot, unit, perpendicular),
    ))
}

fn bounded_face_distance(
    face: &truck_modeling::Face,
    query: Point3,
    foot: Point3,
    normal: UnitVec3,
    perpendicular: f64,
) -> f64 {
    let loops = face_boundary_loops(face);
    let Some((ex, ey)) = face_tangent_basis(normal) else {
        return perpendicular;
    };
    if loops.iter().all(|poly| poly.len() < 2) {
        return perpendicular;
    }
    let project = |p: Point3| (along(p - foot, ex), along(p - foot, ey));
    let segments = || {
        loops.iter().flat_map(|poly| {
            poly.iter()
                .copied()
                .zip(poly.iter().copied().cycle().skip(1))
                .take(poly.len())
        })
    };
    let crossings = segments()
        .filter(|&(a, b)| crosses_positive_x(project(a), project(b)))
        .count();
    if crossings % 2 == 1 {
        perpendicular
    } else {
        segments()
            .map(|(a, b)| point_segment_distance(query, a, b))
            .fold(f64::INFINITY, f64::min)
    }
}

fn face_boundary_loops(face: &truck_modeling::Face) -> Vec<Vec<Point3>> {
    face.boundaries()
        .iter()
        .map(|wire| wire.iter().flat_map(sample_boundary).collect())
        .collect()
}

fn sample_boundary(edge: &Edge) -> Vec<Point3> {
    let curve = edge.curve();
    let (t0, t1) = curve.range_tuple();
    let (_, points) = curve.parameter_division((t0, t1), FACE_BOUNDARY_CHORD_MM);
    points.iter().map(|&p| point_from_truck(p)).collect()
}

fn face_tangent_basis(normal: UnitVec3) -> Option<(UnitVec3, UnitVec3)> {
    let (nx, ny, nz) = normal.components();
    let n = Vec3::from_mm(nx, ny, nz);
    let seed = if nx.abs() <= ny.abs() && nx.abs() <= nz.abs() {
        Vec3::from_mm(1.0, 0.0, 0.0)
    } else if ny.abs() <= nz.abs() {
        Vec3::from_mm(0.0, 1.0, 0.0)
    } else {
        Vec3::from_mm(0.0, 0.0, 1.0)
    };
    let tolerance = Tolerance::new(1.0e-9);
    let ex = seed.cross(n).try_normalize(tolerance).ok()?;
    let (ex_x, ex_y, ex_z) = ex.components();
    let ey = n
        .cross(Vec3::from_mm(ex_x, ex_y, ex_z))
        .try_normalize(tolerance)
        .ok()?;
    Some((ex, ey))
}

fn along(vector: Vec3, axis: UnitVec3) -> f64 {
    let (vx, vy, vz) = vector.coords_mm();
    let (ax, ay, az) = axis.components();
    vx * ax + vy * ay + vz * az
}

fn crosses_positive_x(a: (f64, f64), b: (f64, f64)) -> bool {
    let ((ax, ay), (bx, by)) = (a, b);
    if (ay > 0.0) == (by > 0.0) {
        return false;
    }
    let denom = ay - by;
    denom != 0.0 && ax + (ay / denom) * (bx - ax) > 0.0
}

fn point_segment_distance(p: Point3, a: Point3, b: Point3) -> f64 {
    let ab = b - a;
    let ap = p - a;
    let denom = ab.dot_mm2(ab);
    let t = if denom > 0.0 {
        (ap.dot_mm2(ab) / denom).clamp(0.0, 1.0)
    } else {
        0.0
    };
    let (apx, apy, apz) = ap.coords_mm();
    let (abx, aby, abz) = ab.coords_mm();
    let (dx, dy, dz) = (apx - t * abx, apy - t * aby, apz - t * abz);
    dx.mul_add(dx, dy.mul_add(dy, dz * dz)).sqrt()
}

fn ordered_endpoints(
    incident: &[BrepVertexId],
    curve: &EdgeCurve3,
    vertex_positions: &HashMap<BrepVertexId, Point3>,
) -> [BrepVertexId; 2] {
    let counts = incident.iter().fold(
        HashMap::<BrepVertexId, usize>::new(),
        |mut counts, vertex| {
            *counts.entry(*vertex).or_insert(0) += 1;
            counts
        },
    );
    let mut odd: Vec<BrepVertexId> = counts
        .iter()
        .filter_map(|(vertex, count)| (*count % 2 == 1).then_some(*vertex))
        .collect();
    odd.sort_unstable();
    if let [a, b] = odd.as_slice() {
        let start = curve.evaluate(Parameter::new(0.0));
        let a_dist = vertex_positions
            .get(a)
            .map_or(f64::INFINITY, |pos| squared_distance(*pos, start));
        let b_dist = vertex_positions
            .get(b)
            .map_or(f64::INFINITY, |pos| squared_distance(*pos, start));
        if a_dist <= b_dist { [*a, *b] } else { [*b, *a] }
    } else {
        let representative = counts
            .keys()
            .copied()
            .min()
            .or_else(|| vertex_positions.keys().copied().min())
            .unwrap_or_default();
        [representative, representative]
    }
}

fn squared_distance(a: Point3, b: Point3) -> f64 {
    let (ax, ay, az) = a.coords_mm();
    let (bx, by, bz) = b.coords_mm();
    let dx = ax - bx;
    let dy = ay - by;
    let dz = az - bz;
    dx * dx + dy * dy + dz * dz
}

struct FaceTable {
    faces: SlotMap<BrepFaceId, BrepFace>,
    loops: SlotMap<BrepLoopId, BrepLoop>,
    shells: SlotMap<BrepShellId, BrepShell>,
    face_index: HashMap<FaceID, BrepFaceId>,
    edge_faces: HashMap<BrepEdgeId, Vec<BrepFaceId>>,
}

#[derive(Default)]
struct Gathered {
    adjacency: HashMap<EdgeID, Vec<FaceKey>>,
    patches: HashMap<FaceKey, Vec<Vec<Vec<EdgeID>>>>,
    shell_members: Vec<Vec<FaceKey>>,
    face_to_label: Vec<(FaceID, FaceKey)>,
}

fn build_faces(
    solid: &Solid,
    labeling: &SolidLabeling,
    edges: &EdgeTable,
) -> Result<FaceTable, BrepError> {
    let Gathered {
        adjacency,
        patches,
        shell_members,
        face_to_label,
    } = gather_faces(solid, labeling)?;

    let mut ordered_patches: Vec<(FaceKey, Vec<Vec<Vec<EdgeID>>>)> = patches.into_iter().collect();
    ordered_patches.sort_by_key(|(key, _)| *key);
    let (faces, loops, label_to_face) = ordered_patches.into_iter().fold(
        (
            SlotMap::with_key(),
            SlotMap::with_key(),
            HashMap::<FaceKey, BrepFaceId>::new(),
        ),
        |(mut faces, mut loops, mut label_to_face), (key, patch_list)| {
            let loop_ids = face_loops(&patch_list, key, &adjacency, edges, &mut loops);
            let id = faces.insert_with_key(|id| BrepFace {
                id,
                label: key.0,
                loops: loop_ids,
            });
            label_to_face.insert(key, id);
            (faces, loops, label_to_face)
        },
    );

    let shells = shell_members.into_iter().enumerate().fold(
        SlotMap::with_key(),
        |mut shells, (boundary, members)| {
            let face_ids = dedup_preserving(members)
                .into_iter()
                .map(|key| label_to_face[&key])
                .collect();
            shells.insert_with_key(|id| BrepShell {
                id,
                boundary_index: BoundaryIndex(boundary),
                faces: face_ids,
            });
            shells
        },
    );

    let face_index = face_to_label
        .into_iter()
        .map(|(face_id, key)| (face_id, label_to_face[&key]))
        .collect();

    let edge_faces = edge_face_index(&adjacency, edges, &label_to_face);

    Ok(FaceTable {
        faces,
        loops,
        shells,
        face_index,
        edge_faces,
    })
}

fn edge_face_index(
    adjacency: &HashMap<EdgeID, Vec<FaceKey>>,
    edges: &EdgeTable,
    label_to_face: &HashMap<FaceKey, BrepFaceId>,
) -> HashMap<BrepEdgeId, Vec<BrepFaceId>> {
    adjacency
        .iter()
        .filter_map(|(truck_edge, keys)| {
            edges
                .by_truck
                .get(truck_edge)
                .copied()
                .map(|brep_edge| (brep_edge, keys))
        })
        .fold(
            HashMap::<BrepEdgeId, BTreeSet<FaceKey>>::new(),
            |mut acc, (brep_edge, keys)| {
                let owners = acc.entry(brep_edge).or_default();
                keys.iter().for_each(|key| {
                    owners.insert(*key);
                });
                acc
            },
        )
        .into_iter()
        .map(|(brep_edge, owners)| {
            let faces = owners
                .into_iter()
                .filter_map(|key| label_to_face.get(&key).copied())
                .collect();
            (brep_edge, faces)
        })
        .collect()
}

fn gather_faces(solid: &Solid, labeling: &SolidLabeling) -> Result<Gathered, BrepError> {
    solid
        .boundaries()
        .iter()
        .try_fold(Gathered::default(), |mut gathered, shell| {
            let members = shell
                .face_iter()
                .map(|face| {
                    let label = *labeling
                        .faces
                        .get(&face.id())
                        .ok_or(BrepError::MissingLabel {
                            kind: LabelKind::Face,
                        })?;
                    let fragment = labeling
                        .face_fragments
                        .get(&face.id())
                        .copied()
                        .unwrap_or(FragmentId::WHOLE);
                    let key = (label, fragment);
                    let wires: Vec<Vec<EdgeID>> = face
                        .boundaries()
                        .iter()
                        .map(|wire| wire.edge_iter().map(Edge::id).collect())
                        .collect();
                    wires.iter().flatten().for_each(|edge_id| {
                        gathered.adjacency.entry(*edge_id).or_default().push(key);
                    });
                    gathered.patches.entry(key).or_default().push(wires);
                    gathered.face_to_label.push((face.id(), key));
                    Ok::<FaceKey, BrepError>(key)
                })
                .collect::<Result<Vec<_>, _>>()?;
            gathered.shell_members.push(members);
            Ok(gathered)
        })
}

fn face_loops(
    patches: &[Vec<Vec<EdgeID>>],
    key: FaceKey,
    adjacency: &HashMap<EdgeID, Vec<FaceKey>>,
    edges: &EdgeTable,
    loops: &mut SlotMap<BrepLoopId, BrepLoop>,
) -> Vec<BrepLoopId> {
    if let [single] = patches {
        single
            .iter()
            .filter_map(|wire| {
                let loop_edges =
                    dedup_cycle(wire.iter().map(|edge_id| edges.by_truck[edge_id]).collect());
                (!loop_edges.is_empty()).then(|| {
                    loops.insert_with_key(|id| BrepLoop {
                        id,
                        edges: loop_edges,
                    })
                })
            })
            .collect()
    } else {
        let mut seen = HashSet::new();
        let boundary: Vec<BrepEdgeId> = patches
            .iter()
            .flatten()
            .flatten()
            .filter(|edge_id| !internal_to_group(adjacency, **edge_id, key))
            .filter_map(|edge_id| {
                let brep_id = edges.by_truck[edge_id];
                seen.insert(brep_id).then_some(brep_id)
            })
            .collect();
        chain_loops(boundary, &edges.map)
            .into_iter()
            .map(|loop_edges| {
                loops.insert_with_key(|id| BrepLoop {
                    id,
                    edges: loop_edges,
                })
            })
            .collect()
    }
}

fn internal_to_group(
    adjacency: &HashMap<EdgeID, Vec<FaceKey>>,
    edge_id: EdgeID,
    key: FaceKey,
) -> bool {
    adjacency
        .get(&edge_id)
        .is_some_and(|adjacent| adjacent.len() == 2 && adjacent.iter().all(|other| *other == key))
}

fn chain_loops(
    boundary: Vec<BrepEdgeId>,
    edges: &SlotMap<BrepEdgeId, BrepEdge>,
) -> Vec<Vec<BrepEdgeId>> {
    let mut sorted = boundary;
    sorted.sort_by_key(|edge| edges[*edge].label());
    let (closed, open): (Vec<BrepEdgeId>, Vec<BrepEdgeId>) = sorted
        .into_iter()
        .partition(|edge| edges[*edge].vertices()[0] == edges[*edge].vertices()[1]);
    let closed_loops = closed.into_iter().map(|edge| vec![edge]);
    closed_loops.chain(walk_open(open, edges)).collect()
}

fn walk_open(open: Vec<BrepEdgeId>, edges: &SlotMap<BrepEdgeId, BrepEdge>) -> Vec<Vec<BrepEdgeId>> {
    let mut remaining = open;
    std::iter::from_fn(move || {
        let first = (!remaining.is_empty()).then(|| remaining.remove(0))?;
        let start = edges[first].vertices()[0];
        let next = edges[first].vertices()[1];
        Some(extend_loop(next, start, &mut remaining, edges, vec![first]))
    })
    .collect()
}

fn extend_loop(
    current: BrepVertexId,
    start: BrepVertexId,
    remaining: &mut Vec<BrepEdgeId>,
    edges: &SlotMap<BrepEdgeId, BrepEdge>,
    acc: Vec<BrepEdgeId>,
) -> Vec<BrepEdgeId> {
    if current == start {
        return acc;
    }
    match remaining
        .iter()
        .position(|edge| edges[*edge].vertices().contains(&current))
    {
        None => acc,
        Some(position) => {
            let edge = remaining.remove(position);
            let ends = edges[edge].vertices();
            let next = if ends[0] == current { ends[1] } else { ends[0] };
            extend_loop(
                next,
                start,
                remaining,
                edges,
                acc.into_iter().chain(std::iter::once(edge)).collect(),
            )
        }
    }
}

fn dedup_cycle(ids: Vec<BrepEdgeId>) -> Vec<BrepEdgeId> {
    let mut out = ids.into_iter().fold(Vec::new(), |mut acc, id| {
        if acc.last() != Some(&id) {
            acc.push(id);
        }
        acc
    });
    if out.len() > 1 && out.first() == out.last() {
        out.pop();
    }
    out
}

fn dedup_preserving<T: Eq + std::hash::Hash + Copy>(items: Vec<T>) -> Vec<T> {
    let mut seen = HashSet::new();
    items
        .into_iter()
        .filter(|item| seen.insert(*item))
        .collect()
}

#[cfg(test)]
mod tests {
    use super::revolution_axis;
    use truck_modeling::{Point3 as TruckPoint3, Rad, Surface, Vector3 as TruckVector3, builder};

    #[test]
    fn revolution_axis_recovers_the_turning_axis() {
        let inner = builder::vertex(TruckPoint3::new(2.0, 0.0, 0.0));
        let outer = builder::vertex(TruckPoint3::new(3.0, 0.0, 0.0));
        let edge = builder::line(&inner, &outer);
        let profile = builder::tsweep(&edge, TruckVector3::new(0.0, 0.0, 1.0));
        let tube = builder::rsweep(
            &profile,
            TruckPoint3::new(0.0, 0.0, 0.0),
            TruckVector3::unit_z(),
            Rad(2.0 * std::f64::consts::PI),
        );
        let Some(revolved) = tube
            .face_iter()
            .find(|face| matches!(face.surface(), Surface::RevolutedCurve(_)))
        else {
            panic!("the revolved tube carries a surface of revolution");
        };
        let Some(axis) = revolution_axis(revolved) else {
            panic!("a revolved face yields its turning axis");
        };
        let (dx, dy, dz) = axis.direction().components();
        assert!(dx.abs() < 1.0e-9 && dy.abs() < 1.0e-9 && (dz.abs() - 1.0).abs() < 1.0e-9);
        let (ox, oy, _oz) = axis.origin().coords_mm();
        assert!(ox.abs() < 1.0e-9 && oy.abs() < 1.0e-9);
    }
}
