runmat_geometry_ops/
triangulation.rs

1//! Triangulation helpers shared by MATLAB-facing builtins.
2
3use std::collections::HashMap;
4
5#[derive(Debug, Clone, PartialEq)]
6pub struct Delaunay2d {
7    pub triangles: Vec<[usize; 3]>,
8}
9
10#[derive(Debug, Clone, PartialEq, Eq)]
11pub enum TriangulationError {
12    NonFinitePoint,
13}
14
15impl std::fmt::Display for TriangulationError {
16    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
17        match self {
18            Self::NonFinitePoint => write!(f, "points must be finite"),
19        }
20    }
21}
22
23impl std::error::Error for TriangulationError {}
24
25pub fn delaunay_2d(points: &[[f64; 2]]) -> Result<Delaunay2d, TriangulationError> {
26    let mut unique_points = Vec::with_capacity(points.len());
27    let mut original_indices = Vec::with_capacity(points.len());
28    let mut seen = HashMap::with_capacity(points.len());
29    for (idx, point) in points.iter().enumerate() {
30        if !point[0].is_finite() || !point[1].is_finite() {
31            return Err(TriangulationError::NonFinitePoint);
32        }
33        let key = (stable_float_key(point[0]), stable_float_key(point[1]));
34        if seen.contains_key(&key) {
35            continue;
36        }
37        seen.insert(key, unique_points.len());
38        unique_points.push(delaunator::Point {
39            x: point[0],
40            y: point[1],
41        });
42        original_indices.push(idx);
43    }
44
45    if unique_points.len() < 3 {
46        return Ok(Delaunay2d {
47            triangles: Vec::new(),
48        });
49    }
50
51    let triangulation = delaunator::triangulate(&unique_points);
52    let mut triangles = Vec::with_capacity(triangulation.triangles.len() / 3);
53    for tri in triangulation.triangles.chunks_exact(3) {
54        triangles.push([
55            original_indices[tri[0]],
56            original_indices[tri[1]],
57            original_indices[tri[2]],
58        ]);
59    }
60    Ok(Delaunay2d { triangles })
61}
62
63fn stable_float_key(value: f64) -> u64 {
64    if value == 0.0 {
65        0.0f64.to_bits()
66    } else {
67        value.to_bits()
68    }
69}
70
71pub fn boundary_edges(triangles: &[[usize; 3]]) -> Vec<[usize; 2]> {
72    let mut counts: HashMap<(usize, usize), usize> = HashMap::with_capacity(triangles.len() * 3);
73    for tri in triangles {
74        for [a, b] in [[tri[0], tri[1]], [tri[1], tri[2]], [tri[2], tri[0]]] {
75            let edge = if a <= b { (a, b) } else { (b, a) };
76            *counts.entry(edge).or_insert(0) += 1;
77        }
78    }
79    let mut edges = counts
80        .into_iter()
81        .filter_map(|(edge, count)| (count == 1).then_some([edge.0, edge.1]))
82        .collect::<Vec<_>>();
83    edges.sort_unstable();
84    edges
85}
86
87pub fn nearest_neighbor_indices(points: &[[f64; 2]], queries: &[[f64; 2]]) -> Vec<Option<usize>> {
88    queries
89        .iter()
90        .map(|query| {
91            let mut best = None;
92            let mut best_distance = f64::INFINITY;
93            for (idx, point) in points.iter().enumerate() {
94                let dx = query[0] - point[0];
95                let dy = query[1] - point[1];
96                let distance = dx.mul_add(dx, dy * dy);
97                if distance < best_distance {
98                    best_distance = distance;
99                    best = Some(idx);
100                }
101            }
102            best
103        })
104        .collect()
105}
106
107pub fn point_locations(
108    points: &[[f64; 2]],
109    triangles: &[[usize; 3]],
110    queries: &[[f64; 2]],
111) -> Vec<(Option<usize>, [f64; 3])> {
112    queries
113        .iter()
114        .map(|query| {
115            for (idx, tri) in triangles.iter().enumerate() {
116                let bary = barycentric(points[tri[0]], points[tri[1]], points[tri[2]], *query);
117                if bary
118                    .iter()
119                    .all(|value| *value >= -1.0e-12 && *value <= 1.0 + 1.0e-12)
120                {
121                    return (Some(idx), bary);
122                }
123            }
124            (None, [f64::NAN; 3])
125        })
126        .collect()
127}
128
129fn barycentric(a: [f64; 2], b: [f64; 2], c: [f64; 2], p: [f64; 2]) -> [f64; 3] {
130    let v0 = [b[0] - a[0], b[1] - a[1]];
131    let v1 = [c[0] - a[0], c[1] - a[1]];
132    let v2 = [p[0] - a[0], p[1] - a[1]];
133    let denominator = v0[0] * v1[1] - v1[0] * v0[1];
134    if denominator.abs() <= f64::EPSILON {
135        return [f64::NAN; 3];
136    }
137    let v = (v2[0] * v1[1] - v1[0] * v2[1]) / denominator;
138    let w = (v0[0] * v2[1] - v2[0] * v0[1]) / denominator;
139    [1.0 - v - w, v, w]
140}
141
142#[cfg(test)]
143mod tests {
144    use super::*;
145
146    #[test]
147    fn triangulates_square_as_two_one_based_ready_faces() {
148        let mesh = delaunay_2d(&[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]).unwrap();
149        assert_eq!(mesh.triangles.len(), 2);
150        assert!(mesh
151            .triangles
152            .iter()
153            .all(|tri| tri.iter().all(|idx| *idx < 4)));
154    }
155
156    #[test]
157    fn boundary_edges_ignore_shared_interior_edge() {
158        let edges = boundary_edges(&[[0, 1, 2], [1, 3, 2]]);
159        assert_eq!(edges, vec![[0, 1], [0, 2], [1, 3], [2, 3]]);
160    }
161
162    #[test]
163    fn point_locations_return_triangle_and_barycentric_coordinates() {
164        let locations = point_locations(
165            &[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]],
166            &[[0, 1, 2]],
167            &[[0.25, 0.25], [2.0, 2.0]],
168        );
169        assert_eq!(locations[0].0, Some(0));
170        assert!((locations[0].1[0] - 0.5).abs() < 1.0e-12);
171        assert_eq!(locations[1].0, None);
172        assert!(locations[1].1[0].is_nan());
173    }
174}
runmat_geometry_ops/triangulation.rs

runmat_geometry_ops/
triangulation.rs
