geop_core_geometry/nurb_surface/

split.rs

use geop_core_math::{
    geop_error::{GeopError, GeopResult},
    scalars::Scalar,
    vector::Vector,
};

use super::NurbSurface;
use crate::{aabb::compute_aabb, knot_insertion, nurb_curve::NurbCurve};

/// A parameter direction of a tensor-product surface.
#[derive(Clone, Copy)]
enum Dir {
    U,
    V,
}

impl<const D: usize> NurbSurface<D> {
    fn degree_in(&self, dir: Dir) -> usize {
        match dir {
            Dir::U => self.degree_u,
            Dir::V => self.degree_v,
        }
    }

    fn knots_in(&self, dir: Dir) -> &Vec<Scalar> {
        match dir {
            Dir::U => &self.knot_vector_u,
            Dir::V => &self.knot_vector_v,
        }
    }

    /// The control net as rows running along `dir` (see `knot_insertion`):
    /// the `num_v` columns `P[·][j]` for `U`, the `num_u` rows `P[i][·]` for
    /// `V`.
    fn rows_along(&self, dir: Dir) -> Vec<Vec<Vector<D>>> {
        let (nu, nv) = (self.num_u, self.num_v);
        match dir {
            Dir::U => (0..nv)
                .map(|j| (0..nu).map(|i| self.control_points[i * nv + j]).collect())
                .collect(),
            Dir::V => self.control_points.chunks(nv).map(<[_]>::to_vec).collect(),
        }
    }

    /// This surface with the `dir` knots and rows replaced — the inverse of
    /// [`Self::rows_along`].
    fn with_rows(&self, dir: Dir, knots: Vec<Scalar>, rows: Vec<Vec<Vector<D>>>) -> Self {
        let along = rows[0].len();
        let (num_u, num_v, control_points): (usize, usize, Vec<Vector<D>>) = match dir {
            Dir::U => (
                along,
                rows.len(),
                (0..along)
                    .flat_map(|i| rows.iter().map(move |col| col[i]))
                    .collect(),
            ),
            Dir::V => (rows.len(), along, rows.into_iter().flatten().collect()),
        };
        let (knot_vector_u, knot_vector_v) = match dir {
            Dir::U => (knots, self.knot_vector_v.clone()),
            Dir::V => (self.knot_vector_u.clone(), knots),
        };
        NurbSurface {
            degree_u: self.degree_u,
            degree_v: self.degree_v,
            num_u,
            num_v,
            aabb: compute_aabb(&control_points),
            control_points,
            knot_vector_u,
            knot_vector_v,
        }
    }

    fn split_in(&self, dir: Dir, t: Scalar) -> GeopResult<(Self, Self)> {
        let mut knots = self.knots_in(dir).clone();
        let mut rows = self.rows_along(dir);
        let (right_knots, right_rows) =
            knot_insertion::split(&mut knots, &mut rows, self.degree_in(dir), t)?;
        Ok((
            self.with_rows(dir, knots, rows),
            self.with_rows(dir, right_knots, right_rows),
        ))
    }

    fn restrict_in(&self, dir: Dir, t0: Scalar, t1: Scalar) -> GeopResult<Self> {
        let mut knots = self.knots_in(dir).clone();
        let mut rows = self.rows_along(dir);
        knot_insertion::restrict(&mut knots, &mut rows, self.degree_in(dir), t0, t1)?;
        Ok(self.with_rows(dir, knots, rows))
    }

    /// Split at parameter `t` in the **u** direction.
    ///
    /// `t` must lie strictly inside the u domain.  Returns `(left, right)`.
    ///
    /// `t` is used exactly as given — **not** sharpened here, see `NurbCurve::split`'s own
    /// doc comment for why (this is the same Boehm-insertion construction,
    /// one dimension up: an unsharpened `t` carried into the new knot
    /// vector lets `alpha = (t - e) / (s - e)` blow up over repeated splits
    /// as `s - e` shrinks while `t`'s own width doesn't). Every current caller
    /// therefore sharpens its own midpoint before calling; a caller splitting
    /// at a *located* parameter must validate the sharpened value it is about
    /// to use, not the wide one it started from.
    pub fn split_u(&self, t: Scalar) -> GeopResult<(NurbSurface<D>, NurbSurface<D>)> {
        self.split_in(Dir::U, t)
    }

    /// Split at parameter `t` in the **v** direction.
    ///
    /// `t` must lie strictly inside the v domain.  Returns `(left, right)`.
    /// `t` is used exactly as given — **not** sharpened here, same as [`Self::split_u`].
    pub fn split_v(&self, t: Scalar) -> GeopResult<(NurbSurface<D>, NurbSurface<D>)> {
        self.split_in(Dir::V, t)
    }

    fn split_in_mid(&self, dir: Dir) -> GeopResult<(Self, Self)> {
        let knots = self.knots_in(dir);
        let n = match dir {
            Dir::U => self.num_u,
            Dir::V => self.num_v,
        };
        let (t0, t1) = (knots[self.degree_in(dir)], knots[n]);
        // A self-chosen subdivision point: any value in the interval cuts
        // it equally well, so sharpening loses no accuracy and keeps
        // repeated splits from compounding width (see AGENTS.md).
        let mid = t0.add(t1).div(Scalar::TWO)?.sharpen();
        self.split_in(dir, mid)
    }

