geop_core_geometry/nurb_surface/
revolve.rs

1use geop_core_math::{
2    geop_error::{GeopError, GeopResult},
3    scalars::Scalar,
4    vector::{Vector3, Vector4},
5};
6
7use super::NurbSurface;
8use crate::{nurb_curve::NurbCurve3D, shape::Axis};
9
10impl NurbSurface<4> {
11    /// The surface `profile` sweeps turning about `axis` from the angle
12    /// `from` to the angle `to`, in radians, right-handed about
13    /// `axis.direction`: the profile as it is at angle 0, turned.
14    ///
15    /// `u` is the turn, on `[0, 1]`: rational quadratic pieces of equal
16    /// angle, each at most a quarter turn, the way [`crate::shape::Arc`]
17    /// builds an arc. `v` is the profile's own parameter. So `Su x Sv`
18    /// points away from the axis where the profile runs along
19    /// `axis.direction`.
20    ///
21    /// A profile control point on the axis stays put at every angle: its
22    /// row collapses to one point, a pole of the parametrization.
23    ///
24    /// The angles of the joints between pieces are free choices, taken in
25    /// `f64` as [`crate::shape::Arc::point_at`] does: any angle serves
26    /// equally well, and each joint lies honestly on the circle it is built
27    /// from. Fails unless `to` is more than `from`, by at most a full turn.
28    pub fn revolve(profile: &NurbCurve3D, axis: &Axis, from: f64, to: f64) -> GeopResult<Self> {
29        let sweep = to - from;
30        if !(sweep > 0.0 && sweep <= std::f64::consts::TAU) {
31            return Err(GeopError::new(format!(
32                "NurbSurface::revolve: cannot turn from {from} to {to}: it must be more than none and at most a full turn"
33            )));
34        }
35        let pieces = ((sweep / std::f64::consts::FRAC_PI_2).ceil() as usize).max(1);
36        let step = sweep / pieces as f64;
37        let angles: Vec<f64> = (0..=pieces).map(|k| from + step * k as f64).collect();
38        // Each piece's middle control point is where the tangents at its
39        // ends meet: the ends' radii summed, over `1 + cos(step)`, weighted
40        // `cos(step / 2)`.
41        let one_plus_cos = Scalar::ONE.add(Scalar::from_f64(step).cos());
42        let middle_weight = Scalar::from_f64(step / 2.0).cos();
43
44        let num_v = profile.control_points.len();
45        let mut columns: Vec<Vec<Vector4>> = Vec::with_capacity(num_v);
46        for cp in &profile.control_points {
47            let w = cp[3];
48            let q = cp.dehomogenized()?;
49            let center = axis.project(&q);
50            let x = q.sub(&center);
51            let y = axis.direction.prod_cross(&x);
52            // The radius at `angle`. At angle 0 the profile point itself,
53            // taken as it is: turning it by nothing would only widen it.
54            let radius_at = |angle: f64| {
55                if angle == 0.0 {
56                    x
57                } else {
58                    let a = Scalar::from_f64(angle);
59                    x.prod_scalar(a.cos()).add(&y.prod_scalar(a.sin()))
60                }
61            };
62            let homogeneous = |p: Vector3, weight: Scalar| {
63                Vector4::from_array([p[0].mul(weight), p[1].mul(weight), p[2].mul(weight), weight])
64            };
65            let mut column = Vec::with_capacity(2 * pieces + 1);
66            let mut previous = radius_at(angles[0]);
67            column.push(if angles[0] == 0.0 {
68                *cp
69            } else {
70                homogeneous(center.add(&previous), w)
71            });
72            for &angle in &angles[1..] {
73                let next = radius_at(angle);
74                let middle = center.add(
75                    &previous
76                        .add(&next)
77                        .prod_scalar(Scalar::ONE.div(one_plus_cos)?),
78                );
79                column.push(homogeneous(middle, w.mul(middle_weight)));
80                column.push(homogeneous(center.add(&next), w));
81                previous = next;
82            }
83            columns.push(column);
84        }
85
86        let mut control_points = Vec::with_capacity((2 * pieces + 1) * num_v);
87        for i in 0..2 * pieces + 1 {
88            for column in &columns {
89                control_points.push(column[i]);
90            }
91        }
92        let mut knots_u = vec![Scalar::ZERO; 3];
93        for k in 1..pieces {
94            let knot = Scalar::from_ratio(k as i64, pieces as i64)?;
95            knots_u.extend([knot, knot]);
96        }
97        knots_u.extend([Scalar::ONE; 3]);
98        NurbSurface::try_new(
99            2,
100            profile.degree,
101            control_points,
102            knots_u,
103            profile.knot_vector.clone(),
104        )
105    }
106}
107
108#[cfg(test)]
109mod tests {
110    use crate::test_fixtures::line;
111    use geop_core_math::{scalars::Scalar, vector::Vector3};
112
113    use super::super::NurbSurface3D;
114    use crate::shape::Axis;
115
116    fn v3(x: f64, y: f64, z: f64) -> Vector3 {
117        Vector3::from_array([
118            Scalar::from_f64(x),
119            Scalar::from_f64(y),
120            Scalar::from_f64(z),
121        ])
122    }
123
124    /// A line beside the axis turns into a cylinder: every point at its
125    /// radius, at the height it had, the normal pointing away from the axis.
126    #[test]
127    fn a_line_turns_into_a_cylinder() {
128        let axis = Axis::try_new(v3(1.0, 2.0, 0.0), v3(0.0, 0.0, 1.0)).unwrap();
129        let profile = line([3.0, 2.0, -1.0], [3.0, 2.0, 4.0]);
130        for (from, to) in [(0.0, std::f64::consts::TAU), (0.5, 2.0), (-1.0, 3.0)] {
131            let surface = NurbSurface3D::revolve(&profile, &axis, from, to).unwrap();
132            for (u, v) in [(0.0, 0.0), (0.3, 0.5), (0.77, 1.0), (1.0, 0.25)] {
133                let (u, v) = (Scalar::from_f64(u), Scalar::from_f64(v));
134                let p = surface.evaluate(u, v).unwrap();
135                let r = p.sub(&v3(1.0, 2.0, p[2].to_f64())).norm();
136                assert!(r.could_be_equal(Scalar::from_f64(2.0)), "radius {r:?}");
137                let normal = surface.normal(u, v).unwrap();
138                let outward = p.sub(&axis.project(&p)).normalize().unwrap();
139                assert!(normal.prod_dot(&outward).definitely_greater(Scalar::ZERO));
140            }
141            let start = surface.evaluate(Scalar::ZERO, Scalar::ZERO).unwrap();
142            let (cos, sin) = (Scalar::from_f64(from).cos(), Scalar::from_f64(from).sin());
143            let expected = Vector3::from_array([
144                Scalar::ONE.add(Scalar::TWO.mul(cos)),
145                Scalar::TWO.add(Scalar::TWO.mul(sin)),
146                Scalar::from_f64(-1.0),
147            ]);
148            assert!(
149                start.could_be_equal(&expected),
150                "{from}: {start:?} vs {expected:?}"
151            );
152        }
153    }
154
155    /// A profile touching the axis collapses there: a pole.
156    #[test]
157    fn a_point_on_the_axis_is_a_pole() {
158        let axis = Axis::try_new(v3(0.0, 0.0, 0.0), v3(0.0, 0.0, 1.0)).unwrap();
159        let profile = line([0.0, 0.0, 2.0], [1.0, 0.0, 0.0]);
160        let surface = NurbSurface3D::revolve(&profile, &axis, 0.0, 3.0).unwrap();
161        for u in [0.0, 0.4, 1.0] {
162            let apex = surface.evaluate(Scalar::from_f64(u), Scalar::ZERO).unwrap();
163            assert!(apex.could_be_equal(&v3(0.0, 0.0, 2.0)));
164        }
165        assert!(NurbSurface3D::revolve(&profile, &axis, 1.0, 1.0).is_err());
166        assert!(NurbSurface3D::revolve(&profile, &axis, 0.0, 7.0).is_err());
167    }
168}
geop_core_geometry/nurb_surface/revolve.rs

geop_core_geometry/nurb_surface/
revolve.rs
