geop_core_geometry/nurb_curve/
extrema.rs

1//! Where a curve's coordinate along a direction is stationary: its extremes
2//! along that direction, and the turning points of its projection.
3//!
4//! For a rational piece `C = H / W`, the coordinate `A / W` along a
5//! direction `n` (`A = n · H`) has derivative `(A' W - A W') / W²`. The
6//! weights are positive, so its sign is that of the numerator
7//! `N = A' W - A W'`, a polynomial. On a Bézier piece of degree `p`, `A'`
8//! has degree `p - 1` and `W` degree `p`, and their product's Bernstein
9//! coefficients follow from the factors' by
10//! `B_i^m B_j^n = C(m, i) C(n, j) / C(m + n, i + j) B_{i+j}^{m+n}` — no
11//! evaluation, no sampling. Its roots are then isolated by subdividing the
12//! coefficients (de Casteljau at the middle): a piece whose coefficients all
13//! have one definite sign has none (the convex hull property).
14
15use geop_core_math::{
16    combinatorics::binomial,
17    geop_error::{GeopError, GeopResult},
18    scalars::Scalar,
19    solvers::subdivide::{Order, Searched, Split, subdivide},
20    vector::Vector,
21};
22
23use super::NurbCurve;
24
25/// The Bernstein coefficients (degree `2p - 1`) of `A' W - A W'`, up to the
26/// positive factor `p / (b - a)`, for a Bézier piece of degree `p` with
27/// coordinate coefficients `a` and weights `w`.
28fn derivative_numerator(a: &[Scalar], w: &[Scalar]) -> GeopResult<Vec<Scalar>> {
29    let p = a.len() - 1;
30    let m = 2 * p - 1;
31    let mut n = vec![Scalar::ZERO; m + 1];
32    for i in 0..p {
33        let (da, dw) = (a[i + 1].sub(a[i]), w[i + 1].sub(w[i]));
34        for j in 0..=p {
35            let factor =
36                Scalar::from_ratio(binomial(p - 1, i) * binomial(p, j), binomial(m, i + j))?;
37            let term = da.mul(w[j]).sub(dw.mul(a[j]));
38            n[i + j] = n[i + j].add(factor.mul(term));
39        }
40    }
41    Ok(n)
42}
43
44/// The Bernstein coefficients of the two halves of a polynomial, split at
45/// the middle of its interval (de Casteljau).
46fn halves(coefficients: &[Scalar]) -> GeopResult<(Vec<Scalar>, Vec<Scalar>)> {
47    let mut row = coefficients.to_vec();
48    let mut left = Vec::with_capacity(row.len());
49    let mut right = Vec::with_capacity(row.len());
50    left.push(row[0]);
51    right.push(row[row.len() - 1]);
52    // Each round one coefficient fewer, down to one.
53    for _ in 1..coefficients.len() {
54        row = row
55            .windows(2)
56            .map(|pair| pair[0].add(pair[1]).div(Scalar::TWO))
57            .collect::<GeopResult<_>>()?;
58        left.push(row[0]);
59        right.push(row[row.len() - 1]);
60    }
61    right.reverse();
62    Ok((left, right))
63}
64
65impl<const D: usize> NurbCurve<D> {
66    /// Every parameter at which the curve's coordinate along `direction`
67    /// is stationary — its tangent perpendicular to `direction` — at an
68    /// isolated point, as a parameter box: the curve's extremes along
69    /// `direction`, and where its projection onto a plane containing
70    /// `direction` turns back on itself. Sorted, and each box at most
71    /// `min_subdivision_size` wide unless several merged.
72    ///
73    /// A stretch along which the coordinate is constant (a straight piece
74    /// perpendicular to `direction`) has no isolated stationary point and
75    /// contributes none. A box is a candidate, not a proof: it is where the
76    /// search stopped being able to exclude one. Exhausting `max_nodes` is an
77    /// error, since the result would be incomplete. `C` must be `D - 1`.
78    pub fn stationary_parameters<const C: usize>(
79        &self,
80        direction: &Vector<C>,
81        max_nodes: usize,
82        min_subdivision_size: Scalar,
83    ) -> GeopResult<Vec<Scalar>> {
84        if C + 1 != D {
85            return Err(GeopError::new(format!(
86                "stationary_parameters: a direction of {C} components for a curve of {} \
87                 dimensions",
88                D - 1
89            )));
90        }
91        // Each Bezier piece's derivative numerator, searched depth first,
92        // one piece after another, within one budget.
93        let mut roots = Vec::new();
94        for piece in self.bezier_pieces()? {
95            let p = piece.degree;
96            if p == 0 || piece.control_points.len() != p + 1 {
97                continue;
98            }
99            let a: Vec<Scalar> = piece
100                .control_points
101                .iter()
102                .map(|q| direction.prod_dot(&q.head::<C>()))
103                .collect();
104            let w: Vec<Scalar> = piece.control_points.iter().map(|q| q[D - 1]).collect();
105            let (lo, hi) = piece.domain();
106            roots.push((derivative_numerator(&a, &w)?, lo, hi));
107        }
108        let mut boxes: Vec<(Scalar, Scalar)> = Vec::new();
109        let searched = subdivide(
110            roots,
111            Order::DepthFirst,
112            max_nodes,
113            |(coefficients, lo, hi), _| -> GeopResult<_> {
114                let positive = coefficients
115                    .iter()
116                    .all(|c| c.definitely_greater(Scalar::ZERO));
117                let negative = coefficients.iter().all(|c| c.definitely_less(Scalar::ZERO));
118                let constant = coefficients.iter().all(|c| c.could_be_equal(Scalar::ZERO));
119                if positive || negative || constant {
120                    return Ok(Split::Done);
121                }
122                if !hi.sub(lo).definitely_greater(min_subdivision_size) {
123                    boxes.push((lo, hi));
124                    return Ok(Split::Done);
125                }
126                // Where to cut is the search's own free choice.
127                let mid = lo.add(hi).div(Scalar::TWO)?.sharpen();
128                let (left, right) = halves(&coefficients)?;
129                Ok(Split::Halves((left, lo, mid), (right, mid, hi)))
130            },
131        )?;
132        if let Searched::Exhausted { .. } = searched {
133            return Err(GeopError::new(format!(
134                "stationary_parameters: exhausted max_nodes={max_nodes} along \
135                 {direction:?}; the result would be incomplete"
136            )));
137        }
138        boxes.sort_by(|x, y| x.0.to_f64().total_cmp(&y.0.to_f64()));
139        let mut merged: Vec<(Scalar, Scalar)> = Vec::new();
140        for (lo, hi) in boxes {
141            match merged.last_mut() {
142                Some(last) if !lo.definitely_greater(last.1) => last.1 = last.1.max(hi),
143                _ => merged.push((lo, hi)),
144            }
145        }
146        Ok(merged.into_iter().map(|(lo, hi)| lo.union(hi)).collect())
147    }
148}
149
150impl NurbCurve<4> {
151    /// The box the curve itself spans, per axis: the union of its
152    /// coordinate at both ends and at every stationary parameter along that
153    /// axis ([`NurbCurve::stationary_parameters`]) — its extremes, each
154    /// enclosed, where the control points' box is only a bound.
155    pub fn bounds(
156        &self,
157        max_nodes: usize,
158        min_subdivision_size: Scalar,
159    ) -> GeopResult<[Scalar; 3]> {
160        let (t0, t1) = self.domain();
161        let mut out = [Scalar::ZERO; 3];
162        for (k, axis) in out.iter_mut().enumerate() {
163            let mut direction = Vector::<3>::zero();
164            direction[k] = Scalar::ONE;
165            let mut extent = self.evaluate(t0)?[k].union(self.evaluate(t1)?[k]);
166            for t in self.stationary_parameters(&direction, max_nodes, min_subdivision_size)? {
167                extent = extent.union(self.evaluate(t)?[k]);
168            }
169            *axis = extent;
170        }
171        Ok(out)
172    }
173}
174
175#[cfg(test)]
176mod tests {
177    use crate::{nurb_curve::NurbCurve, shape::Arc, shape::Circle};
178    use geop_core_math::{
179        scalars::Scalar,
180        vector::{Vector3, Vector4},
181    };
182
183    fn v3(x: f64, y: f64, z: f64) -> Vector3 {
184        Vector3::from_array([
185            Scalar::from_f64(x),
186            Scalar::from_f64(y),
187            Scalar::from_f64(z),
188        ])
189    }
190
191    /// A unit circle in the xy plane is extreme along x at `(±1, 0)` and
192    /// along y at `(0, ±1)`, and nowhere along z.
193    #[test]
194    fn circle_extremes() {
195        let circle = Arc {
196            circle: Circle {
197                center: v3(0.0, 0.0, 0.0),
198                normal: v3(0.0, 0.0, 1.0),
199                radius: Scalar::ONE,
200            },
201            start: v3(1.0, 0.0, 0.0),
202            end: v3(1.0, 0.0, 0.0),
203        }
204        .to_curve()
205        .unwrap();
206        let min = Scalar::from_f64(1e-9);
207        for k in 0..2 {
208            let mut direction = v3(0.0, 0.0, 0.0);
209            direction[k] = Scalar::ONE;
210            let ts = circle.stationary_parameters(&direction, 1000, min).unwrap();
211            let values: Vec<Scalar> = ts
212                .iter()
213                .map(|t| circle.evaluate(t.midpoint()).unwrap()[k])
214                .collect();
215            assert!(
216                values.iter().all(|v| v.abs().could_be_equal(Scalar::ONE)),
217                "{values:?}"
218            );
219            assert!(
220                values.iter().any(|v| v.could_be_equal(Scalar::ONE)),
221                "{values:?}"
222            );
223            assert!(
224                values.iter().any(|v| v.could_be_equal(Scalar::ONE.neg())),
225                "{values:?}"
226            );
227        }
228        let ts = circle
229            .stationary_parameters(&v3(0.0, 0.0, 1.0), 1000, min)
230            .unwrap();
231        assert!(ts.is_empty(), "{ts:?}");
232    }
233
234    /// A straight line has no isolated stationary point along any
235    /// direction, including one it is perpendicular to.
236    #[test]
237    fn line_has_none() {
238        let line = NurbCurve::try_new(
239            1,
240            vec![
241                Vector4::from_array([Scalar::ZERO, Scalar::ZERO, Scalar::ZERO, Scalar::ONE]),
242                Vector4::from_array([Scalar::ZERO, Scalar::TWO, Scalar::ZERO, Scalar::ONE]),
243            ],
244            vec![Scalar::ZERO, Scalar::ZERO, Scalar::ONE, Scalar::ONE],
245        )
246        .unwrap();
247        let min = Scalar::from_f64(1e-9);
248        for k in 0..3 {
249            let mut direction = v3(0.0, 0.0, 0.0);
250            direction[k] = Scalar::ONE;
251            assert!(
252                line.stationary_parameters(&direction, 100, min)
253                    .unwrap()
254                    .is_empty()
255            );
256        }
257    }
258}
geop_core_geometry/nurb_curve/extrema.rs

geop_core_geometry/nurb_curve/
extrema.rs
