brep_kernel/construction/sweep_topology/

sweep.rs

use super::*;

/// Sweep a CLOSED PLANAR profile loop along a path curve (§3.3/§5.7).
///
/// Reuses the loft builder rather than a bespoke swept surface: the path is
/// sampled at `STATIONS` uniform stations, a rotation-minimizing frame is
/// propagated along it (double-reflection RMF, Wang et al. 2008 — the frame
/// does NOT spin at inflections the way raw Frenet does), a rigidly
/// transformed copy of the profile is placed at each station (identical
/// degree/knots/weights, only the control points moved — which guarantees
/// loft's per-section compatibility), and the stations are lofted through
/// to produce the tube plus planar end caps.
///
/// Guards return a clear `Err` on an open/non-planar profile, a degenerate
/// path tangent, or a loft failure. A self-intersecting result (path
/// curvature radius smaller than the profile extent) is OUT OF SCOPE — the
/// caller is responsible for keeping the tube from folding onto itself.
pub fn sweep_profile_along_path(
    profile: &[NurbsCurve],
    path: &NurbsCurve,
    name: Option<&str>,
) -> Result<BrepSolid, String> {
    // 32 stations is the original fixed sampling — golden parity pins the
    // emitted geometry to it.  The helix variant raises the count with the
    // turn count instead, hence the shared `_stations` core.
    sweep_profile_along_path_stations(profile, path, name, 32, 0.0, None)
}

/// Twisted path sweep (§5.7): identical to [`sweep_profile_along_path`], but
/// the profile additionally ROTATES about the path tangent, linearly in ARC
/// LENGTH, from 0 at the sweep start to `twist_angle` radians (right-handed
/// about the tangent) at the end.  The arc-length fraction comes from the
/// sampled station polyline, not the raw path parameter, so a non-uniformly
/// parameterized path still twists uniformly in space.
///
/// STATION LAW: 16 stations per quarter turn of twist, floored at the path
/// sweep's 32 and capped at 1024 (the loft's dense interpolation solve is
/// O(stations³) per control column — the same cap the helix uses).  At
/// 16/quarter-turn the inter-station twist step is Δφ ≈ 5.6°, so the cubic
/// v-interpolation error on a profile point circling at radius r is
/// ≈ r·Δφ⁴/384 ≈ 2.4·10⁻¹⁰·r — far below any geometric tolerance.  The cap
/// holds that density up to 16 full turns (|twist| = 32π = 1024/16 quarter
/// turns); a larger twist would silently alias under the cap, so it is
/// REJECTED with an honest error instead.  Everything else the path sweep
/// documents (profile validity, self-intersection being the caller's
/// responsibility) applies unchanged; twisting about the profile's own
/// centroid adds no new radial extent, so no extra collision guard exists
/// to compute here.
pub fn sweep_profile_twisted(
    profile: &[NurbsCurve],
    path: &NurbsCurve,
    twist_angle: f64,
    name: Option<&str>,
) -> Result<BrepSolid, String> {
    use std::f64::consts::{FRAC_PI_2, TAU};

    if !twist_angle.is_finite() {
        return Err("sweep_profile_twisted: twist angle must be finite".into());
    }
    // 16 turns is where the 1024-station cap meets 16 stations/quarter-turn;
    // beyond it the cap would degrade the twist sampling density silently.
    const MAX_TURNS: f64 = 16.0;
    if twist_angle.abs() > MAX_TURNS * TAU {
        return Err(format!(
            "sweep_profile_twisted: twist of {:.3} turns exceeds the {MAX_TURNS}-turn \
             limit the 1024-station cap can resolve at 16 stations per quarter turn; \
             split the sweep or reduce the twist",
            twist_angle.abs() / TAU
        ));
    }
    let quarter_turns = (twist_angle.abs() / FRAC_PI_2).ceil() as usize;
    let stations = (quarter_turns * 16).clamp(32, 1024);
    sweep_profile_along_path_stations(profile, path, name, stations, twist_angle, None)
        .map_err(|error| format!("sweep_profile_twisted: {error}"))
}

