1use 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
53pub trait ShapeSample {
55 type Output;
57
58 fn sample_interior<R: RngExt + ?Sized>(&self, rng: &mut R) -> Self::Output;
72
73 fn sample_boundary<R: RngExt + ?Sized>(&self, rng: &mut R) -> Self::Output;
88
89 fn interior_dist(self) -> impl Distribution<Self::Output>
106 where
107 Self: Sized,
108 {
109 InteriorOf(self)
110 }
111
112 fn boundary_dist(self) -> impl Distribution<Self::Output>
129 where
130 Self: Sized,
131 {
132 BoundaryOf(self)
133 }
134}
135
136#[derive(Clone, Copy)]
137pub struct InteriorOf<T: ShapeSample>(pub T);
140
141#[derive(Clone, Copy)]
142pub 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 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 let theta = FRAC_PI_2 + rng.random_range(-self.half_angle()..self.half_angle());
192 Vec2::from_angle(theta) * self.radius()
193 } else {
194 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#[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 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 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
337fn 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 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
362fn 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 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 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 coords.sort_by(|x, y| x.partial_cmp(y).unwrap());
432
433 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 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 let dist = WeightedIndex::new(areas).unwrap();
459
460 let idx = dist.sample(rng);
462 triangles[idx].sample_interior(rng)
463 } else {
464 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 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 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 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 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 let capsule_vol = cylinder_vol + 4.0 / 3.0 * PI * self.radius * self.radius * self.radius;
562 if capsule_vol > 0.0 {
563 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 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 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 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 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}