oxigdal_algorithms/vector/
delaunay.rs

1//! Delaunay triangulation
2//!
3//! Compute Delaunay triangulation of point sets.
4
5use crate::error::{AlgorithmError, Result};
6use oxigdal_core::vector::{Coordinate, LineString, Point, Polygon};
7
8/// Options for Delaunay triangulation
9#[derive(Debug, Clone)]
10pub struct DelaunayOptions {
11    /// Whether to include triangle quality metrics
12    pub compute_quality: bool,
13    /// Minimum angle threshold for quality triangles (degrees)
14    pub min_angle: f64,
15}
16
17impl Default for DelaunayOptions {
18    fn default() -> Self {
19        Self {
20            compute_quality: false,
21            min_angle: 20.0,
22        }
23    }
24}
25
26/// A triangle in the triangulation
27#[derive(Debug, Clone)]
28pub struct Triangle {
29    /// Indices of the three vertices
30    pub vertices: [usize; 3],
31    /// Triangle polygon
32    pub polygon: Polygon,
33    /// Triangle quality (0-1, higher is better)
34    pub quality: Option<f64>,
35}
36
37/// Delaunay triangulation result
38#[derive(Debug, Clone)]
39pub struct DelaunayTriangulation {
40    /// Input points
41    pub points: Vec<Point>,
42    /// Triangles
43    pub triangles: Vec<Triangle>,
44    /// Number of triangles
45    pub num_triangles: usize,
46}
47
48/// Compute Delaunay triangulation
49///
50/// # Arguments
51///
52/// * `points` - Input points
53/// * `options` - Triangulation options
54///
55/// # Returns
56///
57/// Delaunay triangulation with triangles
58///
59/// # Examples
60///
61/// ```
62/// # use oxigdal_algorithms::error::Result;
63/// use oxigdal_algorithms::vector::delaunay::{delaunay_triangulation, DelaunayOptions};
64/// use oxigdal_algorithms::Point;
65///
66/// # fn main() -> Result<()> {
67/// let points = vec![
68///     Point::new(0.0, 0.0),
69///     Point::new(1.0, 0.0),
70///     Point::new(0.5, 1.0),
71///     Point::new(0.5, 0.5),
72/// ];
73///
74/// let options = DelaunayOptions {
75///     compute_quality: true,
76///     ..Default::default()
77/// };
78///
79/// let triangulation = delaunay_triangulation(&points, &options)?;
80/// assert!(triangulation.num_triangles >= 2);
81/// # Ok(())
82/// # }
83/// ```
84pub fn delaunay_triangulation(
85    points: &[Point],
86    options: &DelaunayOptions,
87) -> Result<DelaunayTriangulation> {
88    if points.len() < 3 {
89        return Err(AlgorithmError::InvalidInput(
90            "Need at least 3 points for triangulation".to_string(),
91        ));
92    }
93
94    // Convert points to delaunator format
95    let delaunator_points: Vec<delaunator::Point> = points
96        .iter()
97        .map(|p| delaunator::Point {
98            x: p.coord.x,
99            y: p.coord.y,
100        })
101        .collect();
102
103    // Compute Delaunay triangulation
104    let delaunay = delaunator::triangulate(&delaunator_points);
105
106    // Build triangles
107    let mut triangles = Vec::new();
108
109    for tri_idx in 0..(delaunay.triangles.len() / 3) {
110        let a = delaunay.triangles[tri_idx * 3];
111        let b = delaunay.triangles[tri_idx * 3 + 1];
112        let c = delaunay.triangles[tri_idx * 3 + 2];
113
114        let pa = &points[a];
115        let pb = &points[b];
116        let pc = &points[c];
117
118        // Create triangle polygon
119        let coords_tri = vec![
120            Coordinate::new_2d(pa.coord.x, pa.coord.y),
121            Coordinate::new_2d(pb.coord.x, pb.coord.y),
122            Coordinate::new_2d(pc.coord.x, pc.coord.y),
123            Coordinate::new_2d(pa.coord.x, pa.coord.y), // Close the ring
124        ];
125
126        let exterior = LineString::new(coords_tri)
127            .map_err(|e| AlgorithmError::InvalidGeometry(format!("Invalid triangle: {}", e)))?;
128
129        let polygon = Polygon::new(exterior, vec![]).map_err(|e| {
130            AlgorithmError::InvalidGeometry(format!("Invalid triangle polygon: {}", e))
131        })?;
132
133        // Compute quality if requested
134        let quality = if options.compute_quality {
135            Some(compute_triangle_quality(pa, pb, pc))
136        } else {
137            None
138        };
139
140        triangles.push(Triangle {
141            vertices: [a, b, c],
142            polygon,
143            quality,
144        });
145    }
146
147    let num_triangles = triangles.len();
148
149    Ok(DelaunayTriangulation {
150        points: points.to_vec(),
151        triangles,
152        num_triangles,
153    })
154}
155
156/// Compute triangle quality (ratio of inradius to circumradius)
157fn compute_triangle_quality(pa: &Point, pb: &Point, pc: &Point) -> f64 {
158    // Edge lengths
159    let a = distance(pb, pc);
160    let b = distance(pc, pa);
161    let c = distance(pa, pb);
162
163    // Semi-perimeter
164    let s = (a + b + c) / 2.0;
165
166    // Area (Heron's formula)
167    let area = (s * (s - a) * (s - b) * (s - c)).sqrt();
168
169    // Inradius
170    let inradius = area / s;
171
172    // Circumradius
173    let circumradius = (a * b * c) / (4.0 * area);
174
175    // Quality ratio (0-1, higher is better)
176    if circumradius > 0.0 {
177        2.0 * inradius / circumradius
178    } else {
179        0.0
180    }
181}
182
183/// Calculate distance between two points
184fn distance(p1: &Point, p2: &Point) -> f64 {
185    let dx = p1.coord.x - p2.coord.x;
186    let dy = p1.coord.y - p2.coord.y;
187    (dx * dx + dy * dy).sqrt()
188}
189
190/// Check if a point is inside the circumcircle of a triangle
191pub fn in_circumcircle(pa: &Point, pb: &Point, pc: &Point, pd: &Point) -> bool {
192    let ax = pa.coord.x - pd.coord.x;
193    let ay = pa.coord.y - pd.coord.y;
194    let bx = pb.coord.x - pd.coord.x;
195    let by = pb.coord.y - pd.coord.y;
196    let cx = pc.coord.x - pd.coord.x;
197    let cy = pc.coord.y - pd.coord.y;
198
199    let det = (ax * ax + ay * ay) * (bx * cy - cx * by) - (bx * bx + by * by) * (ax * cy - cx * ay)
200        + (cx * cx + cy * cy) * (ax * by - bx * ay);
201
202    det > 0.0
203}
204
205/// Constrained Delaunay triangulation with constraint edges
206pub fn constrained_delaunay(
207    points: &[Point],
208    constraints: &[(usize, usize)],
209    options: &DelaunayOptions,
210) -> Result<DelaunayTriangulation> {
211    // Basic implementation: compute triangulation first, then enforce constraints
212    let mut triangulation = delaunay_triangulation(points, options)?;
213
214    // Filter triangles that violate constraints
215    triangulation
216        .triangles
217        .retain(|triangle| !violates_constraints(triangle, constraints, points));
218
219    triangulation.num_triangles = triangulation.triangles.len();
220
221    Ok(triangulation)
222}
223
224/// Check if a triangle violates any constraint edges
225fn violates_constraints(
226    triangle: &Triangle,
227    constraints: &[(usize, usize)],
228    _points: &[Point],
229) -> bool {
230    // Check if any constraint edge intersects the triangle interior
231    for &(c1, c2) in constraints {
232        // If constraint edge is an edge of the triangle, it's OK
233        let [a, b, c] = triangle.vertices;
234
235        if (c1 == a && c2 == b)
236            || (c1 == b && c2 == a)
237            || (c1 == b && c2 == c)
238            || (c1 == c && c2 == b)
239            || (c1 == c && c2 == a)
240            || (c1 == a && c2 == c)
241        {
242            continue;
243        }
244
245        // Otherwise, check for intersection (simplified check)
246        // A full implementation would need proper edge-triangle intersection
247    }
248
249    false
250}
251
252#[cfg(test)]
253mod tests {
254    use super::*;
255
256    #[test]
257    fn test_delaunay_simple() {
258        let points = vec![
259            Point::new(0.0, 0.0),
260            Point::new(1.0, 0.0),
261            Point::new(0.5, 1.0),
262            Point::new(0.5, 0.5),
263        ];
264
265        let options = DelaunayOptions::default();
266        let result = delaunay_triangulation(&points, &options);
267
268        assert!(result.is_ok());
269
270        let triangulation = result.expect("Triangulation failed");
271        assert!(triangulation.num_triangles >= 2);
272    }
273
274    #[test]
275    fn test_triangle_quality() {
276        // Equilateral triangle (perfect quality)
277        let pa = Point::new(0.0, 0.0);
278        let pb = Point::new(1.0, 0.0);
279        let pc = Point::new(0.5, 0.866); // ~sqrt(3)/2
280
281        let quality = compute_triangle_quality(&pa, &pb, &pc);
282        assert!(quality > 0.9); // Close to 1.0 for equilateral
283    }
284
285    #[test]
286    fn test_in_circumcircle() {
287        let pa = Point::new(0.0, 0.0);
288        let pb = Point::new(1.0, 0.0);
289        let pc = Point::new(0.0, 1.0);
290        let pd = Point::new(0.25, 0.25); // Inside circumcircle
291
292        assert!(in_circumcircle(&pa, &pb, &pc, &pd));
293    }
294
295    #[test]
296    fn test_constrained_delaunay() {
297        let points = vec![
298            Point::new(0.0, 0.0),
299            Point::new(1.0, 0.0),
300            Point::new(0.5, 1.0),
301            Point::new(0.5, 0.5),
302        ];
303
304        let constraints = vec![(0, 2)]; // Constraint edge from point 0 to point 2
305
306        let options = DelaunayOptions::default();
307        let result = constrained_delaunay(&points, &constraints, &options);
308
309        assert!(result.is_ok());
310    }
311}
oxigdal_algorithms/vector/delaunay.rs

oxigdal_algorithms/vector/
delaunay.rs