    /// The iso-line of the surface at `u = at` (`u_fixed`), or at `v = at`,
    /// as a curve along the other parameter, on its knots: a row of the
    /// control net, exactly, where `at` is a clamped end of the domain or a
    /// knot of full multiplicity, and otherwise the row knot insertion
    /// makes there. `at` is used exactly as given, as by [`Self::split_u`].
    /// Fails for an `at` outside the domain, or at an end that is not
    /// clamped.
    pub fn iso_curve(&self, u_fixed: bool, at: Scalar) -> GeopResult<NurbCurve<D>> {
        let (dir, other) = if u_fixed {
            (Dir::U, Dir::V)
        } else {
            (Dir::V, Dir::U)
        };
        let (knots, degree) = (self.knots_in(dir), self.degree_in(dir));
        let n = self.rows_along(dir)[0].len();
        let (lo, hi) = (knots[degree], knots[n]);
        let ends = |first: bool| -> GeopResult<Vec<Vector<D>>> {
            if !knot_insertion::is_clamped_at(knots, degree, first) {
                return Err(GeopError::new(format!(
                    "NurbSurface::iso_curve: the domain is not clamped at {at:?}"
                )));
            }
            let i = if first { 0 } else { n - 1 };
            Ok(self.rows_along(dir).iter().map(|row| row[i]).collect())
        };
        let points = if at.is_exactly(lo) {
            ends(true)?
        } else if at.is_exactly(hi) {
            ends(false)?
        } else {
            let (before, _) = self.split_in(dir, at)?;
            let rows = before.rows_along(dir);
            rows.iter().map(|row| row[row.len() - 1]).collect()
        };
        NurbCurve::try_new(self.degree_in(other), points, self.knots_in(other).clone())
    }

    /// Split at the midpoint of the **u** domain.
    pub fn split_u_mid(&self) -> GeopResult<(NurbSurface<D>, NurbSurface<D>)> {
        self.split_in_mid(Dir::U)
    }

    /// Split at the midpoint of the **v** domain.
    pub fn split_v_mid(&self) -> GeopResult<(NurbSurface<D>, NurbSurface<D>)> {
        self.split_in_mid(Dir::V)
    }

    /// The same surface with every span split into `pieces[0]` along `u`
    /// and `pieces[1]` along `v`, at sharp knots: more control points to
    /// fit on.
    pub(crate) fn refined(&self, pieces: [usize; 2]) -> GeopResult<Self> {
        let mut out = self.clone();
        for (dir, pieces) in [Dir::U, Dir::V].into_iter().zip(pieces) {
            let degree = out.degree_in(dir);
            let mut knots = out.knots_in(dir).clone();
            let mut rows = out.rows_along(dir);
            let breaks = crate::spline::breakpoints(&knots, degree, rows[0].len());
            for pair in breaks.windows(2) {
                for k in 1..pieces {
                    let f = Scalar::from_ratio(k as i64, pieces as i64)?;
                    knot_insertion::insert(
                        &mut knots,
                        &mut rows,
                        degree,
                        Scalar::fraction(pair[0], pair[1], f),
                    );
                }
            }
            out = out.with_rows(dir, knots, rows);
        }
        Ok(out)
    }

    /// The same surface with the knots `at` inserted along `u`, each once:
    /// more control points to fit on where they are.
    pub(crate) fn with_knots_u(&self, at: &[Scalar]) -> Self {
        let mut knots = self.knot_vector_u.clone();
        let mut rows = self.rows_along(Dir::U);
        for &t in at {
            knot_insertion::insert(&mut knots, &mut rows, self.degree_u, t);
        }
        self.with_rows(Dir::U, knots, rows)
    }

    /// This surface restricted to `u ∈ [u0, u1]`, `v ∈ [v0, v1]`, with the
    /// same cut rule as `NurbCurve::sub_curve`: only bounds strictly inside
    /// the domain are cut at, so the result covers at least the requested
    /// box ∩ domain. Bounds are used exactly as given and must be narrow.
    pub fn sub_surface(
        &self,
        (u0, u1): (Scalar, Scalar),
        (v0, v1): (Scalar, Scalar),
    ) -> GeopResult<Self> {
        self.restrict_in(Dir::U, u0, u1)?
            .restrict_in(Dir::V, v0, v1)
    }
}

#[cfg(test)]
mod tests {
    use crate::test_fixtures::bilinear;
    use geop_core_math::scalars::Scalar;

    #[test]
    fn split_u_junction_matches() {
        let orig = bilinear();
        let (left, right) = orig.split_u(Scalar::from_f64(0.5)).unwrap();
        for v_val in [0.0_f64, 0.5, 1.0] {
            let v = Scalar::from_f64(v_val);
            let o = orig.evaluate(Scalar::from_f64(0.5), v).unwrap();
            let l = left.evaluate(Scalar::from_f64(0.5), v).unwrap();
            let r = right.evaluate(Scalar::from_f64(0.5), v).unwrap();
            for c in 0..3 {
                assert!(
                    o[c].could_be_equal(l[c]),
                    "split_u left mismatch v={v_val} c={c}"
                );
                assert!(
                    o[c].could_be_equal(r[c]),
                    "split_u right mismatch v={v_val} c={c}"
                );
            }
        }
    }

