cg_math/utils/
mod.rs

1extern crate nalgebra as na;
2extern crate parry2d_f64 as parry;
3
4use crate::geometry::{WeightedEdge, WeightedTriangle, WeightedVertex, Triangulation, TriangleLogic};
5use na::Point2;
6use rand::{distributions::Uniform, rngs::ThreadRng, thread_rng, Rng};
7use std::collections::HashMap;
8
9/// Sample n vertices (i.e. w_i = 0.0), with coordinates drawn form a Uniform distribution
10#[allow(dead_code)]
11pub fn generate_vertices(n: u32) -> Vec<WeightedVertex> {
12    let mut rng: ThreadRng = thread_rng();
13    let distribution: Uniform<f64> = Uniform::new(0.0, 10.0);
14
15    (0..n)
16        .map(|_| WeightedVertex::new(rng.sample(distribution), rng.sample(distribution), 0.0))
17        .collect()
18}
19
20/// Sample n weighted vertices, with coordinates and weights drawn form a Uniform distribution
21#[allow(dead_code)]
22pub fn generate_weighted_vertices(n: u32) -> Vec<WeightedVertex> {
23    let mut rng: ThreadRng = thread_rng();
24    let distribution: Uniform<f64> = Uniform::new(0.0, 10.0);
25
26    (0..n)
27        .map(|_| {
28            WeightedVertex::new(
29                rng.sample(distribution),
30                rng.sample(distribution),
31                rng.sample(distribution),
32            )
33        })
34        .collect()
35}
36
37/// Given an edge of the triangle and, retrieve the point not on the edge
38#[allow(dead_code)]
39pub fn get_fourth_point(triangle: &WeightedTriangle, edge: &WeightedEdge) -> WeightedVertex {
40    if triangle.a != edge.a && triangle.a != edge.b {
41        WeightedVertex::new(triangle.a.x(), triangle.a.y(), triangle.a.w())
42    } else if triangle.b != edge.a && triangle.b != edge.b {
43        WeightedVertex::new(triangle.b.x(), triangle.b.y(), triangle.b.w())
44    } else {
45        WeightedVertex::new(triangle.c.x(), triangle.c.y(), triangle.c.w())
46    }
47}
48
49/// Compute a points barycentrics coordinates w.r.t. the given triangle
50pub fn barycentric_coords(vertex: &WeightedVertex, triangle: &WeightedTriangle) -> (f64, f64, f64) {
51    // adapted from: https://gamedev.stackexchange.com/questions/23743/whats-the-most-efficient-way-to-find-barycentric-coordinates
52    let v0 = triangle.b.as_vec2() - triangle.a.as_vec2();
53    let v1 = triangle.c.as_vec2() - triangle.a.as_vec2();
54    let v2 = vertex.as_vec2() - triangle.a.as_vec2();
55
56    let d00 = v0.dot(&v0);
57    let d01 = v0.dot(&v1);
58    let d11 = v1.dot(&v1);
59    let d20 = v2.dot(&v0);
60    let d21 = v2.dot(&v1);
61
62    let inv_denom = 1. / (d00 * d11 - d01 * d01); // less divisions for optimization
63
64    let v = (d11 * d20 - d01 * d21) * inv_denom;
65    let w = (d00 * d21 - d01 * d20) * inv_denom;
66    let u = 1. - v - w;
67
68    assert!(v + w + u - 1. < 0.0001); // sanity check
69
70    (v, w, u)
71}
72
73/// Check whether the vertex lies inside the given triangle
74pub fn point_in_tri(vertex: &WeightedVertex, triangle: &WeightedTriangle) -> bool {
75    let (v, w, u) = barycentric_coords(vertex, triangle);
76    v >= 0.0 && w >= 0.0 && u >= 0.0
77}
78
79/// Compute the convex hull of the given points in 2D
80pub fn c_hull(points: &Vec<WeightedVertex>) -> Vec<WeightedVertex> {
81    let points2d: Vec<Point2<f64>> = points.iter().map(|v| Point2::new(v.x(), v.y())).collect();
82    let c_hull = parry::transformation::convex_hull(&points2d);
83    
84    c_hull
85        .iter()
86        .map(|p| {
87            // since order is not invariant, we need to map the weights back to the correct vertices.
88            for v in points {
89                if v.x() == p.x && v.y() == p.y {
90                    return WeightedVertex::new(p.x, p.y, v.w());
91                }
92            }
93            panic!("One of the vertices should match the convex hull vertices.")
94        })
95        .collect()
96}
97
98/// A 1-simplex and its two incident 2-simplices yield four points, of which all 2-simplices of the given triangulation that
99/// consist of any permutation of 3 of these vertices form the induced subcomplex.
100/// 
101/// More formally: <br>
102/// Let e be and edge of an arbitrary triangulation T and t1, t2 be it's incident triangles (given they exist). <br>
103/// The induced subcomplex of the four vertices of e, t1, t2, namely V, consists of all simplices in T <br>
104/// spanned by vertices in V.
105/// 
106/// ```
107/// use cg_math::{geometry::{WeightedVertex, WeightedTriangle, Triangulation}, utils::induced_subcomplex};
108/// 
109/// // these 4 vertices span two of the three triangles in the triangulation
110/// let v0 = WeightedVertex::new(6.24, 8.38, 2.24);
111/// let v1 = WeightedVertex::new(9.34, 9.19, 7.46);
112/// let v2 = WeightedVertex::new(8.04, 5.53, 7.53);
113/// let v3 = WeightedVertex::new(8.53, 8.43, 1.34);
114/// 
115/// let mut triangulation = Triangulation::empty();
116/// triangulation.triangles_mut().append(&mut vec![
117///     WeightedTriangle::new(v0, v1, v3),
118///     WeightedTriangle::new(v2 ,v1, v3),
119///     WeightedTriangle::new(v2 ,v1, v0),
120/// ]);
121/// 
122/// let mut ics = induced_subcomplex(&vec![v0, v1, v2, v3], &triangulation);
123/// ics.sort();
124/// triangulation.triangles_mut().sort();
125/// 
126/// // Yet their induced subcomplex consists of all three triangles.
127/// assert!(ics.len() == 3);
128/// assert_eq!(ics, *triangulation.triangles());
129/// ```
130pub fn induced_subcomplex(vertices: &[WeightedVertex], triangulation: &Triangulation) -> Vec<WeightedTriangle> {
131    let subcomplex: Vec<WeightedTriangle> = triangulation.triangles().iter().filter(|tri| {
132        tri.has_vertex(&vertices[0]) && tri.has_vertex(&vertices[1]) && tri.has_vertex(&vertices[2]) 
133        || tri.has_vertex(&vertices[0]) && tri.has_vertex(&vertices[1]) && tri.has_vertex(&vertices[3])
134        || tri.has_vertex(&vertices[1]) && tri.has_vertex(&vertices[2]) && tri.has_vertex(&vertices[3])
135        || tri.has_vertex(&vertices[0]) && tri.has_vertex(&vertices[2]) && tri.has_vertex(&vertices[3])
136    }).copied().collect();
137    subcomplex
138}
139
140
141/// An edge is flippable if the pointwise union of its induced subcomplex equals the convex hull of all the vertices in the induced subcomplex.
142/// ```
143/// use cg_math::{geometry::{WeightedVertex, WeightedTriangle, Triangulation}, utils::{induced_subcomplex, is_flippable2}};
144/// 
145/// // these 4 vertices span two that would need to be 3-to-1 flipped
146/// let v0 = WeightedVertex::new(6.24, 8.38, 2.24); 
147/// let v1 = WeightedVertex::new(9.34, 9.19, 7.46);
148/// let v2 = WeightedVertex::new(8.04, 5.53, 7.53);
149/// let v3 = WeightedVertex::new(8.53, 8.43, 1.34);
150/// 
151/// let mut triangulation = Triangulation::empty();
152/// triangulation.triangles_mut().append(&mut vec![
153///     WeightedTriangle::new(v0, v1, v3),
154///     WeightedTriangle::new(v2 ,v1, v3),
155///     WeightedTriangle::new(v2 ,v1, v0),
156/// ]);
157/// 
158/// let mut ics = induced_subcomplex(&vec![v0, v1, v2, v3], &triangulation);
159/// for tri in ics.iter() {
160///     println!("{}", tri);
161///  }
162/// assert_eq!(is_flippable2(&ics), true);
163/// ```
164pub fn  is_flippable_strict(induced_subcomplex: &[WeightedTriangle]) -> bool {
165    let points = induced_subcomplex
166        .iter()
167        .flat_map(|t| [Point2::new(t.a.x(), t.a.y()), Point2::new(t.b.x(), t.b.y()), Point2::new(t.c.x(), t.c.y())] )
168        .collect::<Vec<Point2<f64>>>();
169
170    let c_hull = parry::transformation::convex_hull(&points);
171    // is even faster for low number of items, than hashing, s.: https://stackoverflow.com/questions/64226562/check-if-vec-contains-all-elements-from-another-vec
172    points.iter().all(|item| c_hull.contains(item))
173}
174
175/// Check whether the edge is flippable in the ploygon described all four points. <br>
176/// "*We call T and e flippable if conv(T) is the underlying space of the induced subcomplex of T.*"
177pub fn is_flippable(point: &WeightedVertex, edge: &WeightedEdge, pk: &WeightedVertex) -> bool {
178    // we could also check if the polygon is convex
179    // hopefully more correct parry version
180    let points: Vec<Point2<f64>> = vec![    
181        Point2::new(point.x(), point.y()),
182        Point2::new(edge.a.x(), edge.a.y()),
183        Point2::new(edge.b.x(), edge.b.y()),
184        Point2::new(pk.x(), pk.y()),
185    ];
186
187    //let points2d: Vec<Point2<f64>> = points.iter().map(|v| Point2::new(v.x, v.y)).collect();
188    let c_hull = parry::transformation::convex_hull(&points);
189
190    //println!("points: {:?}", points);
191    //println!("c_hull: {:?}", c_hull);
192
193    // is even faster for low number of items, than hashing, s.: https://stackoverflow.com/questions/64226562/check-if-vec-contains-all-elements-from-another-vec
194    points.iter().all(|item| c_hull.contains(item))
195}
196
197/// Keep the elements that are unique in the given vector of elements, i.e. that appear exactly once.
198pub fn keep_unique_elements<T: Eq + std::hash::Hash>(elements: Vec<T>) -> Vec<T> {
199    elements
200        .into_iter()
201        .fold(HashMap::new(), |mut map, element| {
202            *map.entry(element).or_insert(0) += 1;
203            map
204        })
205        .into_iter()
206        .filter_map(|(element, count)| (count == 1).then_some(element))
207        .collect()
208}
209
210/// Older, less performant version than `keep_unique_elements`
211/// as found out by Tbolt on discord in rust-questions-and-answers-2
212#[allow(dead_code)]
213pub fn keep_unique_elements2<T: Eq + std::hash::Hash + Clone>(elements: &mut Vec<T>) {
214    let mut freq_map = HashMap::new();
215
216    // count frequency of each element
217    for element in elements.iter() {
218        *freq_map.entry(element.clone()).or_insert(0) += 1;
219    }
220
221    elements.retain(|element| freq_map.get(element) == Some(&1));
222}
223
224// Todo: can we do this more simple / efficient ? -> yes this is O(n^2), O(n) version just above
225/// Deprecated: more efficient and readable version is `keep_unique_elements`; this is kept for learning purposes <br>
226/// Remove duplicates from a list of elements of type T.
227/// So keeping only unique elemtents
228#[allow(warnings)]
229pub fn remove_duplicates<T: std::cmp::PartialEq<T>>(elements: &mut Vec<T>) {
230    let mut indices_to_remove: Vec<usize> = vec![];
231
232    for i in 0..elements.len() {
233        let element: &T = &elements[i];
234        for j in 0..elements.len() {
235            let comp_element: &T = &elements[j];
236            if element == comp_element && i != j {
237                indices_to_remove.push(i);
238                break;
239            }
240        }
241    }
242    // since idx has elements in ascending order, e.g. [2, 7, 13]
243    // we can ierate in reverse and not conflic any array boundaries
244    for idx in indices_to_remove.iter().rev() {
245        elements.remove(*idx);
246    }
247}
248
249#[cfg(test)]
250mod tests {
251    use super::*;
252    use na::{Point2, Vector2};
253
254    #[test]
255    fn test_barycentric_coords() {
256        let tri = WeightedTriangle::new(
257            WeightedVertex::new(0.0, 0.0, 0.0),
258            WeightedVertex::new(1.0, 0.0, 0.0),
259            WeightedVertex::new(0.0, 1.0, 0.0),
260        );
261
262        let point = WeightedVertex::new(0.5, 0.5, 0.0);
263        let (v, w, u) = barycentric_coords(&point, &tri);
264        assert_eq!(v, 0.5);
265        assert_eq!(w, 0.5);
266        assert_eq!(u, 0.0);
267
268        let point = WeightedVertex::new(0.0, 0.0, 0.0);
269        let (v, w, u) = barycentric_coords(&point, &tri);
270        assert_eq!(v, 0.0);
271        assert_eq!(w, 0.0);
272        assert_eq!(u, 1.0);
273
274        let point = WeightedVertex::new(1.0, 0.0, 0.0);
275        let (v, w, u) = barycentric_coords(&point, &tri);
276        assert_eq!(v, 0.0);
277        assert_eq!(w, 1.0);
278        assert_eq!(u, 0.0);
279
280        let point = WeightedVertex::new(0.0, 1.0, 0.0);
281        let (v, w, u) = barycentric_coords(&point, &tri);
282        assert_eq!(v, 1.0);
283        assert_eq!(w, 0.0);
284        assert_eq!(u, 0.0);
285    }
286
287    #[test]
288    fn test_convex_hull() {
289        let points = vec![
290            Vector2::new(0.5, 0.5),
291            Vector2::new(2.5, 2.0),
292            Vector2::new(1.0, 2.5),
293            Vector2::new(2.0, 3.0),
294            Vector2::new(1.0, 1.0),
295        ];
296
297        let true_hull = vec![
298            Point2::new(0.5, 0.5),
299            Point2::new(2.5, 2.0),
300            Point2::new(2.0, 3.0),
301            Point2::new(1.0, 2.5),
302        ];
303
304        let points2d: Vec<Point2<f64>> = points.iter().map(|v| Point2::new(v.x, v.y)).collect();
305        let c_hull = parry::transformation::convex_hull(&points2d);
306        //println!("points: {:?}", points2d);
307        //println!("c_hull: {:?}", c_hull);
308        //println!("true_hull: {:?}", true_hull);
309        assert!(c_hull != points2d);
310        assert!(true_hull.iter().all(|item| true_hull.contains(item)));
311    }
312
313    #[test]
314    fn test_edge_flippable() {
315        pub const LARGE_NUM: f64 = 10000.0; // NOTE: this parameter highly influences the convex hull property of the triangulation
316
317        let p_2_lw: WeightedVertex = WeightedVertex::new(0.0, LARGE_NUM, 0.0);
318
319        let p_1_lw: WeightedVertex = WeightedVertex::new(LARGE_NUM, -LARGE_NUM, 0.0);
320
321        // here the edge in question is the one between (2.5, 2.0) and (1.0, 2.5)
322        // the two triangles are t1 = {(2.5, 2.0), (0.5, 0.5), (1.0, 2.5)} and t2 = {(2.5, 2.0),(1.0, 2.5), (2.0, 3.0)}
323        assert!(is_flippable(
324            &WeightedVertex::new(0.5, 0.5, 0.0),
325            &WeightedEdge::new(
326                WeightedVertex::new(2.5, 2.0, 0.0),
327                WeightedVertex::new(1.0, 2.5, 0.0)
328            ),
329            &WeightedVertex::new(2.0, 3.0, 0.0)
330        ));
331
332        // this is acutally flippable, try with low numbers for L and if necessary make a sketch
333        assert!(is_flippable(
334            &WeightedVertex::new(2.5, 2.0, 0.0),
335            &WeightedEdge::new(WeightedVertex::new(2.0, 3.0, 0.0), p_2_lw),
336            &p_1_lw
337        ));
338    }
339
340    #[test]
341    fn test_edge_not_flippable() {
342        pub const LARGE_NUM: f64 = 10000.0; // NOTE: this parameter highly influences the convex hull property of the triangulation
343
344        let p_2_lw: WeightedVertex = WeightedVertex::new(0.0, LARGE_NUM, 0.0);
345
346        assert!(!is_flippable(
347            &WeightedVertex::new(-5., 1., 0.0),
348            &WeightedEdge::new(
349                WeightedVertex::new(-4., 2., 0.0),
350                WeightedVertex::new(-2.5, 2., 0.0),
351            ),
352            &WeightedVertex::new(-5., 3., 0.0)
353        ));
354
355        assert!(!is_flippable(
356            &WeightedVertex::new(-1.5, 0., 0.0),
357            &WeightedEdge::new(WeightedVertex::new(0., 1.5, 0.0), p_2_lw),
358            &WeightedVertex::new(2., -0.5, 0.0)
359        ));
360    }
361}
cg_math/utils/mod.rs

cg_math/utils/
mod.rs
