truck_geometry/specifieds/
sphere.rs

1use super::*;
2use std::f64::consts::PI;
3
4impl Sphere {
5    /// Creates a sphere
6    #[inline(always)]
7    pub const fn new(center: Point3, radius: f64) -> Sphere { Sphere { center, radius } }
8    /// Returns the center
9    #[inline(always)]
10    pub const fn center(&self) -> Point3 { self.center }
11    /// Returns the radius
12    #[inline(always)]
13    pub const fn radius(&self) -> f64 { self.radius }
14    /// Returns whether the point `pt` is on sphere
15    #[inline(always)]
16    pub fn include(&self, pt: Point3) -> bool { self.center.distance(pt).near(&self.radius) }
17}
18
19impl ParametricSurface for Sphere {
20    type Point = Point3;
21    type Vector = Vector3;
22    #[inline(always)]
23    fn subs(&self, u: f64, v: f64) -> Point3 { self.center() + self.radius * self.normal(u, v) }
24    #[inline(always)]
25    fn uder(&self, u: f64, v: f64) -> Vector3 {
26        self.radius
27            * Vector3::new(
28                f64::cos(u) * f64::cos(v),
29                f64::cos(u) * f64::sin(v),
30                -f64::sin(u),
31            )
32    }
33    #[inline(always)]
34    fn vder(&self, u: f64, v: f64) -> Vector3 {
35        self.radius * f64::sin(u) * Vector3::new(-f64::sin(v), f64::cos(v), 0.0)
36    }
37    #[inline(always)]
38    fn uuder(&self, u: f64, v: f64) -> Vector3 { -self.radius * self.normal(u, v) }
39    #[inline(always)]
40    fn uvder(&self, u: f64, v: f64) -> Vector3 {
41        self.radius * f64::cos(u) * Vector3::new(-f64::sin(v), f64::cos(v), 0.0)
42    }
43    #[inline(always)]
44    fn vvder(&self, u: f64, v: f64) -> Vector3 {
45        -self.radius * f64::sin(u) * Vector3::new(f64::cos(v), f64::sin(v), 0.0)
46    }
47    #[inline(always)]
48    fn parameter_range(&self) -> (ParameterRange, ParameterRange) {
49        (
50            (Bound::Included(0.0), Bound::Included(PI)),
51            (Bound::Included(0.0), Bound::Excluded(2.0 * PI)),
52        )
53    }
54    #[inline(always)]
55    fn v_period(&self) -> Option<f64> { Some(2.0 * PI) }
56}
57
58impl ParametricSurface3D for Sphere {
59    #[inline(always)]
60    fn normal(&self, u: f64, v: f64) -> Vector3 {
61        Vector3::new(
62            f64::sin(u) * f64::cos(v),
63            f64::sin(u) * f64::sin(v),
64            f64::cos(u),
65        )
66    }
67}
68
69#[test]
70fn sphere_derivation_test() {
71    let center = Point3::new(1.0, 2.0, 3.0);
72    let radius = 4.56;
73    let sphere = Sphere::new(center, radius);
74    const N: usize = 100;
75    for i in 0..N {
76        for j in 0..N {
77            let u = PI * i as f64 / N as f64;
78            let v = 2.0 * PI * j as f64 / N as f64;
79            let normal = sphere.normal(u, v);
80            assert!(normal.dot(sphere.uder(u, v)).so_small());
81            assert!(normal.dot(sphere.vder(u, v)).so_small());
82        }
83    }
84}
85
86impl BoundedSurface for Sphere {}
87
88impl IncludeCurve<BSplineCurve<Point3>> for Sphere {
89    #[inline(always)]
90    fn include(&self, curve: &BSplineCurve<Point3>) -> bool {
91        curve.is_const() && self.include(curve.front())
92    }
93}
94
95impl IncludeCurve<NurbsCurve<Vector4>> for Sphere {
96    fn include(&self, curve: &NurbsCurve<Vector4>) -> bool {
97        let (knots, _) = curve.knot_vec().to_single_multi();
98        let degree = curve.degree() * 2;
99        knots
100            .windows(2)
101            .flat_map(move |window| (1..degree).map(move |i| (window, i)))
102            .all(move |(window, i)| {
103                let t = i as f64 / degree as f64;
104                let t = window[0] * (1.0 - t) + window[1] * t;
105                self.include(curve.subs(t))
106            })
107    }
108}
109
110impl ParameterDivision2D for Sphere {
111    #[inline(always)]
112    fn parameter_division(
113        &self,
114        (urange, vrange): ((f64, f64), (f64, f64)),
115        tol: f64,
116    ) -> (Vec<f64>, Vec<f64>) {
117        nonpositive_tolerance!(tol);
118        assert!(
119            tol < self.radius,
120            "Tolerance is larger than the radius of sphere."
121        );
122        let delta = 2.0 * f64::acos(1.0 - tol / self.radius);
123        let u_div = 1 + ((urange.1 - urange.0) / delta).floor() as usize;
124        let v_div = 1 + ((vrange.1 - vrange.0) / delta).floor() as usize;
125        (
126            (0..=u_div)
127                .map(|i| urange.0 + (urange.1 - urange.0) * i as f64 / u_div as f64)
128                .collect(),
129            (0..=v_div)
130                .map(|j| vrange.0 + (vrange.1 - vrange.0) * j as f64 / v_div as f64)
131                .collect(),
132        )
133    }
134}
135
136impl SearchParameter<D2> for Sphere {
137    type Point = Point3;
138    #[inline(always)]
139    fn search_parameter<H: Into<SPHint2D>>(
140        &self,
141        point: Point3,
142        hint: H,
143        _: usize,
144    ) -> Option<(f64, f64)> {
145        let radius = point - self.center;
146        if (self.radius * self.radius).near(&radius.magnitude2()) {
147            let radius = radius.normalize();
148            let u = f64::acos(radius[2]);
149            let sinu = f64::sqrt(1.0 - radius[2] * radius[2]);
150            let cosv = f64::clamp(radius[0] / sinu, -1.0, 1.0);
151            let v = if sinu.so_small() {
152                match hint.into() {
153                    SPHint2D::Parameter(_, hint) => hint,
154                    _ => 0.0,
155                }
156            } else if radius[1] > 0.0 {
157                f64::acos(cosv)
158            } else {
159                2.0 * PI - f64::acos(cosv)
160            };
161            Some((u, v))
162        } else {
163            None
164        }
165    }
166}
167
168impl SearchNearestParameter<D2> for Sphere {
169    type Point = Point3;
170    #[inline(always)]
171    fn search_nearest_parameter<H: Into<SPHint2D>>(
172        &self,
173        point: Point3,
174        _: H,
175        _: usize,
176    ) -> Option<(f64, f64)> {
177        let radius = (point - self.center).normalize();
178        let u = f64::acos(radius[2]);
179        let sinu = f64::sqrt(1.0 - radius[2] * radius[2]);
180        let cosv = radius[0] / sinu;
181        let v = if radius[1] > 0.0 {
182            f64::acos(cosv)
183        } else {
184            2.0 * PI - f64::acos(cosv)
185        };
186        Some((u, v))
187    }
188}
189
190#[cfg(test)]
191fn exec_search_parameter_test() {
192    let center = Point3::new(
193        100.0 * rand::random::<f64>() - 50.0,
194        100.0 * rand::random::<f64>() - 50.0,
195        100.0 * rand::random::<f64>() - 50.0,
196    );
197    let radius = 100.0 * rand::random::<f64>();
198    let sphere = Sphere::new(center, radius);
199    let u = PI * rand::random::<f64>();
200    let v = 2.0 * PI * rand::random::<f64>();
201    let pt = sphere.subs(u, v);
202    let (u0, v0) = sphere.search_parameter(pt, None, 100).unwrap();
203    assert_near!(Vector2::new(u, v), Vector2::new(u0, v0));
204    let pt = pt
205        + Vector3::new(
206            (0.1 * rand::random::<f64>() + 0.01) * f64::signum(rand::random::<f64>() - 0.5),
207            (0.1 * rand::random::<f64>() + 0.01) * f64::signum(rand::random::<f64>() - 0.5),
208            (0.1 * rand::random::<f64>() + 0.01) * f64::signum(rand::random::<f64>() - 0.5),
209        );
210    assert!(sphere.search_parameter(pt, None, 100).is_none());
211    let (u, v) = sphere.search_nearest_parameter(pt, None, 100).unwrap();
212    assert_near!(
213        sphere.subs(u, v),
214        center + (pt - center).normalize() * radius
215    );
216}
217
218#[test]
219fn search_parameter_test() { (0..10).for_each(|_| exec_search_parameter_test()) }
truck_geometry/specifieds/sphere.rs

truck_geometry/specifieds/
sphere.rs
