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