1use core::f32::consts::{FRAC_PI_2, PI, TAU};
42
43use crate::{ops, primitives::*, NormedVectorSpace, ScalarField, Vec2, Vec3};
44use rand::{
45 distr::{
46 uniform::SampleUniform,
47 weighted::{Weight, WeightedIndex},
48 Distribution,
49 },
50 Rng,
51};
52
53pub trait ShapeSample {
55 type Output;
57
58 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output;
71
72 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output;
86
87 fn interior_dist(self) -> impl Distribution<Self::Output>
103 where
104 Self: Sized,
105 {
106 InteriorOf(self)
107 }
108
109 fn boundary_dist(self) -> impl Distribution<Self::Output>
125 where
126 Self: Sized,
127 {
128 BoundaryOf(self)
129 }
130}
131
132#[derive(Clone, Copy)]
133pub struct InteriorOf<T: ShapeSample>(pub T);
136
137#[derive(Clone, Copy)]
138pub struct BoundaryOf<T: ShapeSample>(pub T);
141
142impl<T: ShapeSample> Distribution<<T as ShapeSample>::Output> for InteriorOf<T> {
143 fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> <T as ShapeSample>::Output {
144 self.0.sample_interior(rng)
145 }
146}
147
148impl<T: ShapeSample> Distribution<<T as ShapeSample>::Output> for BoundaryOf<T> {
149 fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> <T as ShapeSample>::Output {
150 self.0.sample_boundary(rng)
151 }
152}
153
154impl ShapeSample for Circle {
155 type Output = Vec2;
156
157 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
158 let theta = rng.random_range(0.0..TAU);
160 let r_squared = rng.random_range(0.0..=(self.radius * self.radius));
161 let r = ops::sqrt(r_squared);
162 let (sin, cos) = ops::sin_cos(theta);
163 Vec2::new(r * cos, r * sin)
164 }
165
166 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
167 let theta = rng.random_range(0.0..TAU);
168 let (sin, cos) = ops::sin_cos(theta);
169 Vec2::new(self.radius * cos, self.radius * sin)
170 }
171}
172
173impl ShapeSample for CircularSector {
174 type Output = Vec2;
175
176 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
177 let theta = rng.random_range(-self.half_angle()..=self.half_angle());
178 let r_squared = rng.random_range(0.0..=(self.radius() * self.radius()));
179 let r = ops::sqrt(r_squared);
180 let (sin, cos) = ops::sin_cos(theta);
181 Vec2::new(r * sin, r * cos)
182 }
183
184 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
185 if rng.random_range(0.0..=1.0) <= self.arc_length() / self.perimeter() {
186 let theta = FRAC_PI_2 + rng.random_range(-self.half_angle()..self.half_angle());
188 Vec2::from_angle(theta) * self.radius()
189 } else {
190 let dir = self.radius() * Vec2::from_angle(FRAC_PI_2 + self.half_angle());
192 let r: f32 = rng.random_range(-1.0..1.0);
193 (-r).clamp(0.0, 1.0) * dir + r.clamp(0.0, 1.0) * dir * Vec2::new(-1.0, 1.0)
194 }
195 }
196}
197
198#[inline]
200fn sample_unit_sphere_boundary<R: Rng + ?Sized>(rng: &mut R) -> Vec3 {
201 let z = rng.random_range(-1f32..=1f32);
202 let (a_sin, a_cos) = ops::sin_cos(rng.random_range(-PI..=PI));
203 let c = ops::sqrt(1f32 - z * z);
204 let x = a_sin * c;
205 let y = a_cos * c;
206
207 Vec3::new(x, y, z)
208}
209
210impl ShapeSample for Sphere {
211 type Output = Vec3;
212
213 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
214 let r_cubed = rng.random_range(0.0..=(self.radius * self.radius * self.radius));
215 let r = ops::cbrt(r_cubed);
216
217 r * sample_unit_sphere_boundary(rng)
218 }
219
220 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
221 self.radius * sample_unit_sphere_boundary(rng)
222 }
223}
224
225impl ShapeSample for Annulus {
226 type Output = Vec2;
227
228 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
229 let inner_radius = self.inner_circle.radius;
230 let outer_radius = self.outer_circle.radius;
231
232 let r_squared =
234 rng.random_range((inner_radius * inner_radius)..(outer_radius * outer_radius));
235 let r = ops::sqrt(r_squared);
236 let theta = rng.random_range(0.0..TAU);
237 let (sin, cos) = ops::sin_cos(theta);
238
239 Vec2::new(r * cos, r * sin)
240 }
241
242 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
243 let total_perimeter = self.inner_circle.perimeter() + self.outer_circle.perimeter();
244 let inner_prob = (self.inner_circle.perimeter() / total_perimeter) as f64;
245
246 let inner = rng.random_bool(inner_prob);
248 if inner {
249 self.inner_circle.sample_boundary(rng)
250 } else {
251 self.outer_circle.sample_boundary(rng)
252 }
253 }
254}
255
256impl ShapeSample for Rhombus {
257 type Output = Vec2;
258
259 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
260 let x: f32 = rng.random_range(0.0..=1.0);
261 let y: f32 = rng.random_range(0.0..=1.0);
262
263 let unit_p = Vec2::NEG_X + x * Vec2::ONE + Vec2::new(y, -y);
264 unit_p * self.half_diagonals
265 }
266
267 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
268 let x: f32 = rng.random_range(-1.0..=1.0);
269 let y_sign = if rng.random() { -1.0 } else { 1.0 };
270
271 let y = (1.0 - ops::abs(x)) * y_sign;
272 Vec2::new(x, y) * self.half_diagonals
273 }
274}
275
276impl ShapeSample for Rectangle {
277 type Output = Vec2;
278
279 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
280 let x = rng.random_range(-self.half_size.x..=self.half_size.x);
281 let y = rng.random_range(-self.half_size.y..=self.half_size.y);
282 Vec2::new(x, y)
283 }
284
285 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
286 let primary_side = rng.random_range(-1.0..1.0);
287 let other_side = if rng.random() { -1.0 } else { 1.0 };
288
289 if self.half_size.x + self.half_size.y > 0.0 {
290 if rng.random_bool((self.half_size.x / (self.half_size.x + self.half_size.y)) as f64) {
291 Vec2::new(primary_side, other_side) * self.half_size
292 } else {
293 Vec2::new(other_side, primary_side) * self.half_size
294 }
295 } else {
296 Vec2::ZERO
297 }
298 }
299}
300
301impl ShapeSample for Cuboid {
302 type Output = Vec3;
303
304 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
305 let x = rng.random_range(-self.half_size.x..=self.half_size.x);
306 let y = rng.random_range(-self.half_size.y..=self.half_size.y);
307 let z = rng.random_range(-self.half_size.z..=self.half_size.z);
308 Vec3::new(x, y, z)
309 }
310
311 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
312 let primary_side1 = rng.random_range(-1.0..1.0);
313 let primary_side2 = rng.random_range(-1.0..1.0);
314 let other_side = if rng.random() { -1.0 } else { 1.0 };
315
316 if let Ok(dist) = WeightedIndex::new([
317 self.half_size.y * self.half_size.z,
318 self.half_size.x * self.half_size.z,
319 self.half_size.x * self.half_size.y,
320 ]) {
321 match dist.sample(rng) {
322 0 => Vec3::new(other_side, primary_side1, primary_side2) * self.half_size,
323 1 => Vec3::new(primary_side1, other_side, primary_side2) * self.half_size,
324 2 => Vec3::new(primary_side1, primary_side2, other_side) * self.half_size,
325 _ => unreachable!(),
326 }
327 } else {
328 Vec3::ZERO
329 }
330 }
331}
332
333fn sample_triangle_interior<P, R>(vertices: [P; 3], rng: &mut R) -> P
335where
336 P: NormedVectorSpace,
337 P::Scalar: SampleUniform + PartialOrd,
338 R: Rng + ?Sized,
339{
340 let [a, b, c] = vertices;
341 let ab = b - a;
342 let ac = c - a;
343
344 let u = rng.random_range(P::Scalar::ZERO..=P::Scalar::ONE);
347 let v = rng.random_range(P::Scalar::ZERO..=P::Scalar::ONE);
348
349 if u + v > P::Scalar::ONE {
350 let u1 = P::Scalar::ONE - v;
351 let v1 = P::Scalar::ONE - u;
352 a + (ab * u1 + ac * v1)
353 } else {
354 a + (ab * u + ac * v)
355 }
356}
357
358fn sample_triangle_boundary<P, R>(vertices: [P; 3], rng: &mut R) -> P
360where
361 P: NormedVectorSpace,
362 P::Scalar: Weight + SampleUniform + PartialOrd + for<'a> ::core::ops::AddAssign<&'a P::Scalar>,
363 R: Rng + ?Sized,
364{
365 let [a, b, c] = vertices;
366 let ab = b - a;
367 let ac = c - a;
368 let bc = c - b;
369
370 let t = rng.random_range(<P::Scalar as ScalarField>::ZERO..=P::Scalar::ONE);
371
372 if let Ok(dist) = WeightedIndex::new([ab.norm(), ac.norm(), bc.norm()]) {
373 match dist.sample(rng) {
374 0 => a.lerp(b, t),
375 1 => a.lerp(c, t),
376 2 => b.lerp(c, t),
377 _ => unreachable!(),
378 }
379 } else {
380 a
383 }
384}
385
386impl ShapeSample for Triangle2d {
387 type Output = Vec2;
388
389 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
390 sample_triangle_interior(self.vertices, rng)
391 }
392
393 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
394 sample_triangle_boundary(self.vertices, rng)
395 }
396}
397
398impl ShapeSample for Triangle3d {
399 type Output = Vec3;
400
401 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
402 sample_triangle_interior(self.vertices, rng)
403 }
404
405 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
406 sample_triangle_boundary(self.vertices, rng)
407 }
408}
409
410impl ShapeSample for Tetrahedron {
411 type Output = Vec3;
412
413 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
414 let [v0, v1, v2, v3] = self.vertices;
415
416 let mut coords: [f32; 3] = [
418 rng.random_range(0.0..1.0),
419 rng.random_range(0.0..1.0),
420 rng.random_range(0.0..1.0),
421 ];
422
423 coords.sort_by(|x, y| x.partial_cmp(y).unwrap());
428
429 let (a, b, c, d) = (
438 coords[0],
439 coords[1] - coords[0],
440 coords[2] - coords[1],
441 1. - coords[2],
442 );
443
444 v0 * a + v1 * b + v2 * c + v3 * d
446 }
447
448 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
449 let triangles = self.faces();
450 let areas = triangles.iter().map(Measured2d::area);
451
452 if areas.clone().sum::<f32>() > 0.0 {
453 let dist = WeightedIndex::new(areas).unwrap();
455
456 let idx = dist.sample(rng);
458 triangles[idx].sample_interior(rng)
459 } else {
460 self.vertices[0]
463 }
464 }
465}
466
467impl ShapeSample for Cylinder {
468 type Output = Vec3;
469
470 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
471 let Vec2 { x, y: z } = self.base().sample_interior(rng);
472 let y = rng.random_range(-self.half_height..=self.half_height);
473 Vec3::new(x, y, z)
474 }
475
476 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
477 if self.radius + 2.0 * self.half_height > 0.0 {
480 if rng.random_bool((self.radius / (self.radius + 2.0 * self.half_height)) as f64) {
481 let Vec2 { x, y: z } = self.base().sample_interior(rng);
482 if rng.random() {
483 Vec3::new(x, self.half_height, z)
484 } else {
485 Vec3::new(x, -self.half_height, z)
486 }
487 } else {
488 let Vec2 { x, y: z } = self.base().sample_boundary(rng);
489 let y = rng.random_range(-self.half_height..=self.half_height);
490 Vec3::new(x, y, z)
491 }
492 } else {
493 Vec3::ZERO
494 }
495 }
496}
497
498impl ShapeSample for Capsule2d {
499 type Output = Vec2;
500
501 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
502 let rectangle_area = self.half_length * self.radius * 4.0;
503 let capsule_area = rectangle_area + PI * self.radius * self.radius;
504 if capsule_area > 0.0 {
505 if rng.random_bool((rectangle_area / capsule_area) as f64) {
507 self.to_inner_rectangle().sample_interior(rng)
508 } else {
509 let circle = Circle::new(self.radius);
510 let point = circle.sample_interior(rng);
511 if point.y > 0.0 {
513 point + Vec2::Y * self.half_length
514 } else {
515 point - Vec2::Y * self.half_length
516 }
517 }
518 } else {
519 Vec2::ZERO
520 }
521 }
522
523 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
524 let rectangle_surface = 4.0 * self.half_length;
525 let capsule_surface = rectangle_surface + TAU * self.radius;
526 if capsule_surface > 0.0 {
527 if rng.random_bool((rectangle_surface / capsule_surface) as f64) {
528 let side_distance =
529 rng.random_range((-2.0 * self.half_length)..=(2.0 * self.half_length));
530 if side_distance < 0.0 {
531 Vec2::new(self.radius, side_distance + self.half_length)
532 } else {
533 Vec2::new(-self.radius, side_distance - self.half_length)
534 }
535 } else {
536 let circle = Circle::new(self.radius);
537 let point = circle.sample_boundary(rng);
538 if point.y > 0.0 {
540 point + Vec2::Y * self.half_length
541 } else {
542 point - Vec2::Y * self.half_length
543 }
544 }
545 } else {
546 Vec2::ZERO
547 }
548 }
549}
550
551impl ShapeSample for Capsule3d {
552 type Output = Vec3;
553
554 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
555 let cylinder_vol = PI * self.radius * self.radius * 2.0 * self.half_length;
556 let capsule_vol = cylinder_vol + 4.0 / 3.0 * PI * self.radius * self.radius * self.radius;
558 if capsule_vol > 0.0 {
559 if rng.random_bool((cylinder_vol / capsule_vol) as f64) {
561 self.to_cylinder().sample_interior(rng)
562 } else {
563 let sphere = Sphere::new(self.radius);
564 let point = sphere.sample_interior(rng);
565 if point.y > 0.0 {
567 point + Vec3::Y * self.half_length
568 } else {
569 point - Vec3::Y * self.half_length
570 }
571 }
572 } else {
573 Vec3::ZERO
574 }
575 }
576
577 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
578 let cylinder_surface = TAU * self.radius * 2.0 * self.half_length;
579 let capsule_surface = cylinder_surface + 4.0 * PI * self.radius * self.radius;
580 if capsule_surface > 0.0 {
581 if rng.random_bool((cylinder_surface / capsule_surface) as f64) {
582 let Vec2 { x, y: z } = Circle::new(self.radius).sample_boundary(rng);
583 let y = rng.random_range(-self.half_length..=self.half_length);
584 Vec3::new(x, y, z)
585 } else {
586 let sphere = Sphere::new(self.radius);
587 let point = sphere.sample_boundary(rng);
588 if point.y > 0.0 {
590 point + Vec3::Y * self.half_length
591 } else {
592 point - Vec3::Y * self.half_length
593 }
594 }
595 } else {
596 Vec3::ZERO
597 }
598 }
599}
600
601impl<P: Primitive2d + Measured2d + ShapeSample<Output = Vec2>> ShapeSample for Extrusion<P> {
602 type Output = Vec3;
603
604 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
605 let base_point = self.base_shape.sample_interior(rng);
606 let depth = rng.random_range(-self.half_depth..self.half_depth);
607 base_point.extend(depth)
608 }
609
610 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
611 let base_area = self.base_shape.area();
612 let total_area = self.area();
613
614 let random = rng.random_range(0.0..total_area);
615 match random {
616 x if x < base_area => self.base_shape.sample_interior(rng).extend(self.half_depth),
617 x if x < 2. * base_area => self
618 .base_shape
619 .sample_interior(rng)
620 .extend(-self.half_depth),
621 _ => self
622 .base_shape
623 .sample_boundary(rng)
624 .extend(rng.random_range(-self.half_depth..self.half_depth)),
625 }
626 }
627}
628
629#[cfg(test)]
630mod tests {
631 use super::*;
632 use rand::SeedableRng;
633 use rand_chacha::ChaCha8Rng;
634
635 #[test]
636 fn circle_interior_sampling() {
637 let mut rng = ChaCha8Rng::from_seed(Default::default());
638 let circle = Circle::new(8.0);
639
640 let boxes = [
641 (-3.0, 3.0),
642 (1.0, 2.0),
643 (-1.0, -2.0),
644 (3.0, -2.0),
645 (1.0, -6.0),
646 (-3.0, -7.0),
647 (-7.0, -3.0),
648 (-6.0, 1.0),
649 ];
650 let mut box_hits = [0; 8];
651
652 for _ in 0..5000 {
654 let point = circle.sample_interior(&mut rng);
655
656 for (i, box_) in boxes.iter().enumerate() {
657 if (point.x > box_.0 && point.x < box_.0 + 4.0)
658 && (point.y > box_.1 && point.y < box_.1 + 4.0)
659 {
660 box_hits[i] += 1;
661 }
662 }
663 }
664
665 assert_eq!(
666 box_hits,
667 [396, 377, 415, 404, 366, 408, 408, 430],
668 "samples will occur across all array items at statistically equal chance"
669 );
670 }
671
672 #[test]
673 fn circle_boundary_sampling() {
674 let mut rng = ChaCha8Rng::from_seed(Default::default());
675 let circle = Circle::new(1.0);
676
677 let mut wedge_hits = [0; 8];
678
679 for _ in 0..5000 {
681 let point = circle.sample_boundary(&mut rng);
682
683 let angle = ops::atan(point.y / point.x) + PI / 2.0;
684 let wedge = ops::floor(angle * 8.0 / PI) as usize;
685 wedge_hits[wedge] += 1;
686 }
687
688 assert_eq!(
689 wedge_hits,
690 [636, 608, 639, 603, 614, 650, 640, 610],
691 "samples will occur across all array items at statistically equal chance"
692 );
693 }
694}