BREP_reconstruction 0.2.1

//! Conservative local BREP reconstruction from the kernel's public mesh segmentation.
//!
//! Plane and safely bounded cylinder/cone regions remain analytic. Unsupported
//! regions, and revolved regions whose interfaces cannot be shared exactly, are
//! represented by maximal edge-connected, provably coplanar planar patches.

use brep_kernel::{
    build_pcurve_on_surface_range, circle_angle_to_parameter, intersect_analytic_pair, make_arc,
    make_line, make_revolution, mesh_to_faceted_brep, project_point_to_curve, segment_mesh_faces,
    solid_signed_volume, ArenaCoedge, ArenaEdge, ArenaFace, ArenaLoop, ArenaShell, ArenaVertex,
    BrepSolid, EdgeId, FaceId, MeshRegion, NurbsCurve, NurbsSurface, RegionCarrier, SegmentOptions,
    TopologyArena, Vec3, Vec4,
};
use std::collections::{BTreeMap, BTreeSet, HashMap, VecDeque};
use std::f64::consts::TAU;

#[derive(Clone, Copy)]
struct PlaneInfo {
    origin: Vec3,
    normal: Vec3,
}

#[derive(Clone, Copy)]
struct CylinderCarrier {
    origin: Vec3,
    axis: Vec3,
    radius: f64,
    sense: i8,
}

#[derive(Clone, Copy)]
struct ConeCarrier {
    apex: Vec3,
    axis: Vec3,
    half_angle: f64,
    sense: i8,
}

struct CylinderInfo {
    carrier: CylinderCarrier,
    base_station: f64,
    height: f64,
    x_axis: Vec3,
    y_axis: Vec3,
    surface: NurbsSurface,
}

struct ConeInfo {
    carrier: ConeCarrier,
    base_station: f64,
    height: f64,
    x_axis: Vec3,
    y_axis: Vec3,
    surface: NurbsSurface,
}

#[derive(Clone)]
struct Chain {
    vertices: Vec<usize>,
    closed: bool,
    adjacent: (u32, u32),
}

#[derive(Clone, Copy)]
struct Traversal {
    chain: usize,
    forward: bool,
}

struct RegionBoundary {
    cycles: Vec<Vec<Traversal>>,
}

struct LocalRegions {
    triangle_regions: Vec<u32>,
    planes: HashMap<u32, PlaneInfo>,
    cylinders: HashMap<u32, CylinderCarrier>,
    cones: HashMap<u32, ConeCarrier>,
    local_to_original: HashMap<u32, u32>,
    count: usize,
    stats: HybridBrepStats,
}

enum HybridBuildFailure {
    Region { original: u32, message: String },
    Global(String),
}

impl From<String> for HybridBuildFailure {
    fn from(message: String) -> Self {
        Self::Global(message)
    }
}

impl From<&str> for HybridBuildFailure {
    fn from(message: &str) -> Self {
        Self::Global(message.to_owned())
    }
}

/// Auditable accounting for an analytic/faceted reconstruction.
#[derive(Clone, Copy, Debug, Default, Eq, PartialEq, serde::Serialize)]
pub(crate) struct HybridBrepStats {
    pub(crate) total_faces: usize,
    pub(crate) analytic_plane_faces: usize,
    pub(crate) analytic_plane_triangles: usize,
    pub(crate) analytic_cylinder_faces: usize,
    pub(crate) analytic_cylinder_triangles: usize,
    pub(crate) analytic_cone_faces: usize,
    pub(crate) analytic_cone_triangles: usize,
    pub(crate) faceted_faces: usize,
    pub(crate) faceted_triangles: usize,
    pub(crate) demoted_regions: usize,
    pub(crate) demoted_triangles: usize,
    pub(crate) segmentation_plane_regions: usize,
    pub(crate) segmentation_cylinder_regions: usize,
    pub(crate) segmentation_cone_regions: usize,
    pub(crate) segmentation_unsupported_regions: usize,
}

/// Validated reconstructed solid and the accounting used to select policy.
pub(crate) struct HybridBrepOutput {
    pub(crate) solid: BrepSolid,
    pub(crate) stats: HybridBrepStats,
}

#[derive(Clone, Copy)]
struct ChainUv {
    region: u32,
    start: [f64; 2],
    end: [f64; 2],
}

struct BuiltEdge {
    edge: EdgeId,
    curve_along_chain: bool,
    revolve_uvs: Vec<ChainUv>,
    custom_revolve_pcurves: HashMap<u32, NurbsCurve>,
}

struct Ids(u64);

impl Ids {
    fn next(&mut self) -> u64 {
        let value = self.0;
        self.0 += 1;
        value
    }
}

/// Rebuild a validated analytic/faceted shell from public segmentation.
///
/// This is intentionally narrower than a general mesh-to-BREP converter.
/// Plane and safely bounded cylinder/cone regions remain analytic. Unsupported
/// regions become maximal edge-connected coplanar planar patches. A cylinder touching
/// those facets is retained only when every shared edge is an exact axial
/// ruling; otherwise the complete cylinder region is recursively demoted to
/// planar triangle faces. Every shared boundary must still be provably a line,
/// circular arc, complete ring, or cylinder ruling.
pub(crate) fn hybrid_plane_cylinder_brep(
    positions: &[f64],
    indices: Option<&[u32]>,
    options: &SegmentOptions,
) -> Result<HybridBrepOutput, String> {
    hybrid_plane_cylinder_brep_retry(positions, indices, options, None)
}

fn hybrid_plane_cylinder_brep_retry(
    positions: &[f64],
    indices: Option<&[u32]>,
    options: &SegmentOptions,
    mut injected_failure: Option<u32>,
) -> Result<HybridBrepOutput, String> {
    let mut forced_demotions = BTreeSet::new();
    loop {
        let attempt = if let Some(original) = injected_failure.take() {
            Err(HybridBuildFailure::Region {
                original,
                message: "injected localized carrier construction failure".into(),
            })
        } else {
            hybrid_plane_cylinder_brep_once(positions, indices, options, &forced_demotions)
        };
        match attempt {
            Ok(output) => return Ok(output),
            Err(HybridBuildFailure::Region { original, message }) => {
                if !forced_demotions.insert(original) {
                    return Err(format!(
                        "hybrid BREP: region {original} remained unconstructible after local demotion: {message}"
                    ));
                }
            }
            Err(HybridBuildFailure::Global(message)) => return Err(message),
        }
    }
}

fn hybrid_plane_cylinder_brep_once(
    positions: &[f64],
    indices: Option<&[u32]>,
    options: &SegmentOptions,
    forced_demotions: &BTreeSet<u32>,
) -> Result<HybridBrepOutput, HybridBuildFailure> {
    let owned_indices;
    let indices = match indices {
        Some(value) => value,
        None => {
            owned_indices = (0..positions.len() as u32 / 3).collect::<Vec<_>>();
            &owned_indices
        }
    };
    if positions.is_empty()
        || !positions.len().is_multiple_of(3)
        || !indices.len().is_multiple_of(3)
    {
        return Err("hybrid BREP: invalid triangle buffers".into());
    }
    let vertices = positions
        .chunks_exact(3)
        .map(|p| Vec3::new(p[0], p[1], p[2]))
        .collect::<Vec<_>>();
    let triangles = indices
        .chunks_exact(3)
        .map(|t| [t[0] as usize, t[1] as usize, t[2] as usize])
        .collect::<Vec<_>>();
    if triangles
        .iter()
        .flatten()
        .any(|&vertex| vertex >= vertices.len())
    {
        return Err("hybrid BREP: triangle index outside vertex buffer".into());
    }

    let segmentation = segment_mesh_faces(positions, indices, options)?;
    if segmentation.welded_vertex_count != vertices.len() {
        return Err(format!(
            "hybrid BREP: segmentation welded {} input vertices to {}; pre-weld the mesh before exact reconstruction",
            vertices.len(), segmentation.welded_vertex_count
        ).into());
    }
    if segmentation.triangle_region_ids.len() != triangles.len() {
        return Err("hybrid BREP: segmentation triangle count drift".into());
    }

    let diagonal = mesh_diagonal(&vertices);
    // Exact shared topology must be substantially tighter than carrier
    // recognition. This still covers deterministic f32 STL quantization.
    let tolerance = (options.fit_tolerance * diagonal)
        .min(2.0e-6 * diagonal)
        .max(1.0e-10);
    let incidence = edge_incidence(&triangles)?;
    if incidence.values().any(|incident| incident.len() != 2) {
        return Err("hybrid BREP: input is not a closed two-manifold".into());
    }
    let local = local_regions(
        &segmentation.regions,
        &segmentation.triangle_region_ids,
        &triangles,
        &vertices,
        &incidence,
        tolerance,
        forced_demotions,
    )?;
    let stats = local.stats;
    let local_to_original = local.local_to_original;
    let planes = local.planes;
    let mut cylinder_carriers = local.cylinders;
    let mut cone_carriers = local.cones;
    let (chains, mut boundaries) =
        region_boundaries(&triangles, &local.triangle_regions, local.count, &incidence)?;

    snap_cylinder_axes(
        &mut cylinder_carriers,
        &planes,
        &chains,
        &vertices,
        tolerance,
    )
    .map_err(|(region, message)| HybridBuildFailure::Region {
        original: local_to_original[&region],
        message,
    })?;
    snap_cone_carriers(
        &mut cone_carriers,
        &cylinder_carriers,
        &planes,
        &chains,
        &vertices,
        tolerance,
    )
    .map_err(|(region, message)| HybridBuildFailure::Region {
        original: local_to_original[&region],
        message,
    })?;
    for (&region, cone) in &cone_carriers {
        let mut low = f64::INFINITY;
        let mut high = f64::NEG_INFINITY;
        for chain in chains
            .iter()
            .filter(|chain| chain.adjacent.0 == region || chain.adjacent.1 == region)
        {
            for point in chain_points(chain, &vertices) {
                let station = cone_station(*cone, point);
                low = low.min(station);
                high = high.max(station);
            }
        }
        validate_cone_span(
            local_to_original[&region],
            low,
            high,
            cone.half_angle,
            tolerance,
        )?;
    }
    let cylinders = build_cylinder_surfaces(&cylinder_carriers, &chains, &vertices, tolerance)
        .map_err(|(region, message)| HybridBuildFailure::Region {
            original: local_to_original[&region],
            message,
        })?;
    let cones = build_cone_surfaces(&cone_carriers, &cylinders, &chains, &vertices, tolerance)
        .map_err(|(region, message)| HybridBuildFailure::Region {
            original: local_to_original[&region],
            message,
        })?;
    sort_region_cycles(
        &mut boundaries,
        &planes,
        &cylinders,
        &cones,
        &chains,
        &vertices,
    );

    // A faceted seed gives us a public, type-safe TopologyArena without
    // depending directly on SlotMap.  The seed is also an independent closed-
    // manifold precondition.  Its topology is cleared before reconstruction.
    let seed = mesh_to_faceted_brep(positions, Some(indices), options.weld_tolerance)?;
    let mut arena = TopologyArena::from_brep(&seed)?;
    arena.vertices.clear();
    arena.edges.clear();
    arena.coedges.clear();
    arena.loops.clear();
    arena.faces.clear();
    arena.shells.clear();
    arena.wire_solid_id = 1;
    arena.genus = 0;
    let mut ids = Ids(2);
    let mut arena_vertices = HashMap::new();

    let mut built_edges = Vec::with_capacity(chains.len());
    for (chain_index, chain) in chains.iter().enumerate() {
        let built = match build_chain_edge(
            chain_index,
            chain,
            &vertices,
            &planes,
            &cylinders,
            &cones,
            tolerance,
            &mut arena,
            &mut arena_vertices,
            &mut ids,
        ) {
            Ok(built) => built,
            Err(message) => {
                let local = [chain.adjacent.0, chain.adjacent.1]
                    .into_iter()
                    .find(|region| cones.contains_key(region))
                    .or_else(|| {
                        [chain.adjacent.0, chain.adjacent.1]
                            .into_iter()
                            .find(|region| cylinders.contains_key(region))
                    });
                if let Some(local) = local {
                    return Err(HybridBuildFailure::Region {
                        original: local_to_original[&local],
                        message,
                    });
                }
                return Err(HybridBuildFailure::Global(message));
            }
        };
        built_edges.push(built);
    }

    let mut faces = Vec::<FaceId>::new();
    for region in 0..local.count as u32 {
        let face = if let Some(plane) = planes.get(&region) {
            build_plane_face(
                region,
                *plane,
                &boundaries[region as usize],
                &chains,
                &built_edges,
                &vertices,
                &mut arena,
                &mut ids,
            )?
        } else if let Some(cylinder) = cylinders.get(&region) {
            build_cylinder_face(
                region,
                cylinder,
                &boundaries[region as usize],
                &built_edges,
                &mut arena,
                &mut ids,
            )
            .map_err(|message| HybridBuildFailure::Region {
                original: local_to_original[&region],
                message,
            })?
        } else {
            build_cone_face(
                region,
                &cones[&region],
                &boundaries[region as usize],
                &built_edges,
                &mut arena,
                &mut ids,
            )
            .map_err(|message| HybridBuildFailure::Region {
                original: local_to_original[&region],
                message,
            })?
        };
        faces.push(face);
    }
    arena.shells.insert(ArenaShell {
        wire_id: ids.next(),
        faces,
    });

    let mut solid = arena.to_brep()?;
    let vertex_count = solid.vertices.len() as i64;
    let edge_count = solid.edges.len() as i64;
    let face_count = solid
        .shells
        .iter()
        .map(|shell| shell.faces.len())
        .sum::<usize>() as i64;
    let extra_loops = solid
        .shells
        .iter()
        .flat_map(|shell| &shell.faces)
        .map(|face| face.loops.len().saturating_sub(1) as i64)
        .sum::<i64>();
    let numerator = 2 - (vertex_count - edge_count + face_count - extra_loops);
    if numerator < 0 || numerator % 2 != 0 {
        return Err(format!(
            "hybrid BREP: non-manifold Euler accounting (2-(V-E+F-H)={numerator})"
        )
        .into());
    }
    solid.genus = numerator / 2;
    let issues = solid.validate();
    if !issues.is_empty() {
        return Err(format!("hybrid BREP: topology validation failed: {issues:?}").into());
    }
    let volume = solid_signed_volume(&solid)?;
    if !volume.is_finite() || volume <= 0.0 {
        return Err(format!(
            "hybrid BREP: reconstructed shell is inverted or empty (signed volume {volume:.6e})"
        )
        .into());
    }
    Ok(HybridBrepOutput { solid, stats })
}

