Skip to main content

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