oxiphysics_geometry/voronoi/
types.rs1use super::functions::{bowyer_watson, is_delaunay};
6use super::functions::{
7 circumcircle, circumsphere_3d, clip_polygon_to_bounds, delaunay_to_voronoi,
8 flip_edge_if_needed, lloyd_relaxation, nearest_site, polygon_area_pt2d, polygon_centroid_pt2d,
9};
10
11#[derive(Debug, Clone)]
13pub struct DelaunayTriangulation {
14 pub points: Vec<Point2D>,
16 pub triangles: Vec<[usize; 3]>,
18}
19impl DelaunayTriangulation {
20 pub fn from_points(pts: &[[f64; 2]]) -> Self {
23 let points: Vec<Point2D> = pts.iter().map(|&p| Point2D::from_array(p)).collect();
24 let raw = bowyer_watson(pts);
25 let triangles: Vec<[usize; 3]> = raw.iter().map(|t| [t.a, t.b, t.c]).collect();
26 Self { points, triangles }
27 }
28 pub fn circumscribed_circle(&self, tri: [usize; 3]) -> ([f64; 2], f64) {
32 let pa = self.points[tri[0]].to_array();
33 let pb = self.points[tri[1]].to_array();
34 let pc = self.points[tri[2]].to_array();
35 match circumcircle(pa, pb, pc) {
36 Some((center, r2)) => (center, r2.sqrt()),
37 None => ([0.0, 0.0], 0.0),
38 }
39 }
40 pub fn is_delaunay(&self) -> bool {
42 let pts: Vec<[f64; 2]> = self.points.iter().map(|p| p.to_array()).collect();
43 let tris: Vec<DelaunayTriangle> = self
44 .triangles
45 .iter()
46 .map(|&t| DelaunayTriangle::new(t[0], t[1], t[2]))
47 .collect();
48 is_delaunay(&pts, &tris)
49 }
50 pub fn flip_edge(&mut self) {
54 let pts: Vec<[f64; 2]> = self.points.iter().map(|p| p.to_array()).collect();
55 let mut tris: Vec<DelaunayTriangle> = self
56 .triangles
57 .iter()
58 .map(|&t| DelaunayTriangle::new(t[0], t[1], t[2]))
59 .collect();
60 let n = tris.len();
61 for i in 0..n {
62 for j in (i + 1)..n {
63 if j < tris.len() {
64 flip_edge_if_needed(&mut tris, &pts, i, j);
65 }
66 }
67 }
68 self.triangles = tris.iter().map(|t| [t.a, t.b, t.c]).collect();
69 }
70 pub fn dual_voronoi(&self) -> VoronoiDiagram {
74 let pts: Vec<[f64; 2]> = self.points.iter().map(|p| p.to_array()).collect();
75 if pts.is_empty() {
76 return VoronoiDiagram {
77 sites: vec![],
78 cells: vec![],
79 bbox: [0.0, 1.0, 0.0, 1.0],
80 };
81 }
82 let (mut xmin, mut xmax, mut ymin, mut ymax) = (f64::MAX, f64::MIN, f64::MAX, f64::MIN);
83 for &p in &pts {
84 if p[0] < xmin {
85 xmin = p[0];
86 }
87 if p[0] > xmax {
88 xmax = p[0];
89 }
90 if p[1] < ymin {
91 ymin = p[1];
92 }
93 if p[1] > ymax {
94 ymax = p[1];
95 }
96 }
97 let margin_x = (xmax - xmin).max(1.0) * 0.2;
98 let margin_y = (ymax - ymin).max(1.0) * 0.2;
99 let bbox = [
100 xmin - margin_x,
101 xmax + margin_x,
102 ymin - margin_y,
103 ymax + margin_y,
104 ];
105 VoronoiDiagram::from_sites(&pts, bbox)
106 }
107}
108#[derive(Debug, Clone)]
112pub struct VoronoiDiagram {
113 pub sites: Vec<VoronoiSite>,
115 pub cells: Vec<VoronoiCell>,
117 pub bbox: [f64; 4],
119}
120impl VoronoiDiagram {
121 pub fn from_sites(pts: &[[f64; 2]], bbox: [f64; 4]) -> Self {
130 let sites: Vec<VoronoiSite> = pts
131 .iter()
132 .enumerate()
133 .map(|(i, &p)| VoronoiSite::new(Point2D::from_array(p), i))
134 .collect();
135 if pts.is_empty() {
136 return Self {
137 sites,
138 cells: vec![],
139 bbox,
140 };
141 }
142 let triangles = bowyer_watson(pts);
143 let legacy_cells = delaunay_to_voronoi(pts, &triangles);
144 let cells: Vec<VoronoiCell> = legacy_cells
145 .iter()
146 .enumerate()
147 .map(|(i, lc)| {
148 let clipped = clip_polygon_to_bounds(&lc.vertices, bbox);
149 let verts: Vec<Point2D> = clipped.iter().map(|&v| Point2D::from_array(v)).collect();
150 VoronoiCell::new(i, verts)
151 })
152 .collect();
153 Self { sites, cells, bbox }
154 }
155 pub fn nearest_site(&self, p: [f64; 2]) -> usize {
157 if self.sites.is_empty() {
158 return 0;
159 }
160 let pts: Vec<[f64; 2]> = self.sites.iter().map(|s| s.pos.to_array()).collect();
161 nearest_site(p, &pts)
162 }
163 pub fn cell_area(&self, i: usize) -> f64 {
165 if i < self.cells.len() {
166 self.cells[i].area
167 } else {
168 0.0
169 }
170 }
171 pub fn centroid(&self, i: usize) -> [f64; 2] {
173 if i < self.cells.len() {
174 let c = self.cells[i].centroid();
175 [c.x, c.y]
176 } else {
177 [0.0, 0.0]
178 }
179 }
180 pub fn lloyd_relaxation(&self, n_iter: usize) -> Self {
183 let pts_arr: Vec<[f64; 2]> = self.sites.iter().map(|s| s.pos.to_array()).collect();
184 let relaxed = lloyd_relaxation(&pts_arr, self.bbox, n_iter);
185 Self::from_sites(&relaxed, self.bbox)
186 }
187}
188#[derive(Debug, Clone, Copy, PartialEq)]
190pub struct Point2D {
191 pub x: f64,
193 pub y: f64,
195}
196impl Point2D {
197 pub fn new(x: f64, y: f64) -> Self {
199 Self { x, y }
200 }
201 #[inline]
203 pub fn dist2(&self, other: &Point2D) -> f64 {
204 (self.x - other.x).powi(2) + (self.y - other.y).powi(2)
205 }
206 #[inline]
208 pub fn dist(&self, other: &Point2D) -> f64 {
209 self.dist2(other).sqrt()
210 }
211 #[inline]
213 pub fn to_array(&self) -> [f64; 2] {
214 [self.x, self.y]
215 }
216 #[inline]
218 pub fn from_array(a: [f64; 2]) -> Self {
219 Self { x: a[0], y: a[1] }
220 }
221}
222#[derive(Debug, Clone, Copy, PartialEq, Eq)]
224pub struct DelaunayTriangle {
225 pub a: usize,
227 pub b: usize,
229 pub c: usize,
231}
232impl DelaunayTriangle {
233 pub fn new(a: usize, b: usize, c: usize) -> Self {
235 Self { a, b, c }
236 }
237 pub fn indices(&self) -> [usize; 3] {
239 [self.a, self.b, self.c]
240 }
241 pub fn contains_super_vertex(&self, offset: usize) -> bool {
243 self.a >= offset || self.b >= offset || self.c >= offset
244 }
245}
246#[derive(Debug, Clone)]
250pub struct VoronoiCell {
251 pub site_index: usize,
253 pub vertices: Vec<Point2D>,
255 pub area: f64,
257}
258impl VoronoiCell {
259 pub fn new(site_index: usize, vertices: Vec<Point2D>) -> Self {
261 let area = polygon_area_pt2d(&vertices);
262 Self {
263 site_index,
264 vertices,
265 area,
266 }
267 }
268 pub fn centroid(&self) -> Point2D {
270 polygon_centroid_pt2d(&self.vertices)
271 }
272 pub fn perimeter(&self) -> f64 {
274 let n = self.vertices.len();
275 if n < 2 {
276 return 0.0;
277 }
278 (0..n)
279 .map(|i| self.vertices[i].dist(&self.vertices[(i + 1) % n]))
280 .sum()
281 }
282}
283#[derive(Debug, Clone)]
287pub struct LegacyVoronoiCell {
288 pub site: [f64; 2],
290 pub vertices: Vec<[f64; 2]>,
292}
293impl LegacyVoronoiCell {
294 pub fn new(site: [f64; 2], vertices: Vec<[f64; 2]>) -> Self {
296 Self { site, vertices }
297 }
298}
299#[derive(Debug, Clone, Copy, PartialEq, Eq)]
301pub struct DelaunayTetrahedron {
302 pub v: [usize; 4],
304}
305impl DelaunayTetrahedron {
306 pub fn new(a: usize, b: usize, c: usize, d: usize) -> Self {
308 Self { v: [a, b, c, d] }
309 }
310 pub fn contains_super(&self, offset: usize) -> bool {
312 self.v.iter().any(|&i| i >= offset)
313 }
314 pub fn circumsphere(&self, pts: &[[f64; 3]]) -> Option<([f64; 3], f64)> {
316 circumsphere_3d(
317 pts[self.v[0]],
318 pts[self.v[1]],
319 pts[self.v[2]],
320 pts[self.v[3]],
321 )
322 }
323}
324#[derive(Debug, Clone, Copy)]
329pub struct WeightedSite {
330 pub position: [f64; 2],
332 pub weight: f64,
334}
335impl WeightedSite {
336 pub fn new(position: [f64; 2], weight: f64) -> Self {
338 Self { position, weight }
339 }
340 pub fn power_distance(&self, q: [f64; 2]) -> f64 {
342 let dx = q[0] - self.position[0];
343 let dy = q[1] - self.position[1];
344 dx * dx + dy * dy - self.weight
345 }
346}
347#[derive(Debug, Clone, Copy, PartialEq)]
349pub struct VoronoiSite {
350 pub pos: Point2D,
352 pub index: usize,
354}
355impl VoronoiSite {
356 pub fn new(pos: Point2D, index: usize) -> Self {
358 Self { pos, index }
359 }
360}