fn validate_cone_span(
    original: u32,
    low: f64,
    high: f64,
    half_angle: f64,
    tolerance: f64,
) -> Result<(), HybridBuildFailure> {
    if !low.is_finite() || high - low <= tolerance || low * half_angle.tan() <= tolerance {
        return Err(HybridBuildFailure::Region {
            original,
            message: "pointed or zero-span cone requires unsupported pole topology".into(),
        });
    }
    Ok(())
}

fn mesh_diagonal(vertices: &[Vec3]) -> f64 {
    let mut low = vertices[0];
    let mut high = vertices[0];
    for &point in &vertices[1..] {
        low.x = low.x.min(point.x);
        low.y = low.y.min(point.y);
        low.z = low.z.min(point.z);
        high.x = high.x.max(point.x);
        high.y = high.y.max(point.y);
        high.z = high.z.max(point.z);
    }
    high.sub(low).length()
}

type EdgeKey = (usize, usize);

fn edge_key(a: usize, b: usize) -> EdgeKey {
    if a < b {
        (a, b)
    } else {
        (b, a)
    }
}

fn edge_incidence(triangles: &[[usize; 3]]) -> Result<BTreeMap<EdgeKey, Vec<usize>>, String> {
    let mut incidence = BTreeMap::<EdgeKey, Vec<usize>>::new();
    for (triangle_id, triangle) in triangles.iter().enumerate() {
        if triangle[0] == triangle[1] || triangle[1] == triangle[2] || triangle[2] == triangle[0] {
            return Err("hybrid BREP: repeated vertex in triangle".into());
        }
        for corner in 0..3 {
            incidence
                .entry(edge_key(triangle[corner], triangle[(corner + 1) % 3]))
                .or_default()
                .push(triangle_id);
        }
    }
    Ok(incidence)
}

fn carrier_cylinder(region: &MeshRegion) -> Result<Option<CylinderCarrier>, String> {
    let RegionCarrier::Cylinder {
        axis_point,
        axis_dir,
        radius,
        sense,
    } = region.carrier
    else {
        return Ok(None);
    };
    if !radius.is_finite() || radius <= 0.0 || !matches!(sense, -1 | 1) {
        return Ok(None);
    }
    Ok(Some(CylinderCarrier {
        origin: axis_point,
        axis: axis_dir.normalized()?,
        radius,
        sense,
    }))
}

fn carrier_cone(region: &MeshRegion) -> Result<Option<ConeCarrier>, String> {
    let RegionCarrier::Cone {
        apex,
        axis_dir,
        half_angle_rad,
        sense,
    } = region.carrier
    else {
        return Ok(None);
    };
    if !half_angle_rad.is_finite()
        || !(1.0e-6..std::f64::consts::FRAC_PI_2 - 1.0e-6).contains(&half_angle_rad)
        || !matches!(sense, -1 | 1)
    {
        return Ok(None);
    }
    Ok(Some(ConeCarrier {
        apex,
        axis: axis_dir.normalized()?,
        half_angle: half_angle_rad,
        sense,
    }))
}

fn edge_is_cylinder_ruling(cylinder: CylinderCarrier, a: Vec3, b: Vec3, tolerance: f64) -> bool {
    let ra = radial(cylinder, a);
    let rb = radial(cylinder, b);
    let Ok(ua) = ra.normalized() else {
        return false;
    };
    let Ok(ub) = rb.normalized() else {
        return false;
    };
    let chord = b.sub(a);
    let Ok(direction) = chord.normalized() else {
        return false;
    };
    (ra.length() - cylinder.radius).abs() <= tolerance
        && (rb.length() - cylinder.radius).abs() <= tolerance
        && ua.cross(ub).length() <= 0.01
        && direction.cross(cylinder.axis).length() <= 0.01
}

fn cylinder_plane_edge_supported(
    cylinder: CylinderCarrier,
    plane: PlaneInfo,
    a: Vec3,
    b: Vec3,
    tolerance: f64,
) -> bool {
    let Ok(normal) = plane.normal.normalized() else {
        return false;
    };
    if a.sub(plane.origin).dot(normal).abs() > tolerance
        || b.sub(plane.origin).dot(normal).abs() > tolerance
    {
        return false;
    }
    let ra = radial(cylinder, a);
    let rb = radial(cylinder, b);
    let looks_like_ruling = (station(cylinder, a) - station(cylinder, b)).abs() > tolerance
        && ra
            .normalized()
            .ok()
            .zip(rb.normalized().ok())
            .is_some_and(|(ua, ub)| ua.cross(ub).length() <= 0.05);
    if looks_like_ruling {
        // The exact axial line is contained by the plane only when its
        // direction is tangent to that plane.
        return normal.dot(cylinder.axis).abs() <= 0.01;
    }
    let on_ring = (station(cylinder, a) - station(cylinder, b)).abs() <= tolerance * 4.0;
    // A constant-station circle lies in the neighboring plane only when that
    // plane is perpendicular to the cylinder axis.
    on_ring && normal.cross(cylinder.axis).length() <= 0.01
}

fn cone_station(cone: ConeCarrier, point: Vec3) -> f64 {
    point.sub(cone.apex).dot(cone.axis)
}

fn cone_radial(cone: ConeCarrier, point: Vec3) -> Vec3 {
    let delta = point.sub(cone.apex);
    delta.sub(cone.axis.scale(delta.dot(cone.axis)))
}

fn point_on_cone(cone: ConeCarrier, point: Vec3, tolerance: f64) -> bool {
    let station = cone_station(cone, point);
    station >= -tolerance
        && (cone_radial(cone, point).length() - station * cone.half_angle.tan()).abs() <= tolerance
}

fn edge_is_cone_ruling(cone: ConeCarrier, a: Vec3, b: Vec3, tolerance: f64) -> bool {
    if !point_on_cone(cone, a, tolerance) || !point_on_cone(cone, b, tolerance) {
        return false;
    }
    let ra = cone_radial(cone, a);
    let rb = cone_radial(cone, b);
    let (Ok(ua), Ok(ub), Ok(direction)) = (ra.normalized(), rb.normalized(), b.sub(a).normalized())
    else {
        return false;
    };
    let expected = cone
        .axis
        .add(ua.scale(cone.half_angle.tan()))
        .normalized()
        .unwrap_or(cone.axis);
    ua.cross(ub).length() <= 0.01
        && direction
            .cross(expected)
            .length()
            .min(direction.cross(expected.scale(-1.0)).length())
            <= 0.01
}

fn cone_plane_edge_supported(
    cone: ConeCarrier,
    plane: PlaneInfo,
    a: Vec3,
    b: Vec3,
    tolerance: f64,
) -> bool {
    let Ok(normal) = plane.normal.normalized() else {
        return false;
    };
    let da = a.sub(plane.origin).dot(normal).abs();
    let db = b.sub(plane.origin).dot(normal).abs();
    if da > tolerance || db > tolerance {
        return false;
    }
    let ra = cone_radial(cone, a);
    let rb = cone_radial(cone, b);
    let looks_like_ruling = (cone_station(cone, a) - cone_station(cone, b)).abs() > tolerance
        && ra
            .normalized()
            .ok()
            .zip(rb.normalized().ok())
            .is_some_and(|(ua, ub)| ua.cross(ub).length() <= 0.05);
    if looks_like_ruling {
        let apex_in_plane = cone.apex.sub(plane.origin).dot(normal).abs() <= tolerance * 4.0;
        if apex_in_plane && normal.dot(b.sub(a)).abs() <= tolerance {
            return true;
        }
    }
    let on_ring = (cone_station(cone, a) - cone_station(cone, b)).abs() <= tolerance * 4.0;
    if on_ring && normal.cross(cone.axis).length() <= 0.01 {
        return true;
    }
    // A finite oblique trim is represented later as the exact projective map
    // of only its reference-circle arc. This also covers local hyperbolic and
    // parabolic branches without the mixed weights of a full-circle map.
    true
}