    #[test]
    fn split_u_left_start_matches() {
        let orig = bilinear();
        let (left, _) = orig.split_u(Scalar::from_f64(0.5)).unwrap();
        for v_val in [0.0_f64, 0.5, 1.0] {
            let v = Scalar::from_f64(v_val);
            let o = orig.evaluate(Scalar::ZERO, v).unwrap();
            let l = left.evaluate(Scalar::ZERO, v).unwrap();
            for c in 0..3 {
                assert!(o[c].could_be_equal(l[c]), "c={c} v={v_val}");
            }
        }
    }

    #[test]
    fn split_u_right_end_matches() {
        let orig = bilinear();
        let (_, right) = orig.split_u(Scalar::from_f64(0.5)).unwrap();
        for v_val in [0.0_f64, 0.5, 1.0] {
            let v = Scalar::from_f64(v_val);
            let o = orig.evaluate(Scalar::ONE, v).unwrap();
            let r = right.evaluate(Scalar::ONE, v).unwrap();
            for c in 0..3 {
                assert!(o[c].could_be_equal(r[c]), "c={c} v={v_val}");
            }
        }
    }

    #[test]
    fn split_u_at_boundary_returns_err() {
        let s = bilinear();
        assert!(s.split_u(Scalar::ZERO).is_err());
        assert!(s.split_u(Scalar::ONE).is_err());
    }

    #[test]
    fn split_v_junction_matches() {
        let orig = bilinear();
        let (left, right) = orig.split_v(Scalar::from_f64(0.5)).unwrap();
        for u_val in [0.0_f64, 0.5, 1.0] {
            let u = Scalar::from_f64(u_val);
            let o = orig.evaluate(u, Scalar::from_f64(0.5)).unwrap();
            let l = left.evaluate(u, Scalar::from_f64(0.5)).unwrap();
            let r = right.evaluate(u, Scalar::from_f64(0.5)).unwrap();
            for c in 0..3 {
                assert!(
                    o[c].could_be_equal(l[c]),
                    "split_v left mismatch u={u_val} c={c}"
                );
                assert!(
                    o[c].could_be_equal(r[c]),
                    "split_v right mismatch u={u_val} c={c}"
                );
            }
        }
    }

    #[test]
    fn split_v_left_start_matches() {
        let orig = bilinear();
        let (left, _) = orig.split_v(Scalar::from_f64(0.5)).unwrap();
        for u_val in [0.0_f64, 0.5, 1.0] {
            let u = Scalar::from_f64(u_val);
            let o = orig.evaluate(u, Scalar::ZERO).unwrap();
            let l = left.evaluate(u, Scalar::ZERO).unwrap();
            for c in 0..3 {
                assert!(o[c].could_be_equal(l[c]), "c={c} u={u_val}");
            }
        }
    }

    #[test]
    fn split_v_right_end_matches() {
        let orig = bilinear();
        let (_, right) = orig.split_v(Scalar::from_f64(0.5)).unwrap();
        for u_val in [0.0_f64, 0.5, 1.0] {
            let u = Scalar::from_f64(u_val);
            let o = orig.evaluate(u, Scalar::ONE).unwrap();
            let r = right.evaluate(u, Scalar::ONE).unwrap();
            for c in 0..3 {
                assert!(o[c].could_be_equal(r[c]), "c={c} u={u_val}");
            }
        }
    }

    #[test]
    fn split_v_at_boundary_returns_err() {
        let s = bilinear();
        assert!(s.split_v(Scalar::ZERO).is_err());
        assert!(s.split_v(Scalar::ONE).is_err());
    }

    #[test]
    fn split_u_domains() {
        let orig = bilinear();
        let (left, right) = orig.split_u(Scalar::from_f64(0.5)).unwrap();
        let (lu_min, lu_max) = left.domain_u();
        let (ru_min, ru_max) = right.domain_u();
        assert!(lu_min.could_be_equal(Scalar::ZERO));
        assert!(lu_max.could_be_equal(Scalar::from_f64(0.5)));
        assert!(ru_min.could_be_equal(Scalar::from_f64(0.5)));
        assert!(ru_max.could_be_equal(Scalar::ONE));
    }

    #[test]
    fn split_v_domains() {
        let orig = bilinear();
        let (left, right) = orig.split_v(Scalar::from_f64(0.5)).unwrap();
        let (lv_min, lv_max) = left.domain_v();
        let (rv_min, rv_max) = right.domain_v();
        assert!(lv_min.could_be_equal(Scalar::ZERO));
        assert!(lv_max.could_be_equal(Scalar::from_f64(0.5)));
        assert!(rv_min.could_be_equal(Scalar::from_f64(0.5)));
        assert!(rv_max.could_be_equal(Scalar::ONE));
    }
}