Skip to main content

oxigaf_flame/
sampler.rs

1//! Area-weighted random point sampling on a triangle mesh surface.
2
3use nalgebra as na;
4use rand::{Rng, RngExt};
5
6use crate::mesh::Mesh;
7
8/// A point sampled on a mesh surface, with its face binding information.
9#[derive(Debug, Clone)]
10pub struct SurfacePoint {
11    /// World-space position.
12    pub position: na::Point3<f32>,
13    /// Interpolated unit normal.
14    pub normal: na::Vector3<f32>,
15    /// Index of the face this point lies on.
16    pub face_index: u32,
17    /// Barycentric coordinates `[u, v, w]` with `u + v + w ≈ 1`.
18    pub barycentric: [f32; 3],
19}
20
21/// Sample `count` points uniformly (area-weighted) on the surface of `mesh`.
22///
23/// Uses the standard method:
24/// 1. Build a CDF over triangle areas.
25/// 2. For each sample, pick a random triangle and generate random barycentric
26///    coordinates.
27#[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    // 1. Compute face areas and build CDF
34    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    // Ensure last entry is exactly 1.0
48    if let Some(last) = cdf.last_mut() {
49        *last = 1.0;
50    }
51
52    // 2. Sample points
53    let mut points = Vec::with_capacity(count);
54
55    for _ in 0..count {
56        // Pick a face via the CDF
57        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        // Random barycentric coordinates (uniform on triangle)
70        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    /// Build a minimal unit triangle mesh for testing.
96    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            // Point should be in the triangle: x >= 0, y >= 0, x + y <= 1
116            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            // z should be 0 (planar triangle on XY plane)
124            assert!(p.position.z.abs() < 1e-6);
125            // Barycentric coords should sum to ~1
126            let sum: f32 = p.barycentric.iter().sum();
127            assert!((sum - 1.0).abs() < 1e-5, "bary sum = {sum}");
128        }
129    }
130}