fn cylinder_cone_edge_supported(
    cylinder: CylinderCarrier,
    cone: ConeCarrier,
    a: Vec3,
    b: Vec3,
    tolerance: f64,
) -> bool {
    if cylinder.axis.cross(cone.axis).length() > 0.01 {
        return false;
    }
    let cylinder_station_spread = (station(cylinder, a) - station(cylinder, b)).abs();
    let cone_station_spread = (cone_station(cone, a) - cone_station(cone, b)).abs();
    cylinder_station_spread <= tolerance * 4.0
        && cone_station_spread <= tolerance * 4.0
        && (radial(cylinder, a).length() - cylinder.radius).abs() <= tolerance * 4.0
        && (radial(cylinder, b).length() - cylinder.radius).abs() <= tolerance * 4.0
        && point_on_cone(cone, a, tolerance * 4.0)
        && point_on_cone(cone, b, tolerance * 4.0)
}

fn facet_plane(triangle: [usize; 3], vertices: &[Vec3]) -> Result<PlaneInfo, String> {
    let a = vertices[triangle[0]];
    let b = vertices[triangle[1]];
    let c = vertices[triangle[2]];
    let normal = b.sub(a).cross(c.sub(a)).normalized().map_err(|_| {
        "hybrid BREP: cannot create a planar face from a degenerate source triangle".to_owned()
    })?;
    Ok(PlaneInfo { origin: a, normal })
}

fn facets_are_coplanar(
    first: [usize; 3],
    second: [usize; 3],
    vertices: &[Vec3],
    tolerance: f64,
) -> Result<bool, String> {
    let first_plane = facet_plane(first, vertices)?;
    let second_plane = facet_plane(second, vertices)?;
    // Do not turn a shallow crease into a plane merely because a distance
    // tolerance happens to cover its short triangles. The distance checks
    // below cover f32 STL coordinate noise; this angular check remains
    // deliberately close to exact coplanarity.
    if first_plane.normal.dot(second_plane.normal) <= 0.0
        || first_plane.normal.cross(second_plane.normal).length() > 1.0e-8
    {
        return Ok(false);
    }
    Ok(first.iter().chain(second.iter()).all(|&vertex| {
        let point = vertices[vertex];
        point.sub(first_plane.origin).dot(first_plane.normal).abs() <= tolerance
            && point
                .sub(second_plane.origin)
                .dot(second_plane.normal)
                .abs()
                <= tolerance
    }))
}

/// Replace unsupported regions with safely merged coplanar patches and
/// recursively demote revolved regions that cannot share exact edges with them.
fn local_regions(
    regions: &[MeshRegion],
    triangle_regions: &[u32],
    triangles: &[[usize; 3]],
    vertices: &[Vec3],
    incidence: &BTreeMap<EdgeKey, Vec<usize>>,
    tolerance: f64,
    forced_demotions: &BTreeSet<u32>,
) -> Result<LocalRegions, String> {
    let mut original_cylinders = HashMap::<u32, CylinderCarrier>::new();
    let mut original_cones = HashMap::<u32, ConeCarrier>::new();
    let mut original_planes = HashMap::<u32, PlaneInfo>::new();
    let mut retained = BTreeSet::<u32>::new();
    for region in regions {
        if let RegionCarrier::Plane { origin, normal } = region.carrier {
            original_planes.insert(region.id, PlaneInfo { origin, normal });
        }
        if let Some(cylinder) = carrier_cylinder(region)? {
            original_cylinders.insert(region.id, cylinder);
            if !forced_demotions.contains(&region.id) {
                retained.insert(region.id);
            }
        }
        if let Some(cone) = carrier_cone(region)? {
            original_cones.insert(region.id, cone);
            if !forced_demotions.contains(&region.id) {
                retained.insert(region.id);
            }
        }
    }

    loop {
        let mut demote = BTreeSet::new();
        for (&(a, b), incident) in incidence {
            let first = triangle_regions[incident[0]];
            let second = triangle_regions[incident[1]];
            if first == second {
                continue;
            }
            for (candidate, neighbor) in [(first, second), (second, first)] {
                if !retained.contains(&candidate) {
                    continue;
                }
                if retained.contains(&neighbor) {
                    match (
                        original_cylinders.get(&candidate),
                        original_cones.get(&candidate),
                        original_cylinders.get(&neighbor),
                        original_cones.get(&neighbor),
                    ) {
                        (Some(&cylinder), None, None, Some(&cone))
                        | (None, Some(&cone), Some(&cylinder), None) => {
                            if !cylinder_cone_edge_supported(
                                cylinder,
                                cone,
                                vertices[a],
                                vertices[b],
                                tolerance,
                            ) {
                                demote.insert(candidate);
                                demote.insert(neighbor);
                            }
                        }
                        _ => {
                            // Same-kind revolve intersections and cone/cone
                            // intersections are not yet proven exact.
                            demote.insert(candidate);
                            demote.insert(neighbor);
                        }
                    }
                    continue;
                }
                // A cone may have been locally demoted because one of its
                // other trims is an unsupported hyperbolic/parabolic section.
                // Its shared circle with a neighboring cylinder remains exact,
                // so that cylinder need not be demoted with it.
                if let (Some(&cylinder), Some(&cone)) = (
                    original_cylinders.get(&candidate),
                    original_cones.get(&neighbor),
                ) {
                    if cylinder_cone_edge_supported(
                        cylinder,
                        cone,
                        vertices[a],
                        vertices[b],
                        tolerance,
                    ) {
                        continue;
                    }
                }
                if let Some(&plane) = original_planes.get(&neighbor) {
                    let supported = if let Some(&cylinder) = original_cylinders.get(&candidate) {
                        cylinder_plane_edge_supported(
                            cylinder,
                            plane,
                            vertices[a],
                            vertices[b],
                            tolerance,
                        )
                    } else {
                        cone_plane_edge_supported(
                            original_cones[&candidate],
                            plane,
                            vertices[a],
                            vertices[b],
                            tolerance,
                        )
                    };
                    if !supported {
                        demote.insert(candidate);
                    }
                    continue;
                }
                let supported = if let Some(&cylinder) = original_cylinders.get(&candidate) {
                    edge_is_cylinder_ruling(cylinder, vertices[a], vertices[b], tolerance)
                } else {
                    edge_is_cone_ruling(
                        original_cones[&candidate],
                        vertices[a],
                        vertices[b],
                        tolerance,
                    )
                };
                if !supported {
                    demote.insert(candidate);
                }
            }
        }
        if demote.is_empty() {
            break;
        }
        for region in demote {
            retained.remove(&region);
        }
    }

    let demoted = original_cylinders
        .keys()
        .chain(original_cones.keys())
        .copied()
        .filter(|id| !retained.contains(id))
        .collect::<BTreeSet<_>>();
    let mut planes = HashMap::new();
    let mut cylinders = HashMap::new();
    let mut cones = HashMap::new();
    let mut original_to_local = HashMap::<u32, u32>::new();
    let mut local_to_original = HashMap::<u32, u32>::new();
    let mut next = 0_u32;
    for region in regions {
        match region.carrier {
            RegionCarrier::Plane { origin, normal } => {
                original_to_local.insert(region.id, next);
                local_to_original.insert(next, region.id);
                planes.insert(next, PlaneInfo { origin, normal });
                next += 1;
            }
            RegionCarrier::Cylinder { .. } if retained.contains(&region.id) => {
                original_to_local.insert(region.id, next);
                local_to_original.insert(next, region.id);
                cylinders.insert(next, original_cylinders[&region.id]);
                next += 1;
            }
            RegionCarrier::Cone { .. } if retained.contains(&region.id) => {
                original_to_local.insert(region.id, next);
                local_to_original.insert(next, region.id);
                cones.insert(next, original_cones[&region.id]);
                next += 1;
            }
            _ => {}
        }
    }

    let fallback = triangle_regions
        .iter()
        .map(|original| !original_to_local.contains_key(original))
        .collect::<Vec<_>>();
    let mut neighbors = vec![Vec::<usize>::new(); triangles.len()];
    for incident in incidence.values() {
        let first = incident[0];
        let second = incident[1];
        if fallback[first] && fallback[second] {
            neighbors[first].push(second);
            neighbors[second].push(first);
        }
    }

    // Check every member against one fixed seed carrier. Pairwise-only
    // transitive merging could otherwise accumulate a long sequence of tiny
    // creases into a patch that is not actually planar.
    let mut fallback_component = vec![None; triangles.len()];
    let mut fallback_components = 0_usize;
    for seed in 0..triangles.len() {
        if !fallback[seed] || fallback_component[seed].is_some() {
            continue;
        }
        let component = fallback_components;
        fallback_components += 1;
        fallback_component[seed] = Some(component);
        let mut queue = VecDeque::from([seed]);
        while let Some(current) = queue.pop_front() {
            for &neighbor in &neighbors[current] {
                if fallback_component[neighbor].is_none()
                    && facets_are_coplanar(
                        triangles[seed],
                        triangles[neighbor],
                        vertices,
                        tolerance,
                    )?
                {
                    fallback_component[neighbor] = Some(component);
                    queue.push_back(neighbor);
                }
            }
        }
    }

    let mut fallback_to_local = HashMap::<usize, u32>::new();
    let mut local_triangle_regions = Vec::with_capacity(triangles.len());
    for (triangle_id, &original) in triangle_regions.iter().enumerate() {
        if let Some(&local) = original_to_local.get(&original) {
            local_triangle_regions.push(local);
        } else {
            let component = fallback_component[triangle_id]
                .ok_or_else(|| "hybrid BREP: fallback component was not assigned".to_owned())?;
            let local = if let Some(&local) = fallback_to_local.get(&component) {
                local
            } else {
                let local = next;
                next += 1;
                planes.insert(local, facet_plane(triangles[triangle_id], vertices)?);
                fallback_to_local.insert(component, local);
                local
            };
            local_triangle_regions.push(local);
        }
    }
    let analytic_plane_triangles = triangle_regions
        .iter()
        .filter(|region| original_planes.contains_key(region))
        .count();
    let analytic_cylinder_triangles = triangle_regions
        .iter()
        .filter(|region| retained.contains(region) && original_cylinders.contains_key(region))
        .count();
    let analytic_cone_triangles = triangle_regions
        .iter()
        .filter(|region| retained.contains(region) && original_cones.contains_key(region))
        .count();
    let demoted_triangles = triangle_regions
        .iter()
        .filter(|region| demoted.contains(region))
        .count();
    let faceted_triangles = triangles.len()
        - analytic_plane_triangles
        - analytic_cylinder_triangles
        - analytic_cone_triangles;
    let stats = HybridBrepStats {
        total_faces: next as usize,
        analytic_plane_faces: original_planes.len(),
        analytic_plane_triangles,
        analytic_cylinder_faces: original_cylinders
            .keys()
            .filter(|id| retained.contains(id))
            .count(),
        analytic_cylinder_triangles,
        analytic_cone_faces: original_cones
            .keys()
            .filter(|id| retained.contains(id))
            .count(),
        analytic_cone_triangles,
        faceted_faces: fallback_to_local.len(),
        faceted_triangles,
        demoted_regions: demoted.len(),
        demoted_triangles,
        segmentation_plane_regions: original_planes.len(),
        segmentation_cylinder_regions: original_cylinders.len(),
        segmentation_cone_regions: original_cones.len(),
        segmentation_unsupported_regions: regions.len()
            - original_planes.len()
            - original_cylinders.len()
            - original_cones.len(),
    };
    Ok(LocalRegions {
        triangle_regions: local_triangle_regions,
        planes,
        cylinders,
        cones,
        local_to_original,
        count: next as usize,
        stats,
    })
}

