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