Uniform point sampling methods for some primitive shapes.

crates/bevy_math/Cargo.toml

@@ -14,11 +14,17 @@ thiserror = "1.0"
 serde = { version = "1", features = ["derive"], optional = true }
 libm = { version = "0.2", optional = true }
 approx = { version = "0.5", optional = true }
+rand = { version = "0.8", features = [
+ "alloc",
+], default-features = false, optional = true }
 
 [dev-dependencies]
 approx = "0.5"
+# Enable the use of thread_rng in examples
+rand = "0.8"
 
 [features]
+default = ["rand"]
 serialize = ["dep:serde", "glam/serde"]
 # Enable approx for glam types to approximate floating point equality comparisons and assertions
 approx = ["dep:approx", "glam/approx"]
@@ -31,6 +37,8 @@ libm = ["dep:libm", "glam/libm"]
 glam_assert = ["glam/glam-assert"]
 # Enable assertions in debug builds to check the validity of parameters passed to glam
 debug_glam_assert = ["glam/debug-glam-assert"]
+# Enable the rand dependency for shape_sampling
+rand = ["dep:rand"]
 
 [lints]
 workspace = true

crates/bevy_math/src/lib.rs

@@ -14,16 +14,23 @@ pub mod primitives;
 mod ray;
 mod rects;
 mod rotation2d;
+#[cfg(feature = "rand")]
+mod shape_sampling;
 
 pub use affine3::*;
 pub use aspect_ratio::AspectRatio;
 pub use direction::*;
 pub use ray::{Ray2d, Ray3d};
 pub use rects::*;
 pub use rotation2d::Rotation2d;
