geop_core_geometry/nurb_curve/

helix.rs

//! Helices: [`NurbCurve3D::helix`], and the control rows it is built from
//! ([`helix_rows`]), which a screw sweep carries a whole profile along.
//!
//! A helix is not a rational curve: its height grows linearly with the angle
//! it has turned, and no rational parametrization of a circle turns at a
//! constant rate. So it is approximated, the same way every circle of the
//! kernel is built: one rational quadratic per span of at most a quarter
//! turn, its middle control point where the tangents of the circle meet,
//! weighted `cos(step / 2)`. That part is exact, so **the curve lies exactly
//! on its cylinder**: seen down the axis it is the kernel's own circle. Only
//! the height is approximate. It rises from span end to span end exactly as
//! a helix does, the middle control row sits at the height halfway between
//! (where the curve, by symmetry, passes at the span's middle angle, exactly
//! on the helix), and in between the height runs ahead of or behind the true
//! helix by at most **0.53 % of the pitch** for quarter-turn spans (`5.3e-3`
//! times the pitch, measured over the span; see the test). Any other middle height
//! is worse: the error is odd about the middle of the span, and moving the
//! middle row adds an even term.
//!
//! For an ISO metric thread that is 5 µm on an M6 (pitch 1 mm), well inside
//! the tolerance of the thread itself.
//!
//! A helix may also be **tapered** ([`Helix::taper`]), its radius growing in
//! step with its height, as a pipe thread's (NPT) does: on a cone instead of
//! a cylinder. The radius is then built as the height is, exactly at every
//! span's end and its middle row at the mean of the two, so the curve leaves
//! the cone between them as the height leaves the helix, by a fraction of
//! how much the radius grows in a quarter turn (see the test).

use geop_core_math::{
    geop_error::{GeopError, GeopResult},
    primitives::CoordinateSystem,
    scalars::Scalar,
    vector::{Vector2, Vector3},
};

use super::{NurbCurve, NurbCurve3D};

/// Which way a helix turns as it rises: right-handed, like an ordinary
/// screw — counter-clockwise seen from where it rises to — or left-handed.
#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
pub enum Handedness {
    #[default]
    Right,
    Left,
}

/// The screw motion a helix follows: rising `pitch` per turn through
/// `turns` turns, turning as `handedness` says, and moving out from the axis
/// by `taper` for every unit it rises — none for a helix on a cylinder.
#[derive(Clone, Copy, Debug)]
pub struct Helix {
    pub pitch: Scalar,
    pub turns: f64,
    pub handedness: Handedness,
    pub taper: Scalar,
}

impl Helix {
    /// The motion along a helix on a cylinder.
    pub fn new(pitch: Scalar, turns: f64, handedness: Handedness) -> Self {
        Self {
            pitch,
            turns,
            handedness,
            taper: Scalar::ZERO,
        }
    }

    /// This motion, moving out from the axis by `taper` for every unit it
    /// rises: on a cone whose radius grows by `taper` per unit of height,
    /// half a pipe thread's taper on the diameter.
    pub fn tapered(self, taper: Scalar) -> Self {
        Self { taper, ..self }
    }
}

/// One control row of a helix around `frame.w()` through `frame.origin()`:
/// a point at distance `r` from the axis, raised by `z` above the row, lies
/// at `origin + (r + out) radial + (rise + z) w`, weighted `weight` — none
/// for a row at a span's end, whose weight is one.
#[derive(Clone, Debug)]
pub struct HelixRow {
    /// Unit length at a span's end; at a span's middle, as far out as the
    /// tangents of the unit circle at the span's ends meet.
    pub radial: Vector3,
    pub rise: Scalar,
    /// How far a tapered helix has moved out from its axis here: its taper
    /// times the rise. None on a cylinder.
    pub out: Option<Scalar>,
    pub weight: Option<Scalar>,
}

