1use std::f32::consts::{PI, TAU};
42
43use crate::{primitives::*, NormedVectorSpace, Vec2, Vec3};
44use rand::{
45 distributions::{Distribution, WeightedIndex},
46 Rng,
47};
48
49pub trait ShapeSample {
51 type Output;
53
54 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output;
67
68 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output;
82
83 fn interior_dist(self) -> impl Distribution<Self::Output>
99 where
100 Self: Sized,
101 {
102 InteriorOf(self)
103 }
104
105 fn boundary_dist(self) -> impl Distribution<Self::Output>
121 where
122 Self: Sized,
123 {
124 BoundaryOf(self)
125 }
126}
127
128#[derive(Clone, Copy)]
129pub struct InteriorOf<T: ShapeSample>(pub T);
132
133#[derive(Clone, Copy)]
134pub 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 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 Vec2::new(r * theta.cos(), r * theta.sin())
159 }
160
161 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
162 let theta = rng.gen_range(0.0..TAU);
163 Vec2::new(self.radius * theta.cos(), self.radius * theta.sin())
164 }
165}
166
167impl ShapeSample for Sphere {
168 type Output = Vec3;
169
170 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
171 let theta = rng.gen_range(0.0..TAU);
173 let phi = rng.gen_range(-1.0_f32..1.0).acos();
174 let r_cubed = rng.gen_range(0.0..=(self.radius * self.radius * self.radius));
175 let r = r_cubed.cbrt();
176 Vec3 {
177 x: r * phi.sin() * theta.cos(),
178 y: r * phi.sin() * theta.sin(),
179 z: r * phi.cos(),
180 }
181 }
182
183 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
184 let theta = rng.gen_range(0.0..TAU);
185 let phi = rng.gen_range(-1.0_f32..1.0).acos();
186 Vec3 {
187 x: self.radius * phi.sin() * theta.cos(),
188 y: self.radius * phi.sin() * theta.sin(),
189 z: self.radius * phi.cos(),
190 }
191 }
192}
193
194impl ShapeSample for Annulus {
195 type Output = Vec2;
196
197 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
198 let inner_radius = self.inner_circle.radius;
199 let outer_radius = self.outer_circle.radius;
200
201 let r_squared = rng.gen_range((inner_radius * inner_radius)..(outer_radius * outer_radius));
203 let r = r_squared.sqrt();
204 let theta = rng.gen_range(0.0..TAU);
205
206 Vec2::new(r * theta.cos(), r * theta.sin())
207 }
208
209 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
210 let total_perimeter = self.inner_circle.perimeter() + self.outer_circle.perimeter();
211 let inner_prob = (self.inner_circle.perimeter() / total_perimeter) as f64;
212
213 let inner = rng.gen_bool(inner_prob);
215 if inner {
216 self.inner_circle.sample_boundary(rng)
217 } else {
218 self.outer_circle.sample_boundary(rng)
219 }
220 }
221}
222
223impl ShapeSample for Rectangle {
224 type Output = Vec2;
225
226 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
227 let x = rng.gen_range(-self.half_size.x..=self.half_size.x);
228 let y = rng.gen_range(-self.half_size.y..=self.half_size.y);
229 Vec2::new(x, y)
230 }
231
232 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
233 let primary_side = rng.gen_range(-1.0..1.0);
234 let other_side = if rng.gen() { -1.0 } else { 1.0 };
235
236 if self.half_size.x + self.half_size.y > 0.0 {
237 if rng.gen_bool((self.half_size.x / (self.half_size.x + self.half_size.y)) as f64) {
238 Vec2::new(primary_side, other_side) * self.half_size
239 } else {
240 Vec2::new(other_side, primary_side) * self.half_size
241 }
242 } else {
243 Vec2::ZERO
244 }
245 }
246}
247
248impl ShapeSample for Cuboid {
249 type Output = Vec3;
250
251 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
252 let x = rng.gen_range(-self.half_size.x..=self.half_size.x);
253 let y = rng.gen_range(-self.half_size.y..=self.half_size.y);
254 let z = rng.gen_range(-self.half_size.z..=self.half_size.z);
255 Vec3::new(x, y, z)
256 }
257
258 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
259 let primary_side1 = rng.gen_range(-1.0..1.0);
260 let primary_side2 = rng.gen_range(-1.0..1.0);
261 let other_side = if rng.gen() { -1.0 } else { 1.0 };
262
263 if let Ok(dist) = WeightedIndex::new([
264 self.half_size.y * self.half_size.z,
265 self.half_size.x * self.half_size.z,
266 self.half_size.x * self.half_size.y,
267 ]) {
268 match dist.sample(rng) {
269 0 => Vec3::new(other_side, primary_side1, primary_side2) * self.half_size,
270 1 => Vec3::new(primary_side1, other_side, primary_side2) * self.half_size,
271 2 => Vec3::new(primary_side1, primary_side2, other_side) * self.half_size,
272 _ => unreachable!(),
273 }
274 } else {
275 Vec3::ZERO
276 }
277 }
278}
279
280fn sample_triangle_interior<P: NormedVectorSpace, R: Rng + ?Sized>(
282 vertices: [P; 3],
283 rng: &mut R,
284) -> P {
285 let [a, b, c] = vertices;
286 let ab = b - a;
287 let ac = c - a;
288
289 let u = rng.gen_range(0.0..=1.0);
292 let v = rng.gen_range(0.0..=1.0);
293
294 if u + v > 1. {
295 let u1 = 1. - v;
296 let v1 = 1. - u;
297 a + (ab * u1 + ac * v1)
298 } else {
299 a + (ab * u + ac * v)
300 }
301}
302
303fn sample_triangle_boundary<P: NormedVectorSpace, R: Rng + ?Sized>(
305 vertices: [P; 3],
306 rng: &mut R,
307) -> P {
308 let [a, b, c] = vertices;
309 let ab = b - a;
310 let ac = c - a;
311 let bc = c - b;
312
313 let t = rng.gen_range(0.0..=1.0);
314
315 if let Ok(dist) = WeightedIndex::new([ab.norm(), ac.norm(), bc.norm()]) {
316 match dist.sample(rng) {
317 0 => a.lerp(b, t),
318 1 => a.lerp(c, t),
319 2 => b.lerp(c, t),
320 _ => unreachable!(),
321 }
322 } else {
323 a
326 }
327}
328
329impl ShapeSample for Triangle2d {
330 type Output = Vec2;
331
332 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
333 sample_triangle_interior(self.vertices, rng)
334 }
335
336 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
337 sample_triangle_boundary(self.vertices, rng)
338 }
339}
340
341impl ShapeSample for Triangle3d {
342 type Output = Vec3;
343
344 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
345 sample_triangle_interior(self.vertices, rng)
346 }
347
348 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
349 sample_triangle_boundary(self.vertices, rng)
350 }
351}
352
353impl ShapeSample for Tetrahedron {
354 type Output = Vec3;
355
356 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
357 let [v0, v1, v2, v3] = self.vertices;
358
359 let mut coords: [f32; 3] = [
361 rng.gen_range(0.0..1.0),
362 rng.gen_range(0.0..1.0),
363 rng.gen_range(0.0..1.0),
364 ];
365
366 coords.sort_by(|x, y| x.partial_cmp(y).unwrap());
371
372 let (a, b, c, d) = (
381 coords[0],
382 coords[1] - coords[0],
383 coords[2] - coords[1],
384 1. - coords[2],
385 );
386
387 v0 * a + v1 * b + v2 * c + v3 * d
389 }
390
391 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
392 let triangles = self.faces();
393 let areas = triangles.iter().map(|t| t.area());
394
395 if areas.clone().sum::<f32>() > 0.0 {
396 let dist = WeightedIndex::new(areas).unwrap();
398
399 let idx = dist.sample(rng);
401 triangles[idx].sample_interior(rng)
402 } else {
403 self.vertices[0]
406 }
407 }
408}
409
410impl ShapeSample for Cylinder {
411 type Output = Vec3;
412
413 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
414 let Vec2 { x, y: z } = self.base().sample_interior(rng);
415 let y = rng.gen_range(-self.half_height..=self.half_height);
416 Vec3::new(x, y, z)
417 }
418
419 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
420 if self.radius + 2.0 * self.half_height > 0.0 {
423 if rng.gen_bool((self.radius / (self.radius + 2.0 * self.half_height)) as f64) {
424 let Vec2 { x, y: z } = self.base().sample_interior(rng);
425 if rng.gen() {
426 Vec3::new(x, self.half_height, z)
427 } else {
428 Vec3::new(x, -self.half_height, z)
429 }
430 } else {
431 let Vec2 { x, y: z } = self.base().sample_boundary(rng);
432 let y = rng.gen_range(-self.half_height..=self.half_height);
433 Vec3::new(x, y, z)
434 }
435 } else {
436 Vec3::ZERO
437 }
438 }
439}
440
441impl ShapeSample for Capsule2d {
442 type Output = Vec2;
443
444 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
445 let rectangle_area = self.half_length * self.radius * 4.0;
446 let capsule_area = rectangle_area + PI * self.radius * self.radius;
447 if capsule_area > 0.0 {
448 if rng.gen_bool((rectangle_area / capsule_area) as f64) {
450 let rectangle = Rectangle::new(self.radius, self.half_length * 2.0);
451 rectangle.sample_interior(rng)
452 } else {
453 let circle = Circle::new(self.radius);
454 let point = circle.sample_interior(rng);
455 if point.y > 0.0 {
457 point + Vec2::Y * self.half_length
458 } else {
459 point - Vec2::Y * self.half_length
460 }
461 }
462 } else {
463 Vec2::ZERO
464 }
465 }
466
467 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec2 {
468 let rectangle_surface = 4.0 * self.half_length;
469 let capsule_surface = rectangle_surface + TAU * self.radius;
470 if capsule_surface > 0.0 {
471 if rng.gen_bool((rectangle_surface / capsule_surface) as f64) {
472 let side_distance =
473 rng.gen_range((-2.0 * self.half_length)..=(2.0 * self.half_length));
474 if side_distance < 0.0 {
475 Vec2::new(self.radius, side_distance + self.half_length)
476 } else {
477 Vec2::new(-self.radius, side_distance - self.half_length)
478 }
479 } else {
480 let circle = Circle::new(self.radius);
481 let point = circle.sample_boundary(rng);
482 if point.y > 0.0 {
484 point + Vec2::Y * self.half_length
485 } else {
486 point - Vec2::Y * self.half_length
487 }
488 }
489 } else {
490 Vec2::ZERO
491 }
492 }
493}
494
495impl ShapeSample for Capsule3d {
496 type Output = Vec3;
497
498 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
499 let cylinder_vol = PI * self.radius * self.radius * 2.0 * self.half_length;
500 let capsule_vol = cylinder_vol + 4.0 / 3.0 * PI * self.radius * self.radius * self.radius;
502 if capsule_vol > 0.0 {
503 if rng.gen_bool((cylinder_vol / capsule_vol) as f64) {
505 self.to_cylinder().sample_interior(rng)
506 } else {
507 let sphere = Sphere::new(self.radius);
508 let point = sphere.sample_interior(rng);
509 if point.y > 0.0 {
511 point + Vec3::Y * self.half_length
512 } else {
513 point - Vec3::Y * self.half_length
514 }
515 }
516 } else {
517 Vec3::ZERO
518 }
519 }
520
521 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Vec3 {
522 let cylinder_surface = TAU * self.radius * 2.0 * self.half_length;
523 let capsule_surface = cylinder_surface + 4.0 * PI * self.radius * self.radius;
524 if capsule_surface > 0.0 {
525 if rng.gen_bool((cylinder_surface / capsule_surface) as f64) {
526 let Vec2 { x, y: z } = Circle::new(self.radius).sample_boundary(rng);
527 let y = rng.gen_range(-self.half_length..=self.half_length);
528 Vec3::new(x, y, z)
529 } else {
530 let sphere = Sphere::new(self.radius);
531 let point = sphere.sample_boundary(rng);
532 if point.y > 0.0 {
534 point + Vec3::Y * self.half_length
535 } else {
536 point - Vec3::Y * self.half_length
537 }
538 }
539 } else {
540 Vec3::ZERO
541 }
542 }
543}
544
545impl<P: Primitive2d + Measured2d + ShapeSample<Output = Vec2>> ShapeSample for Extrusion<P> {
546 type Output = Vec3;
547
548 fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
549 let base_point = self.base_shape.sample_interior(rng);
550 let depth = rng.gen_range(-self.half_depth..self.half_depth);
551 base_point.extend(depth)
552 }
553
554 fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
555 let base_area = self.base_shape.area();
556 let total_area = self.area();
557
558 let random = rng.gen_range(0.0..total_area);
559 match random {
560 x if x < base_area => self.base_shape.sample_interior(rng).extend(self.half_depth),
561 x if x < 2. * base_area => self
562 .base_shape
563 .sample_interior(rng)
564 .extend(-self.half_depth),
565 _ => self
566 .base_shape
567 .sample_boundary(rng)
568 .extend(rng.gen_range(-self.half_depth..self.half_depth)),
569 }
570 }
571}
572
573#[cfg(test)]
574mod tests {
575 use super::*;
576 use rand::SeedableRng;
577 use rand_chacha::ChaCha8Rng;
578
579 #[test]
580 fn circle_interior_sampling() {
581 let mut rng = ChaCha8Rng::from_seed(Default::default());
582 let circle = Circle::new(8.0);
583
584 let boxes = [
585 (-3.0, 3.0),
586 (1.0, 2.0),
587 (-1.0, -2.0),
588 (3.0, -2.0),
589 (1.0, -6.0),
590 (-3.0, -7.0),
591 (-7.0, -3.0),
592 (-6.0, 1.0),
593 ];
594 let mut box_hits = [0; 8];
595
596 for _ in 0..5000 {
598 let point = circle.sample_interior(&mut rng);
599
600 for (i, box_) in boxes.iter().enumerate() {
601 if (point.x > box_.0 && point.x < box_.0 + 4.0)
602 && (point.y > box_.1 && point.y < box_.1 + 4.0)
603 {
604 box_hits[i] += 1;
605 }
606 }
607 }
608
609 assert_eq!(
610 box_hits,
611 [396, 377, 415, 404, 366, 408, 408, 430],
612 "samples will occur across all array items at statistically equal chance"
613 );
614 }
615
616 #[test]
617 fn circle_boundary_sampling() {
618 let mut rng = ChaCha8Rng::from_seed(Default::default());
619 let circle = Circle::new(1.0);
620
621 let mut wedge_hits = [0; 8];
622
623 for _ in 0..5000 {
625 let point = circle.sample_boundary(&mut rng);
626
627 let angle = f32::atan(point.y / point.x) + PI / 2.0;
628 let wedge = (angle * 8.0 / PI).floor() as usize;
629 wedge_hits[wedge] += 1;
630 }
631
632 assert_eq!(
633 wedge_hits,
634 [636, 608, 639, 603, 614, 650, 640, 610],
635 "samples will occur across all array items at statistically equal chance"
636 );
637 }
638}