+#[cfg(feature = "rand")]
+pub use shape_sampling::ShapeSample;
 
 /// The `bevy_math` prelude.
 pub mod prelude {
+ #[doc(hidden)]
+ #[cfg(feature = "rand")]
+ pub use crate::shape_sampling::ShapeSample;
  #[doc(hidden)]
  pub use crate::{
  cubic_splines::{

crates/bevy_math/src/shape_sampling.rs

@@ -0,0 +1,312 @@
+use std::f32::consts::{PI, TAU};
+
+use crate::{primitives::*, Vec2, Vec3};
+use rand::{
+ distributions::{Distribution, WeightedIndex},
+ Rng,
+};
+
+/// Exposes methods to uniformly sample a variety of primitive shapes.
+pub trait ShapeSample {
+ /// The type of vector returned by the sample methods, [`Vec2`] for 2D shapes and [`Vec3`] for 3D shapes.
+ type Output;
+
+ /// Uniformly sample a point from inside the area/volume of this shape, centered on 0.
+ ///
+ /// Shapes like [`Cylinder`], [`Capsule2d`] and [`Capsule3d`] are oriented along the y-axis.
+ ///
+ /// # Example
+ /// ```
+ /// # use bevy_math::prelude::*;
+ /// let square = Rectangle::new(2.0, 2.0);
+ ///
+ /// // Returns a Vec2 with both x and y between -1 and 1.
+ /// println!("{:?}", square.sample_interior(&mut rand::thread_rng()));
+ /// ```
+ fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output;
+
+ /// Uniformly sample a point from the surface of this shape, centered on 0.
+ ///
+ /// Shapes like [`Cylinder`], [`Capsule2d`] and [`Capsule3d`] are oriented along the y-axis.
+ ///
+ /// # Example
+ /// ```
+ /// # use bevy_math::prelude::*;
+ /// let square = Rectangle::new(2.0, 2.0);
+ ///
+ /// // Returns a Vec2 where one of the coordinates is at ±1,
+ /// // and the other is somewhere between -1 and 1.
+ /// println!("{:?}", square.sample_boundary(&mut rand::thread_rng()));
+ /// ```
+ fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output;
+}
+
+impl ShapeSample for Circle {
+ type Output = Vec2;
+
+ fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
+ // https://mathworld.wolfram.com/DiskPointPicking.html
+ let theta = rng.gen_range(0.0..TAU);
+ let r_squared = rng.gen_range(0.0..=(self.radius * self.radius));
+ let r = r_squared.sqrt();
+ Vec2::new(r * theta.cos(), r * theta.sin())
+ }
+
+ fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
+ let theta = rng.gen_range(0.0..TAU);
+ Vec2::new(self.radius * theta.cos(), self.radius * theta.sin())
+ }
+}
+
+impl ShapeSample for Sphere {
+ type Output = Vec3;
+
+ fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
+ // https://mathworld.wolfram.com/SpherePointPicking.html
+ let theta = rng.gen_range(0.0..TAU);
+ let phi = rng.gen_range(-1.0_f32..1.0).acos();
+ let r_cubed = rng.gen_range(0.0..=(self.radius * self.radius * self.radius));
+ let r = r_cubed.cbrt();
+ Vec3 {
+ x: r * phi.sin() * theta.cos(),
+ y: r * phi.sin() * theta.sin(),
+ z: r * phi.cos(),
+ }
+ }
+
+ fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
+ let theta = rng.gen_range(0.0..TAU);
+ let phi = rng.gen_range(-1.0_f32..1.0).acos();
+ Vec3 {
+ x: self.radius * phi.sin() * theta.cos(),
+ y: self.radius * phi.sin() * theta.sin(),
+ z: self.radius * phi.cos(),
+ }
+ }
+}
+
+impl ShapeSample for Rectangle {
+ type Output = Vec2;
+
+ fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
+ let x = rng.gen_range(-self.half_size.x..=self.half_size.x);
+ let y = rng.gen_range(-self.half_size.y..=self.half_size.y);
+ Vec2::new(x, y)
+ }
+
+ fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
+ let primary_side = rng.gen_range(-1.0..1.0);
+ let other_side = if rng.gen() { -1.0 } else { 1.0 };
+
+ if rng.gen_bool((self.half_size.x / (self.half_size.x + self.half_size.y)) as f64) {
+ Vec2::new(primary_side, other_side) * self.half_size
+ } else {
+ Vec2::new(other_side, primary_side) * self.half_size
+ }
+ }
+}
+
+impl ShapeSample for Cuboid {
+ type Output = Vec3;
+
+ fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
+ let x = rng.gen_range(-self.half_size.x..=self.half_size.x);
+ let y = rng.gen_range(-self.half_size.y..=self.half_size.y);
+ let z = rng.gen_range(-self.half_size.z..=self.half_size.z);
+ Vec3::new(x, y, z)
+ }
+
+ fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
+ let primary_side = rng.gen_range(-1.0..1.0);
+ let other_side1 = if rng.gen() { -1.0 } else { 1.0 };
+ let other_side2 = if rng.gen() { -1.0 } else { 1.0 };
+
+ let dist = WeightedIndex::new(self.half_size.to_array()).unwrap();
+ match dist.sample(rng) {
+ 0 => Vec3::new(primary_side, other_side1, other_side2) * self.half_size,
+ 1 => Vec3::new(other_side1, primary_side, other_side2) * self.half_size,
+ 2 => Vec3::new(other_side1, other_side2, primary_side) * self.half_size,
+ _ => unreachable!(),
+ }
+ }
+}
+
+impl ShapeSample for Cylinder {
+ type Output = Vec3;
+
+ fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
+ let Vec2 { x, y: z } = self.base().sample_interior(rng);
+ let y = rng.gen_range(-self.half_height..=self.half_height);
+ Vec3::new(x, y, z)
+ }
+
+ fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
+ // This uses the area of the ends divided by the overall surface area (optimised)
+ // [2 (\pi r^2)]/[2 (\pi r^2) + 2 \pi r h] = r/(r + h)
+ if rng.gen_bool((self.radius / (self.radius + 2.0 * self.half_height)) as f64) {
+ let Vec2 { x, y: z } = self.base().sample_interior(rng);
+ if rng.gen() {
+ Vec3::new(x, self.half_height, z)
+ } else {
+ Vec3::new(x, -self.half_height, z)
+ }
+ } else {
+ let Vec2 { x, y: z } = self.base().sample_boundary(rng);
+ let y = rng.gen_range(-self.half_height..=self.half_height);
+ Vec3::new(x, y, z)
+ }
+ }
+}
+
+impl ShapeSample for Capsule2d {
+ type Output = Vec2;
+
+ fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
+ let rectangle_area = self.half_length * self.radius * 4.0;
+ let capsule_area = rectangle_area + PI * self.radius * self.radius;
+ // Check if the random point should be inside the rectangle
+ if rng.gen_bool((rectangle_area / capsule_area) as f64) {
+ let rectangle = Rectangle::new(self.radius, self.half_length * 2.0);
+ rectangle.sample_interior(rng)
+ } else {
+ let circle = Circle::new(self.radius);
+ let point = circle.sample_interior(rng);
+ // Add half length if it is the top semi-circle, otherwise subtract half
+ if point.y > 0.0 {
+ point + Vec2::Y * self.half_length
+ } else {
+ point - Vec2::Y * self.half_length
+ }
+ }
+ }
+
+ fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
+ let rectangle_surface = 4.0 * self.half_length;
+ let capsule_surface = rectangle_surface + TAU * self.radius;
+ if rng.gen_bool((rectangle_surface / capsule_surface) as f64) {
+ let side_distance = rng.gen_range((-2.0 * self.half_length)..=(2.0 * self.half_length));
+ if side_distance < 0.0 {
+ Vec2::new(self.radius, side_distance + self.half_length)
+ } else {
+ Vec2::new(-self.radius, side_distance - self.half_length)
+ }
+ } else {
+ let circle = Circle::new(self.radius);
+ let point = circle.sample_boundary(rng);
+ // Add half length if it is the top semi-circle, otherwise subtract half
+ if point.y > 0.0 {
+ point + Vec2::Y * self.half_length
+ } else {
+ point - Vec2::Y * self.half_length
+ }
+ }
+ }
+}
+
+impl ShapeSample for Capsule3d {
+ type Output = Vec3;
+
+ fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
+ let cylinder_vol = PI * self.radius * self.radius * 2.0 * self.half_length;
+ // Add 4/3 pi r^3
+ let capsule_vol = cylinder_vol + 4.0 / 3.0 * PI * self.radius * self.radius * self.radius;
+ // Check if the random point should be inside the cylinder
+ if rng.gen_bool((cylinder_vol / capsule_vol) as f64) {
+ self.to_cylinder().sample_interior(rng)
+ } else {
+ let sphere = Sphere::new(self.radius);
+ let point = sphere.sample_interior(rng);
+ // Add half length if it is the top semi-sphere, otherwise subtract half
+ if point.y > 0.0 {
+ point + Vec3::Y * self.half_length
+ } else {
+ point - Vec3::Y * self.half_length
+ }
+ }
+ }
+
+ fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
+ let cylinder_surface = TAU * self.radius * 2.0 * self.half_length;
+ let capsule_surface = cylinder_surface + 4.0 * PI * self.radius * self.radius;
+ if rng.gen_bool((cylinder_surface / capsule_surface) as f64) {
+ let Vec2 { x, y: z } = Circle::new(self.radius).sample_boundary(rng);
+ let y = rng.gen_range(-self.half_length..=self.half_length);
+ Vec3::new(x, y, z)
+ } else {
+ let sphere = Sphere::new(self.radius);
+ let point = sphere.sample_boundary(rng);
+ // Add half length if it is the top semi-sphere, otherwise subtract half
+ if point.y > 0.0 {
+ point + Vec3::Y * self.half_length
+ } else {
+ point - Vec3::Y * self.half_length
+ }
+ }
+ }
+}
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+ use rand::{rngs::StdRng, SeedableRng};
+
+ #[test]
+ fn circle_interior_sampling() {
+ let mut rng = StdRng::from_seed(Default::default());
+ let circle = Circle::new(8.0);
+
+ let boxes = [
+ (-3.0, 3.0),
+ (1.0, 2.0),
+ (-1.0, -2.0),
+ (3.0, -2.0),
+ (1.0, -6.0),
+ (-3.0, -7.0),
+ (-7.0, -3.0),
+ (-6.0, 1.0),
+ ];
+ let mut box_hits = [0; 8];
+
+ // Checks which boxes (if any) the sampled points are in
+ for _ in 0..5000 {
+ let point = circle.sample_interior(&mut rng);
+
+ for (i, box_) in boxes.iter().enumerate() {
+ if (point.x > box_.0 && point.x < box_.0 + 4.0)
+ && (point.y > box_.1 && point.y < box_.1 + 4.0)
+ {
+ box_hits[i] += 1;
+ }
+ }
+ }
+
+ assert_eq!(
+ box_hits,
+ [400, 367, 365, 394, 445, 417, 405, 382],
+ "samples will occur across all array items at statistically equal chance"
+ );
+ }
+
+ #[test]
+ fn circle_boundary_sampling() {
+ let mut rng = StdRng::from_seed(Default::default());
+ let circle = Circle::new(1.0);
+
+ let mut wedge_hits = [0; 8];
+
+ // Checks in which eighth of the circle each sampled point is in
+ for _ in 0..5000 {
+ let point = circle.sample_boundary(&mut rng);
+
+ let angle = f32::atan(point.y / point.x) + PI / 2.0;
+ let wedge = (angle * 8.0 / PI).floor() as usize;
+ wedge_hits[wedge] += 1;
+ }
+
+ assert_eq!(
+ wedge_hits,
+ [677, 638, 649, 625, 594, 624, 600, 593],
+ "samples will occur across all array items at statistically equal chance"
+ );
+ }
+}