cg_math/
geometry.rs

1extern crate nalgebra as na;
2
3use core::panic;
4use egui::plot::PlotPoints;
5use na::{Matrix3, Matrix4, Vector2, Vector3};
6use ordered_float::NotNan;
7use std::fmt;
8
9/// Describes common behaviour between regular and weighted triangles <br>
10// V can be Vertex or WeightedVertex
11// E can be Edge or WeightedEdge
12// T can be Triangle or WeightedTriangle
13pub trait TriangleLogic<V, E, T> {
14    fn area(&self) -> f64;
15    fn area3d(&self) -> f64;
16    fn edges(&self) -> Vec<E>;
17    fn has_edge(&self, edge: &E) -> bool;
18    fn has_vertex(&self, vertex: &V) -> bool;
19    fn is_neighbor(&self, other: &T) -> Option<E>;
20    fn regularity(&self, other: &V) -> f64;
21    /// Check whether the triangle is regular, w.r.t to some vertex; by tolerance of epsilon. <br>
22    fn is_regular(&self, vertex: &V) -> bool;
23    fn is_eps_regular(&self, vertex: &V, epsilon: f64) -> bool;
24    fn orientation(&self) -> f64;
25    fn vertices(&self) -> Vec<WeightedVertex>;
26}
27
28/// A weighted 0-simplex in 2-dimensional space.
29#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, PartialOrd, Ord)]
30pub struct WeightedVertex {
31    x: NotNan<f64>,
32    y: NotNan<f64>,
33    w: NotNan<f64>,
34}
35
36impl WeightedVertex {
37    pub fn new(x: f64, y: f64, weight: f64) -> Self {
38        Self {
39            x: NotNan::new(x).unwrap(),
40            y: NotNan::new(y).unwrap(),
41            w: NotNan::new(weight).unwrap(),
42        }
43    }
44
45    pub fn as_vec2(&self) -> Vector2<f64> {
46        Vector2::new(self.x(), self.y())
47    }
48
49    // Lifts the vertex to 3d such that $z = xˆ2 + yˆ2 - w$
50    pub fn lift(&self) -> Vector3<f64> {
51        Vector3::new(self.x(), self.y(), self.x().powi(2) + self.y().powi(2) - self.w())
52    }
53
54    pub fn w(&self) -> f64 {
55        self.w.into_inner()
56    }
57
58    pub fn x(&self) -> f64 {
59        self.x.into_inner()
60    }
61
62    pub fn y(&self) -> f64 {
63        self.y.into_inner()
64    }
65}
66
67impl fmt::Display for WeightedVertex {
68    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
69        write!(f, "({}, {}, w: {})", self.x(), self.y(), self.w())
70    }
71}
72
73/// A weighted 1-simplex in 2-dimensional space.
74#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, PartialOrd, Ord)]
75pub struct WeightedEdge {
76    pub a: WeightedVertex,
77    pub b: WeightedVertex,
78    triangles: [Option<WeightedTriangle>; 2], // store the potential 1, 2 incident triangles
79}
80
81impl WeightedEdge {
82    pub fn new(a: WeightedVertex, b: WeightedVertex) -> Self {
83        let mut vertices = vec![a, b];
84        vertices.sort_unstable();
85        Self {
86            a: vertices[0],
87            b: vertices[1],
88            triangles: [None, None],
89        }
90    }
91
92    /// Create a new edge with the given triangles as incident triangles
93    pub fn with_triangles(
94        a: WeightedVertex,
95        b: WeightedVertex,
96        triangles: [Option<WeightedTriangle>; 2],
97    ) -> Self {
98        let mut vertices = vec![a, b];
99        vertices.sort_unstable();
100        Self {
101            a: vertices[0],
102            b: vertices[1],
103            triangles,
104        }
105    }
106
107    /// Set the incident triangles of the edge
108    pub fn set_triangles(&mut self, triangles: [Option<WeightedTriangle>; 2]) {
109        self.triangles = triangles;
110    }
111
112    pub fn vertices(&self) -> Vec<WeightedVertex> {
113        vec![self.a, self.b]
114    }
115}
116
117impl fmt::Display for WeightedEdge {
118    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
119        write!(f, "Edge {{ {}, {} }}", self.a, self.b)
120    }
121}
122
123/// A weighted 2-simplex in 2-dimensional space.
124#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, PartialOrd, Ord)]
125pub struct WeightedTriangle {
126    pub a: WeightedVertex,
127    pub b: WeightedVertex,
128    pub c: WeightedVertex,
129}
130
131impl WeightedTriangle {
132    pub fn new(a: WeightedVertex, b: WeightedVertex, c: WeightedVertex) -> Self {
133        let mut vertices = vec![a, b, c];
134        vertices.sort_unstable();
135        Self {
136            a: vertices[0],
137            b: vertices[1],
138            c: vertices[2],
139        }
140    }
141
142    /// Given an edge of the triangle retrieve the point not on the edge
143    pub fn get_third_point(&self, edge: &WeightedEdge) -> WeightedVertex {
144        if self.a != edge.a && self.a != edge.b {
145            self.a
146        } else if self.b != edge.a && self.b != edge.b {
147            self.b
148        } else {
149            self.c
150        }
151    }
152}
153
154impl TriangleLogic<WeightedVertex, WeightedEdge, WeightedTriangle> for WeightedTriangle {
155    /// The area of the triangle.
156    /// ### Example
157    /// 
158    /// ```
159    /// use cg_math::geometry::{WeightedVertex, WeightedTriangle, TriangleLogic};
160    /// 
161    /// let v0 = WeightedVertex::new(2.5, 2.0, 0.0);
162    /// let v1 = WeightedVertex::new(1.0, 2.5, 0.0);
163    /// let v2 = WeightedVertex::new(2.0, 3.0, 0.0);
164    /// 
165    /// let triangle = WeightedTriangle::new(v0, v1, v2);
166    /// 
167    /// assert_eq!(triangle.area(), 0.625);
168    /// ```
169    fn area(&self) -> f64 {
170        ((self.a.x() * (self.b.y() - self.c.y())
171            + self.b.x() * (self.c.y() - self.a.y())
172            + self.c.x() * (self.a.y() - self.b.y()))
173            / 2.0)
174            .abs()
175    }
176
177
178    /// The area of the triangle lifted to 3d with $z = x^2 + y^2 - w$.
179    /// ### Example 
180    /// 
181    /// ```
182    /// use cg_math::geometry::{WeightedVertex, WeightedTriangle, TriangleLogic};
183    /// 
184    /// let v0 = WeightedVertex::new(2.5, 2.0, 1.0);
185    /// let v1 = WeightedVertex::new(1.0, 2.5, 2.0);
186    /// let v2 = WeightedVertex::new(2.0, 3.0, 3.0);
187    /// 
188    /// let triangle = WeightedTriangle::new(v0, v1, v2);
189    /// 
190    /// //assert_eq!(triangle.area3d(), 1.3975424859373686);
191    /// ```
192    fn area3d(&self) -> f64 {
193        let ab = self.b.lift() - self.a.lift();
194        let ac = self.c.lift() - self.a.lift();
195        let cross = ab.cross(&ac);
196        let cross_norm = cross.norm();
197        cross_norm / 2.0
198    }
199
200    fn edges(&self) -> Vec<WeightedEdge> {
201        vec![
202            WeightedEdge::new(self.a, self.b),
203            WeightedEdge::new(self.b, self.c),
204            WeightedEdge::new(self.c, self.a),
205        ]
206    }
207
208    /// Check if the given edge is part of the triangle.
209    fn has_edge(&self, edge: &WeightedEdge) -> bool {
210        for e in self.edges() {
211            if e == *edge {
212                return true;
213            }
214        }
215        false
216    }
217
218    /// Check if the given vertex is part of the triangle.
219    fn has_vertex(&self, vertex: &WeightedVertex) -> bool {
220        self.a == *vertex || self.b == *vertex || self.c == *vertex
221    }
222
223    /// Check whether the triangle is regular w.r.t to the given vertex.
224    /// 
225    /// ```
226    /// use cg_math::{geometry::{WeightedVertex, WeightedTriangle, TriangleLogic}, utils::c_hull};
227    /// 
228    /// // these 4 vertices span two of the three triangles in the triangulation
229    /// let v0 = WeightedVertex::new(6.24, 8.38, 2.24);
230    /// let v1 = WeightedVertex::new(9.34, 9.19, 7.46);  // assume this one was inserted last
231    /// let v2 = WeightedVertex::new(8.04, 5.53, 7.53);
232    /// let v3 = WeightedVertex::new(8.53, 8.43, 1.34);
233    /// 
234    /// assert_eq!(WeightedTriangle::new(v0, v2, v3).is_regular(&v1), false);
235    /// ```
236    fn is_regular(&self, vertex: &WeightedVertex) -> bool {
237        self.regularity(vertex) < 0.0
238    }
239
240    /// Check whether the triangle is regular w.r.t to the given vertex; by tolerance of epsilon.
241    fn is_eps_regular(&self, vertex: &WeightedVertex, epsilon: f64) -> bool {
242        self.regularity(vertex) - epsilon < 0.0
243    }
244
245    /// The orientation of the triangle; clockwise (cw) or counter-cw.
246    /// # Panics
247    /// If the triangle consists of three collinear points.
248    fn orientation(&self) -> f64 {
249        let matrix = Matrix3::new(
250            self.a.x(),
251            self.a.y(),
252            1.0,
253            self.b.x(),
254            self.b.y(),
255            1.0,
256            self.c.x(),
257            self.c.y(),
258            1.0,
259        );
260        match matrix.determinant() {
261            d if d > 0.0 => {
262                1.0 //println!("CounterClockwise");
263            }
264            d if d < 0.0 => {
265                -1.0 //println!("Clockwise");
266            }
267            _ => panic!("Triangle consists of three collinear points"),
268        }
269    }
270
271    /// Compute the level of local regularity w.r.t to the given vertex
272    fn regularity(&self, vertex: &WeightedVertex) -> f64 {
273        let orientation = self.orientation();
274
275        // TODO: we could also discuss to shift each weight by epsilon
276        let matrix = Matrix4::new(
277            self.a.x(),
278            self.a.y(),
279            self.a.x().powi(2) + self.a.y().powi(2) - self.a.w(),
280            1.0,
281            self.b.x(),
282            self.b.y(),
283            self.b.x().powi(2) + self.b.y().powi(2) - self.b.w(),
284            1.0,
285            self.c.x(),
286            self.c.y(),
287            self.c.x().powi(2) + self.c.y().powi(2) - self.c.w(),
288            1.0,
289            vertex.x(),
290            vertex.y(),
291            vertex.x().powi(2) + vertex.y().powi(2) - vertex.w(),
292            1.0,
293        );
294
295        //println!("matrix det * orientation: {}", matrix.determinant() * orientation);
296        matrix.determinant() * orientation
297    }
298
299    /// Check if this and the other triangle share an edge, i.e. are neighbors
300    /// and return it if existing
301    fn is_neighbor(&self, other: &WeightedTriangle) -> Option<WeightedEdge> {
302        let self_edges = self.edges();
303        let other_edges = other.edges();
304
305        for self_edge in &self_edges {
306            for other_edge in &other_edges {
307                if *self_edge == *other_edge {
308                    return Some(*self_edge);
309                }
310            }
311        }
312        None
313    }
314
315    fn vertices(&self) -> Vec<WeightedVertex> {
316        vec![self.a, self.b, self.c]
317    }
318}
319
320impl From<&WeightedTriangle> for PlotPoints {
321    fn from(triangle: &WeightedTriangle) -> Self {
322        let points: Vec<[f64; 2]> = triangle
323            .vertices()
324            .iter()
325            .map(|v| [v.x(), v.y()])
326            .collect();
327        PlotPoints::from(points)
328    }
329}
330
331impl fmt::Display for WeightedTriangle {
332    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
333        write!(f, "{{ {}, {}, {} }}", self.a, self.b, self.c)
334    }
335}
336
337/// A triangulation of a set of vertices in 2-dimensional space.
338pub struct WeightedTriangulation<'a>(pub &'a Vec<WeightedTriangle>);
339
340impl fmt::Display for WeightedTriangulation<'_> {
341    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
342        for (index, triangle) in self.0.iter().enumerate() {
343            writeln!(f, "{}: {}", index, triangle)?;
344        }
345        Ok(())
346    }
347}
348
349/// A triangulation of a set of vertices in 2-dimensional space.
350// TODO merge with WeightedTriangulation
351#[derive(PartialEq)]
352pub struct Triangulation {
353    pub triangles: Vec<WeightedTriangle>,
354    pub used_vertices: Vec<WeightedVertex>,
355    pub ignored_vertices: Vec<WeightedVertex>,
356    pub cache: Vec<Vec<WeightedTriangle>>
357}
358
359impl Triangulation {
360    /// Create a new empty trangulation.
361    pub fn empty() -> Self {
362        Self {
363            triangles: vec![],
364            used_vertices: vec![],
365            ignored_vertices: vec![],
366            cache: vec![],
367        }
368    }
369    
370    /// Create a new trangulation with the given starting triangle.
371    pub fn new(starting_triangle: WeightedTriangle) -> Self {
372        Self {
373            triangles: vec![starting_triangle],
374            used_vertices: vec![],
375            ignored_vertices: vec![],
376            cache: vec![],
377        }
378    }
379
380    pub fn triangles(&self) -> &Vec<WeightedTriangle> {
381        &self.triangles
382    }
383
384    pub fn triangles_mut(&mut self) -> &mut Vec<WeightedTriangle> {
385        &mut self.triangles
386    }
387
388    pub fn used_vertices(&self) -> &Vec<WeightedVertex> {
389        &self.used_vertices
390    }
391
392    pub fn used_vertices_mut(&mut self) -> &mut Vec<WeightedVertex> {
393        &mut self.used_vertices
394    }
395
396    pub fn ignored_vertices(&self) -> &Vec<WeightedVertex> {
397        &self.ignored_vertices
398    }
399
400    pub fn ignored_vertices_mut(&mut self) -> &mut Vec<WeightedVertex> {
401        &mut self.ignored_vertices
402    }
403}
404
405#[cfg(test)]
406mod tests {
407    use super::*;
408
409    #[test]
410    fn edges_equal() {
411        let e0 = WeightedEdge::new(
412            WeightedVertex::new(0.0, 1.0, 0.0),
413            WeightedVertex::new(3.0, 4.0, 0.0),
414        );
415        let e1 = WeightedEdge::new(
416            WeightedVertex::new(3.0, 4.0, 0.0),
417            WeightedVertex::new(0.0, 1.0, 0.0),
418        );
419        let e2 = WeightedEdge::new(
420            WeightedVertex::new(0.0, 0.0, 0.0),
421            WeightedVertex::new(1.0, 1.0, 0.0),
422        );
423        assert_eq!(e0, e1);
424        assert_ne!(e0, e2);
425        assert_ne!(e1, e2);
426    }
427
428    #[test]
429    fn triangles_equal() {
430        let t0 = WeightedTriangle::new(
431            WeightedVertex::new(0.0, 1.0, 0.0),
432            WeightedVertex::new(1.0, 0.0, 0.0),
433            WeightedVertex::new(-1.0, 0.0, 0.0),
434        );
435        let t1 = WeightedTriangle::new(
436            WeightedVertex::new(1.0, 0.0, 0.0),
437            WeightedVertex::new(-1.0, 0.0, 0.0),
438            WeightedVertex::new(0.0, 1.0, 0.0),
439        );
440        assert_eq!(t0, t1);
441    }
442
443    #[test]
444    fn triangle_and_point_regular() {
445        let vertex = WeightedVertex::new(0.6984975495419887, 0.6787937052516924, 0.0);
446
447        let tri = WeightedTriangle::new(
448            WeightedVertex::new(0.5222943585280901, 0.9101950968457111, 0.0),
449            WeightedVertex::new(0.7504630819229956, 0.42525042963771664, 0.0),
450            WeightedVertex::new(0.9030435560273877, 0.6666048060052796, 0.0),
451        );
452
453        assert!(!tri.is_regular(&vertex));
454        assert!(tri.regularity(&vertex) - 0.1 < 0.0);
455    }
456
457    #[test]
458    fn neighbors() {
459        let t0 = WeightedTriangle::new(
460            WeightedVertex::new(2.0, 0.0, 0.0),
461            WeightedVertex::new(2.8, 2.8, 0.0),
462            WeightedVertex::new(2.4, 2.4, 0.0),
463        );
464
465        let t1 = WeightedTriangle::new(
466            WeightedVertex::new(3.0, 2.0, 0.0),
467            WeightedVertex::new(2.8, 2.8, 0.0),
468            WeightedVertex::new(2.4, 2.4, 0.0),
469        );
470
471        assert_eq!(
472            t0.is_neighbor(&t1),
473            Some(WeightedEdge::new(
474                WeightedVertex::new(2.8, 2.8, 0.0),
475                WeightedVertex::new(2.4, 2.4, 0.0)
476            ))
477        );
478        assert_eq!(
479            t1.is_neighbor(&t0),
480            Some(WeightedEdge::new(
481                WeightedVertex::new(2.8, 2.8, 0.0),
482                WeightedVertex::new(2.4, 2.4, 0.0)
483            ))
484        );
485    }
486}
cg_math/geometry.rs

cg_math/
geometry.rs
