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