fn region_boundaries(
    triangles: &[[usize; 3]],
    regions: &[u32],
    region_count: usize,
    incidence: &BTreeMap<EdgeKey, Vec<usize>>,
) -> Result<(Vec<Chain>, Vec<RegionBoundary>), String> {
    let mut chains = Vec::<Chain>::new();
    let mut chain_by_edges = BTreeMap::<Vec<EdgeKey>, usize>::new();
    let mut boundaries = Vec::with_capacity(region_count);
    for region in 0..region_count as u32 {
        let mut outgoing = BTreeMap::<usize, (usize, u32)>::new();
        for (triangle_id, triangle) in triangles.iter().enumerate() {
            if regions[triangle_id] != region {
                continue;
            }
            for corner in 0..3 {
                let a = triangle[corner];
                let b = triangle[(corner + 1) % 3];
                let incident = &incidence[&edge_key(a, b)];
                let other = incident.iter().copied().find(|&id| id != triangle_id);
                let neighbor = other.map(|id| regions[id]).unwrap_or(u32::MAX);
                if neighbor != region && outgoing.insert(a, (b, neighbor)).is_some() {
                    return Err(format!(
                        "hybrid BREP: region {region} boundary pinches at vertex {a}"
                    ));
                }
            }
        }
        let starts = outgoing.keys().copied().collect::<Vec<_>>();
        let mut visited = BTreeSet::new();
        let mut cycles = Vec::new();
        for start in starts {
            if visited.contains(&start) {
                continue;
            }
            let mut steps = Vec::<(usize, usize, u32)>::new();
            let mut current = start;
            loop {
                if !visited.insert(current) {
                    return Err(format!(
                        "hybrid BREP: region {region} boundary self-crosses"
                    ));
                }
                let &(next, neighbor) = outgoing
                    .get(&current)
                    .ok_or_else(|| format!("hybrid BREP: region {region} boundary walk is open"))?;
                steps.push((current, next, neighbor));
                current = next;
                if current == start {
                    break;
                }
            }
            let all_one_neighbor = steps.iter().all(|step| step.2 == steps[0].2);
            let mut groups = Vec::<Vec<(usize, usize, u32)>>::new();
            if all_one_neighbor {
                groups.push(steps);
            } else {
                let cut = (0..steps.len())
                    .find(|&index| {
                        steps[index].2 != steps[(index + steps.len() - 1) % steps.len()].2
                    })
                    .unwrap();
                steps.rotate_left(cut);
                let mut begin = 0;
                while begin < steps.len() {
                    let neighbor = steps[begin].2;
                    let mut end = begin + 1;
                    while end < steps.len() && steps[end].2 == neighbor {
                        end += 1;
                    }
                    groups.push(steps[begin..end].to_vec());
                    begin = end;
                }
            }
            let mut traversals = Vec::new();
            for group in groups {
                // `steps` is one boundary cycle. It may be only one of several
                // cycles owned by this region, so comparing against the
                // region-wide outgoing-edge count would misclassify a hole ring.
                let closed = all_one_neighbor;
                let mut chain_vertices = vec![group[0].0];
                chain_vertices.extend(group.iter().map(|step| step.1));
                if closed {
                    chain_vertices.pop();
                }
                let mut key = group
                    .iter()
                    .map(|step| edge_key(step.0, step.1))
                    .collect::<Vec<_>>();
                key.sort_unstable();
                let chain_id = if let Some(&existing) = chain_by_edges.get(&key) {
                    existing
                } else {
                    let id = chains.len();
                    let adjacent = (region.min(group[0].2), region.max(group[0].2));
                    chains.push(Chain {
                        vertices: chain_vertices.clone(),
                        closed,
                        adjacent,
                    });
                    chain_by_edges.insert(key, id);
                    id
                };
                let stored = &chains[chain_id];
                let a = chain_vertices[0];
                let b = chain_vertices[1];
                let position = stored
                    .vertices
                    .iter()
                    .position(|&vertex| vertex == a)
                    .ok_or("hybrid BREP: shared chain start is missing")?;
                let forward = if stored.closed {
                    stored.vertices[(position + 1) % stored.vertices.len()] == b
                } else if position + 1 < stored.vertices.len() && stored.vertices[position + 1] == b
                {
                    true
                } else if position > 0 && stored.vertices[position - 1] == b {
                    false
                } else {
                    return Err("hybrid BREP: shared chain directions disagree".into());
                };
                traversals.push(Traversal {
                    chain: chain_id,
                    forward,
                });
            }
            cycles.push(traversals);
        }
        boundaries.push(RegionBoundary { cycles });
    }
    Ok((chains, boundaries))
}

fn chain_points(chain: &Chain, vertices: &[Vec3]) -> Vec<Vec3> {
    chain.vertices.iter().map(|&id| vertices[id]).collect()
}

fn snap_cylinder_axes(
    cylinders: &mut HashMap<u32, CylinderCarrier>,
    planes: &HashMap<u32, PlaneInfo>,
    chains: &[Chain],
    vertices: &[Vec3],
    tolerance: f64,
) -> Result<(), (u32, String)> {
    for (&region, cylinder) in cylinders.iter_mut() {
        let mut candidates = Vec::new();
        for chain in chains.iter().filter(|chain| {
            chain.closed && (chain.adjacent.0 == region || chain.adjacent.1 == region)
        }) {
            let neighbor = if chain.adjacent.0 == region {
                chain.adjacent.1
            } else {
                chain.adjacent.0
            };
            let Some(plane) = planes.get(&neighbor) else {
                continue;
            };
            if plane.normal.cross(cylinder.axis).length() > 0.01 {
                continue;
            }
            let points = chain_points(chain, vertices);
            let spread = points
                .iter()
                .map(|point| point.sub(plane.origin).dot(plane.normal).abs())
                .fold(0.0_f64, f64::max);
            if spread <= tolerance {
                let axis = if plane.normal.dot(cylinder.axis) >= 0.0 {
                    plane.normal
                } else {
                    plane.normal.scale(-1.0)
                };
                candidates.push(axis.normalized().map_err(|error| (region, error))?);
            }
        }
        if let Some(&axis) = candidates.first() {
            // Independent plane fits to rounded STL rings need not have
            // bit-identical normals. Bound their disagreement by the
            // positional error it induces around this cylinder's rim,
            // using the same tolerance as the shared topology checks.
            if candidates
                .iter()
                .any(|other| cylinder.radius * other.cross(axis).length() > tolerance)
            {
                return Err((
                    region,
                    format!("hybrid BREP: cylinder {region} has inconsistent ring normals"),
                ));
            }
            cylinder.axis = axis;
        }
    }
    Ok(())
}

