procedural_modelling/mesh/
triangulation.rs

1use crate::math::{IndexType, LineSegment2D, Polygon, Scalar, ScalarIteratorExt, Vector2D};
2use std::collections::{HashMap, HashSet};
3
4/// A vertex with its index in the global structure
5#[derive(Debug, Clone, Copy, PartialEq, Eq)]
6pub struct IndexedVertex2D<V: IndexType, Vec2: Vector2D> {
7    /// Position of the point
8    pub vec: Vec2,
9    /// Index in the global structure
10    pub index: V,
11}
12
13impl<V: IndexType, Vec2: Vector2D> IndexedVertex2D<V, Vec2> {
14    /// Create a new indexed vertex
15    pub fn new(vec: Vec2, index: V) -> Self {
16        IndexedVertex2D { vec, index }
17    }
18
19    /// Convert a vector of Vector2Ds to a vector of indexed vertices
20    pub fn from_vector(vec: Vec<Vec2>) -> Vec<Self> {
21        vec.into_iter()
22            .enumerate()
23            .map(|(i, v)| IndexedVertex2D::new(v, V::new(i)))
24            .collect()
25    }
26}
27
28/// A triangulation of a polygon.
29/// Will borrow the index buffer and append new triangles to it.
30/// Most methods will only look at the indices that are added after the borrow startet.
31/// It's fine to add triangles to the index buffer directly while it is borrowed.
32pub struct Triangulation<'a, V: IndexType> {
33    /// The index buffer
34    indices: &'a mut Vec<V>,
35
36    /// The position of the index where _this_ `Triangulation` begins
37    start: usize,
38}
39
40impl<V: IndexType> std::fmt::Debug for Triangulation<'_, V> {
41    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
42        write!(
43            f,
44            "Triangulation({} triangles; start {})",
45            self.len(),
46            self.start
47        )?;
48        for i in 0..self.len() {
49            let (a, b, c) = self.get_triangle(i);
50            write!(f, "\n{} {} {}", a, b, c)?;
51        }
52        Ok(())
53    }
54}
55
56impl<'a, V: IndexType> Triangulation<'a, V> {
57    /// Create a new triangulation
58    pub fn new(indices: &'a mut Vec<V>) -> Self {
59        Triangulation {
60            start: indices.len(),
61            indices,
62        }
63    }
64
65    /// Insert a triangle into the triangulation using global indices
66    pub fn insert_triangle(&mut self, a: V, b: V, c: V) {
67        self.indices.extend([a, b, c]);
68    }
69
70    /// Get the ith index that was added to the triangulation
71    pub fn get_index(&self, i: usize) -> V {
72        self.indices[self.start + i]
73    }
74
75    /// Get a triangle from the triangulation using the number of the triangle in the triangulation
76    pub fn get_triangle(&self, i: usize) -> (V, V, V) {
77        (
78            self.indices[self.start + 3 * i],
79            self.indices[self.start + 3 * i + 1],
80            self.indices[self.start + 3 * i + 2],
81        )
82    }
83
84    /// Get the area of a triangle in the triangulation
85    pub fn get_triangle_area<Vec2: Vector2D>(
86        &self,
87        i: usize,
88        vec_hm: &HashMap<V, Vec2>,
89    ) -> Vec2::S {
90        let (i1, i2, i3) = self.get_triangle(i);
91        let v0 = vec_hm[&i1];
92        let v1 = vec_hm[&i2];
93        let v2 = vec_hm[&i3];
94
95        // Use the determinant to calculate the area of the triangle
96        (v1.x() - v0.x()) * (v2.y() - v0.y()) - (v1.y() - v0.y()) * (v2.x() - v0.x())
97    }
98
99    /// Insert a triangle into the triangulation using local indices
100    pub fn insert_triangle_local<Vec2: Vector2D>(
101        &mut self,
102        a: usize,
103        b: usize,
104        c: usize,
105        vec2s: &Vec<IndexedVertex2D<V, Vec2>>,
106    ) {
107        self.indices
108            .extend([vec2s[a].index, vec2s[b].index, vec2s[c].index]);
109    }
110
111    /// Map the indices in the triangulation using a hashmap
112    pub fn map_indices(&mut self, id_map: &HashMap<V, V>) {
113        for i in self.start..self.indices.len() {
114            self.indices[i] = id_map[&self.indices[i]];
115        }
116    }
117
118    /// Get the number of triangles inserted into the index buffer since the triangulation was created
119    pub fn len(&self) -> usize {
120        let n = self.indices.len() - self.start;
121        assert!(n % 3 == 0, "Invalid number of indices in triangulation");
122        n / 3
123    }
124
125    /// Get the next index that will be added to the index buffer
126    pub fn next_pos(&self) -> usize {
127        self.indices.len()
128    }
129
130    /// Flip the edge of the two triangles
131    pub fn flip_edge(
132        &mut self,
133        a: V,
134        b: V,
135        triangle_ab: usize,
136        triangle_ba: usize,
137    ) -> Result<(), ()> {
138        let offset_ab = if self.indices[triangle_ab + 0] == a {
139            0
140        } else if self.indices[triangle_ab + 1] == a {
141            1
142        } else {
143            2
144        };
145        if self.indices[triangle_ab + offset_ab] != a
146            || self.indices[triangle_ab + ((offset_ab + 1) % 3)] != b
147        {
148            return Err(());
149        }
150
151        let offset_ba = if self.indices[triangle_ba + 0] == a {
152            0
153        } else if self.indices[triangle_ba + 1] == a {
154            1
155        } else {
156            2
157        };
158        if self.indices[triangle_ba + offset_ba] != a
159            || self.indices[triangle_ba + ((offset_ba + 2) % 3)] != b
160        {
161            return Err(());
162        }
163
164        let c = self.indices[triangle_ab + ((offset_ab + 2) % 3)];
165        let d = self.indices[triangle_ba + ((offset_ba + 1) % 3)];
166
167        // Apply the flip
168        // abc -> adc
169        // adb -> dbc
170        self.indices[triangle_ab + 0] = a;
171        self.indices[triangle_ab + 1] = d;
172        self.indices[triangle_ab + 2] = c;
173        self.indices[triangle_ba + 0] = d;
174        self.indices[triangle_ba + 1] = b;
175        self.indices[triangle_ba + 2] = c;
176
177        Ok(())
178    }
179
180    /// Check for non-degenerate triangles (no zero-area triangles)
181    pub fn verify_non_degenerate_triangle<Vec2: Vector2D>(&self, vec_hm: &HashMap<V, Vec2>) {
182        for i in self.start..self.len() {
183            let area = self.get_triangle_area(i, vec_hm);
184            /*assert!(
185                area.abs() > Vec2::S::ZERO,
186                "Triangle has zero or negative area"
187            );*/
188            // degenerate triangles are ok. But not negative ones!
189            if !(area >= -Vec2::S::EPS.sqrt()) {
190                println!("Triangle area: {}", area);
191                assert!(area >= -Vec2::S::EPS.sqrt(), "Triangle has negative area");
192            }
193        }
194    }
195
196    /// Check for valid indices (i.e., they should be within the bounds of the vertices)
197    pub fn verify_indices<Vec2: Vector2D>(&self, vec_hm: &HashMap<V, Vec2>) {
198        // Check that the triangulation returns the correct number of triangles
199        let num_vertices = vec_hm.len();
200        let num_triangles = self.len();
201
202        assert!(
203            num_triangles == num_vertices - 2,
204            "Expected {} triangles but found {}",
205            num_vertices - 2,
206            num_triangles
207        );
208        for i in self.start..self.indices.len() {
209            assert!(
210                vec_hm.get(&self.indices[i]).is_some(),
211                "Index {} out of bounds in triangulation",
212                self.indices[i]
213            );
214        }
215    }
216
217    /// Check that no two triangles have intersecting edges
218    pub fn verify_no_intersections<Vec2: Vector2D>(&self, vec_hm: &HashMap<V, Vec2>) {
219        let num_vertices = self.indices.len() - self.start;
220        /*for i in (0..num_vertices).step_by(3) {
221            println!(
222                "tri: {:?}",
223                (
224                    self.get_index(i),
225                    self.get_index(i + 1),
226                    self.get_index(i + 2)
227                ),
228            );
229        }*/
230
231        for i in (self.start..num_vertices).step_by(3) {
232            for j in (self.start..num_vertices).step_by(3) {
233                if i == j {
234                    continue;
235                }
236                for k in 0..3 {
237                    for l in 0..3 {
238                        let i0 = self.get_index(i + k);
239                        let v0 = vec_hm[&i0];
240                        let i1 = self.get_index(i + (k + 1) % 3);
241                        let v1 = vec_hm[&i1];
242
243                        let i2 = self.get_index(j + l);
244                        let v2 = vec_hm[&i2];
245                        let i3 = self.get_index(j + (l + 1) % 3);
246                        let v3 = vec_hm[&i3];
247
248                        // If they share a vertex, they can't intersect
249                        if i0 == i2 || i0 == i3 || i1 == i2 || i1 == i3 {
250                            continue;
251                        }
252
253                        let l1 = LineSegment2D::new(v0, v1);
254                        let l2 = LineSegment2D::new(v2, v3);
255                        let length = l1.length() + l2.length();
256                        let inter = l1.intersect_line(
257                            &l2,
258                            Vec2::S::EPS.sqrt(), // be strict about parallel edges
259                            -Vec2::S::EPS.sqrt() * length, // Allow intersections/touching at the endpoints up to a portion of sqrt(eps), i.e., 0.0345% for f32
260                        );
261                        assert!(
262                            inter.is_none(),
263                            "Edges: \n{} {:?} -> {} {:?}\n{} {:?} -> {} {:?}\nintersect in {:?} (shortest distance: {} * sqrt(eps))\nTriangles {:?} and {:?}",
264                            i0,
265                            v0,
266                            i1,
267                            v1,
268                            i2,
269                            v2,
270                            i3,
271                            v3,
272                            inter.unwrap(),
273                            [v0,v1,v2,v3].iter().map(|v| inter.unwrap().distance(&v)).min_by(|a, b| a.partial_cmp(b).unwrap()).unwrap() / Vec2::S::EPS.sqrt(),
274                            (self.get_index(i), self.get_index(i+1), self.get_index(i+2)),
275                            (self.get_index(j), self.get_index(j+1), self.get_index(j+2)),
276                        );
277                    }
278                }
279            }
280        }
281    }
282
283    /// Sum the area of all triangles added to the index buffer since the triangulation was created
284    pub fn get_area<Vec2: Vector2D>(&self, vec_hm: &HashMap<V, Vec2>) -> Vec2::S {
285        Vec2::S::HALF
286            * (0..self.len())
287                .into_iter()
288                .map(|i| self.get_triangle_area(i, vec_hm).abs())
289                .stable_sum()
290    }
291
292    /// Calculate the total edge weight of the triangulation
293    pub fn total_edge_weight<Vec2: Vector2D>(&self, vec_hm: &HashMap<V, Vec2>) -> Vec2::S {
294        let mut total = Vec2::S::ZERO;
295        for i in self.start..self.len() {
296            let (i1, i2, i3) = self.get_triangle(i);
297            let v0 = vec_hm[&i1];
298            let v1 = vec_hm[&i2];
299            let v2 = vec_hm[&i3];
300            total += v1.distance(&v0) + v2.distance(&v1) + v0.distance(&v2);
301        }
302        total
303    }
304
305    /// Calculate the area of the polygon and check if it is the same as the sum of the areas of the triangles
306    pub fn verify_area<Vec2: Vector2D, Poly: Polygon<Vec2>>(
307        &self,
308        vec2s: &Vec<IndexedVertex2D<V, Vec2>>,
309        vec_hm: &HashMap<V, Vec2>,
310    ) {
311        let area = self.get_area(vec_hm);
312        let reference = Poly::from_iter(vec2s.iter().map(|v| v.vec)).area();
313
314        // Check if the area of the polygon is the same as the sum of the areas of the triangles
315        assert!(
316            (Vec2::S::ONE - area / reference).abs()
317                <= (Vec2::S::ONE + Vec2::S::from_usize(5) * Vec2::S::EPS),
318            "Area of the polygon is not equal to the sum of the areas of the triangles ({} != {})",
319            area,
320            reference
321        );
322    }
323
324    /// Check that the set of used indices exactly matches the set of indices in the triangulation
325    pub fn verify_all_indices_used<Vec2: Vector2D>(&self, vec2s: &Vec<IndexedVertex2D<V, Vec2>>) {
326        let mut seen = HashSet::new();
327        for i in self.start..self.indices.len() {
328            seen.insert(self.indices[i]);
329        }
330
331        for vertex in vec2s.iter() {
332            assert!(
333                seen.remove(&vertex.index),
334                "Vertex not used in triangulation {}",
335                vertex.index.index()
336            );
337        }
338
339        assert!(
340            seen.is_empty(),
341            "Foreign indices used in triangulation: {:?}",
342            seen.iter().map(|i| i.index()).collect::<Vec<_>>()
343        );
344    }
345
346    /// Runs a large number of tests on the triangulation to verify that it is well-formed
347    pub fn verify_full<Vec2: Vector2D, Poly: Polygon<Vec2>>(
348        &self,
349        vec2s: &Vec<IndexedVertex2D<V, Vec2>>,
350    ) {
351        let vec_hm: HashMap<V, Vec2> = vec2s.iter().map(|v| (v.index, v.vec)).collect();
352
353        self.verify_indices(&vec_hm);
354        self.verify_all_indices_used(&vec2s);
355        self.verify_no_intersections(&vec_hm);
356        self.verify_non_degenerate_triangle(&vec_hm);
357        self.verify_area::<Vec2, Poly>(&vec2s, &vec_hm);
358    }
359}
procedural_modelling/mesh/triangulation.rs

procedural_modelling/mesh/
triangulation.rs