/// `(cos, sin)` of `degrees` — exact at every quarter turn, so that a
/// revolve through a right angle lands exactly on the plane it should, and
/// a helix's spans of a quarter turn sit on exactly the circle a revolve's
/// do.
pub fn cos_sin(degrees: f64) -> (f64, f64) {
    match degrees.rem_euclid(360.0) {
        0.0 => (1.0, 0.0),
        90.0 => (0.0, 1.0),
        180.0 => (-1.0, 0.0),
        270.0 => (0.0, -1.0),
        d => (d.to_radians().cos(), d.to_radians().sin()),
    }
}

/// `sqrt(2) / 2`, the weight of the middle control point that makes a
/// rational quadratic a quarter of a circle exactly.
pub fn sqrt2_over_2() -> Scalar {
    Scalar::from_f64(std::f64::consts::SQRT_2 / 2.0)
}

/// The middle control point of an arc through `step` degrees, at most a
/// half turn either way: `(weight, outwards)`. It lies where the tangents
/// at the arc's ends meet, `outwards = 1 / (1 + cos step)` along the sum of
/// the unit directions to the ends, weighted `cos(step / 2)` —
/// [`sqrt2_over_2`] for a quarter turn. Revolves and helices both build
/// their spans with it, so a helix's quarter lies on the very circle a
/// revolved quarter does.
pub fn arc_span(step: f64) -> (Scalar, Scalar) {
    let (cos_step, _) = cos_sin(step);
    let weight = if step.abs() == 90.0 {
        sqrt2_over_2()
    } else {
        Scalar::from_f64((step / 2.0).to_radians().cos())
    };
    (weight, Scalar::from_f64(1.0 / (1.0 + cos_step)))
}

/// The control rows of a helix around `frame.w()`, from `frame.u()`, moving
/// as `helix` says along `w`: `2 n + 1` rows for `n` equal spans of at most
/// a quarter turn — a span's ends at even indices, its middle in between.
pub fn helix_rows(frame: &CoordinateSystem, helix: &Helix) -> GeopResult<Vec<HelixRow>> {
    let Helix {
        pitch,
        turns,
        handedness,
        taper,
    } = *helix;
    if !(turns.is_finite() && turns > 0.0) {
        return Err(GeopError::new(format!(
            "helix: cannot turn {turns} times: it must be more than none"
        )));
    }
    if !pitch.definitely_greater(Scalar::ZERO) {
        return Err(GeopError::new(format!(
            "helix: the pitch {pitch:?} must be definitely positive"
        )));
    }
    let spans = (turns * 4.0).ceil() as usize;
    let sign = match handedness {
        Handedness::Right => 1.0,
        Handedness::Left => -1.0,
    };
    let step = sign * 360.0 * turns / spans as f64;
    let direction = |degrees: f64| {
        let (c, s) = cos_sin(degrees);
        frame.uv_direction(&Vector2::from_array([
            Scalar::from_f64(c),
            Scalar::from_f64(s),
        ]))
    };
    let (weight, outwards) = arc_span(step);
    let total = pitch.mul(Scalar::from_f64(turns));
    // The height at the end of span `k`: none at the start and all of it at
    // the end, taken as they are rather than multiplied by 0 or 1.
    let rise = |k: usize| -> GeopResult<Scalar> {
        Ok(match k {
            0 => Scalar::ZERO,
            k if k == spans => total,
            k => total.mul(Scalar::from_ratio(k as i64, spans as i64)?),
        })
    };
    // Out by the taper times the rise, built as the rise is: none where
    // there is no taper, rather than multiplied by zero.
    let out = |rise: Scalar| (!taper.is_zero()).then(|| taper.mul(rise));
    let mut rows = Vec::with_capacity(2 * spans + 1);
    let mut start = HelixRow {
        radial: direction(0.0),
        rise: Scalar::ZERO,
        out: out(Scalar::ZERO),
        weight: None,
    };
    for k in 0..spans {
        let rise_end = rise(k + 1)?;
        let end = HelixRow {
            radial: direction((k + 1) as f64 * step),
            rise: rise_end,
            out: out(rise_end),
            weight: None,
        };
        let rise_middle = start.rise.add(end.rise).div(Scalar::TWO)?;
        let middle = HelixRow {
            radial: start.radial.add(&end.radial).prod_scalar(outwards),
            rise: rise_middle,
            out: out(rise_middle),
            weight: Some(weight),
        };
        rows.push(start);
        rows.push(middle);
        start = end;
    }
    rows.push(start);
    Ok(rows)
}