fn snap_cone_carriers(
    cones: &mut HashMap<u32, ConeCarrier>,
    cylinders: &HashMap<u32, CylinderCarrier>,
    planes: &HashMap<u32, PlaneInfo>,
    chains: &[Chain],
    vertices: &[Vec3],
    tolerance: f64,
) -> Result<(), (u32, String)> {
    for (&region, cone) in cones.iter_mut() {
        let relevant = chains
            .iter()
            .filter(|chain| chain.adjacent.0 == region || chain.adjacent.1 == region)
            .collect::<Vec<_>>();
        let mut snapped = None::<(Vec3, Vec3)>;

        // A cylinder/cone join is a shared circle. The cylinder fixes the
        // coaxial frame and the cone apex then follows from its opening angle.
        for chain in &relevant {
            let neighbor = if chain.adjacent.0 == region {
                chain.adjacent.1
            } else {
                chain.adjacent.0
            };
            let Some(cylinder) = cylinders.get(&neighbor).copied() else {
                continue;
            };
            let mut axis = cylinder.axis;
            if axis.dot(cone.axis) < 0.0 {
                axis = axis.scale(-1.0);
            }
            let points = chain_points(chain, vertices);
            let center_station = points
                .iter()
                .map(|&point| point.sub(cylinder.origin).dot(axis))
                .sum::<f64>()
                / points.len() as f64;
            let center = cylinder.origin.add(axis.scale(center_station));

            let mut rings = Vec::<(f64, f64)>::new();
            for boundary_chain in &relevant {
                let boundary_points = chain_points(boundary_chain, vertices);
                let stations = boundary_points
                    .iter()
                    .map(|&point| point.sub(cylinder.origin).dot(axis))
                    .collect::<Vec<_>>();
                let low_station = stations.iter().copied().fold(f64::INFINITY, f64::min);
                let high_station = stations.iter().copied().fold(f64::NEG_INFINITY, f64::max);
                if high_station - low_station > tolerance * 4.0 {
                    continue;
                }
                let station = stations.iter().sum::<f64>() / stations.len() as f64;
                let neighbor = if boundary_chain.adjacent.0 == region {
                    boundary_chain.adjacent.1
                } else {
                    boundary_chain.adjacent.0
                };
                let radius = if cylinders.contains_key(&neighbor) {
                    cylinder.radius
                } else {
                    boundary_points
                        .iter()
                        .map(|&point| {
                            let delta = point.sub(cylinder.origin);
                            delta.sub(axis.scale(delta.dot(axis))).length()
                        })
                        .sum::<f64>()
                        / boundary_points.len() as f64
                };
                rings.push((station, radius));
            }
            rings.sort_by(|a, b| a.0.total_cmp(&b.0));
            if let (Some(&(s0, r0)), Some(&(s1, r1))) = (rings.first(), rings.last()) {
                let slope = (r1 - r0) / (s1 - s0);
                if slope.is_finite() && slope > 1.0e-6 {
                    let apex_station = s0 - r0 / slope;
                    cone.apex = cylinder.origin.add(axis.scale(apex_station));
                    cone.axis = axis;
                    cone.half_angle = slope.atan();
                    snapped = Some((cone.apex, axis));
                    break;
                }
            }

            // Short cone patches fit their apex poorly from tessellated
            // normals. Their two side planes, however, contain the exact cone
            // apex. Intersect those planes with the coaxial line fixed by the
            // adjacent cylinder, then derive the half-angle from the common
            // cylinder/cone ring.
            let mut apex_stations = Vec::new();
            for side_chain in &relevant {
                let side_neighbor = if side_chain.adjacent.0 == region {
                    side_chain.adjacent.1
                } else {
                    side_chain.adjacent.0
                };
                let Some(side_plane) = planes.get(&side_neighbor) else {
                    continue;
                };
                let side_points = chain_points(side_chain, vertices);
                let side_spread = side_points
                    .iter()
                    .map(|&point| point.sub(cylinder.origin).dot(axis))
                    .fold((f64::INFINITY, f64::NEG_INFINITY), |(low, high), value| {
                        (low.min(value), high.max(value))
                    });
                if side_spread.1 - side_spread.0 <= tolerance * 4.0 {
                    continue;
                }
                let normal = side_plane
                    .normal
                    .normalized()
                    .map_err(|error| (region, error))?;
                let denominator = normal.dot(axis);
                if denominator.abs() <= 1.0e-8 {
                    continue;
                }
                apex_stations
                    .push(normal.dot(side_plane.origin.sub(cylinder.origin)) / denominator);
            }
            let side_solution = if !apex_stations.is_empty() {
                let apex_station = apex_stations.iter().sum::<f64>() / apex_stations.len() as f64;
                let axial_radius = center_station - apex_station;
                if axial_radius <= tolerance {
                    None
                } else {
                    Some((
                        cylinder.origin.add(axis.scale(apex_station)),
                        (cylinder.radius / axial_radius).atan(),
                    ))
                }
            } else {
                None
            };
            // Accept a side-plane apex only when it corroborates the fitted
            // opening angle. Otherwise those planes cut the cone in genuine
            // oblique conics and do not contain its apex.
            let (apex, half_angle) = match side_solution {
                Some((apex, angle)) if (angle - cone.half_angle).abs() <= 0.01 => (apex, angle),
                _ => (
                    center.sub(axis.scale(cylinder.radius / cone.half_angle.tan())),
                    cone.half_angle,
                ),
            };
            cone.half_angle = half_angle;
            snapped = Some((apex, axis));
            break;
        }

        // Otherwise a ring in an axis-normal plane supplies an equally exact
        // frame. This covers truncated cones without an adjacent cylinder.
        if snapped.is_none() {
            for chain in relevant.iter().filter(|chain| chain.closed) {
                let neighbor = if chain.adjacent.0 == region {
                    chain.adjacent.1
                } else {
                    chain.adjacent.0
                };
                let Some(plane) = planes.get(&neighbor).copied() else {
                    continue;
                };
                let mut axis = plane.normal.normalized().map_err(|error| (region, error))?;
                if axis.dot(cone.axis) < 0.0 {
                    axis = axis.scale(-1.0);
                }
                if axis.cross(cone.axis).length() > 0.01 {
                    continue;
                }
                let points = chain_points(chain, vertices);
                let mean = points
                    .iter()
                    .copied()
                    .fold(Vec3::new(0.0, 0.0, 0.0), Vec3::add)
                    .scale(1.0 / points.len() as f64);
                let center = mean.sub(axis.scale(mean.sub(plane.origin).dot(axis)));
                let radius = points
                    .iter()
                    .map(|&point| {
                        let delta = point.sub(center);
                        delta.sub(axis.scale(delta.dot(axis))).length()
                    })
                    .sum::<f64>()
                    / points.len() as f64;
                if radius > tolerance {
                    let apex = center.sub(axis.scale(radius / cone.half_angle.tan()));
                    snapped = Some((apex, axis));
                    break;
                }
            }
        }

        if let Some((apex, axis)) = snapped {
            cone.apex = apex;
            cone.axis = axis;
        }
    }
    Ok(())
}

fn station(carrier: CylinderCarrier, point: Vec3) -> f64 {
    point.sub(carrier.origin).dot(carrier.axis)
}

fn radial(carrier: CylinderCarrier, point: Vec3) -> Vec3 {
    let delta = point.sub(carrier.origin);
    delta.sub(carrier.axis.scale(delta.dot(carrier.axis)))
}

fn build_cylinder_surfaces(
    carriers: &HashMap<u32, CylinderCarrier>,
    chains: &[Chain],
    vertices: &[Vec3],
    tolerance: f64,
) -> Result<HashMap<u32, CylinderInfo>, (u32, String)> {
    let mut result = HashMap::new();
    for (&region, &carrier) in carriers {
        let relevant = chains
            .iter()
            .filter(|chain| chain.adjacent.0 == region || chain.adjacent.1 == region)
            .collect::<Vec<_>>();
        let mut low = f64::INFINITY;
        let mut high = f64::NEG_INFINITY;
        for chain in &relevant {
            for point in chain_points(chain, vertices) {
                let value = station(carrier, point);
                low = low.min(value);
                high = high.max(value);
            }
        }
        if !low.is_finite() || high - low <= tolerance {
            return Err((
                region,
                format!("hybrid BREP: cylinder {region} has no finite axial span"),
            ));
        }
        let mut seam = None;
        for chain in &relevant {
            if chain.closed {
                continue;
            }
            let points = chain_points(chain, vertices);
            let s0 = station(carrier, points[0]);
            let s1 = station(carrier, *points.last().unwrap());
            let r0 = radial(carrier, points[0]);
            let r1 = radial(carrier, *points.last().unwrap());
            if (s1 - s0).abs() > tolerance && r0.cross(r1).length() <= tolerance * carrier.radius {
                seam = Some(r0.normalized().map_err(|error| (region, error))?);
                break;
            }
        }
        let x_axis = match seam {
            Some(axis) => axis,
            None => carrier
                .axis
                .perpendicular()
                .map_err(|error| (region, error))?,
        };
        let y_axis = carrier
            .axis
            .cross(x_axis)
            .normalized()
            .map_err(|error| (region, error))?;
        let base = carrier.origin.add(carrier.axis.scale(low));
        let profile = make_line(
            base.add(x_axis.scale(carrier.radius)),
            base.add(carrier.axis.scale(high - low))
                .add(x_axis.scale(carrier.radius)),
        )
        .map_err(|error| (region, error))?;
        let surface =
            make_revolution(base, carrier.axis, &profile, TAU).map_err(|error| (region, error))?;
        result.insert(
            region,
            CylinderInfo {
                carrier,
                base_station: low,
                height: high - low,
                x_axis,
                y_axis,
                surface,
            },
        );
    }
    Ok(result)
}

fn build_cone_surfaces(
    carriers: &HashMap<u32, ConeCarrier>,
    cylinders: &HashMap<u32, CylinderInfo>,
    chains: &[Chain],
    vertices: &[Vec3],
    tolerance: f64,
) -> Result<HashMap<u32, ConeInfo>, (u32, String)> {
    let mut result = HashMap::new();
    for (&region, &carrier) in carriers {
        let relevant = chains
            .iter()
            .filter(|chain| chain.adjacent.0 == region || chain.adjacent.1 == region)
            .collect::<Vec<_>>();
        let mut low = f64::INFINITY;
        let mut high = f64::NEG_INFINITY;
        for chain in &relevant {
            for point in chain_points(chain, vertices) {
                let value = cone_station(carrier, point);
                low = low.min(value);
                high = high.max(value);
            }
        }
        // Pointed cones need a degenerate pole edge, which this conservative
        // mixed builder deliberately does not synthesize.
        if !low.is_finite()
            || low * carrier.half_angle.tan() <= tolerance
            || high - low <= tolerance
        {
            return Err((
                region,
                format!("hybrid BREP: cone {region} is pointed or has no finite axial span"),
            ));
        }
        // Keep fitted STL boundary samples away from the NURBS domain ends.
        // The trim loops still define the face; this small carrier extension
        // merely permits exact plane/cone intersections through quantized end
        // vertices without evaluating outside the surface domain.
        let margin = ((high - low) * 1.0e-4).max(tolerance * 2.0);
        low = (low - margin).max(tolerance / carrier.half_angle.tan());
        high += margin;
        let mut seam = relevant.iter().find_map(|chain| {
            let neighbor = if chain.adjacent.0 == region {
                chain.adjacent.1
            } else {
                chain.adjacent.0
            };
            cylinders.get(&neighbor).map(|cylinder| cylinder.x_axis)
        });
        for chain in &relevant {
            if chain.closed {
                continue;
            }
            let points = chain_points(chain, vertices);
            let s0 = cone_station(carrier, points[0]);
            let s1 = cone_station(carrier, *points.last().unwrap());
            let r0 = cone_radial(carrier, points[0]);
            let r1 = cone_radial(carrier, *points.last().unwrap());
            if (s1 - s0).abs() > tolerance
                && r0.cross(r1).length() <= tolerance * r0.length().max(r1.length())
            {
                seam = Some(r0.normalized().map_err(|error| (region, error))?);
                break;
            }
        }
        let x_axis = match seam {
            Some(axis) => axis,
            None => carrier
                .axis
                .perpendicular()
                .map_err(|error| (region, error))?,
        };
        let y_axis = carrier
            .axis
            .cross(x_axis)
            .normalized()
            .map_err(|error| (region, error))?;
        let radius_low = low * carrier.half_angle.tan();
        let radius_high = high * carrier.half_angle.tan();
        let base = carrier.apex.add(carrier.axis.scale(low));
        let profile = make_line(
            base.add(x_axis.scale(radius_low)),
            base.add(carrier.axis.scale(high - low))
                .add(x_axis.scale(radius_high)),
        )
        .map_err(|error| (region, error))?;
        let surface =
            make_revolution(base, carrier.axis, &profile, TAU).map_err(|error| (region, error))?;
        result.insert(
            region,
            ConeInfo {
                carrier,
                base_station: low,
                height: high - low,
                x_axis,
                y_axis,
                surface,
            },
        );
    }
    Ok(result)
}

fn max_line_deviation(points: &[Vec3]) -> f64 {
    if points.len() < 3 {
        return 0.0;
    }
    let start = points[0];
    let chord = points[points.len() - 1].sub(start);
    let denominator = chord.dot(chord);
    points[1..points.len() - 1]
        .iter()
        .map(|point| {
            let delta = point.sub(start);
            let t = if denominator > 0.0 {
                (delta.dot(chord) / denominator).clamp(0.0, 1.0)
            } else {
                0.0
            };
            delta.sub(chord.scale(t)).length()
        })
        .fold(0.0, f64::max)
}

