bevy_math/sampling/
shape_sampling.rs

1//! The [`ShapeSample`] trait, allowing random sampling from geometric shapes.
2//!
3//! At the most basic level, this allows sampling random points from the interior and boundary of
4//! geometric primitives. For example:
5//! ```
6//! # use bevy_math::primitives::*;
7//! # use bevy_math::ShapeSample;
8//! # use rand::SeedableRng;
9//! # use rand::rngs::StdRng;
10//! // Get some `Rng`:
11//! let rng = &mut StdRng::from_os_rng();
12//! // Make a circle of radius 2:
13//! let circle = Circle::new(2.0);
14//! // Get a point inside this circle uniformly at random:
15//! let interior_pt = circle.sample_interior(rng);
16//! // Get a point on the circle's boundary uniformly at random:
17//! let boundary_pt = circle.sample_boundary(rng);
18//! ```
19//!
20//! For repeated sampling, `ShapeSample` also includes methods for accessing a [`Distribution`]:
21//! ```
22//! # use bevy_math::primitives::*;
23//! # use bevy_math::{Vec2, ShapeSample};
24//! # use rand::SeedableRng;
25//! # use rand::rngs::StdRng;
26//! # use rand::distr::Distribution;
27//! # let rng1 = StdRng::from_os_rng();
28//! # let rng2 = StdRng::from_os_rng();
29//! // Use a rectangle this time:
30//! let rectangle = Rectangle::new(1.0, 2.0);
31//! // Get an iterator that spits out random interior points:
32//! let interior_iter = rectangle.interior_dist().sample_iter(rng1);
33//! // Collect random interior points from the iterator:
34//! let interior_pts: Vec<Vec2> = interior_iter.take(1000).collect();
35//! // Similarly, get an iterator over many random boundary points and collect them:
36//! let boundary_pts: Vec<Vec2> = rectangle.boundary_dist().sample_iter(rng2).take(1000).collect();
37//! ```
38//!
39//! In any case, the [`Rng`] used as the source of randomness must be provided explicitly.
40
41use core::f32::consts::{FRAC_PI_2, PI, TAU};
42
43use crate::{ops, primitives::*, NormedVectorSpace, ScalarField, Vec2, Vec3};
44use rand::{
45    distr::{
46        uniform::SampleUniform,
47        weighted::{Weight, WeightedIndex},
48        Distribution,
49    },
50    Rng,
51};
52
53/// Exposes methods to uniformly sample a variety of primitive shapes.
54pub trait ShapeSample {
55    /// The type of vector returned by the sample methods, [`Vec2`] for 2D shapes and [`Vec3`] for 3D shapes.
56    type Output;
57
58    /// Uniformly sample a point from inside the area/volume of this shape, centered on 0.
59    ///
60    /// Shapes like [`Cylinder`], [`Capsule2d`] and [`Capsule3d`] are oriented along the y-axis.
61    ///
62    /// # Example
63    /// ```
64    /// # use bevy_math::prelude::*;
65    /// let square = Rectangle::new(2.0, 2.0);
66    ///
67    /// // Returns a Vec2 with both x and y between -1 and 1.
68    /// println!("{}", square.sample_interior(&mut rand::rng()));
69    /// ```
70    fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output;
71
72    /// Uniformly sample a point from the surface of this shape, centered on 0.
73    ///
74    /// Shapes like [`Cylinder`], [`Capsule2d`] and [`Capsule3d`] are oriented along the y-axis.
75    ///
76    /// # Example
77    /// ```
78    /// # use bevy_math::prelude::*;
79    /// let square = Rectangle::new(2.0, 2.0);
80    ///
81    /// // Returns a Vec2 where one of the coordinates is at ±1,
82    /// //  and the other is somewhere between -1 and 1.
83    /// println!("{}", square.sample_boundary(&mut rand::rng()));
84    /// ```
85    fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output;
86
87    /// Extract a [`Distribution`] whose samples are points of this shape's interior, taken uniformly.
88    ///
89    /// # Example
90    ///
91    /// ```
92    /// # use bevy_math::prelude::*;
93    /// # use rand::distr::Distribution;
94    /// let square = Rectangle::new(2.0, 2.0);
95    /// let rng = rand::rng();
96    ///
97    /// // Iterate over points randomly drawn from `square`'s interior:
98    /// for random_val in square.interior_dist().sample_iter(rng).take(5) {
99    ///     println!("{}", random_val);
100    /// }
101    /// ```
102    fn interior_dist(self) -> impl Distribution<Self::Output>
103    where
104        Self: Sized,
105    {
106        InteriorOf(self)
107    }
108
109    /// Extract a [`Distribution`] whose samples are points of this shape's boundary, taken uniformly.
110    ///
111    /// # Example
112    ///
113    /// ```
114    /// # use bevy_math::prelude::*;
115    /// # use rand::distr::Distribution;
116    /// let square = Rectangle::new(2.0, 2.0);
117    /// let rng = rand::rng();
118    ///
119    /// // Iterate over points randomly drawn from `square`'s boundary:
120    /// for random_val in square.boundary_dist().sample_iter(rng).take(5) {
121    ///     println!("{}", random_val);
122    /// }
123    /// ```
124    fn boundary_dist(self) -> impl Distribution<Self::Output>
125    where
126        Self: Sized,
127    {
128        BoundaryOf(self)
129    }
130}
131
132#[derive(Clone, Copy)]
133/// A wrapper struct that allows interior sampling from a [`ShapeSample`] type directly as
134/// a [`Distribution`].
135pub struct InteriorOf<T: ShapeSample>(pub T);
136
137#[derive(Clone, Copy)]
138/// A wrapper struct that allows boundary sampling from a [`ShapeSample`] type directly as
139/// a [`Distribution`].
140pub struct BoundaryOf<T: ShapeSample>(pub T);
141
142impl<T: ShapeSample> Distribution<<T as ShapeSample>::Output> for InteriorOf<T> {
143    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> <T as ShapeSample>::Output {
144        self.0.sample_interior(rng)
145    }
146}
147
148impl<T: ShapeSample> Distribution<<T as ShapeSample>::Output> for BoundaryOf<T> {
149    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> <T as ShapeSample>::Output {
150        self.0.sample_boundary(rng)
151    }
152}
153
154impl ShapeSample for Circle {
155    type Output = Vec2;
156
157    fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
158        // https://mathworld.wolfram.com/DiskPointPicking.html
159        let theta = rng.random_range(0.0..TAU);
160        let r_squared = rng.random_range(0.0..=(self.radius * self.radius));
161        let r = ops::sqrt(r_squared);
162        let (sin, cos) = ops::sin_cos(theta);
163        Vec2::new(r * cos, r * sin)
164    }
165
166    fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
167        let theta = rng.random_range(0.0..TAU);
168        let (sin, cos) = ops::sin_cos(theta);
169        Vec2::new(self.radius * cos, self.radius * sin)
170    }
171}
172
173impl ShapeSample for CircularSector {
174    type Output = Vec2;
175
176    fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
177        let theta = rng.random_range(-self.half_angle()..=self.half_angle());
178        let r_squared = rng.random_range(0.0..=(self.radius() * self.radius()));
179        let r = ops::sqrt(r_squared);
180        let (sin, cos) = ops::sin_cos(theta);
181        Vec2::new(r * sin, r * cos)
182    }
183
184    fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
185        if rng.random_range(0.0..=1.0) <= self.arc_length() / self.perimeter() {
186            // Sample on the arc
187            let theta = FRAC_PI_2 + rng.random_range(-self.half_angle()..self.half_angle());
188            Vec2::from_angle(theta) * self.radius()
189        } else {
190            // Sample on the "inner" straight lines
191            let dir = self.radius() * Vec2::from_angle(FRAC_PI_2 + self.half_angle());
192            let r: f32 = rng.random_range(-1.0..1.0);
193            (-r).clamp(0.0, 1.0) * dir + r.clamp(0.0, 1.0) * dir * Vec2::new(-1.0, 1.0)
194        }
195    }
196}
197
198/// Boundary sampling for unit-spheres
199#[inline]
200fn sample_unit_sphere_boundary<R: Rng + ?Sized>(rng: &mut R) -> Vec3 {
201    let z = rng.random_range(-1f32..=1f32);
202    let (a_sin, a_cos) = ops::sin_cos(rng.random_range(-PI..=PI));
203    let c = ops::sqrt(1f32 - z * z);
204    let x = a_sin * c;
205    let y = a_cos * c;
206
207    Vec3::new(x, y, z)
208}
209
210impl ShapeSample for Sphere {
211    type Output = Vec3;
212
213    fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
214        let r_cubed = rng.random_range(0.0..=(self.radius * self.radius * self.radius));
215        let r = ops::cbrt(r_cubed);
216
217        r * sample_unit_sphere_boundary(rng)
218    }
219
220    fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
221        self.radius * sample_unit_sphere_boundary(rng)
222    }
223}
224
225impl ShapeSample for Annulus {
226    type Output = Vec2;
227
228    fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
229        let inner_radius = self.inner_circle.radius;
230        let outer_radius = self.outer_circle.radius;
231
232        // Like random sampling for a circle, radius is weighted by the square.
233        let r_squared =
234            rng.random_range((inner_radius * inner_radius)..(outer_radius * outer_radius));
235        let r = ops::sqrt(r_squared);
236        let theta = rng.random_range(0.0..TAU);
237        let (sin, cos) = ops::sin_cos(theta);
238
239        Vec2::new(r * cos, r * sin)
240    }
241
242    fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
243        let total_perimeter = self.inner_circle.perimeter() + self.outer_circle.perimeter();
244        let inner_prob = (self.inner_circle.perimeter() / total_perimeter) as f64;
245
246        // Sample from boundary circles, choosing which one by weighting by perimeter:
247        let inner = rng.random_bool(inner_prob);
248        if inner {
249            self.inner_circle.sample_boundary(rng)
250        } else {
251            self.outer_circle.sample_boundary(rng)
252        }
253    }
254}
255
256impl ShapeSample for Rhombus {
257    type Output = Vec2;
258
259    fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
260        let x: f32 = rng.random_range(0.0..=1.0);
261        let y: f32 = rng.random_range(0.0..=1.0);
262
263        let unit_p = Vec2::NEG_X + x * Vec2::ONE + Vec2::new(y, -y);
264        unit_p * self.half_diagonals
265    }
266
267    fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
268        let x: f32 = rng.random_range(-1.0..=1.0);
269        let y_sign = if rng.random() { -1.0 } else { 1.0 };
270
271        let y = (1.0 - ops::abs(x)) * y_sign;
272        Vec2::new(x, y) * self.half_diagonals
273    }
274}
275
276impl ShapeSample for Rectangle {
277    type Output = Vec2;
278
279    fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
280        let x = rng.random_range(-self.half_size.x..=self.half_size.x);
281        let y = rng.random_range(-self.half_size.y..=self.half_size.y);
282        Vec2::new(x, y)
283    }
284
285    fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
286        let primary_side = rng.random_range(-1.0..1.0);
287        let other_side = if rng.random() { -1.0 } else { 1.0 };
288
289        if self.half_size.x + self.half_size.y > 0.0 {
290            if rng.random_bool((self.half_size.x / (self.half_size.x + self.half_size.y)) as f64) {
291                Vec2::new(primary_side, other_side) * self.half_size
292            } else {
293                Vec2::new(other_side, primary_side) * self.half_size
294            }
295        } else {
296            Vec2::ZERO
297        }
298    }
299}
300
301impl ShapeSample for Cuboid {
302    type Output = Vec3;
303
304    fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
305        let x = rng.random_range(-self.half_size.x..=self.half_size.x);
306        let y = rng.random_range(-self.half_size.y..=self.half_size.y);
307        let z = rng.random_range(-self.half_size.z..=self.half_size.z);
308        Vec3::new(x, y, z)
309    }
310
311    fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
312        let primary_side1 = rng.random_range(-1.0..1.0);
313        let primary_side2 = rng.random_range(-1.0..1.0);
314        let other_side = if rng.random() { -1.0 } else { 1.0 };
315
316        if let Ok(dist) = WeightedIndex::new([
317            self.half_size.y * self.half_size.z,
318            self.half_size.x * self.half_size.z,
319            self.half_size.x * self.half_size.y,
320        ]) {
321            match dist.sample(rng) {
322                0 => Vec3::new(other_side, primary_side1, primary_side2) * self.half_size,
323                1 => Vec3::new(primary_side1, other_side, primary_side2) * self.half_size,
324                2 => Vec3::new(primary_side1, primary_side2, other_side) * self.half_size,
325                _ => unreachable!(),
326            }
327        } else {
328            Vec3::ZERO
329        }
330    }
331}
332
333/// Interior sampling for triangles which doesn't depend on the ambient dimension.
334fn sample_triangle_interior<P, R>(vertices: [P; 3], rng: &mut R) -> P
335where
336    P: NormedVectorSpace,
337    P::Scalar: SampleUniform + PartialOrd,
338    R: Rng + ?Sized,
339{
340    let [a, b, c] = vertices;
341    let ab = b - a;
342    let ac = c - a;
343
344    // Generate random points on a parallelepiped and reflect so that
345    // we can use the points that lie outside the triangle
346    let u = rng.random_range(P::Scalar::ZERO..=P::Scalar::ONE);
347    let v = rng.random_range(P::Scalar::ZERO..=P::Scalar::ONE);
348
349    if u + v > P::Scalar::ONE {
350        let u1 = P::Scalar::ONE - v;
351        let v1 = P::Scalar::ONE - u;
352        a + (ab * u1 + ac * v1)
353    } else {
354        a + (ab * u + ac * v)
355    }
356}
357
358/// Boundary sampling for triangles which doesn't depend on the ambient dimension.
359fn sample_triangle_boundary<P, R>(vertices: [P; 3], rng: &mut R) -> P
360where
361    P: NormedVectorSpace,
362    P::Scalar: Weight + SampleUniform + PartialOrd + for<'a> ::core::ops::AddAssign<&'a P::Scalar>,
363    R: Rng + ?Sized,
364{
365    let [a, b, c] = vertices;
366    let ab = b - a;
367    let ac = c - a;
368    let bc = c - b;
369
370    let t = rng.random_range(<P::Scalar as ScalarField>::ZERO..=P::Scalar::ONE);
371
372    if let Ok(dist) = WeightedIndex::new([ab.norm(), ac.norm(), bc.norm()]) {
373        match dist.sample(rng) {
374            0 => a.lerp(b, t),
375            1 => a.lerp(c, t),
376            2 => b.lerp(c, t),
377            _ => unreachable!(),
378        }
379    } else {
380        // This should only occur when the triangle is 0-dimensional degenerate
381        // so this is actually the correct result.
382        a
383    }
384}
385
386impl ShapeSample for Triangle2d {
387    type Output = Vec2;
388
389    fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
390        sample_triangle_interior(self.vertices, rng)
391    }
392
393    fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
394        sample_triangle_boundary(self.vertices, rng)
395    }
396}
397
398impl ShapeSample for Triangle3d {
399    type Output = Vec3;
400
401    fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
402        sample_triangle_interior(self.vertices, rng)
403    }
404
405    fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
406        sample_triangle_boundary(self.vertices, rng)
407    }
408}
409
410impl ShapeSample for Tetrahedron {
411    type Output = Vec3;
412
413    fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
414        let [v0, v1, v2, v3] = self.vertices;
415
416        // Generate a random point in a cube:
417        let mut coords: [f32; 3] = [
418            rng.random_range(0.0..1.0),
419            rng.random_range(0.0..1.0),
420            rng.random_range(0.0..1.0),
421        ];
422
423        // The cube is broken into six tetrahedra of the form 0 <= c_0 <= c_1 <= c_2 <= 1,
424        // where c_i are the three euclidean coordinates in some permutation. (Since 3! = 6,
425        // there are six of them). Sorting the coordinates folds these six tetrahedra into the
426        // tetrahedron 0 <= x <= y <= z <= 1 (i.e. a fundamental domain of the permutation action).
427        coords.sort_by(|x, y| x.partial_cmp(y).unwrap());
428
429        // Now, convert a point from the fundamental tetrahedron into barycentric coordinates by
430        // taking the four successive differences of coordinates; note that these telescope to sum
431        // to 1, and this transformation is linear, hence preserves the probability density, since
432        // the latter comes from the Lebesgue measure.
433        //
434        // (See https://en.wikipedia.org/wiki/Lebesgue_measure#Properties — specifically, that
435        // Lebesgue measure of a linearly transformed set is its original measure times the
436        // determinant.)
437        let (a, b, c, d) = (
438            coords[0],
439            coords[1] - coords[0],
440            coords[2] - coords[1],
441            1. - coords[2],
442        );
443
444        // This is also a linear mapping, so probability density is still preserved.
445        v0 * a + v1 * b + v2 * c + v3 * d
446    }
447
448    fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
449        let triangles = self.faces();
450        let areas = triangles.iter().map(Measured2d::area);
451
452        if areas.clone().sum::<f32>() > 0.0 {
453            // There is at least one triangle with nonzero area, so this unwrap succeeds.
454            let dist = WeightedIndex::new(areas).unwrap();
455
456            // Get a random index, then sample the interior of the associated triangle.
457            let idx = dist.sample(rng);
458            triangles[idx].sample_interior(rng)
459        } else {
460            // In this branch the tetrahedron has zero surface area; just return a point that's on
461            // the tetrahedron.
462            self.vertices[0]
463        }
464    }
465}
466
467impl ShapeSample for Cylinder {
468    type Output = Vec3;
469
470    fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
471        let Vec2 { x, y: z } = self.base().sample_interior(rng);
472        let y = rng.random_range(-self.half_height..=self.half_height);
473        Vec3::new(x, y, z)
474    }
475
476    fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
477        // This uses the area of the ends divided by the overall surface area (optimized)
478        // [2 (\pi r^2)]/[2 (\pi r^2) + 2 \pi r h] = r/(r + h)
479        if self.radius + 2.0 * self.half_height > 0.0 {
480            if rng.random_bool((self.radius / (self.radius + 2.0 * self.half_height)) as f64) {
481                let Vec2 { x, y: z } = self.base().sample_interior(rng);
482                if rng.random() {
483                    Vec3::new(x, self.half_height, z)
484                } else {
485                    Vec3::new(x, -self.half_height, z)
486                }
487            } else {
488                let Vec2 { x, y: z } = self.base().sample_boundary(rng);
489                let y = rng.random_range(-self.half_height..=self.half_height);
490                Vec3::new(x, y, z)
491            }
492        } else {
493            Vec3::ZERO
494        }
495    }
496}
497
498impl ShapeSample for Capsule2d {
499    type Output = Vec2;
500
501    fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
502        let rectangle_area = self.half_length * self.radius * 4.0;
503        let capsule_area = rectangle_area + PI * self.radius * self.radius;
504        if capsule_area > 0.0 {
505            // Check if the random point should be inside the rectangle
506            if rng.random_bool((rectangle_area / capsule_area) as f64) {
507                self.to_inner_rectangle().sample_interior(rng)
508            } else {
509                let circle = Circle::new(self.radius);
510                let point = circle.sample_interior(rng);
511                // Add half length if it is the top semi-circle, otherwise subtract half
512                if point.y > 0.0 {
513                    point + Vec2::Y * self.half_length
514                } else {
515                    point - Vec2::Y * self.half_length
516                }
517            }
518        } else {
519            Vec2::ZERO
520        }
521    }
522
523    fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
524        let rectangle_surface = 4.0 * self.half_length;
525        let capsule_surface = rectangle_surface + TAU * self.radius;
526        if capsule_surface > 0.0 {
527            if rng.random_bool((rectangle_surface / capsule_surface) as f64) {
528                let side_distance =
529                    rng.random_range((-2.0 * self.half_length)..=(2.0 * self.half_length));
530                if side_distance < 0.0 {
531                    Vec2::new(self.radius, side_distance + self.half_length)
532                } else {
533                    Vec2::new(-self.radius, side_distance - self.half_length)
534                }
535            } else {
536                let circle = Circle::new(self.radius);
537                let point = circle.sample_boundary(rng);
538                // Add half length if it is the top semi-circle, otherwise subtract half
539                if point.y > 0.0 {
540                    point + Vec2::Y * self.half_length
541                } else {
542                    point - Vec2::Y * self.half_length
543                }
544            }
545        } else {
546            Vec2::ZERO
547        }
548    }
549}
550
551impl ShapeSample for Capsule3d {
552    type Output = Vec3;
553
554    fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
555        let cylinder_vol = PI * self.radius * self.radius * 2.0 * self.half_length;
556        // Add 4/3 pi r^3
557        let capsule_vol = cylinder_vol + 4.0 / 3.0 * PI * self.radius * self.radius * self.radius;
558        if capsule_vol > 0.0 {
559            // Check if the random point should be inside the cylinder
560            if rng.random_bool((cylinder_vol / capsule_vol) as f64) {
561                self.to_cylinder().sample_interior(rng)
562            } else {
563                let sphere = Sphere::new(self.radius);
564                let point = sphere.sample_interior(rng);
565                // Add half length if it is the top semi-sphere, otherwise subtract half
566                if point.y > 0.0 {
567                    point + Vec3::Y * self.half_length
568                } else {
569                    point - Vec3::Y * self.half_length
570                }
571            }
572        } else {
573            Vec3::ZERO
574        }
575    }
576
577    fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
578        let cylinder_surface = TAU * self.radius * 2.0 * self.half_length;
579        let capsule_surface = cylinder_surface + 4.0 * PI * self.radius * self.radius;
580        if capsule_surface > 0.0 {
581            if rng.random_bool((cylinder_surface / capsule_surface) as f64) {
582                let Vec2 { x, y: z } = Circle::new(self.radius).sample_boundary(rng);
583                let y = rng.random_range(-self.half_length..=self.half_length);
584                Vec3::new(x, y, z)
585            } else {
586                let sphere = Sphere::new(self.radius);
587                let point = sphere.sample_boundary(rng);
588                // Add half length if it is the top semi-sphere, otherwise subtract half
589                if point.y > 0.0 {
590                    point + Vec3::Y * self.half_length
591                } else {
592                    point - Vec3::Y * self.half_length
593                }
594            }
595        } else {
596            Vec3::ZERO
597        }
598    }
599}
600
601impl<P: Primitive2d + Measured2d + ShapeSample<Output = Vec2>> ShapeSample for Extrusion<P> {
602    type Output = Vec3;
603
604    fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
605        let base_point = self.base_shape.sample_interior(rng);
606        let depth = rng.random_range(-self.half_depth..self.half_depth);
607        base_point.extend(depth)
608    }
609
610    fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
611        let base_area = self.base_shape.area();
612        let total_area = self.area();
613
614        let random = rng.random_range(0.0..total_area);
615        match random {
616            x if x < base_area => self.base_shape.sample_interior(rng).extend(self.half_depth),
617            x if x < 2. * base_area => self
618                .base_shape
619                .sample_interior(rng)
620                .extend(-self.half_depth),
621            _ => self
622                .base_shape
623                .sample_boundary(rng)
624                .extend(rng.random_range(-self.half_depth..self.half_depth)),
625        }
626    }
627}
628
629#[cfg(test)]
630mod tests {
631    use super::*;
632    use rand::SeedableRng;
633    use rand_chacha::ChaCha8Rng;
634
635    #[test]
636    fn circle_interior_sampling() {
637        let mut rng = ChaCha8Rng::from_seed(Default::default());
638        let circle = Circle::new(8.0);
639
640        let boxes = [
641            (-3.0, 3.0),
642            (1.0, 2.0),
643            (-1.0, -2.0),
644            (3.0, -2.0),
645            (1.0, -6.0),
646            (-3.0, -7.0),
647            (-7.0, -3.0),
648            (-6.0, 1.0),
649        ];
650        let mut box_hits = [0; 8];
651
652        // Checks which boxes (if any) the sampled points are in
653        for _ in 0..5000 {
654            let point = circle.sample_interior(&mut rng);
655
656            for (i, box_) in boxes.iter().enumerate() {
657                if (point.x > box_.0 && point.x < box_.0 + 4.0)
658                    && (point.y > box_.1 && point.y < box_.1 + 4.0)
659                {
660                    box_hits[i] += 1;
661                }
662            }
663        }
664
665        assert_eq!(
666            box_hits,
667            [396, 377, 415, 404, 366, 408, 408, 430],
668            "samples will occur across all array items at statistically equal chance"
669        );
670    }
671
672    #[test]
673    fn circle_boundary_sampling() {
674        let mut rng = ChaCha8Rng::from_seed(Default::default());
675        let circle = Circle::new(1.0);
676
677        let mut wedge_hits = [0; 8];
678
679        // Checks in which eighth of the circle each sampled point is in
680        for _ in 0..5000 {
681            let point = circle.sample_boundary(&mut rng);
682
683            let angle = ops::atan(point.y / point.x) + PI / 2.0;
684            let wedge = ops::floor(angle * 8.0 / PI) as usize;
685            wedge_hits[wedge] += 1;
686        }
687
688        assert_eq!(
689            wedge_hits,
690            [636, 608, 639, 603, 614, 650, 640, 610],
691            "samples will occur across all array items at statistically equal chance"
692        );
693    }
694}