1use nalgebra as na;
4use rand::{Rng, RngExt};
5
6use crate::mesh::Mesh;
7
8#[derive(Debug, Clone)]
10pub struct SurfacePoint {
11 pub position: na::Point3<f32>,
13 pub normal: na::Vector3<f32>,
15 pub face_index: u32,
17 pub barycentric: [f32; 3],
19}
20
21#[must_use]
28pub fn sample_mesh_surface(mesh: &Mesh, count: usize, rng: &mut impl Rng) -> Vec<SurfacePoint> {
29 if mesh.faces.is_empty() || count == 0 {
30 return Vec::new();
31 }
32
33 let areas: Vec<f32> = mesh.faces.iter().map(|f| mesh.face_area(f)).collect();
35 let total_area: f32 = areas.iter().sum();
36
37 if total_area < 1e-12 {
38 return Vec::new();
39 }
40
41 let mut cdf = Vec::with_capacity(areas.len());
42 let mut cumulative = 0.0f32;
43 for &a in &areas {
44 cumulative += a / total_area;
45 cdf.push(cumulative);
46 }
47 if let Some(last) = cdf.last_mut() {
49 *last = 1.0;
50 }
51
52 let mut points = Vec::with_capacity(count);
54
55 for _ in 0..count {
56 let r: f32 = rng.random();
58 let face_idx = cdf.partition_point(|&x| x < r).min(mesh.faces.len() - 1);
59 let face = &mesh.faces[face_idx];
60
61 let v0 = &mesh.vertices[face[0] as usize];
62 let v1 = &mesh.vertices[face[1] as usize];
63 let v2 = &mesh.vertices[face[2] as usize];
64
65 let n0 = &mesh.normals[face[0] as usize];
66 let n1 = &mesh.normals[face[1] as usize];
67 let n2 = &mesh.normals[face[2] as usize];
68
69 let r1: f32 = rng.random::<f32>().sqrt();
71 let r2: f32 = rng.random();
72 let u = 1.0 - r1;
73 let v = r2 * r1;
74 let w = 1.0 - u - v;
75
76 let position = na::Point3::from(v0.coords * u + v1.coords * v + v2.coords * w);
77 let normal = (n0 * u + n1 * v + n2 * w).normalize();
78
79 points.push(SurfacePoint {
80 position,
81 normal,
82 face_index: face_idx as u32,
83 barycentric: [u, v, w],
84 });
85 }
86
87 points
88}
89
90#[cfg(test)]
91mod tests {
92 use super::*;
93 use rand::SeedableRng;
94
95 fn unit_triangle() -> Mesh {
97 Mesh::new(
98 vec![
99 na::Point3::new(0.0, 0.0, 0.0),
100 na::Point3::new(1.0, 0.0, 0.0),
101 na::Point3::new(0.0, 1.0, 0.0),
102 ],
103 vec![[0, 1, 2]],
104 )
105 }
106
107 #[test]
108 fn sampled_points_are_on_surface() {
109 let mesh = unit_triangle();
110 let mut rng = rand::rngs::StdRng::seed_from_u64(42);
111 let points = sample_mesh_surface(&mesh, 100, &mut rng);
112
113 assert_eq!(points.len(), 100);
114 for p in &points {
115 assert!(p.position.x >= -1e-6, "x = {}", p.position.x);
117 assert!(p.position.y >= -1e-6, "y = {}", p.position.y);
118 assert!(
119 p.position.x + p.position.y <= 1.0 + 1e-6,
120 "x+y = {}",
121 p.position.x + p.position.y
122 );
123 assert!(p.position.z.abs() < 1e-6);
125 let sum: f32 = p.barycentric.iter().sum();
127 assert!((sum - 1.0).abs() < 1e-5, "bary sum = {sum}");
128 }
129 }
130}