fn trim_curve(curve: &NurbsCurve, start: f64, end: f64) -> Result<NurbsCurve, String> {
    let domain = curve.domain()?;
    let mut trimmed = if end < domain[1] - 1.0e-12 {
        curve.split(end)?.0
    } else {
        curve.clone()
    };
    if start > domain[0] + 1.0e-12 {
        trimmed = trimmed.split(start)?.1;
    }
    Ok(trimmed)
}

fn cylinder_coordinates(info: &CylinderInfo, point: Vec3) -> [f64; 2] {
    let r = radial(info.carrier, point);
    let angle = r.dot(info.y_axis).atan2(r.dot(info.x_axis)).rem_euclid(TAU);
    // A rational quadratic circle is not angle-linear inside a knot span.
    // Use the public exact angle-to-NURBS mapping used by the kernel carrier.
    let u = circle_angle_to_parameter(4, TAU, angle);
    let v = (station(info.carrier, point) - info.base_station) / info.height;
    [u, v]
}

fn cone_coordinates(info: &ConeInfo, point: Vec3) -> [f64; 2] {
    let r = cone_radial(info.carrier, point);
    let angle = r.dot(info.y_axis).atan2(r.dot(info.x_axis)).rem_euclid(TAU);
    let u = circle_angle_to_parameter(4, TAU, angle);
    let v = (cone_station(info.carrier, point) - info.base_station) / info.height;
    [u, v]
}

fn exact_local_cone_plane_conic(
    info: &ConeInfo,
    plane: PlaneInfo,
    points: &[Vec3],
    tolerance: f64,
) -> Result<NurbsCurve, String> {
    let mut angles = points
        .iter()
        .map(|&point| {
            let radial = cone_radial(info.carrier, point);
            radial.dot(info.y_axis).atan2(radial.dot(info.x_axis))
        })
        .collect::<Vec<_>>();
    for index in 1..angles.len() {
        while angles[index] - angles[index - 1] > std::f64::consts::PI {
            angles[index] -= TAU;
        }
        while angles[index] - angles[index - 1] < -std::f64::consts::PI {
            angles[index] += TAU;
        }
    }
    if !monotone_parameter_walk(&angles, 1.0e-5) {
        return Err("hybrid BREP: cone/plane source angles are not monotone".into());
    }
    let low = angles[0].min(*angles.last().unwrap());
    let high = angles[0].max(*angles.last().unwrap());
    if high - low <= 1.0e-10 || high - low > TAU + 1.0e-10 {
        return Err("hybrid BREP: cone/plane source angle span is invalid".into());
    }

    // Map a finite exact reference-circle arc projectively through the cone
    // apex onto the cutting plane. Unlike mapping a full circle, this keeps
    // hyperbolic/parabolic branches one-signed over the actual trim interval.
    let reference_station = info.base_station + 0.5 * info.height;
    let reference_radius = reference_station * info.carrier.half_angle.tan();
    let center = info
        .carrier
        .apex
        .add(info.carrier.axis.scale(reference_station));
    let circle = make_arc(
        center,
        info.x_axis,
        info.y_axis,
        reference_radius,
        low,
        high,
    )?;
    let normal = plane.normal.normalized()?;
    let apex = info.carrier.apex;
    let scale = plane.origin.sub(apex).dot(normal);
    if scale.abs() <= tolerance {
        return Err("hybrid BREP: cone/plane section passes through the apex".into());
    }
    let mut controls = circle
        .control_points
        .iter()
        .map(|point| {
            let relative = Vec3::new(
                point.x - point.w * apex.x,
                point.y - point.w * apex.y,
                point.z - point.w * apex.z,
            );
            let weight = relative.dot(normal);
            let scaled = relative.scale(scale);
            Vec4 {
                x: apex.x * weight + scaled.x,
                y: apex.y * weight + scaled.y,
                z: apex.z * weight + scaled.z,
                w: weight,
            }
        })
        .collect::<Vec<_>>();
    if controls.iter().any(|point| point.w.abs() <= 1.0e-12) {
        return Err("hybrid BREP: cone/plane local conic crosses an asymptote".into());
    }
    let sign = controls[0].w.signum();
    if controls.iter().any(|point| point.w.signum() != sign) {
        return Err("hybrid BREP: cone/plane local conic has mixed weight signs".into());
    }
    if sign < 0.0 {
        for point in &mut controls {
            point.x = -point.x;
            point.y = -point.y;
            point.z = -point.z;
            point.w = -point.w;
        }
    }
    NurbsCurve::new(circle.degree, circle.knots, controls)
}

fn unwrap(values: &mut [f64]) {
    for index in 1..values.len() {
        while values[index] - values[index - 1] > 0.5 {
            values[index] -= 1.0;
        }
        while values[index] - values[index - 1] < -0.5 {
            values[index] += 1.0;
        }
    }
}

fn monotone_parameter_walk(values: &[f64], tolerance: f64) -> bool {
    let net = values.last().unwrap() - values[0];
    if net.abs() <= tolerance {
        return false;
    }
    let sign = net.signum();
    let variation = values
        .windows(2)
        .map(|pair| pair[1] - pair[0])
        .try_fold(0.0, |sum, step| {
            (step * sign >= -tolerance).then_some(sum + step.abs())
        });
    variation.is_some_and(|total| (total - net.abs()).abs() <= tolerance * values.len() as f64)
}

fn exact_arena_vertex(
    mesh_vertex: usize,
    exact_point: Vec3,
    tolerance: f64,
    arena: &mut TopologyArena,
    arena_vertices: &mut HashMap<usize, brep_kernel::VertexId>,
    ids: &mut Ids,
) -> Result<brep_kernel::VertexId, String> {
    if let Some(&vertex) = arena_vertices.get(&mesh_vertex) {
        let existing = arena
            .vertices
            .get(vertex)
            .ok_or_else(|| "hybrid BREP: reused arena vertex disappeared".to_owned())?;
        let disagreement = existing.point.sub(exact_point).length();
        if !disagreement.is_finite() || disagreement > tolerance {
            return Err(format!(
                "hybrid BREP: exact curves disagree by {disagreement:.6e} at shared mesh vertex {mesh_vertex}"
            ));
        }
        return Ok(vertex);
    }
    let vertex = arena.vertices.insert(ArenaVertex {
        wire_id: ids.next(),
        point: exact_point,
    });
    arena_vertices.insert(mesh_vertex, vertex);
    Ok(vertex)
}

#[derive(Clone, Copy)]
enum RevolveInfoRef<'a> {
    Cylinder(&'a CylinderInfo),
    Cone(&'a ConeInfo),
}

impl<'a> RevolveInfoRef<'a> {
    fn surface(self) -> &'a NurbsSurface {
        match self {
            Self::Cylinder(info) => &info.surface,
            Self::Cone(info) => &info.surface,
        }
    }

    fn coordinates(self, point: Vec3) -> [f64; 2] {
        match self {
            Self::Cylinder(info) => cylinder_coordinates(info, point),
            Self::Cone(info) => cone_coordinates(info, point),
        }
    }

    fn station(self, point: Vec3) -> f64 {
        match self {
            Self::Cylinder(info) => station(info.carrier, point),
            Self::Cone(info) => cone_station(info.carrier, point),
        }
    }

    fn radial_error(self, point: Vec3) -> f64 {
        match self {
            Self::Cylinder(info) => {
                (radial(info.carrier, point).length() - info.carrier.radius).abs()
            }
            Self::Cone(info) => {
                let s = cone_station(info.carrier, point);
                (cone_radial(info.carrier, point).length() - s * info.carrier.half_angle.tan())
                    .abs()
                    * info.carrier.half_angle.cos()
            }
        }
    }
}

fn revolve_info<'a>(
    region: u32,
    cylinders: &'a HashMap<u32, CylinderInfo>,
    cones: &'a HashMap<u32, ConeInfo>,
) -> Option<RevolveInfoRef<'a>> {
    cylinders
        .get(&region)
        .map(RevolveInfoRef::Cylinder)
        .or_else(|| cones.get(&region).map(RevolveInfoRef::Cone))
}