impl NurbCurve3D {
    /// The helix of `radius` around `frame.w()` through `frame.origin()`:
    /// starting at `origin + radius u`, moving as `helix` says along `w` —
    /// right-handed turning from `u` towards `v`, left-handed the other
    /// way. One rational quadratic per span of at most a quarter turn,
    /// exactly on the cylinder, its height within 0.53 % of the pitch of
    /// the true helix (see the module docs); tapered, on the cone as
    /// closely. The parameter runs over `[0, 1]`, a span's end at every
    /// `k / n`. Refused where the radius could come to nothing.
    pub fn helix(frame: &CoordinateSystem, radius: Scalar, helix: &Helix) -> GeopResult<Self> {
        let rows = helix_rows(frame, helix)?;
        let radius_at = |row: &HelixRow| row.out.map_or(radius, |out| radius.add(out));
        if let Some(row) = rows
            .iter()
            .find(|row| !radius_at(row).definitely_greater(Scalar::ZERO))
        {
            return Err(GeopError::new(format!(
                "helix: the radius {:?} must be definitely positive, starting at {radius:?}",
                radius_at(row)
            )));
        }
        let spans = rows.len() / 2;
        let (origin, w) = (frame.origin(), frame.w());
        let control_points = rows
            .iter()
            .map(|row| {
                origin
                    .add(&row.radial.prod_scalar(radius_at(row)))
                    .add(&w.prod_scalar(row.rise))
                    .homogeneous(row.weight.unwrap_or(Scalar::ONE))
            })
            .collect();
        let mut knots = vec![Scalar::ZERO; 3];
        for k in 1..spans {
            let knot = Scalar::from_ratio(k as i64, spans as i64)?;
            knots.extend([knot, knot]);
        }
        knots.extend([Scalar::ONE; 3]);
        NurbCurve::try_new(2, control_points, knots)
    }
}

#[cfg(test)]
mod tests {
    use geop_core_math::{primitives::CoordinateSystem, scalars::Scalar, vector::Vector3};

    use super::*;

    fn frame() -> CoordinateSystem {
        CoordinateSystem::world_at(Vector3::from_array([1.0, 2.0, 3.0].map(Scalar::from_f64)))
    }

    /// Every point of the helix lies on its cylinder, it ends where it
    /// should, and its height stays within 0.53 % of the pitch of the true
    /// helix's at the same angle — turning the right way round.
    #[test]
    fn helix_is_on_its_cylinder_and_close_to_the_true_helix() {
        let (radius, pitch, turns) = (2.5, 0.8, 2.3);
        for handedness in [Handedness::Right, Handedness::Left] {
            let curve = NurbCurve3D::helix(
                &frame(),
                Scalar::from_f64(radius),
                &Helix::new(Scalar::from_f64(pitch), turns, handedness),
            )
            .unwrap();
            let sign = if handedness == Handedness::Right {
                1.0
            } else {
                -1.0
            };
            let samples = 997;
            let mut unwrapped = 0.0f64;
            let mut last = 0.0f64;
            let mut worst = 0.0f64;
            for i in 0..=samples {
                let t = Scalar::from_ratio(i, samples).unwrap();
                let p = curve.evaluate(t).unwrap();
                let centre = Vector3::from_array([1.0, 2.0, p[2].to_f64()].map(Scalar::from_f64));
                let off = p.sub(&centre);
                let r2 = off[0].mul(off[0]).add(off[1].mul(off[1]));
                assert!(
                    r2.could_be_equal(Scalar::from_f64(radius * radius)),
                    "{p:?} is off the cylinder"
                );
                let (x, y, z) = (
                    p[0].to_f64() - 1.0,
                    p[1].to_f64() - 2.0,
                    p[2].to_f64() - 3.0,
                );
                let angle = sign * y.atan2(x);
                let mut d = angle - last;
                if d < -std::f64::consts::PI {
                    d += 2.0 * std::f64::consts::PI;
                }
                if d > std::f64::consts::PI {
                    d -= 2.0 * std::f64::consts::PI;
                }
                // It never turns back.
                assert!(d >= -1e-12, "{i}: {d}");
                unwrapped += d;
                last = angle;
                let expected = pitch * unwrapped / (2.0 * std::f64::consts::PI);
                worst = worst.max((z - expected).abs());
            }
            assert!(
                (unwrapped - turns * 2.0 * std::f64::consts::PI).abs() < 1e-9,
                "{unwrapped}"
            );
            assert!(worst <= 5.3e-3 * pitch, "{worst}");
            // Not merely a circle lifted: it does rise with the angle.
            assert!(worst > 1e-4 * pitch, "{worst}");
            let end = curve.evaluate(Scalar::ONE).unwrap();
            assert!(
                end[2].could_be_equal(
                    Scalar::from_f64(3.0).add(Scalar::from_f64(pitch).mul(Scalar::from_f64(turns)))
                ),
                "{end:?}"
            );
        }
    }

