Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Uniform point sampling methods for some primitive shapes. #12484

Merged
merged 18 commits into from
Mar 17, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
d02b7db
Add an optional (default) dependency for rand
13ros27 Mar 14, 2024
5745659
Add rand to dev-dependencies so that it will have thread_rng
13ros27 Mar 14, 2024
b5516f3
Added shape_sampling.rs with support for Circle, Sphere, Rectangle, C…
13ros27 Mar 14, 2024
efef007
Added shape_sampling to bevy_math lib.rs
13ros27 Mar 14, 2024
26b1e5e
Added ShapeSample for Capsule3d
13ros27 Mar 14, 2024
2b4b744
Merge branch 'main' into shape_sampling
13ros27 Mar 14, 2024
118a6fc
Added a ?Sized to all the rng arguments
13ros27 Mar 14, 2024
2f9b4ca
Changed from the weird `gen_range(0..=1) as f32 * 2.0 - 1.0` to `if g…
13ros27 Mar 15, 2024
ff44767
Change `gen_ratio(1, 2)` to `gen()`
13ros27 Mar 15, 2024
ac91a7a
Add the `alloc` feature to rand so that we get access to distributions
13ros27 Mar 15, 2024
5c9537b
Fixed rectangle and cuboid sampling
13ros27 Mar 15, 2024
0ac8d06
Appease taplo (this seems considerably worse formatted to me?)
13ros27 Mar 15, 2024
b314ca8
Added documentation specifying the direction ambiguous shapes point in
13ros27 Mar 15, 2024
4679fbe
Rename `sample_volume` to `sample_interior` and `sample_surface` to `…
13ros27 Mar 15, 2024
2f8a5a5
Added a test for how uniform Circle.sample_interior and Circle.sample…
13ros27 Mar 15, 2024
7b19264
Switch to `ChaCha8Rng` from `StdRng` for stability
13ros27 Mar 15, 2024
54890a0
Fixed the implementation for Cuboid::sample_boundary again
13ros27 Mar 15, 2024
f74d0a4
Added checks to make it return VecN::ZERO rather than panicking with …
13ros27 Mar 15, 2024
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 9 additions & 0 deletions crates/bevy_math/Cargo.toml
Original file line number Diff line number Diff line change
Expand Up @@ -14,11 +14,18 @@ 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"
# Supply rngs for examples and tests
rand = "0.8"
13ros27 marked this conversation as resolved.
Show resolved Hide resolved
rand_chacha = "0.3"

[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"]
Expand All @@ -31,6 +38,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
Expand Down
7 changes: 7 additions & 0 deletions crates/bevy_math/src/lib.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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::{
Expand Down
345 changes: 345 additions & 0 deletions crates/bevy_math/src/shape_sampling.rs
Original file line number Diff line number Diff line change
@@ -0,0 +1,345 @@
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 self.half_size.x + self.half_size.y > 0.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
}
} else {
Vec2::ZERO
}
13ros27 marked this conversation as resolved.
Show resolved Hide resolved
}
}

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 {
13ros27 marked this conversation as resolved.
Show resolved Hide resolved
let primary_side1 = rng.gen_range(-1.0..1.0);
let primary_side2 = rng.gen_range(-1.0..1.0);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

primary_side2... sounds almost like a secondary_side. :)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, my naming there is a little questionable, its really the sides that are fully varied in comparison to the one which is just -1 or 1

let other_side = if rng.gen() { -1.0 } else { 1.0 };

if let Ok(dist) = WeightedIndex::new([
self.half_size.y * self.half_size.z,
self.half_size.x * self.half_size.z,
self.half_size.x * self.half_size.y,
]) {
match dist.sample(rng) {
0 => Vec3::new(other_side, primary_side1, primary_side2) * self.half_size,
1 => Vec3::new(primary_side1, other_side, primary_side2) * self.half_size,
2 => Vec3::new(primary_side1, primary_side2, other_side) * self.half_size,
_ => unreachable!(),
}
} else {
Vec3::ZERO
}
13ros27 marked this conversation as resolved.
Show resolved Hide resolved
}
}

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 self.radius + 2.0 * self.half_height > 0.0 {
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)
}
} else {
Vec3::ZERO
}
}
}

13ros27 marked this conversation as resolved.
Show resolved Hide resolved
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;
if capsule_area > 0.0 {
// 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
}
}
} else {
Vec2::ZERO
}
}

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 capsule_surface > 0.0 {
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
}
}
} else {
Vec2::ZERO
}
}
}

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;
if capsule_vol > 0.0 {
// 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
}
}
} else {
Vec3::ZERO
}
}

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 capsule_surface > 0.0 {
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
}
}
} else {
Vec3::ZERO
}
}
}
13ros27 marked this conversation as resolved.
Show resolved Hide resolved

#[cfg(test)]
mod tests {
use super::*;
use rand::SeedableRng;
use rand_chacha::ChaCha8Rng;

#[test]
fn circle_interior_sampling() {
let mut rng = ChaCha8Rng::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,
[396, 377, 415, 404, 366, 408, 408, 430],
"samples will occur across all array items at statistically equal chance"
);
}

#[test]
fn circle_boundary_sampling() {
let mut rng = ChaCha8Rng::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,
[636, 608, 639, 603, 614, 650, 640, 610],
"samples will occur across all array items at statistically equal chance"
);
}
}