brep_kernel/construction/loft_topology/

basic.rs

use super::*;

#[derive(Clone)]
struct SectionFrame {
    normal: Vec3,
    centroid: Vec3,
    samples: Vec<Vec3>,
    planar: bool,
}

pub(super) fn closed_points(curves: &[NurbsCurve], tolerance: f64) -> Result<Vec<Vec3>, String> {
    if curves.len() < 2 {
        return Err("loftSolid: a section needs at least 2 curves".into());
    }
    let mut points = Vec::with_capacity(curves.len());
    for (index, curve) in curves.iter().enumerate() {
        let [start, end] = curve.domain()?;
        let next = &curves[(index + 1) % curves.len()];
        let next_start = next.domain()?[0];
        if curve
            .evaluate(end)?
            .sub(next.evaluate(next_start)?)
            .length()
            > tolerance
        {
            return Err(format!("loftSolid: section is open at curve {index}"));
        }
        points.push(curve.evaluate(start)?);
    }
    Ok(points)
}

fn section_frame(curves: &[NurbsCurve], tolerance: f64) -> Result<SectionFrame, String> {
    let mut samples = Vec::new();
    for curve in curves {
        let [start, end] = curve.domain()?;
        for index in 0..16 {
            samples.push(curve.evaluate(start + (end - start) * index as f64 / 16.0)?);
        }
    }
    let mut normal = Vec3::default();
    let mut centroid = Vec3::default();
    for index in 0..samples.len() {
        let point = samples[index];
        let next = samples[(index + 1) % samples.len()];
        normal.x += (point.y - next.y) * (point.z + next.z);
        normal.y += (point.z - next.z) * (point.x + next.x);
        normal.z += (point.x - next.x) * (point.y + next.y);
        centroid = centroid.add(point);
    }
    normal = normal.normalized()?;
    centroid = centroid.scale(1.0 / samples.len() as f64);
    let planar = samples
        .iter()
        .all(|point| point.sub(samples[0]).dot(normal).abs() <= tolerance * 100.0);
    Ok(SectionFrame {
        normal,
        centroid,
        samples,
        planar,
    })
}

fn reverse_section(curves: &[NurbsCurve]) -> Result<Vec<NurbsCurve>, String> {
    curves.iter().rev().map(NurbsCurve::reversed).collect()
}

fn cap_face(
    curves: &[NurbsCurve],
    edge_ids: &[u64],
    frame: &SectionFrame,
    axis: Vec3,
    outward: bool,
    next_id: &mut u64,
) -> Result<FaceRecord, String> {
    let normal = if frame.normal.dot(axis) >= 0.0 {
        frame.normal
    } else {
        frame.normal.scale(-1.0)
    };
    let x_axis = normal.perpendicular()?;
    let y_axis = normal.cross(x_axis).normalized()?;
    let (mut min_x, mut min_y) = (f64::INFINITY, f64::INFINITY);
    let (mut max_x, mut max_y) = (f64::NEG_INFINITY, f64::NEG_INFINITY);
    for point in &frame.samples {
        let delta = point.sub(frame.centroid);
        min_x = min_x.min(delta.dot(x_axis));
        max_x = max_x.max(delta.dot(x_axis));
        min_y = min_y.min(delta.dot(y_axis));
        max_y = max_y.max(delta.dot(y_axis));
    }
    let padding = (max_x - min_x).max(max_y - min_y) * 0.05 + 1e-6;
    let origin = frame
        .centroid
        .add(x_axis.scale(min_x - padding))
        .add(y_axis.scale(min_y - padding));
    let surface = make_plane(
        origin,
        x_axis,
        y_axis,
        max_x - min_x + 2.0 * padding,
        max_y - min_y + 2.0 * padding,
    )?;
    let mut coedges = Vec::with_capacity(curves.len());
    if outward {
        for (index, curve) in curves.iter().enumerate() {
            coedges.push(CoedgeRecord {
                id: *next_id,
                edge_id: edge_ids[index],
                forward: true,
                pcurve: curve_to_plane_parameters(curve, origin, x_axis, y_axis)?,
            });
            *next_id += 1;
        }
    } else {
        for index in (0..curves.len()).rev() {
            coedges.push(CoedgeRecord {
                id: *next_id,
                edge_id: edge_ids[index],
                forward: false,
                pcurve: curve_to_plane_parameters(&curves[index], origin, x_axis, y_axis)?
                    .reversed()?,
            });
            *next_id += 1;
        }
    }
    let loop_id = *next_id;
    *next_id += 1;
    let face_id = *next_id;
    *next_id += 1;
    Ok(FaceRecord {
        id: face_id,
        surface,
        same_sense: outward,
        loops: vec![LoopRecord {
            id: loop_id,
            coedges,
        }],
        name: None,
    })
}