/// The placement anchor a path sweep transplants its profile with: the plane
/// frame `(origin = boundary-sample centroid, normal, pu, pv)` the station
/// loop maps profile points through (`local = p − origin` → `station + ri·(local
/// ·pu) + si·(local·pv)`). Extracted as data so a HOLE loop can sweep with its
/// OUTER loop's anchor — sweeping each loop with its OWN centroid would
/// re-center every loop onto the path and lose the hole's in-plane offset.
#[derive(Debug, Clone, Copy)]
pub struct ProfileAnchor {
    pub origin: Vec3,
    pub normal: Vec3,
    pub pu: Vec3,
    pub pv: Vec3,
}

/// [`sweep_profile_along_path`] with an explicit placement anchor (see
/// [`ProfileAnchor`]) — the hole-loop cutter path: the swept loop is validated
/// as usual but PLACED in its outer loop's frame.
pub fn sweep_profile_along_path_anchored(
    profile: &[NurbsCurve],
    path: &NurbsCurve,
    name: Option<&str>,
    anchor: ProfileAnchor,
) -> Result<BrepSolid, String> {
    sweep_profile_along_path_stations(profile, path, name, 32, 0.0, Some(anchor))
}

/// [`sweep_profile_twisted`] with an explicit placement anchor: the hole loop
/// twists about the SAME path axis as its outer loop (shared anchor), so the
/// cutter stays registered with the outer wall through the whole twist.
pub fn sweep_profile_twisted_anchored(
    profile: &[NurbsCurve],
    path: &NurbsCurve,
    twist_angle: f64,
    name: Option<&str>,
    anchor: ProfileAnchor,
) -> Result<BrepSolid, String> {
    use std::f64::consts::{FRAC_PI_2, TAU};

    if !twist_angle.is_finite() {
        return Err("sweep_profile_twisted: twist angle must be finite".into());
    }
    const MAX_TURNS: f64 = 16.0;
    if twist_angle.abs() > MAX_TURNS * TAU {
        return Err(format!(
            "sweep_profile_twisted: twist of {:.3} turns exceeds the {MAX_TURNS}-turn \
             limit the 1024-station cap can resolve at 16 stations per quarter turn; \
             split the sweep or reduce the twist",
            twist_angle.abs() / TAU
        ));
    }
    let quarter_turns = (twist_angle.abs() / FRAC_PI_2).ceil() as usize;
    let stations = (quarter_turns * 16).clamp(32, 1024);
    sweep_profile_along_path_stations(profile, path, name, stations, twist_angle, Some(anchor))
        .map_err(|error| format!("sweep_profile_twisted: {error}"))
}

