honeycomb_kernels/grisubal/routines/
pre_processing.rs

1//! Step 0 implementation
2
3use std::collections::HashSet;
4
5use honeycomb_core::geometry::{CoordsFloat, Vertex2};
6
7use crate::grisubal::{
8    GrisubalError,
9    model::{Geometry2, GridCellId},
10};
11
12/// Check for orientation issue **per boundary**.
13///
14/// This function check for the most obvious orientation issue; given a boundary, are all segments making it up
15/// oriented consistently. If it is not the case, then there is at least one of:
16///
17/// - a vertex being the origin of two segment
18/// - a vertex being the end-point of two segment
19///
20/// This does not cover consistent orientation across distinct boundaries (e.g. a geometry with a hole in it).
21pub fn detect_orientation_issue<T: CoordsFloat>(
22    geometry: &Geometry2<T>,
23) -> Result<(), GrisubalError> {
24    let mut origins = HashSet::new();
25    let mut endpoints = HashSet::new();
26
27    for (orig, endp) in &geometry.segments {
28        if !origins.insert(orig) || !endpoints.insert(endp) {
29            return Err(GrisubalError::InconsistentOrientation(
30                "in-boundary inconsistency",
31            ));
32        }
33    }
34
35    Ok(())
36}
37
38#[allow(clippy::cast_precision_loss)]
39pub fn compute_overlapping_grid<T: CoordsFloat>(
40    geometry: &Geometry2<T>,
41    [len_cell_x, len_cell_y]: [T; 2],
42) -> Result<([usize; 2], Vertex2<T>), GrisubalError> {
43    // compute the minimum bounding box
44    let (mut min_x, mut max_x, mut min_y, mut max_y): (T, T, T, T) = {
45        let Some(tmp) = geometry.vertices.first() else {
46            return Err(GrisubalError::InvalidShape("no vertex in shape"));
47        };
48
49        (tmp.x(), tmp.x(), tmp.y(), tmp.y())
50    };
51
52    geometry.vertices.iter().for_each(|v| {
53        min_x = min_x.min(v.x());
54        max_x = max_x.max(v.x()); // may not be optimal
55        min_y = min_y.min(v.y()); // don't care
56        max_y = max_y.max(v.y());
57    });
58
59    if max_x <= min_x {
60        return Err(GrisubalError::InvalidShape(
61            "bounding values along X axis are equal",
62        ));
63    }
64    if max_y <= min_y {
65        return Err(GrisubalError::InvalidShape(
66            "bounding values along Y axis are equal",
67        ));
68    }
69
70    // compute characteristics of the overlapping Cartesian grid
71
72    // create a ~one-and-a-half cell buffer to contain the geometry
73    // this, along with the `+1` below, guarantees that
74    // dart at the boundary of the grid are not intersected by the geometry
75    let mut og_x = min_x - len_cell_x * T::from(1.5).unwrap();
76    let mut og_y = min_y - len_cell_y * T::from(1.5).unwrap();
77    // we check for some extremely annoying cases here
78    // if some are detected, the origin is incrementally shifted
79    let (mut on_corner, mut reflect) =
80        detect_overlaps(geometry, [len_cell_x, len_cell_y], Vertex2(og_x, og_y));
81    let mut i = 1;
82
83    while on_corner | reflect {
84        eprintln!(
85            "W: land on corner: {on_corner} - reflect on an axis: {reflect}, shifting origin"
86        );
87        og_x += len_cell_x * T::from(1. / (2_i32.pow(i + 1) as f32)).unwrap();
88        og_y += len_cell_y * T::from(1. / (2_i32.pow(i + 1) as f32)).unwrap();
89        (on_corner, reflect) =
90            detect_overlaps(geometry, [len_cell_x, len_cell_y], Vertex2(og_x, og_y));
91        i += 1;
92    }
93
94    let n_cells_x = ((max_x - og_x) / len_cell_x).ceil().to_usize().unwrap() + 1;
95    let n_cells_y = ((max_y - og_y) / len_cell_y).ceil().to_usize().unwrap() + 1;
96
97    Ok(([n_cells_x, n_cells_y], Vertex2(og_x, og_y)))
98}
99
100/// Remove from their geometry points of interest that intersect with a grid of specified dimension.
101///
102/// This function works under the assumption that the grid is Cartesian & has its origin on `(0.0, 0.0)`.
103pub fn remove_redundant_poi<T: CoordsFloat>(
104    geometry: &mut Geometry2<T>,
105    [cx, cy]: [T; 2],
106    origin: Vertex2<T>,
107) {
108    // PoI that land on the grid create a number of issues; removing them is ok since we're intersecting the grid
109    // at their coordinates, so the shape will be captured via intersection anyway
110    geometry.poi.retain(|idx| {
111        let v = geometry.vertices[*idx];
112        // origin is assumed to be (0.0, 0.0)
113        let on_x_axis = ((v.x() - origin.x()) % cx).is_zero();
114        let on_y_axis = ((v.y() - origin.y()) % cy).is_zero();
115        !(on_x_axis | on_y_axis)
116    });
117}
118
119#[allow(clippy::map_all_any_identity)]
120pub fn detect_overlaps<T: CoordsFloat>(
121    geometry: &Geometry2<T>,
122    [cx, cy]: [T; 2],
123    origin: Vertex2<T>,
124) -> (bool, bool) {
125    let on_corner = geometry
126        .vertices
127        .iter()
128        .map(|v| {
129            let on_x_axis = ((v.x() - origin.x()) % cx).is_zero();
130            let on_y_axis = ((v.y() - origin.y()) % cy).is_zero();
131            on_x_axis && on_y_axis
132        })
133        .any(|a| a);
134
135    let bad_reflection = geometry
136        .vertices
137        .iter()
138        .enumerate()
139        .filter_map(|(id, v)| {
140            let on_x_axis = ((v.x() - origin.x()) % cx).is_zero();
141            let on_y_axis = ((v.y() - origin.y()) % cy).is_zero();
142            if on_x_axis | on_y_axis {
143                return Some(id);
144            }
145            None
146        })
147        // skip vertices that do not belong to the boundary
148        .filter(|id| {
149            geometry
150                .segments
151                .iter()
152                .any(|(v1, v2)| (id == v1) || (id == v2))
153        })
154        .map(|id| {
155            // if a vertex appear in the boundary, there should be both a segment landing and a
156            // segment starting on the vertex; hence `.expect()`
157            let vid_in = geometry
158                .segments
159                .iter()
160                .find_map(|(vin, ref_id)| {
161                    if id == *ref_id {
162                        return Some(*vin);
163                    }
164                    None
165                })
166                .expect("E: found a vertex with no incident segment - is the geometry open?");
167            // same
168            let vid_out = geometry
169                .segments
170                .iter()
171                .find_map(|(ref_id, vout)| {
172                    if id == *ref_id {
173                        return Some(*vout);
174                    }
175                    None
176                })
177                .expect("E: found a vertex with no incident segment - is the geometry open?");
178            let v_in = geometry.vertices[vid_in];
179            let v_out = geometry.vertices[vid_out];
180            let Vertex2(ox, oy) = origin;
181            let (c_in, c_out) = (
182                GridCellId(
183                    ((v_in.x() - ox) / cx).floor().to_usize().unwrap(),
184                    ((v_in.y() - oy) / cy).floor().to_usize().unwrap(),
185                ),
186                GridCellId(
187                    ((v_out.x() - ox) / cx).floor().to_usize().unwrap(),
188                    ((v_out.y() - oy) / cy).floor().to_usize().unwrap(),
189                ),
190            );
191            // if v_in and v_out belong to the same grid cell, there was a "reflection" on one
192            // of the grid's axis
193            c_in == c_out
194        })
195        .any(|a| a);
196
197    (on_corner, bad_reflection)
198}
honeycomb_kernels/grisubal/routines/pre_processing.rs

honeycomb_kernels/grisubal/routines/
pre_processing.rs
