egml_core/model/geometry/primitives/
polygon.rs

1use crate::model::common::{
2    ApplyTransform, ComputeEnvelope, IterGeometries, Triangulate, Triangulation,
3};
4use crate::model::geometry::primitives::{
5    AbstractRingKind, AbstractSurface, AsAbstractSurface, AsAbstractSurfaceMut,
6};
7use crate::model::geometry::refs::AbstractGeometryKindRef;
8use crate::model::geometry::{DirectPosition, Envelope};
9use crate::util::plane::Plane;
10use crate::util::triangulate::triangulate;
11use crate::{
12    Error, impl_abstract_surface_mut_traits, impl_abstract_surface_traits, impl_has_geometry_type,
13};
14use nalgebra::{Isometry3, Rotation3, Scale3, Transform3, Vector3};
15use rayon::prelude::*;
16
17#[derive(Debug, Clone, PartialEq)]
18pub struct Polygon {
19    pub abstract_surface: AbstractSurface,
20    exterior: Option<AbstractRingKind>,
21    interior: Vec<AbstractRingKind>,
22}
23
24impl Polygon {
25    pub fn new(
26        exterior: Option<AbstractRingKind>,
27        interior: impl IntoIterator<Item = AbstractRingKind>,
28    ) -> Result<Self, Error> {
29        Ok(Self {
30            abstract_surface: AbstractSurface::default(),
31            exterior,
32            interior: interior.into_iter().collect(),
33        })
34    }
35
36    pub fn from_abstract_surface(
37        abstract_surface: AbstractSurface,
38        exterior: Option<AbstractRingKind>,
39        interior: impl IntoIterator<Item = AbstractRingKind>,
40    ) -> Self {
41        Self {
42            abstract_surface,
43            exterior,
44            interior: interior.into_iter().collect(),
45        }
46    }
47
48    pub fn exterior(&self) -> Option<&AbstractRingKind> {
49        self.exterior.as_ref()
50    }
51
52    pub fn set_exterior(&mut self, exterior: AbstractRingKind) {
53        self.exterior = Some(exterior);
54    }
55
56    pub fn set_exterior_opt(&mut self, exterior: Option<AbstractRingKind>) {
57        self.exterior = exterior;
58    }
59
60    pub fn clear_exterior(&mut self) {
61        self.exterior = None;
62    }
63
64    pub fn interior(&self) -> &[AbstractRingKind] {
65        &self.interior
66    }
67
68    pub fn set_interior(&mut self, interior: Vec<AbstractRingKind>) {
69        self.interior = interior;
70    }
71
72    pub fn push_interior(&mut self, ring: AbstractRingKind) {
73        self.interior.push(ring);
74    }
75
76    pub fn extend_interiors(&mut self, rings: impl IntoIterator<Item = AbstractRingKind>) {
77        self.interior.extend(rings);
78    }
79}
80
81impl AsAbstractSurface for Polygon {
82    fn abstract_surface(&self) -> &AbstractSurface {
83        &self.abstract_surface
84    }
85}
86
87impl AsAbstractSurfaceMut for Polygon {
88    fn abstract_surface_mut(&mut self) -> &mut AbstractSurface {
89        &mut self.abstract_surface
90    }
91}
92
93impl_abstract_surface_traits!(Polygon);
94impl_abstract_surface_mut_traits!(Polygon);
95impl_has_geometry_type!(Polygon, Polygon);
96
97impl Polygon {
98    ///
99    /// See also <https://www.khronos.org/opengl/wiki/Calculating_a_Surface_Normal#Newell.27s_Method>
100    fn normal(&self) -> Vector3<f64> {
101        let mut enclosed_boundary_points =
102            self.exterior().expect("should be there").points().to_vec();
103        let first = enclosed_boundary_points
104            .first()
105            .copied()
106            .expect("should be there");
107        enclosed_boundary_points.push(first);
108
109        let mut normal = Vector3::new(0.0, 0.0, 0.0);
110        for current_point_pair in enclosed_boundary_points.windows(2) {
111            let current_first_point: Vector3<f64> = current_point_pair[0].into();
112            let current_second_point: Vector3<f64> = current_point_pair[1].into();
113
114            normal += (current_first_point - current_second_point)
115                .cross(&(current_first_point + current_second_point));
116        }
117
118        normal.normalize()
119    }
120
121    pub fn plane_equation(&self) -> Plane {
122        let envelope = self.compute_envelope().expect("should have envelope");
123        Plane::new(*envelope.lower_corner(), self.normal())
124    }
125
126    /// Returns the net 3D area_3d of this polygon: exterior area_3d minus the sum of all interior hole area_3ds.
127    ///
128    /// # Errors
129    ///
130    /// Returns [`Error::MissingExteriorRing`] if the polygon has no exterior ring.
131    pub fn area_3d(&self) -> Result<f64, Error> {
132        let exterior = self
133            .exterior
134            .as_ref()
135            .ok_or(Error::MissingExteriorRing)?
136            .area_3d();
137        let holes = self.interior.iter().map(|r| r.area_3d()).sum::<f64>();
138
139        Ok(exterior - holes)
140    }
141
142    pub fn points(&self) -> Vec<&DirectPosition> {
143        let mut all_points = Vec::new();
144        if let Some(exterior) = &self.exterior {
145            all_points.extend(exterior.points());
146        }
147        for ring in &self.interior {
148            all_points.extend(ring.points());
149        }
150
151        all_points
152    }
153}
154
155impl ApplyTransform for Polygon {
156    fn apply_transform(&mut self, transform: Transform3<f64>) {
157        if let Some(exterior) = &mut self.exterior {
158            exterior.apply_transform(transform);
159        }
160        self.interior
161            .par_iter_mut()
162            .for_each(|p| p.apply_transform(transform));
163    }
164
165    fn apply_isometry(&mut self, isometry: Isometry3<f64>) {
166        if let Some(exterior) = &mut self.exterior {
167            exterior.apply_isometry(isometry);
168        }
169        self.interior
170            .par_iter_mut()
171            .for_each(|p| p.apply_isometry(isometry));
172    }
173
174    fn apply_translation(&mut self, vector: Vector3<f64>) {
175        if let Some(exterior) = &mut self.exterior {
176            exterior.apply_translation(vector);
177        }
178        self.interior
179            .par_iter_mut()
180            .for_each(|p| p.apply_translation(vector));
181    }
182
183    fn apply_rotation(&mut self, rotation: Rotation3<f64>) {
184        if let Some(exterior) = &mut self.exterior {
185            exterior.apply_rotation(rotation);
186        }
187        self.interior
188            .par_iter_mut()
189            .for_each(|p| p.apply_rotation(rotation));
190    }
191
192    fn apply_scale(&mut self, scale: Scale3<f64>) {
193        if let Some(exterior) = &mut self.exterior {
194            exterior.apply_scale(scale);
195        }
196        self.interior
197            .par_iter_mut()
198            .for_each(|p| p.apply_scale(scale));
199    }
200}
201
202impl ComputeEnvelope for Polygon {
203    fn compute_envelope(&self) -> Option<Envelope> {
204        if let Some(exterior) = &self.exterior
205            && let Some(e) = exterior.compute_envelope()
206        {
207            return Some(e);
208        }
209
210        let envelopes = self
211            .interior
212            .iter()
213            .filter_map(|x| x.compute_envelope())
214            .collect::<Vec<_>>();
215
216        Envelope::from_envelopes(&envelopes)
217    }
218}
219
220impl Triangulate for Polygon {
221    fn triangulate(&self) -> Result<Triangulation, Error> {
222        let surface = triangulate(self.exterior.clone(), self.interior.to_vec())?;
223        Ok(Triangulation::new(surface, Vec::new()))
224    }
225}
226
227impl IterGeometries for Polygon {
228    fn iter_geometries(&self) -> Box<dyn Iterator<Item = AbstractGeometryKindRef<'_>> + '_> {
229        Box::new(
230            std::iter::once(self.into())
231                .chain(
232                    self.exterior
233                        .as_ref()
234                        .into_iter()
235                        .flat_map(|x| x.iter_geometries()),
236                )
237                .chain(self.interior.iter().flat_map(|x| x.iter_geometries())),
238        )
239    }
240}
241
242#[cfg(test)]
243mod test {
244    use super::*;
245    use crate::model::geometry::DirectPosition;
246    use crate::model::geometry::primitives::{AbstractRingKind, AsSurface, LinearRing};
247    use nalgebra::Vector3;
248
249    #[test]
250    fn area_3d_unit_square() {
251        let ring = LinearRing::new([
252            DirectPosition::new(0.0, 0.0, 1.0).unwrap(),
253            DirectPosition::new(1.0, 0.0, 1.0).unwrap(),
254            DirectPosition::new(1.0, 1.0, 1.0).unwrap(),
255            DirectPosition::new(0.0, 1.0, 1.0).unwrap(),
256        ])
257        .unwrap();
258        let polygon = Polygon::new(Some(AbstractRingKind::LinearRing(ring)), []).unwrap();
259        assert!((polygon.area_3d().expect("has exterior ring") - 1.0).abs() < 1e-10);
260    }
261
262    #[test]
263    fn area_3d_with_hole() {
264        // 4×4 outer square with a 1×1 hole — net area_3d should be 15.
265        let exterior = LinearRing::new([
266            DirectPosition::new(0.0, 0.0, 0.0).unwrap(),
267            DirectPosition::new(4.0, 0.0, 0.0).unwrap(),
268            DirectPosition::new(4.0, 4.0, 0.0).unwrap(),
269            DirectPosition::new(0.0, 4.0, 0.0).unwrap(),
270        ])
271        .unwrap();
272        let hole = LinearRing::new([
273            DirectPosition::new(1.0, 1.0, 0.0).unwrap(),
274            DirectPosition::new(2.0, 1.0, 0.0).unwrap(),
275            DirectPosition::new(2.0, 2.0, 0.0).unwrap(),
276            DirectPosition::new(1.0, 2.0, 0.0).unwrap(),
277        ])
278        .unwrap();
279        let polygon = Polygon::new(
280            Some(AbstractRingKind::LinearRing(exterior)),
281            vec![AbstractRingKind::LinearRing(hole)],
282        )
283        .unwrap();
284        assert!((polygon.area_3d().expect("has exterior ring") - 15.0).abs() < 1e-10);
285    }
286
287    #[test]
288    fn area_3d_no_exterior_ring() {
289        let polygon = Polygon::new(None, []).unwrap();
290        assert_eq!(polygon.area_3d(), Err(Error::MissingExteriorRing));
291    }
292
293    #[test]
294    fn triangulate_no_exterior_ring() {
295        let polygon = Polygon::new(None, []).unwrap();
296        assert_eq!(
297            polygon.triangulate().err(),
298            Some(Error::MissingExteriorRing)
299        );
300    }
301
302    #[test]
303    fn basic_normal_vector() {
304        let point_a = DirectPosition::new(0.0, 0.0, 1.0).unwrap();
305        let point_b = DirectPosition::new(1.0, 0.0, 1.0).unwrap();
306        let point_c = DirectPosition::new(1.0, 1.0, 1.0).unwrap();
307        let point_d = DirectPosition::new(0.0, 1.0, 1.0).unwrap();
308        let linear_ring = LinearRing::new([point_a, point_b, point_c, point_d]).unwrap();
309        let linear_ring = AbstractRingKind::LinearRing(linear_ring);
310        let polygon = Polygon::new(Some(linear_ring), []).unwrap();
311        let normal = polygon.normal();
312
313        assert_eq!(normal, Vector3::new(0.0, 0.0, 1.0));
314    }
315
316    #[test]
317    fn basic_plane_equation() {
318        let point_a = DirectPosition::new(0.0, 0.0, 1.0).unwrap();
319        let point_b = DirectPosition::new(1.0, 0.0, 1.0).unwrap();
320        let point_c = DirectPosition::new(1.0, 1.0, 1.0).unwrap();
321        let point_d = DirectPosition::new(0.0, 1.0, 1.0).unwrap();
322        let linear_ring = LinearRing::new([point_a, point_b, point_c, point_d]).unwrap();
323        let linear_ring = AbstractRingKind::LinearRing(linear_ring);
324        let polygon = Polygon::new(Some(linear_ring), []).unwrap();
325        let plane_equation = polygon.plane_equation();
326
327        assert_eq!(
328            plane_equation.point,
329            DirectPosition::new(0.0, 0.0, 1.0).unwrap()
330        );
331        assert_eq!(plane_equation.normal(), Vector3::new(0.0, 0.0, 1.0));
332    }
333
334    #[test]
335    fn test_polygon_triangulation() {
336        let linear_ring_exterior = LinearRing::new([
337            DirectPosition::new(0.0, 0.0, 0.0).expect("should work"),
338            DirectPosition::new(1.0, 0.0, 0.0).expect("should work"),
339            DirectPosition::new(1.0, 1.0, 2.0).expect("should work"),
340            DirectPosition::new(0.0, 1.0, 2.0).expect("should work"),
341        ])
342        .expect("should work");
343        let linear_ring_exterior = AbstractRingKind::LinearRing(linear_ring_exterior);
344
345        let polygon = Polygon::new(Some(linear_ring_exterior), vec![]).expect("should work");
346        let triangulation = polygon.triangulate().expect("should work");
347        assert_eq!(triangulation.surface().patches_len(), 2);
348    }
349
350    #[test]
351    fn test_polygon_with_interior_triangulation() {
352        // 2×2 rectangle in the tilted plane z = 2y (true size 2 × 2√5) with a
353        // 1×1 hole (true size 1 × √5) — net area 4√5 − √5 = 3√5.
354        let linear_ring_exterior = LinearRing::new([
355            DirectPosition::new(0.0, 0.0, 0.0).expect("should work"),
356            DirectPosition::new(2.0, 0.0, 0.0).expect("should work"),
357            DirectPosition::new(2.0, 2.0, 4.0).expect("should work"),
358            DirectPosition::new(0.0, 2.0, 4.0).expect("should work"),
359        ])
360        .expect("should work");
361        let linear_ring_interior = LinearRing::new([
362            DirectPosition::new(0.5, 0.5, 1.0).expect("should work"),
363            DirectPosition::new(0.5, 1.5, 3.0).expect("should work"),
364            DirectPosition::new(1.5, 1.5, 3.0).expect("should work"),
365            DirectPosition::new(1.5, 0.5, 1.0).expect("should work"),
366        ])
367        .expect("should work");
368
369        let polygon = Polygon::new(
370            Some(AbstractRingKind::LinearRing(linear_ring_exterior)),
371            vec![AbstractRingKind::LinearRing(linear_ring_interior)],
372        )
373        .expect("should work");
374        let triangulation = polygon.triangulate().expect("should work");
375
376        let expected = 3.0 * 5.0_f64.sqrt();
377        assert!(triangulation.skipped().is_empty());
378        assert!((triangulation.surface().area_3d().expect("should work") - expected).abs() < 1e-10);
379        assert!((polygon.area_3d().expect("should work") - expected).abs() < 1e-10);
380    }
381}
egml_core/model/geometry/primitives/polygon.rs

egml_core/model/geometry/primitives/
polygon.rs