pub fn loft_profile_brep(input_sections: &[Vec<NurbsCurve>]) -> Result<BrepSolid, String> {
    loft_profile_brep_core(input_sections, None)
}

/// §5.8 loft with END TANGENCY: the skin leaves the first section along
/// `start_direction` and arrives at the last along `end_direction` (unit
/// directions; each interpolation column scales them by its own chord length,
/// the standard magnitude that keeps the v-parametrization well conditioned).
/// Exact by construction — the column interpolant reproduces the prescribed
/// end derivatives.
pub fn loft_profile_brep_tangent(
    input_sections: &[Vec<NurbsCurve>],
    start_direction: Vec3,
    end_direction: Vec3,
) -> Result<BrepSolid, String> {
    let start = start_direction
        .normalized()
        .map_err(|_| "loftSolid: start tangent must be a nonzero direction".to_string())?;
    let end = end_direction
        .normalized()
        .map_err(|_| "loftSolid: end tangent must be a nonzero direction".to_string())?;
    loft_profile_brep_core(input_sections, Some((start, end)))
}

fn loft_profile_brep_core(
    input_sections: &[Vec<NurbsCurve>],
    end_tangents: Option<(Vec3, Vec3)>,
) -> Result<BrepSolid, String> {
    let tolerance = 1e-6;
    let section_count = input_sections.len();
    if section_count < 2 {
        return Err("loftSolid: need at least 2 sections".into());
    }
    let mut sections = input_sections.to_vec();
    let curve_count = sections[0].len();
    if sections.iter().any(|section| section.len() != curve_count) {
        return Err("loftSolid: sections must have the same curve count".into());
    }
    for section in &sections {
        closed_points(section, tolerance)?;
    }
    for curve_index in 0..curve_count {
        let reference = &sections[0][curve_index];
        for (section_index, section) in sections.iter().enumerate().skip(1) {
            let curve = &section[curve_index];
            if curve.degree != reference.degree
                || curve.control_points.len() != reference.control_points.len()
            {
                return Err(format!(
                    "loftSolid: section {section_index} curve {curve_index} incompatible with section 0"
                ));
            }
            if curve.knots.len() != reference.knots.len()
                || curve
                    .knots
                    .iter()
                    .zip(&reference.knots)
                    .any(|(a, b)| (a - b).abs() > 1e-9)
            {
                return Err(format!(
                    "loftSolid: section {section_index} curve {curve_index} has different knots"
                ));
            }
            if curve
                .control_points
                .iter()
                .zip(&reference.control_points)
                .any(|(a, b)| (a.w - b.w).abs() > 1e-9)
            {
                return Err(format!(
                    "loftSolid: section {section_index} curve {curve_index} has different weights"
                ));
            }
        }
    }

    let first_frame = section_frame(&sections[0], tolerance)?;
    let last_frame = section_frame(&sections[section_count - 1], tolerance)?;
    if !first_frame.planar || !last_frame.planar {
        return Err("loftSolid: end sections must be planar".into());
    }
    let axis = last_frame
        .centroid
        .sub(first_frame.centroid)
        .normalized()
        .map_err(|_| "loftSolid: end sections coincide".to_string())?;
    let first_normal = if first_frame.normal.dot(axis) >= 0.0 {
        first_frame.normal
    } else {
        first_frame.normal.scale(-1.0)
    };
    if first_normal.dot(axis).abs() < 0.1 {
        return Err("loftSolid: loft direction nearly parallel to end section plane".into());
    }
    let x_axis = first_normal.perpendicular()?;
    let y_axis = first_normal.cross(x_axis).normalized()?;
    // Winding is normalized PER SECTION against this one frame, not decided once
    // from section 0 and applied to every section.
    //
    // The skin rails control point `k` of curve `i` of one section to the same
    // `(i, k)` of the next, so two sections traversed OPPOSITE ways rail corner
    // to opposite corner and the skin twists into a bow-tie. Deciding one flip
    // from section 0 and applying it to all of them cannot see that: it reverses
    // them together, which preserves the disagreement exactly.
    //
    // Measured before this changed (`examples/feature_refusal_census_probe.rs`,
    // `GATE loft_profile_brep.opposed_winding`): two 2x2 squares 5 apart, wound
    // opposite ways, lofted to a solid of **volume 0 that passes `validate()`**
    // — a silently-wrong result, not a refusal. The shipping LOFT feature
    // reached it with two ordinary sketches on opposed planes.
    //
    // The comparison is CHAIN-RELATIVE — each section against its PREDECESSOR,
    // not against section 0.
    //
    // Measuring every section against section 0's frame looks equivalent and is
    // not: a frame-guided loft rotates each station's plane with the guide, so
    // on a bend past 90 degrees a station's projected area goes negative purely
    // because its plane has turned (cos > 90 deg < 0), and reversing it un-does
    // the guide's own rotation. Measured on a 150-degree arc: volume 17.5685
    // against a Pappus-exact 104.7198, a 83 % error, where the predecessor rule
    // is exact.
    //
    // Consecutive stations of any sane loft turn far less than 90 degrees, so
    // the predecessor comparison leaves a bend alone, while a genuinely opposed
    // pair still flips. It also removes the knife edge at exactly 90 degrees,
    // where a fixed-frame projection is float noise.
    //
    // `section_frame(..).normal` is the Newell normal, whose ORIENTATION already
    // encodes the traversal direction; only its sign is read. Section 0 keeps
    // the original rule verbatim, so its behaviour — and the side-face
    // permutation below, which keys off section 0 alone — is byte-identical.
    let mut section_reversed = Vec::with_capacity(section_count);
    section_reversed.push(profile_area(&sections[0], first_frame.centroid, x_axis, y_axis)? < 0.0);
    let mut previous_normal = if section_reversed[0] {
        first_frame.normal.scale(-1.0)
    } else {
        first_frame.normal
    };
    for section in sections.iter().skip(1) {
        let normal = section_frame(section, tolerance)?.normal;
        let flip = normal.dot(previous_normal) < 0.0;
        section_reversed.push(flip);
        previous_normal = if flip { normal.scale(-1.0) } else { normal };
    }
    // Side faces are named by the FIRST section's input-curve order, so the
    // face permutation below still keys off section 0 alone.
    let reversed_winding = section_reversed[0];
    for (index, flip) in section_reversed.into_iter().enumerate() {
        if flip {
            sections[index] = reverse_section(&sections[index])?;
        }
    }

    let mut parameters = vec![0.0; section_count];
    let mut accumulated = vec![0.0; section_count];
    let mut columns = 0usize;
    for curve_index in 0..curve_count {
        for control_index in 0..sections[0][curve_index].control_points.len() {
            let mut total = 0.0;
            let mut chords = vec![0.0; section_count];
            for section_index in 1..section_count {
                let previous = sections[section_index - 1][curve_index].control_points
                    [control_index]
                    .point()?;
                let current =
                    sections[section_index][curve_index].control_points[control_index].point()?;
                total += current.sub(previous).length();
                chords[section_index] = total;
            }
            if total <= tolerance {
                continue;
            }
            for section_index in 0..section_count {
                accumulated[section_index] += chords[section_index] / total;
            }
            columns += 1;
        }
    }
    if columns == 0 {
        return Err("loftSolid: sections coincide".into());
    }
    for index in 0..section_count {
        parameters[index] = accumulated[index] / columns as f64;
    }
    parameters[0] = 0.0;
    parameters[section_count - 1] = 1.0;
    if parameters.windows(2).any(|pair| pair[1] <= pair[0] + 1e-9) {
        return Err("loftSolid: sections are not strictly ordered".into());
    }
    // Tangent lofts always interpolate cubically — the end-derivative rows
    // need the two extra control points even for a 2-section Hermite loft.
    let degree_v = if end_tangents.is_some() {
        3
    } else {
        3usize.min(section_count - 1)
    };
    let mut skins = Vec::with_capacity(curve_count);
    for curve_index in 0..curve_count {
        let reference = &sections[0][curve_index];
        let mut grid = Vec::with_capacity(reference.control_points.len());
        let mut knots_v = Vec::new();
        for control_index in 0..reference.control_points.len() {
            let points = sections
                .iter()
                .map(|section| section[curve_index].control_points[control_index].point())
                .collect::<Result<Vec<_>, _>>()?;
            let interpolated = match end_tangents {
                None => interpolate_curve(&points, degree_v, &parameters)?,
                Some((start_direction, end_direction)) => {
                    let chord: f64 = points
                        .windows(2)
                        .map(|pair| pair[1].sub(pair[0]).length())
                        .sum();
                    // A column whose stations coincide still needs a usable
                    // tangent magnitude; fall back to the section spacing.
                    let magnitude = if chord > tolerance { chord } else { 1.0 };
                    crate::interpolate_curve_with_end_tangents(
                        &points,
                        &parameters,
                        start_direction.scale(magnitude),
                        end_direction.scale(magnitude),
                    )?
                }
            };
            knots_v = interpolated.knots.clone();
            let weight = reference.control_points[control_index].w;
            grid.push(
                interpolated
                    .control_points
                    .iter()
                    .map(|point| Vec4::from_point(point.point().unwrap(), weight))
                    .collect(),
            );
        }
        skins.push(NurbsSurface::new(
            reference.degree,
            degree_v,
            reference.knots.clone(),
            knots_v,
            grid,
        )?);
    }

    let bottom = &sections[0];
    let top = &sections[section_count - 1];
    let bottom_points = closed_points(bottom, tolerance)?;
    let top_points = closed_points(top, tolerance)?;
    let mut vertices = Vec::with_capacity(2 * curve_count);
    for (index, point) in bottom_points.iter().chain(&top_points).enumerate() {
        vertices.push(VertexRecord {
            id: index as u64 + 1,
            point: *point,
        });
    }
    let mut edges = Vec::with_capacity(3 * curve_count);
    let mut bottom_edge_ids = Vec::new();
    let mut top_edge_ids = Vec::new();
    let mut vertical_edge_ids = Vec::new();
    for index in 0..curve_count {
        let [start, end] = bottom[index].domain()?;
        let bottom_id = 10 + index as u64;
        let top_id = 10 + curve_count as u64 + index as u64;
        let vertical_id = 10 + 2 * curve_count as u64 + index as u64;
        bottom_edge_ids.push(bottom_id);
        top_edge_ids.push(top_id);
        vertical_edge_ids.push(vertical_id);
        edges.push(EdgeRecord {
            id: bottom_id,
            curve: bottom[index].clone(),
            t0: start,
            t1: end,
            start_vertex_id: index as u64 + 1,
            end_vertex_id: ((index + 1) % curve_count) as u64 + 1,
            degenerate: false,
            name: None,
        });
        edges.push(EdgeRecord {
            id: top_id,
            curve: top[index].clone(),
            t0: start,
            t1: end,
            start_vertex_id: (curve_count + index) as u64 + 1,
            end_vertex_id: (curve_count + (index + 1) % curve_count) as u64 + 1,
            degenerate: false,
            name: None,
        });
        let vertical = skins[index].iso_curve_u(start)?;
        let [v_start, v_end] = vertical.domain()?;
        edges.push(EdgeRecord {
            id: vertical_id,
            curve: vertical,
            t0: v_start,
            t1: v_end,
            start_vertex_id: index as u64 + 1,
            end_vertex_id: (curve_count + index) as u64 + 1,
            degenerate: false,
            name: None,
        });
    }

    let mut next_id = 1000u64;
    let mut faces = Vec::with_capacity(curve_count + 2);
    for index in 0..curve_count {
        let [start, end] = bottom[index].domain()?;
        let coedges = vec![
            CoedgeRecord {
                id: next_id,
                edge_id: bottom_edge_ids[index],
                forward: true,
                pcurve: parameter_line(start, 0.0, end, 0.0)?,
            },
            CoedgeRecord {
                id: next_id + 1,
                edge_id: vertical_edge_ids[(index + 1) % curve_count],
                forward: true,
                pcurve: parameter_line(end, 0.0, end, 1.0)?,
            },
            CoedgeRecord {
                id: next_id + 2,
                edge_id: top_edge_ids[index],
                forward: false,
                pcurve: parameter_line(end, 1.0, start, 1.0)?,
            },
            CoedgeRecord {
                id: next_id + 3,
                edge_id: vertical_edge_ids[index],
                forward: false,
                pcurve: parameter_line(start, 1.0, start, 0.0)?,
            },
        ];
        next_id += 4;
        let loop_id = next_id;
        next_id += 1;
        let face_id = next_id;
        next_id += 1;
        faces.push(FaceRecord {
            id: face_id,
            surface: skins[index].clone(),
            same_sense: true,
            loops: vec![LoopRecord {
                id: loop_id,
                coedges,
            }],
            name: None,
        });
    }
    if reversed_winding {
        // The winding normalization reversed the section curve order above.
        // Callers stamp side-face names by INPUT-curve order (the extrude
        // and revolve wall-naming permutation, loft edition) — emit side
        // faces in input order.
        faces.reverse();
    }
    faces.push(cap_face(
        bottom,
        &bottom_edge_ids,
        &first_frame,
        axis,
        false,
        &mut next_id,
    )?);
    faces.push(cap_face(
        top,
        &top_edge_ids,
        &last_frame,
        axis,
        true,
        &mut next_id,
    )?);
    let solid = BrepSolid {
        id: next_id + 1,
        vertices,
        edges,
        shells: vec![ShellRecord { id: next_id, faces }],
        genus: 0,
    };
    let issues = solid.validate();
    if issues.is_empty() {
        Ok(solid)
    } else {
        Err(format!(
            "Rust loft builder produced invalid topology: {issues:?}"
        ))
    }
}