    /// Nothing to turn, no pitch, no radius: refused.
    #[test]
    fn degenerate_helices_are_refused() {
        use geop_core_math::scalars::Scalar;
        let f = frame();
        let one = Scalar::from_f64(1.0);
        let right = |pitch: Scalar, turns: f64| Helix::new(pitch, turns, Handedness::Right);
        assert!(NurbCurve3D::helix(&f, one, &right(one, 0.0)).is_err());
        assert!(NurbCurve3D::helix(&f, one, &right(Scalar::ZERO, 1.0)).is_err());
        assert!(NurbCurve3D::helix(&f, Scalar::ZERO, &right(one, 1.0)).is_err());
        // Tapered in so far that it would reach its axis.
        let inwards = right(one, 2.0).tapered(Scalar::from_f64(-0.6));
        assert!(NurbCurve3D::helix(&f, one, &inwards).is_err());
    }

    /// A pipe thread's helix, tapered 1 in 32 on the radius: it starts and
    /// ends where the cone has it, at every quarter turn it is on the cone,
    /// and in between it stays within a hundredth of what the radius grows
    /// in a quarter turn.
    #[test]
    fn tapered_helix_is_on_its_cone_at_every_quarter_and_close_between() {
        let (radius, pitch, turns, taper) = (10.0, 1.8, 3.5, 1.0 / 32.0);
        let helix = Helix::new(Scalar::from_f64(pitch), turns, Handedness::Right)
            .tapered(Scalar::from_f64(taper));
        let curve = NurbCurve3D::helix(&frame(), Scalar::from_f64(radius), &helix).unwrap();
        let off_cone = |p: &Vector3| {
            let (x, y, z) = (
                p[0].to_f64() - 1.0,
                p[1].to_f64() - 2.0,
                p[2].to_f64() - 3.0,
            );
            x.hypot(y) - (radius + taper * z)
        };
        let spans = (turns * 4.0).ceil() as i64;
        for k in 0..=spans {
            let p = curve
                .evaluate(Scalar::from_ratio(k, spans).unwrap())
                .unwrap();
            assert!(off_cone(&p).abs() < 1e-12, "quarter {k}: {p:?}");
        }
        let quarter = taper * pitch / 4.0;
        let samples = 997;
        let worst = (0..=samples)
            .map(|i| {
                off_cone(
                    &curve
                        .evaluate(Scalar::from_ratio(i, samples).unwrap())
                        .unwrap(),
                )
            })
            .fold(0.0f64, |w, d| w.max(d.abs()));
        assert!(worst < 1e-2 * quarter, "{worst} against {quarter}");
        assert!(worst > 0.0);
        let end = curve.evaluate(Scalar::ONE).unwrap();
        assert!(
            off_cone(&end).abs() < 1e-12 && (end[2].to_f64() - 3.0 - pitch * turns).abs() < 1e-12,
            "{end:?}"
        );
    }
}