#[allow(clippy::too_many_arguments)]
fn build_chain_edge(
    chain_index: usize,
    chain: &Chain,
    vertices: &[Vec3],
    planes: &HashMap<u32, PlaneInfo>,
    cylinders: &HashMap<u32, CylinderInfo>,
    cones: &HashMap<u32, ConeInfo>,
    tolerance: f64,
    arena: &mut TopologyArena,
    arena_vertices: &mut HashMap<usize, brep_kernel::VertexId>,
    ids: &mut Ids,
) -> Result<BuiltEdge, String> {
    let points = chain_points(chain, vertices);
    let mut custom_revolve_pcurves = HashMap::new();
    let revolve_ids = [chain.adjacent.0, chain.adjacent.1]
        .into_iter()
        .filter(|id| cylinders.contains_key(id) || cones.contains_key(id))
        .collect::<Vec<_>>();
    if revolve_ids.len() > 2 {
        return Err(format!(
            "hybrid BREP: chain {chain_index} has too many revolved neighbors"
        ));
    }
    let (curve, t0, t1, curve_along_chain, primary_uv) = if let Some(&region) = revolve_ids.first()
    {
        let info = revolve_info(region, cylinders, cones).unwrap();
        let stations = points.iter().map(|&p| info.station(p)).collect::<Vec<_>>();
        let station_spread = stations.iter().copied().fold(f64::NEG_INFINITY, f64::max)
            - stations.iter().copied().fold(f64::INFINITY, f64::min);
        let radius_error = points
            .iter()
            .map(|&point| info.radial_error(point))
            .fold(0.0, f64::max);
        let mut coordinates = points
            .iter()
            .map(|&p| info.coordinates(p))
            .collect::<Vec<_>>();
        let mut us = coordinates.iter().map(|uv| uv[0]).collect::<Vec<_>>();
        unwrap(&mut us);
        for (uv, u) in coordinates.iter_mut().zip(&us) {
            uv[0] = *u;
        }
        if chain.closed {
            let first = points[0];
            let mut ring_points = points.clone();
            ring_points.push(first);
            let mut ring_us = ring_points
                .iter()
                .map(|&p| info.coordinates(p)[0])
                .collect::<Vec<_>>();
            unwrap(&mut ring_us);
            let sweep = ring_us[ring_us.len() - 1] - ring_us[0];
            if station_spread > tolerance
                || radius_error > tolerance
                || (sweep.abs() - 1.0).abs() > 0.02
                || !monotone_parameter_walk(&ring_us, 1.0e-5)
            {
                return Err(format!(
                    "hybrid BREP: closed chain {chain_index} is not a cylinder ring"
                ));
            }
            let v = coordinates.iter().map(|uv| uv[1]).sum::<f64>() / coordinates.len() as f64;
            let curve = info.surface().iso_curve_v(v)?;
            (
                curve,
                0.0,
                1.0,
                sweep > 0.0,
                Some(ChainUv {
                    region,
                    start: [0.0, v],
                    end: [if sweep > 0.0 { 1.0 } else { -1.0 }, v],
                }),
            )
        } else if station_spread <= tolerance && radius_error <= tolerance {
            if !monotone_parameter_walk(&us, 1.0e-5) {
                return Err(format!(
                    "hybrid BREP: cylinder arc chain {chain_index} backtracks"
                ));
            }
            let mut a = us[0];
            let mut b = *us.last().unwrap();
            let shift = (-2..=2)
                .map(f64::from)
                .find(|shift| {
                    a + shift >= -1.0e-8
                        && b + shift >= -1.0e-8
                        && a + shift <= 1.0 + 1.0e-8
                        && b + shift <= 1.0 + 1.0e-8
                })
                .ok_or_else(|| {
                    format!("hybrid BREP: cylinder arc chain {chain_index} crosses the chosen seam")
                })?;
            a += shift;
            b += shift;
            if a.abs() <= 1.0e-8 {
                a = 0.0;
            } else if (a - 1.0).abs() <= 1.0e-8 {
                a = 1.0;
            }
            if b.abs() <= 1.0e-8 {
                b = 0.0;
            } else if (b - 1.0).abs() <= 1.0e-8 {
                b = 1.0;
            }
            if (b - a).abs() <= 1.0e-9 || (b - a).abs() >= 1.0 - 1.0e-9 {
                return Err(format!(
                    "hybrid BREP: cylinder arc chain {chain_index} has invalid sweep"
                ));
            }
            let v = coordinates.iter().map(|uv| uv[1]).sum::<f64>() / coordinates.len() as f64;
            let low = a.min(b);
            let high = a.max(b);
            let curve = trim_curve(&info.surface().iso_curve_v(v)?, low, high)?;
            let domain = curve.domain()?;
            (
                curve,
                domain[0],
                domain[1],
                b > a,
                Some(ChainUv {
                    region,
                    start: [a, v],
                    end: [b, v],
                }),
            )
        } else if max_line_deviation(&points) <= tolerance
            || us.iter().copied().fold(f64::NEG_INFINITY, f64::max)
                - us.iter().copied().fold(f64::INFINITY, f64::min)
                <= 1.0e-4
        {
            let a = coordinates[0];
            let mut b = *coordinates.last().unwrap();
            while b[0] - a[0] > 0.5 {
                b[0] -= 1.0;
            }
            while b[0] - a[0] < -0.5 {
                b[0] += 1.0;
            }
            if (b[0] - a[0]).abs() > 1.0e-4 {
                return Err(format!(
                    "hybrid BREP: chain {chain_index} is not an axial ruling"
                ));
            }
            let u = (a[0] + b[0]) * 0.5;
            let low = a[1].min(b[1]);
            let high = a[1].max(b[1]);
            let curve = trim_curve(&info.surface().iso_curve_u(u.rem_euclid(1.0))?, low, high)?;
            let domain = curve.domain()?;
            (
                curve,
                domain[0],
                domain[1],
                b[1] > a[1],
                Some(ChainUv {
                    region,
                    start: [u, a[1]],
                    end: [u, b[1]],
                }),
            )
        } else if let (Some(cone), Some(plane)) = (
            cones.get(&region),
            [chain.adjacent.0, chain.adjacent.1]
                .into_iter()
                .find_map(|neighbor| planes.get(&neighbor)),
        ) {
            let normal = plane.normal.normalized()?;
            let x = normal.perpendicular()?;
            let y = normal.cross(x).normalized()?;
            let extent = points
                .iter()
                .map(|point| point.sub(plane.origin).length())
                .fold(1.0_f64, f64::max)
                * 4.0;
            let plane_surface = brep_kernel::make_plane(
                plane.origin.sub(x.scale(extent)).sub(y.scale(extent)),
                x,
                y,
                extent * 2.0,
                extent * 2.0,
            )?;
            let branches = intersect_analytic_pair(&plane_surface, &cone.surface, tolerance)
                .filter(|branches| !branches.is_empty())
                .unwrap_or_else(|| {
                    exact_local_cone_plane_conic(cone, *plane, &points, tolerance)
                        .into_iter()
                        .collect()
                });
            let mut candidates = branches
                .into_iter()
                .map(|curve| {
                    let projections = points
                        .iter()
                        .map(|&point| project_point_to_curve(&curve, point))
                        .collect::<Result<Vec<_>, String>>()?;
                    let score: f64 = projections
                        .iter()
                        .map(|projection| projection.distance)
                        .sum();
                    Ok((score, curve, projections))
                })
                .collect::<Result<Vec<_>, String>>()?;
            candidates.sort_by(|a, b| a.0.total_cmp(&b.0));
            let (_, curve, projections) = candidates.into_iter().next().ok_or_else(|| {
                format!("hybrid BREP: exact cone/plane conic is empty for chain {chain_index}")
            })?;
            let source_error = projections
                .iter()
                .map(|projection| projection.distance)
                .fold(0.0_f64, f64::max);
            if source_error > tolerance * 2.0 {
                return Err(format!(
                    "hybrid BREP: exact cone/plane conic chain {chain_index} misses source by {source_error:.6e}"
                ));
            }
            let p0 = projections[0];
            let p1 = *projections.last().unwrap();
            let low = p0.u.min(p1.u);
            let high = p0.u.max(p1.u);
            let direction = (p1.u - p0.u).signum();
            let monotone = direction != 0.0
                && projections.windows(2).all(|pair| {
                    (pair[1].u - pair[0].u) * direction >= -1.0e-8
                        && pair[1].u >= low - 1.0e-8
                        && pair[1].u <= high + 1.0e-8
                });
            if !monotone || high - low <= 1.0e-10 {
                return Err(format!(
                    "hybrid BREP: exact cone/plane conic chain {chain_index} crosses its branch seam"
                ));
            }
            let pcurve = build_pcurve_on_surface_range(
                &cone.surface,
                &curve,
                low,
                high,
                true,
                (tolerance * 0.05).max(1.0e-9),
            )?;
            let pcurve_domain = pcurve.domain()?;
            for sample in 0..=16 {
                let fraction = sample as f64 / 16.0;
                let t = low + (high - low) * fraction;
                let q = pcurve_domain[0] + (pcurve_domain[1] - pcurve_domain[0]) * fraction;
                let edge_point = curve.evaluate(t)?;
                let uv = pcurve.evaluate(q)?;
                let surface_point = cone.surface.evaluate(uv.x, uv.y)?;
                if edge_point.sub(surface_point).length() > tolerance * 0.1
                    || edge_point.sub(plane.origin).dot(normal).abs() > tolerance * 0.1
                {
                    return Err(format!(
                        "hybrid BREP: cone/plane conic chain {chain_index} pcurve failed agreement"
                    ));
                }
            }
            let edge_uv_start = pcurve.evaluate(pcurve_domain[0])?;
            let edge_uv_end = pcurve.evaluate(pcurve_domain[1])?;
            let curve_along_chain = p1.u > p0.u;
            let (chain_uv_start, chain_uv_end) = if curve_along_chain {
                (edge_uv_start, edge_uv_end)
            } else {
                (edge_uv_end, edge_uv_start)
            };
            custom_revolve_pcurves.insert(region, pcurve);
            (
                curve,
                low,
                high,
                curve_along_chain,
                Some(ChainUv {
                    region,
                    start: [chain_uv_start.x, chain_uv_start.y],
                    end: [chain_uv_end.x, chain_uv_end.y],
                }),
            )
        } else {
            return Err(format!(
                "hybrid BREP: revolved chain {chain_index} (region {region}, adjacent {:?}, closed={}, station spread {station_spread:.6e}, radial error {radius_error:.6e}, line deviation {:.6e}, uv {:?}->{:?}, points {:?}) is neither arc, ring, nor ruling",
                chain.adjacent,
                chain.closed,
                max_line_deviation(&points),
                coordinates.first(),
                coordinates.last(),
                points,
                // Source samples are included to make conservative demotions
                // diagnosable in focused regression tests.
            ));
        }
    } else {
        if chain.closed || max_line_deviation(&points) > tolerance {
            return Err(format!(
                "hybrid BREP: non-cylinder chain {chain_index} is not a line"
            ));
        }
        (
            make_line(points[0], *points.last().unwrap())?,
            0.0,
            1.0,
            true,
            None,
        )
    };

    let first = chain.vertices[0];
    let last = if chain.closed {
        first
    } else {
        *chain.vertices.last().unwrap()
    };
    let (start_mesh, end_mesh) = if curve_along_chain {
        (first, last)
    } else {
        (last, first)
    };
    let exact_start = curve.evaluate(t0)?;
    let exact_end = curve.evaluate(t1)?;
    let start = exact_arena_vertex(
        start_mesh,
        exact_start,
        tolerance,
        arena,
        arena_vertices,
        ids,
    )?;
    let end = exact_arena_vertex(end_mesh, exact_end, tolerance, arena, arena_vertices, ids)?;
    let edge = arena.edges.insert(ArenaEdge {
        wire_id: ids.next(),
        curve,
        t0,
        t1,
        start,
        end,
        degenerate: false,
        name: None,
    });
    let mut revolve_uvs = primary_uv.into_iter().collect::<Vec<_>>();
    if revolve_ids.len() == 2 {
        let secondary_region = revolve_ids[1];
        let secondary = revolve_info(secondary_region, cylinders, cones).unwrap();
        let mut coordinates = points
            .iter()
            .map(|&point| secondary.coordinates(point))
            .collect::<Vec<_>>();
        let mut us = coordinates.iter().map(|uv| uv[0]).collect::<Vec<_>>();
        unwrap(&mut us);
        for (uv, u) in coordinates.iter_mut().zip(&us) {
            uv[0] = *u;
        }
        let uv = if chain.closed {
            let mut ring_points = points.clone();
            ring_points.push(points[0]);
            let mut ring_us = ring_points
                .iter()
                .map(|&point| secondary.coordinates(point)[0])
                .collect::<Vec<_>>();
            unwrap(&mut ring_us);
            let sweep = ring_us.last().unwrap() - ring_us[0];
            let v = coordinates.iter().map(|uv| uv[1]).sum::<f64>() / coordinates.len() as f64;
            ChainUv {
                region: secondary_region,
                start: [0.0, v],
                end: [if sweep > 0.0 { 1.0 } else { -1.0 }, v],
            }
        } else {
            let start = coordinates[0];
            let end = *coordinates.last().unwrap();
            if (start[1] - end[1]).abs() > 1.0e-5 {
                return Err(format!(
                    "hybrid BREP: shared revolve chain {chain_index} is not a common ring arc"
                ));
            }
            ChainUv {
                region: secondary_region,
                start,
                end,
            }
        };
        revolve_uvs.push(uv);
    }
    Ok(BuiltEdge {
        edge,
        curve_along_chain,
        revolve_uvs,
        custom_revolve_pcurves,
    })
}