/// Validate a closed planar profile loop and derive its placement anchor —
/// the path sweep's §1 block, extracted bit-identically: 16 samples per curve,
/// closure at `tolerance`, Newell normal, boundary-sample-mean origin,
/// planarity at `tolerance * 100`, `pu = np.perpendicular()`, `pv = np × pu`.
pub fn profile_anchor(profile: &[NurbsCurve]) -> Result<ProfileAnchor, String> {
    let tolerance = 1e-6;
    if profile.len() < 2 {
        return Err("sweepSolid: profile needs at least 2 curves forming a closed loop".into());
    }
    let mut samples = Vec::new();
    for (index, curve) in profile.iter().enumerate() {
        let [start, end] = curve.domain()?;
        let next = &profile[(index + 1) % profile.len()];
        let next_start = next.domain()?[0];
        if curve
            .evaluate(end)?
            .sub(next.evaluate(next_start)?)
            .length()
            > tolerance
        {
            return Err(format!(
                "sweepSolid: profile is not closed at curve {index}"
            ));
        }
        for sample in 0..16 {
            samples.push(curve.evaluate(start + (end - start) * sample 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);
    }
    let np = normal
        .normalized()
        .map_err(|_| "sweepSolid: profile is degenerate (zero enclosed area)".to_string())?;
    let origin = centroid.scale(1.0 / samples.len() as f64);
    if samples
        .iter()
        .any(|point| point.sub(origin).dot(np).abs() > tolerance * 100.0)
    {
        return Err("sweepSolid: profile is not planar".into());
    }
    let pu = np.perpendicular()?;
    let pv = np.cross(pu).normalized()?;
    Ok(ProfileAnchor {
        origin,
        normal: np,
        pu,
        pv,
    })
}

/// Core of the path sweep with an explicit station count.  Every consumer of
/// the station count (path sampling, RMF propagation, section placement)
/// derives from the one `stations` argument so the density scales as a unit.
/// `twist_angle` (radians; 0 for the untwisted variants) rotates the placed
/// profile about the path tangent linearly in sampled arc length — the
/// `twist_angle == 0.0` fast path leaves the RMF axes bit-identical, so the
/// untwisted callers keep golden parity.
fn sweep_profile_along_path_stations(
    profile: &[NurbsCurve],
    path: &NurbsCurve,
    name: Option<&str>,
    stations: usize,
    twist_angle: f64,
    anchor: Option<ProfileAnchor>,
) -> Result<BrepSolid, String> {
    // Loft carries no face names; the app stamps them onto the emitted face
    // order.  Accept `name` for ABI symmetry with the other builders.
    let _ = name;
    let tolerance = 1e-6;
    if stations < 2 {
        return Err("sweepSolid: need at least 2 stations".into());
    }

    // --- 1. Validate the profile: closed + planar; derive (origin, np, pu, pv).
    //        A caller-supplied anchor OVERRIDES the placement frame (the swept
    //        loop is still validated against its own plane), so a hole loop
    //        rides the path in its outer loop's frame instead of re-centering.
    let computed = profile_anchor(profile)?;
    let ProfileAnchor {
        origin, pu, pv, ..
    } = anchor.unwrap_or(computed);

    // --- 2. Sample the path; require a non-degenerate tangent at every station.
    let [t0, t1] = path.domain()?;
    if (t1 - t0).abs() <= tolerance {
        return Err("sweepSolid: path domain is degenerate".into());
    }
    let mut points = Vec::with_capacity(stations);
    let mut tangents = Vec::with_capacity(stations);
    for index in 0..stations {
        let t = t0 + (t1 - t0) * index as f64 / (stations - 1) as f64;
        let derivatives = path.derivatives(t, 1)?;
        let tangent = derivatives[1]
            .normalized()
            .map_err(|_| format!("sweepSolid: path tangent is degenerate at station {index}"))?;
        points.push(derivatives[0]);
        tangents.push(tangent);
    }

    // --- 2b. Twist distribution: cumulative ARC-LENGTH fractions over the
    //         sampled station polyline (chord sums), so the twist advances
    //         uniformly in space even on a non-uniformly parameterized path.
    //         Only computed when a twist is actually requested — the
    //         `twist_angle == 0.0` path must stay bit-identical to the
    //         pre-twist builder.
    let twist_fractions: Option<Vec<f64>> = if twist_angle != 0.0 {
        let mut cumulative = vec![0.0; stations];
        let mut total = 0.0;
        for index in 1..stations {
            total += points[index].sub(points[index - 1]).length();
            cumulative[index] = total;
        }
        if total <= tolerance {
            return Err("sweepSolid: path has zero length; cannot distribute the twist".into());
        }
        for length in &mut cumulative {
            *length /= total;
        }
        Some(cumulative)
    } else {
        None
    };

    // --- 3. Rotation-minimizing frames via the double-reflection method.
    let mut r_axes = Vec::with_capacity(stations);
    let mut s_axes = Vec::with_capacity(stations);
    let r0 = tangents[0].perpendicular()?; // any unit vector ⟂ T0
    s_axes.push(tangents[0].cross(r0).normalized()?);
    r_axes.push(r0);
    for index in 0..stations - 1 {
        let t_next = tangents[index + 1];
        let v1 = points[index + 1].sub(points[index]);
        let c1 = v1.dot(v1);
        let r_candidate = if c1 <= 1e-18 {
            // Coincident stations: carry the reference axis forward unchanged.
            r_axes[index]
        } else {
            // First reflection across the plane bisecting the step vector.
            let reflected_r = r_axes[index].sub(v1.scale(2.0 / c1 * v1.dot(r_axes[index])));
            let reflected_t = tangents[index].sub(v1.scale(2.0 / c1 * v1.dot(tangents[index])));
            // Second reflection across the plane bisecting the tangents.
            let v2 = t_next.sub(reflected_t);
            let c2 = v2.dot(v2);
            if c2 <= 1e-18 {
                reflected_r
            } else {
                reflected_r.sub(v2.scale(2.0 / c2 * v2.dot(reflected_r)))
            }
        };
        // Re-orthogonalize against the new tangent to shed floating drift.
        let r_next = r_candidate
            .sub(t_next.scale(r_candidate.dot(t_next)))
            .normalized()
            .map_err(|_| format!("sweepSolid: frame degenerated at station {index}"))?;
        s_axes.push(t_next.cross(r_next).normalized()?);
        r_axes.push(r_next);
    }

    // --- 4. Place a rigidly transformed copy of the profile at each station.
    let mut sections: Vec<Vec<NurbsCurve>> = Vec::with_capacity(stations);
    for station in 0..stations {
        let station_origin = points[station];
        // Rotate the RMF axes about the tangent by the station's twist angle
        // (Rodrigues on vectors ⟂ the tangent: r' = r·cosφ + s·sinφ,
        // s' = s·cosφ − r·sinφ, since s = t × r and t × s = −r).
        let (ri, si) = match &twist_fractions {
            Some(fractions) => {
                let phi = twist_angle * fractions[station];
                let (sin_phi, cos_phi) = phi.sin_cos();
                let r = r_axes[station];
                let s = s_axes[station];
                (
                    r.scale(cos_phi).add(s.scale(sin_phi)),
                    s.scale(cos_phi).sub(r.scale(sin_phi)),
                )
            }
            None => (r_axes[station], s_axes[station]),
        };
        let mut section = Vec::with_capacity(profile.len());
        for curve in profile {
            let control_points = curve
                .control_points
                .iter()
                .map(|point| {
                    let weight = point.w;
                    let euclidean = Vec3::new(point.x / weight, point.y / weight, point.z / weight);
                    let local = euclidean.sub(origin);
                    let world = station_origin
                        .add(ri.scale(local.dot(pu)))
                        .add(si.scale(local.dot(pv)));
                    Vec4 {
                        x: world.x * weight,
                        y: world.y * weight,
                        z: world.z * weight,
                        w: weight,
                    }
                })
                .collect();
            section.push(NurbsCurve::new(
                curve.degree,
                curve.knots.clone(),
                control_points,
            )?);
        }
        sections.push(section);
    }

    // --- 5. Loft through the swept sections (side walls + planar end caps).
    loft_profile_brep(&sections)
        .map_err(|error| format!("sweepSolid: loft through swept sections failed: {error}"))
}

/// Helical sweep (§5.7): sweep a CLOSED PLANAR profile loop along a helix of
/// `helix_radius` about the axis through `axis_origin` with direction
/// `axis_direction`, rising `pitch` per revolution for `turns` revolutions.
///
/// The helix is transcendental, not exactly NURBS-representable, so the path
/// is FIT: points (R·cosθ, R·sinθ, pitch·θ/2π) in the axis frame are sampled
/// uniformly in θ and globally interpolated with a cubic (`interpolate_curve`
/// — the same machinery every other fitted path here uses).  DENSITY: 64
/// samples per turn, capped at 1025 total nodes.  Cubic interpolation error
/// on a circle of radius R with node spacing Δθ is ≈ R·Δθ⁴/384: ~2·10⁻⁷·R at
/// 64/turn, and still ~6·10⁻⁵·R at the cap's worst case (16/turn at the
/// 64-turn limit) — orders below any profile a caller could sweep without
/// self-intersecting.  Uniform-in-θ parameters are exact chord-length for a
/// helix (constant speed), which is what the averaged-knot interpolation
/// assumes.
///
/// The fitted path then drives the EXISTING path-sweep core with 32 stations
/// per turn (min 32, capped at 1024 — the loft's dense interpolation solve is
/// O(stations³) per control column, so the cap trades per-turn density, never
/// correctness, at high turn counts).
///
/// GUARDS beyond the path sweep's own: the path sweep documents
/// self-intersection as caller responsibility, so the helix variant — which
/// knows its curvature analytically — rejects the two garbage modes itself:
///   • fold-over: the helix curvature radius (R² + c²)/R (c = pitch/2π) must
///     exceed the profile's max extent about its centroid, or the tube folds
///     through itself on the inner side (the torus tube-radius > major-radius
///     failure, pitch-relaxed);
///   • coil collision (turns ≥ 1): the normal gap between consecutive coils,
///     pitch·2πR/√((2πR)² + pitch²), must exceed the profile diameter.
/// Both use the profile's max sample distance from its centroid — conservative
/// for asymmetric profiles (extent in a harmless direction still counts), but
/// a false reject beats silent garbage.
///
/// KNOWN LIMITATION: whole-turn shallow helixes (pitch ≲ 0.63·R) are rejected
/// by the loft's cap-plane guard (`|n·axis| ≥ 0.1` — the end-to-end axis is
/// purely axial while the cap normal is nearly tangential); the error
/// propagates honestly rather than being worked around.
pub fn sweep_profile_helix(
    profile: &[NurbsCurve],
    axis_origin: Vec3,
    axis_direction: Vec3,
    helix_radius: f64,
    pitch: f64,
    turns: f64,
    name: Option<&str>,
) -> Result<BrepSolid, String> {
    use std::f64::consts::TAU;

    // --- 1. Validate the helix parameters with honest errors.
    let w = axis_direction
        .normalized()
        .map_err(|_| "sweep_profile_helix: axis direction is degenerate".to_string())?;
    if !(helix_radius.is_finite() && helix_radius > 0.0) {
        return Err("sweep_profile_helix: helix radius must be positive".into());
    }
    if !(pitch.is_finite() && pitch > 0.0) {
        return Err("sweep_profile_helix: pitch must be positive".into());
    }
    if !(turns.is_finite() && turns > 0.0) {
        return Err("sweep_profile_helix: turns must be positive".into());
    }
    // 64 turns bounds the loft's O(stations³) interpolation solve; beyond it
    // the station cap would silently degrade per-turn density anyway.
    if turns > 64.0 {
        return Err("sweep_profile_helix: turns must be at most 64".into());
    }

    // --- 2. Profile extent about its centroid, sampled exactly like the path
    //        sweep derives its placement origin (boundary-sample mean), so the
    //        extent is measured about the point that actually rides the path.
    let mut samples = Vec::new();
    for curve in profile {
        let [start, end] = curve.domain()?;
        for sample in 0..16 {
            samples.push(curve.evaluate(start + (end - start) * sample as f64 / 16.0)?);
        }
    }
    if !samples.is_empty() {
        let mut centroid = Vec3::default();
        for point in &samples {
            centroid = centroid.add(*point);
        }
        let origin = centroid.scale(1.0 / samples.len() as f64);
        let extent = samples
            .iter()
            .map(|point| point.sub(origin).length())
            .fold(0.0, f64::max);
        let c = pitch / TAU; // axial rise per radian
                             // Fold-over: profile reaches past the helix's center of curvature.
        let curvature_radius = (helix_radius * helix_radius + c * c) / helix_radius;
        if extent >= curvature_radius {
            return Err(format!(
                "sweep_profile_helix: profile extent {extent:.6} reaches the helix \
                 curvature radius {curvature_radius:.6}; the tube would fold through \
                 itself — increase the helix radius or pitch, or shrink the profile"
            ));
        }
        // Coil collision: only possible once the sweep spans a full revolution.
        if turns >= 1.0 {
            let circumference = TAU * helix_radius;
            let turn_length = (circumference * circumference + pitch * pitch).sqrt();
            let coil_gap = pitch * circumference / turn_length;
            if coil_gap <= 2.0 * extent {
                return Err(format!(
                    "sweep_profile_helix: consecutive turns would self-intersect — \
                     coil gap {coil_gap:.6} does not clear the profile diameter {:.6}; \
                     increase the pitch or shrink the profile",
                    2.0 * extent
                ));
            }
        }
    }

    // --- 3. Fit the helical path (see the density rationale in the fn docs).
    let u = w.perpendicular()?;
    let v = w.cross(u).normalized()?;
    let total_angle = turns * TAU;
    let rise = pitch / TAU;
    let count = ((turns * 64.0).ceil() as usize + 1).clamp(9, 1025);
    let mut points = Vec::with_capacity(count);
    let mut parameters = Vec::with_capacity(count);
    for index in 0..count {
        let s = index as f64 / (count - 1) as f64;
        let theta = total_angle * s;
        points.push(
            axis_origin
                .add(u.scale(helix_radius * theta.cos()))
                .add(v.scale(helix_radius * theta.sin()))
                .add(w.scale(rise * theta)),
        );
        parameters.push(s);
    }
    let path = interpolate_curve(&points, 3, &parameters)
        .map_err(|error| format!("sweep_profile_helix: helix path fit failed: {error}"))?;

    // --- 4. Drive the existing sweep core; its (or the loft's) failures
    //        propagate with helix context prepended.
    let stations = ((turns * 32.0).ceil() as usize).clamp(32, 1024);
    sweep_profile_along_path_stations(profile, &path, name, stations, 0.0, None)
        .map_err(|error| format!("sweep_profile_helix: {error}"))
}