geop_core_geometry/nurb_surface/
carry.rs

1//! One surface carried on over another: whether a patch lies on another
2//! patch's surface, and where ([`NurbSurface3D::carrying`]); a patch
3//! carried on until it reaches given points ([`NurbSurface3D::covering`]).
4
5use geop_core_math::{
6    geop_error::{GeopError, GeopResult},
7    scalars::Scalar,
8    vector::{Vector2, Vector3},
9};
10
11use super::NurbSurface3D;
12use crate::{nurb_curve::NurbCurve2D, search::PROJECT_ITERATIONS};
13
14/// How often [`NurbSurface3D::covering`] carries a surface on, at most: each
15/// time by as much as its domain is long, so it grows past any reach in a
16/// few.
17const MOST_COVERINGS: usize = 8;
18
19/// An affine map of `(u, v)` that takes each parameter to one parameter:
20/// `x' = m x + shift`, `m` diagonal, or — the parameters swapped —
21/// anti-diagonal.
22#[derive(Clone, Copy, Debug)]
23pub struct UvMap {
24    pub m: [[Scalar; 2]; 2],
25    pub shift: Vector2,
26}
27
28impl UvMap {
29    /// Where `uv` goes.
30    pub fn apply(&self, uv: Vector2) -> Vector2 {
31        Vector2::from_array(std::array::from_fn(|k| {
32            self.m[k][0]
33                .mul(uv[0])
34                .add(self.m[k][1].mul(uv[1]))
35                .add(self.shift[k])
36        }))
37    }
38
39    /// Which parameter the `k`-th maps to.
40    fn image(&self, k: usize) -> usize {
41        if self.m[0][k].could_be_equal(Scalar::ZERO) {
42            1
43        } else {
44            0
45        }
46    }
47
48    /// The pcurve `curve` mapped: each homogeneous control point `(w x, w)`
49    /// to `(m (w x) + w shift, w)`, linear in it, so the mapped curve is
50    /// exactly the curve mapped, rational or not.
51    pub fn apply_to(&self, curve: &NurbCurve2D) -> GeopResult<NurbCurve2D> {
52        // The identity multiplies nothing: rounding would only widen.
53        let exactly = |s: Scalar, x: f64| s.is_sharp() && s.to_f64() == x;
54        if exactly(self.m[0][0], 1.0)
55            && exactly(self.m[1][1], 1.0)
56            && exactly(self.m[0][1], 0.0)
57            && exactly(self.m[1][0], 0.0)
58            && exactly(self.shift[0], 0.0)
59            && exactly(self.shift[1], 0.0)
60        {
61            return Ok(curve.clone());
62        }
63        let points = (curve.control_points.iter())
64            .map(|cp| {
65                let w = cp[2];
66                let mut mapped = *cp;
67                for k in 0..2 {
68                    mapped[k] = self.m[k][0]
69                        .mul(cp[0])
70                        .add(self.m[k][1].mul(cp[1]))
71                        .add(self.shift[k].mul(w));
72                }
73                mapped
74            })
75            .collect();
76        NurbCurve2D::try_new(curve.degree, points, curve.knot_vector.clone())
77    }
78}
79
80/// Where `point` is on `surface`, from the middle of its domain: its foot
81/// point, if it lies on the surface there.
82fn on(surface: &NurbSurface3D, point: Vector3) -> GeopResult<Option<Vector2>> {
83    let ((u0, u1), (v0, v1)) = (surface.domain_u(), surface.domain_v());
84    let half = Scalar::from_f64(0.5);
85    let (u, v) = surface.project_in_domain(
86        point,
87        Scalar::interpolate(u0, u1, half).sharpen(),
88        Scalar::interpolate(v0, v1, half).sharpen(),
89        PROJECT_ITERATIONS,
90    )?;
91    Ok((surface.evaluate(u, v)?.could_be_equal(&point)).then(|| Vector2::from_array([u, v])))
92}
93
94impl NurbSurface3D {
95    /// This surface carried on past its sides ([`NurbSurface3D::extended`])
96    /// until the foot point of each of `points` on it lies inside its
97    /// domain, not on a side: a surface something is to lie on has to reach
98    /// it. Carried on past a side each time by as much as the domain is
99    /// long that way — how much is a free choice — and refused, naming the
100    /// points, where it cannot be carried on far enough.
101    pub fn covering(&self, points: &[Vector3]) -> GeopResult<Self> {
102        let mut surface = self.clone();
103        for _ in 0..MOST_COVERINGS {
104            let mut sides = [[false; 2]; 2];
105            for point in points {
106                let ((u0, u1), (v0, v1)) = (surface.domain_u(), surface.domain_v());
107                let half = Scalar::from_f64(0.5);
108                let (u, v) = surface.project_in_domain(
109                    *point,
110                    Scalar::interpolate(u0, u1, half).sharpen(),
111                    Scalar::interpolate(v0, v1, half).sharpen(),
112                    PROJECT_ITERATIONS,
113                )?;
114                for (k, (t, lo, hi)) in [(u, u0, u1), (v, v0, v1)].into_iter().enumerate() {
115                    sides[k][0] |= t.could_be_equal(lo);
116                    sides[k][1] |= t.could_be_equal(hi);
117                }
118            }
119            if sides.iter().flatten().all(|&side| !side) {
120                return Ok(surface);
121            }
122            for (k, ends) in sides.iter().enumerate() {
123                for (at_end, &short) in ends.iter().enumerate() {
124                    if !short {
125                        continue;
126                    }
127                    let (lo, hi) = if k == 0 {
128                        surface.domain_u()
129                    } else {
130                        surface.domain_v()
131                    };
132                    surface = surface.extended(k == 0, at_end == 1, hi.sub(lo).sharpen())?;
133                }
134            }
135        }
136        Err(GeopError::new(format!(
137            "the surface, carried on {MOST_COVERINGS} times, does not reach the points {points:?}"
138        )))
139    }
140
141    /// Where the patch `other` lies on this surface, if it does: this
142    /// surface carried on until it covers `other` ([`Self::covering`]), and
143    /// the map from `other`'s `(u, v)` to its own — the very same surface,
144    /// each parameter of `other` one of its own shifted and scaled, the way
145    /// patches of one surface cut apart, mirrored or built again are, its
146    /// normal the same way. `None` where `other` is another surface, faces
147    /// the other way, or is parametrized otherwise.
148    ///
149    /// The map comes from where `other`'s corners are on this surface. That
150    /// it holds everywhere is a polynomial identity: on each cell between
151    /// the two patches' knots, both are rational of the same degree `p` in
152    /// each parameter, and two such that agree at `2p + 1` parameters in
153    /// each direction agree everywhere — their difference, cleared of
154    /// denominators, is a polynomial of degree `2p` with more roots than
155    /// its degree. So they are compared there, and only there; within the
156    /// enclosures' widths, as everything is.
157    pub fn carrying(&self, other: &Self) -> GeopResult<Option<(Self, UvMap)>> {
158        // The very same patch — faces cut from one face — needs no map.
159        if self.could_be_equal(other) {
160            let identity = UvMap {
161                m: [[Scalar::ONE, Scalar::ZERO], [Scalar::ZERO, Scalar::ONE]],
162                shift: Vector2::zero(),
163            };
164            return Ok(Some((self.clone(), identity)));
165        }
166        let ((u0, u1), (v0, v1)) = (other.domain_u(), other.domain_v());
167        let corners = [
168            other.evaluate(u0, v0)?,
169            other.evaluate(u1, v0)?,
170            other.evaluate(u0, v1)?,
171            other.evaluate(u1, v1)?,
172        ];
173        let Ok(carried) = self.covering(&corners) else {
174            return Ok(None);
175        };
176        let mut at = Vec::with_capacity(3);
177        for &corner in &corners[..3] {
178            match on(&carried, corner)? {
179                Some(uv) => at.push(uv),
180                None => return Ok(None),
181            }
182        }
183        let (p00, p10, p01) = (at[0], at[1], at[2]);
184        let (du, dv) = (p10.sub(&p00), p01.sub(&p00));
185        let zero = |s: Scalar| s.could_be_equal(Scalar::ZERO);
186        let swapped = if zero(du[1]) && zero(dv[0]) {
187            false
188        } else if zero(du[0]) && zero(dv[1]) {
189            true
190        } else {
191            return Ok(None);
192        };
193        let (su, sv) = (
194            du.prod_scalar(Scalar::ONE.div(u1.sub(u0))?),
195            dv.prod_scalar(Scalar::ONE.div(v1.sub(v0))?),
196        );
197        let m = if swapped {
198            [[Scalar::ZERO, sv[0]], [su[1], Scalar::ZERO]]
199        } else {
200            [[su[0], Scalar::ZERO], [Scalar::ZERO, sv[1]]]
201        };
202        // The same way round: the normal of this surface along the map is
203        // its own times the map's determinant.
204        let det = m[0][0].mul(m[1][1]).sub(m[0][1].mul(m[1][0]));
205        if !det.definitely_greater(Scalar::ZERO) {
206            return Ok(None);
207        }
208        let degrees = if swapped {
209            (self.degree_v, self.degree_u)
210        } else {
211            (self.degree_u, self.degree_v)
212        };
213        if degrees != (other.degree_u, other.degree_v) {
214            return Ok(None);
215        }
216        let shift = Vector2::from_array(std::array::from_fn(|k| {
217            p00[k].sub(m[k][0].mul(u0)).sub(m[k][1].mul(v0))
218        }));
219        let map = UvMap { m, shift };
220        // The cells: between `other`'s knots, and the carried surface's
221        // taken back into `other`'s parameters.
222        let cells = |k: usize| -> GeopResult<Vec<Scalar>> {
223            let own = if k == 0 {
224                other.breakpoints_u()
225            } else {
226                other.breakpoints_v()
227            };
228            let image = map.image(k);
229            let theirs = if image == 0 {
230                carried.breakpoints_u()
231            } else {
232                carried.breakpoints_v()
233            };
234            let (lo, hi) = (own[0], *own.last().expect("a domain"));
235            let mut breaks = own;
236            for b in theirs {
237                let back = b.sub(map.shift[image]).div(map.m[image][k])?;
238                if back.definitely_greater(lo) && back.definitely_less(hi) {
239                    breaks.push(back);
240                }
241            }
242            breaks.sort_by(|a, b| a.to_f64().total_cmp(&b.to_f64()));
243            Ok(breaks)
244        };
245        let samples = |breaks: &[Scalar], degree: usize| -> GeopResult<Vec<Scalar>> {
246            let n = 2 * degree + 1;
247            let mut at = Vec::new();
248            for pair in breaks.windows(2) {
249                for i in 0..=n {
250                    let f = Scalar::from_ratio(i as i64, n as i64)?;
251                    at.push(Scalar::interpolate(pair[0], pair[1], f));
252                }
253            }
254            Ok(at)
255        };
256        let (us, vs) = (cells(0)?, cells(1)?);
257        for &u in &samples(&us, other.degree_u)? {
258            for &v in &samples(&vs, other.degree_v)? {
259                let mapped = map.apply(Vector2::from_array([u, v]));
260                let theirs = other.evaluate(u, v)?;
261                let ours = carried.evaluate(mapped[0], mapped[1])?;
262                if !theirs.could_be_equal(&ours) {
263                    return Ok(None);
264                }
265            }
266        }
267        Ok(Some((carried, map)))
268    }
269}
270
271#[cfg(test)]
272mod tests {
273    use geop_core_math::{
274        scalars::Scalar,
275        vector::{Vector2, Vector3},
276    };
277
278    use crate::nurb_surface::NurbSurface3D;
279
280    fn p(x: f64, y: f64, z: f64) -> Vector3 {
281        Vector3::from_array([x, y, z].map(Scalar::from_f64))
282    }
283
284    fn square(x0: f64, x1: f64, z: f64) -> NurbSurface3D {
285        NurbSurface3D::bilinear(p(x0, 0., z), p(x1, 0., z), p(x1, 1., z), p(x0, 1., z)).unwrap()
286    }
287
288    /// Two squares side by side in one plane, parametrized alike: the
289    /// second lies on the first carried on, one further along `u`.
290    #[test]
291    fn a_square_beside_another_lies_on_it() {
292        let (_, map) = square(0., 1., 1.)
293            .carrying(&square(1., 3., 1.))
294            .unwrap()
295            .expect("one plane");
296        let corner = map.apply(Vector2::from_array([Scalar::ONE, Scalar::ONE]));
297        assert!(
298            corner[0].could_be_equal(Scalar::from_f64(3.0)),
299            "{corner:?}"
300        );
301    }
302
303    /// Parametrized the other way round in both directions, as a mirrored
304    /// copy may be, it still lies on it; facing the other way, it does not.
305    #[test]
306    fn a_square_turned_round_lies_on_it_facing_alike() {
307        let turned =
308            NurbSurface3D::bilinear(p(3., 1., 1.), p(1., 1., 1.), p(1., 0., 1.), p(3., 0., 1.))
309                .unwrap();
310        assert!(square(0., 1., 1.).carrying(&turned).unwrap().is_some());
311        let facing_away = turned.reverse_u();
312        assert!(square(0., 1., 1.).carrying(&facing_away).unwrap().is_none());
313    }
314
315    /// A square a hair above the plane of another is another surface.
316    #[test]
317    fn a_square_above_another_does_not() {
318        assert!(
319            square(0., 1., 1.)
320                .carrying(&square(1., 3., 1.001))
321                .unwrap()
322                .is_none()
323        );
324    }
325}
geop_core_geometry/nurb_surface/carry.rs

geop_core_geometry/nurb_surface/
carry.rs