fn plane_frame(
    plane: PlaneInfo,
    boundary: &RegionBoundary,
    chains: &[Chain],
    vertices: &[Vec3],
) -> Result<(Vec3, Vec3, Vec3, f64, f64), String> {
    let normal = plane.normal.normalized()?;
    let x = normal.perpendicular()?;
    let y = normal.cross(x).normalized()?;
    let mut min_u = f64::INFINITY;
    let mut max_u = f64::NEG_INFINITY;
    let mut min_v = f64::INFINITY;
    let mut max_v = f64::NEG_INFINITY;
    for traversal in boundary.cycles.iter().flatten() {
        for &vertex in &chains[traversal.chain].vertices {
            let delta = vertices[vertex].sub(plane.origin);
            let u = delta.dot(x);
            let v = delta.dot(y);
            min_u = min_u.min(u);
            max_u = max_u.max(u);
            min_v = min_v.min(v);
            max_v = max_v.max(v);
        }
    }
    if !min_u.is_finite() || max_u - min_u <= 1.0e-12 || max_v - min_v <= 1.0e-12 {
        return Err("hybrid BREP: degenerate planar face bounds".into());
    }
    Ok((
        plane.origin.add(x.scale(min_u)).add(y.scale(min_v)),
        x,
        y,
        max_u - min_u,
        max_v - min_v,
    ))
}

fn affine_pcurve(curve: &NurbsCurve, origin: Vec3, x: Vec3, y: Vec3) -> Result<NurbsCurve, String> {
    let controls = curve
        .control_points
        .iter()
        .map(|control| {
            if !control.w.is_finite() || control.w.abs() <= f64::EPSILON {
                return Err("hybrid BREP: pcurve source has an invalid homogeneous weight".into());
            }
            let point = Vec3::new(
                control.x / control.w,
                control.y / control.w,
                control.z / control.w,
            );
            let delta = point.sub(origin);
            Ok(Vec4 {
                x: delta.dot(x) * control.w,
                y: delta.dot(y) * control.w,
                z: 0.0,
                w: control.w,
            })
        })
        .collect::<Result<Vec<_>, String>>()?;
    NurbsCurve::new(curve.degree, curve.knots.clone(), controls)
}

#[allow(clippy::too_many_arguments)]
fn build_plane_face(
    _region: u32,
    plane: PlaneInfo,
    boundary: &RegionBoundary,
    chains: &[Chain],
    built_edges: &[BuiltEdge],
    vertices: &[Vec3],
    arena: &mut TopologyArena,
    ids: &mut Ids,
) -> Result<FaceId, String> {
    let (origin, x, y, u_extent, v_extent) = plane_frame(plane, boundary, chains, vertices)?;
    let surface = brep_kernel::make_plane(origin, x, y, u_extent, v_extent)?;
    let mut loops = Vec::new();
    for cycle in &boundary.cycles {
        let mut coedges = Vec::new();
        for traversal in cycle {
            let built = &built_edges[traversal.chain];
            let edge = &arena.edges[built.edge];
            let mut pcurve = affine_pcurve(&edge.curve, origin, x, y)?;
            let forward = traversal.forward == built.curve_along_chain;
            if !forward {
                pcurve = pcurve.reversed()?;
            }
            coedges.push(arena.coedges.insert(ArenaCoedge {
                wire_id: ids.next(),
                edge: built.edge,
                forward,
                pcurve,
            }));
        }
        loops.push(arena.loops.insert(ArenaLoop {
            wire_id: ids.next(),
            coedges,
        }));
    }
    Ok(arena.faces.insert(ArenaFace {
        wire_id: ids.next(),
        surface,
        same_sense: true,
        loops,
        name: None,
    }))
}

fn parameter_line(start: [f64; 2], end: [f64; 2]) -> Result<NurbsCurve, String> {
    make_line(
        Vec3::new(start[0], start[1], 0.0),
        Vec3::new(end[0], end[1], 0.0),
    )
}

fn build_cylinder_face(
    region: u32,
    cylinder: &CylinderInfo,
    boundary: &RegionBoundary,
    built_edges: &[BuiltEdge],
    arena: &mut TopologyArena,
    ids: &mut Ids,
) -> Result<FaceId, String> {
    let mut loops = Vec::new();
    for cycle in &boundary.cycles {
        let mut coedges = Vec::new();
        let mut previous_end: Option<[f64; 2]> = None;
        for traversal in cycle {
            let built = &built_edges[traversal.chain];
            let uv = built
                .revolve_uvs
                .iter()
                .find(|uv| uv.region == region)
                .copied()
                .ok_or_else(|| format!("hybrid BREP: cylinder {region} has non-parametric edge"))?;
            let (mut start, mut end) = if traversal.forward {
                (uv.start, uv.end)
            } else {
                (uv.end, uv.start)
            };
            if let Some(previous) = previous_end {
                let shift = (previous[0] - start[0]).round();
                start[0] += shift;
                end[0] += shift;
            }
            previous_end = Some(end);
            let forward = traversal.forward == built.curve_along_chain;
            let mut pcurve = built
                .custom_revolve_pcurves
                .get(&region)
                .cloned()
                .unwrap_or(parameter_line(start, end)?);
            if !forward && built.custom_revolve_pcurves.contains_key(&region) {
                pcurve = pcurve.reversed()?;
            }
            coedges.push(arena.coedges.insert(ArenaCoedge {
                wire_id: ids.next(),
                edge: built.edge,
                forward,
                pcurve,
            }));
        }
        loops.push(arena.loops.insert(ArenaLoop {
            wire_id: ids.next(),
            coedges,
        }));
    }
    Ok(arena.faces.insert(ArenaFace {
        wire_id: ids.next(),
        surface: cylinder.surface.clone(),
        same_sense: cylinder.carrier.sense == 1,
        loops,
        name: None,
    }))
}

fn build_cone_face(
    region: u32,
    cone: &ConeInfo,
    boundary: &RegionBoundary,
    built_edges: &[BuiltEdge],
    arena: &mut TopologyArena,
    ids: &mut Ids,
) -> Result<FaceId, String> {
    let mut loops = Vec::new();
    for cycle in &boundary.cycles {
        let mut coedges = Vec::new();
        let mut previous_end: Option<[f64; 2]> = None;
        for traversal in cycle {
            let built = &built_edges[traversal.chain];
            let uv = built
                .revolve_uvs
                .iter()
                .find(|uv| uv.region == region)
                .copied()
                .ok_or_else(|| format!("hybrid BREP: cone {region} has non-parametric edge"))?;
            let (mut start, mut end) = if traversal.forward {
                (uv.start, uv.end)
            } else {
                (uv.end, uv.start)
            };
            if let Some(previous) = previous_end {
                let shift = (previous[0] - start[0]).round();
                start[0] += shift;
                end[0] += shift;
            }
            previous_end = Some(end);
            let forward = traversal.forward == built.curve_along_chain;
            let mut pcurve = built
                .custom_revolve_pcurves
                .get(&region)
                .cloned()
                .unwrap_or(parameter_line(start, end)?);
            if !forward && built.custom_revolve_pcurves.contains_key(&region) {
                pcurve = pcurve.reversed()?;
            }
            coedges.push(arena.coedges.insert(ArenaCoedge {
                wire_id: ids.next(),
                edge: built.edge,
                forward,
                pcurve,
            }));
        }
        loops.push(arena.loops.insert(ArenaLoop {
            wire_id: ids.next(),
            coedges,
        }));
    }
    Ok(arena.faces.insert(ArenaFace {
        wire_id: ids.next(),
        surface: cone.surface.clone(),
        same_sense: cone.carrier.sense == 1,
        loops,
        name: None,
    }))
}

fn sort_region_cycles(
    boundaries: &mut [RegionBoundary],
    planes: &HashMap<u32, PlaneInfo>,
    cylinders: &HashMap<u32, CylinderInfo>,
    cones: &HashMap<u32, ConeInfo>,
    chains: &[Chain],
    vertices: &[Vec3],
) {
    for (region, boundary) in boundaries.iter_mut().enumerate() {
        let mut scored = boundary
            .cycles
            .drain(..)
            .map(|cycle| {
                let mut points = Vec::new();
                for traversal in &cycle {
                    let chain = &chains[traversal.chain];
                    let iter: Box<dyn Iterator<Item = &usize>> = if traversal.forward {
                        Box::new(chain.vertices.iter())
                    } else {
                        Box::new(chain.vertices.iter().rev())
                    };
                    for vertex in iter {
                        points.push(vertices[*vertex]);
                    }
                }
                let area = if let Some(plane) = planes.get(&(region as u32)) {
                    let normal = plane.normal.normalized().unwrap_or(plane.normal);
                    let x = normal.perpendicular().unwrap_or(Vec3::new(1.0, 0.0, 0.0));
                    let y = normal.cross(x);
                    polygon_area(
                        &points
                            .iter()
                            .map(|p| {
                                let d = p.sub(plane.origin);
                                [d.dot(x), d.dot(y)]
                            })
                            .collect::<Vec<_>>(),
                    )
                } else if let Some(cylinder) = cylinders.get(&(region as u32)) {
                    let mut uv = points
                        .iter()
                        .map(|&p| cylinder_coordinates(cylinder, p))
                        .collect::<Vec<_>>();
                    let mut us = uv.iter().map(|p| p[0]).collect::<Vec<_>>();
                    unwrap(&mut us);
                    for (p, u) in uv.iter_mut().zip(us) {
                        p[0] = u;
                    }
                    polygon_area(&uv)
                } else if let Some(cone) = cones.get(&(region as u32)) {
                    let mut uv = points
                        .iter()
                        .map(|&p| cone_coordinates(cone, p))
                        .collect::<Vec<_>>();
                    let mut us = uv.iter().map(|p| p[0]).collect::<Vec<_>>();
                    unwrap(&mut us);
                    for (p, u) in uv.iter_mut().zip(us) {
                        p[0] = u;
                    }
                    polygon_area(&uv)
                } else {
                    0.0
                };
                let same = planes.contains_key(&(region as u32))
                    || cylinders
                        .get(&(region as u32))
                        .is_some_and(|c| c.carrier.sense == 1)
                    || cones
                        .get(&(region as u32))
                        .is_some_and(|c| c.carrier.sense == 1);
                let score = if same { area } else { -area };
                (score, cycle)
            })
            .collect::<Vec<_>>();
        scored.sort_by(|a, b| b.0.total_cmp(&a.0));
        boundary.cycles = scored.into_iter().map(|(_, cycle)| cycle).collect();
    }
}

fn polygon_area(points: &[[f64; 2]]) -> f64 {
    if points.len() < 3 {
        return 0.0;
    }
    (0..points.len())
        .map(|i| {
            let a = points[i];
            let b = points[(i + 1) % points.len()];
            a[0] * b[1] - b[0] * a[1]
        })
        .sum::<f64>()
        * 0.5
}

// BREP private tests: fcd92d1fb084babc