Skip to main content

egml_core/model/geometry/primitives/
polygon.rs

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