gistools/geometry/tools/
clip.rs

1use alloc::{vec, vec::Vec};
2use s2json::{
3    Axis, BBox3D, VectorFeature, VectorGeometry, VectorGeometryType, VectorLineString,
4    VectorLineStringGeometry, VectorMultiLineOffset, VectorMultiLineString,
5    VectorMultiLineStringGeometry, VectorMultiPointGeometry, VectorMultiPolygon,
6    VectorMultiPolygonGeometry, VectorMultiPolygonOffset, VectorPoint, VectorPointGeometry,
7    VectorPolygonGeometry,
8};
9
10// TODO: Cases of `to_vec` clones large swathes of data. Can we optimize this?
11
12/// Internal clip function for a collection of VectorFeatures
13pub fn clip_features<M, P: Clone + Default, D: Clone + Default>(
14    features: &[VectorFeature<M, P, D>],
15    scale: f64,
16    k1: f64,
17    k2: f64,
18    axis: Axis,
19    base_buffer: f64,
20) -> Option<Vec<VectorFeature<M, P, D>>>
21where
22    M: Clone,
23{
24    // scale
25    let k1 = k1 / scale;
26    let k2 = k2 / scale;
27    // prep buffer and result container
28    let buffer = base_buffer / scale;
29    let k1b = k1 - buffer;
30    let k2b = k2 + buffer;
31    let mut clipped: Vec<VectorFeature<M, P, D>> = vec![];
32    let axis_x = axis == Axis::X;
33
34    for feature in features {
35        let geometry = &feature.geometry;
36        // trivial accept and reject cases
37        if let Some(vec_bbox) = geometry.vec_bbox() {
38            let min = if axis_x { vec_bbox.left } else { vec_bbox.bottom };
39            let max = if axis_x { vec_bbox.right } else { vec_bbox.top };
40            if min >= k1 && max < k2 {
41                clipped.push(feature.clone());
42                continue;
43            } else if max < k1 || min >= k2 {
44                continue;
45            }
46        }
47        // build the new clipped geometry
48        let new_geometry: Option<VectorGeometry<D>> = match geometry {
49            VectorGeometry::Point(geo) => clip_point(geo, axis, k1, k2),
50            VectorGeometry::MultiPoint(geo) => clip_multi_point(geo, axis, k1, k2),
51            VectorGeometry::LineString(geo) => clip_line_string(geo, axis, k1b, k2b),
52            VectorGeometry::MultiLineString(geo) => {
53                clip_multi_line_string(geo, axis, k1b, k2b, false)
54            }
55            VectorGeometry::Polygon(geo) => clip_polygon(geo, axis, k1b, k2b),
56            VectorGeometry::MultiPolygon(geo) => clip_multi_polygon(geo, axis, k1b, k2b),
57        };
58        // store if the geometry was inside the range
59        if let Some(new_geometry) = new_geometry {
60            clipped.push(VectorFeature::from_vector_feature(feature, Some(new_geometry)));
61        }
62    }
63
64    if clipped.is_empty() { None } else { Some(clipped) }
65}
66
67/// Clip a point to an axis and range
68pub fn clip_point<M: Clone + Default>(
69    geometry: &VectorPointGeometry<M>,
70    axis: Axis,
71    k1: f64,
72    k2: f64,
73) -> Option<VectorGeometry<M>> {
74    let coords = &geometry.coordinates;
75    let value = if axis == Axis::X { coords.x } else { coords.y };
76    if value >= k1 && value < k2 { Some(VectorGeometry::Point(geometry.clone())) } else { None }
77}
78
79/// Clip a MultiPoint to an axis and range
80pub fn clip_multi_point<M: Clone + Default>(
81    geometry: &VectorMultiPointGeometry<M>,
82    axis: Axis,
83    k1: f64,
84    k2: f64,
85) -> Option<VectorGeometry<M>> {
86    let mut new_geo = geometry.clone();
87    new_geo.coordinates = geometry
88        .coordinates
89        .iter()
90        .filter(|point| {
91            let value = if axis == Axis::X { point.x } else { point.y };
92            value >= k1 && value < k2
93        })
94        .cloned()
95        .collect();
96
97    if new_geo.coordinates.is_empty() { None } else { Some(VectorGeometry::MultiPoint(new_geo)) }
98}
99
100/// Clip a LineString to an axis and range
101pub fn clip_line_string<M: Clone + Default>(
102    geometry: &VectorLineStringGeometry<M>,
103    axis: Axis,
104    k1: f64,
105    k2: f64,
106) -> Option<VectorGeometry<M>> {
107    let VectorLineStringGeometry { is_3d, coordinates: line, bbox, vec_bbox, .. } = geometry;
108    let init_o = geometry.offset.unwrap_or(0.);
109    let mut new_offsets: VectorMultiLineOffset = vec![];
110    let mut new_lines: VectorMultiLineString<M> = vec![];
111    for clip in
112        _clip_line(ClipLineResult { line: line.to_vec(), offset: init_o }, k1, k2, axis, false)
113    {
114        new_offsets.push(clip.offset);
115        new_lines.push(clip.line);
116    }
117    if new_lines.is_empty() {
118        None
119    } else {
120        Some(VectorGeometry::MultiLineString(VectorMultiLineStringGeometry {
121            _type: VectorGeometryType::MultiLineString,
122            is_3d: *is_3d,
123            coordinates: new_lines,
124            bbox: *bbox,
125            offset: Some(new_offsets),
126            vec_bbox: Some(vec_bbox.unwrap_or_default().clip(axis, k1, k2)),
127            ..Default::default()
128        }))
129    }
130}
131
132/// Clip a MultiLineString geometry to an axis and range
133pub fn clip_multi_line_string<M: Clone + Default>(
134    geometry: &VectorMultiLineStringGeometry<M>,
135    axis: Axis,
136    k1: f64,
137    k2: f64,
138    is_polygon: bool,
139) -> Option<VectorGeometry<M>> {
140    let VectorMultiLineStringGeometry::<M> { is_3d, coordinates, bbox, vec_bbox, .. } = geometry;
141    let init_o =
142        geometry.offset.clone().unwrap_or_else(|| coordinates.iter().map(|_| 0.).collect());
143    let mut new_offsets: VectorMultiLineOffset = vec![];
144    let mut new_lines: VectorMultiLineString<M> = vec![];
145    let vec_bbox = vec_bbox.unwrap_or_default().clip(axis, k1, k2);
146    coordinates.iter().enumerate().for_each(|(i, line)| {
147        for clip in _clip_line(
148            ClipLineResult { line: line.to_vec(), offset: init_o[i] },
149            k1,
150            k2,
151            axis,
152            is_polygon,
153        ) {
154            new_offsets.push(clip.offset);
155            new_lines.push(clip.line);
156        }
157    });
158    if new_lines.is_empty() || (is_polygon && new_lines[0].len() < 4) {
159        None
160    } else if !is_polygon {
161        Some(VectorGeometry::MultiLineString(VectorMultiLineStringGeometry {
162            _type: VectorGeometryType::MultiLineString,
163            is_3d: *is_3d,
164            coordinates: new_lines,
165            bbox: *bbox,
166            offset: Some(new_offsets),
167            vec_bbox: Some(vec_bbox),
168            ..Default::default()
169        }))
170    } else {
171        Some(VectorGeometry::Polygon(VectorPolygonGeometry {
172            _type: VectorGeometryType::Polygon,
173            is_3d: *is_3d,
174            coordinates: new_lines,
175            bbox: *bbox,
176            offset: Some(new_offsets),
177            vec_bbox: Some(vec_bbox),
178            ..Default::default()
179        }))
180    }
181}
182
183/// Clip a Polygon geometry to an axis and range
184pub fn clip_polygon<M: Clone + Default>(
185    geometry: &VectorPolygonGeometry<M>,
186    axis: Axis,
187    k1: f64,
188    k2: f64,
189) -> Option<VectorGeometry<M>> {
190    clip_multi_line_string(geometry, axis, k1, k2, true)
191}
192
193/// Clip a MultiPolygon geometry to an axis and range
194pub fn clip_multi_polygon<M: Clone + Default>(
195    geometry: &VectorMultiPolygonGeometry<M>,
196    axis: Axis,
197    k1: f64,
198    k2: f64,
199) -> Option<VectorGeometry<M>> {
200    let VectorMultiPolygonGeometry { is_3d, coordinates, bbox, vec_bbox, .. } = geometry;
201    let init_o = geometry
202        .offset
203        .clone()
204        .unwrap_or_else(|| coordinates.iter().map(|l| l.iter().map(|_| 0.).collect()).collect());
205    let mut new_coordinates: VectorMultiPolygon<M> = vec![];
206    let mut new_offsets: VectorMultiPolygonOffset = vec![];
207    coordinates.iter().enumerate().for_each(|(p, polygon)| {
208        let new_polygon = clip_polygon(
209            &VectorPolygonGeometry {
210                _type: VectorGeometryType::Polygon,
211                is_3d: *is_3d,
212                coordinates: polygon.to_vec(),
213                offset: Some(init_o[p].clone()),
214                ..Default::default()
215            },
216            axis,
217            k1,
218            k2,
219        );
220        if let Some(VectorGeometry::Polygon(new_polygon)) = new_polygon {
221            new_coordinates.push(new_polygon.coordinates);
222            if let Some(new_offset) = new_polygon.offset {
223                new_offsets.push(new_offset);
224            }
225        }
226    });
227
228    if new_coordinates.is_empty() {
229        None
230    } else {
231        Some(VectorGeometry::MultiPolygon(VectorMultiPolygonGeometry {
232            _type: VectorGeometryType::MultiPolygon,
233            is_3d: *is_3d,
234            coordinates: new_coordinates,
235            bbox: *bbox,
236            offset: Some(new_offsets),
237            vec_bbox: Some(vec_bbox.unwrap_or_default().clip(axis, k1, k2)),
238            ..Default::default()
239        }))
240    }
241}
242
243/// After clipping a line, return the altered line,
244/// the offset the new line starts at,
245/// and if the line is ccw
246#[derive(Debug)]
247pub struct ClipLineResult<M: Clone + Default> {
248    /// The clipped line
249    pub line: VectorLineString<M>,
250    /// The offset the new line starts at
251    pub offset: f64,
252}
253/// Ensuring `vec_bbox` exists
254#[derive(Debug)]
255pub struct ClipLineResultWithBBox<M: Clone + Default> {
256    /// The clipped line
257    pub line: VectorLineString<M>,
258    /// The offset the new line starts at
259    pub offset: f64,
260    /// The new vector bounding box
261    pub vec_bbox: BBox3D,
262}
263
264/// clip an input line to a bounding box
265/// Data should always be in a 0->1 coordinate system to use this clip function
266pub fn clip_line<M: Clone + Default>(
267    geom: &VectorLineString<M>,
268    bbox: BBox3D,
269    is_polygon: bool,
270    offset: Option<f64>,
271    buffer: Option<f64>, /* default for a full size tile. Assuming 1024 extent and a 64 point buffer */
272) -> Vec<ClipLineResultWithBBox<M>> {
273    let offset = offset.unwrap_or(0.);
274    let buffer = buffer.unwrap_or(0.0625);
275    let mut res: Vec<ClipLineResult<M>> = vec![];
276    let BBox3D { left, bottom, right, top, .. } = bbox;
277    // clip horizontally
278    let horizontal_clips = _clip_line(
279        ClipLineResult { line: geom.clone(), offset },
280        left - buffer,
281        right + buffer,
282        Axis::X,
283        is_polygon,
284    );
285    for clip in horizontal_clips {
286        // clip vertically
287        let mut vertical_clips =
288            _clip_line(clip, bottom - buffer, top + buffer, Axis::Y, is_polygon);
289        res.append(&mut vertical_clips);
290    }
291    res.iter_mut()
292        .map(|clip| {
293            let mut vec_bbox: Option<BBox3D> = None;
294            for p in clip.line.iter() {
295                match &mut vec_bbox {
296                    Some(bbox) => bbox.extend_from_point(p),
297                    None => vec_bbox = Some(BBox3D::from_point(p)),
298                }
299            }
300            ClipLineResultWithBBox {
301                line: core::mem::take(&mut clip.line),
302                offset: clip.offset,
303                vec_bbox: vec_bbox.unwrap(),
304            }
305        })
306        .collect()
307}
308
309/// Interal clip tool
310fn _clip_line<M: Clone + Default>(
311    input: ClipLineResult<M>,
312    k1: f64,
313    k2: f64,
314    axis: Axis,
315    is_polygon: bool,
316) -> Vec<ClipLineResult<M>> {
317    //   let { line: geom, offset: startOffset } = input;
318    let geom = &input.line;
319    let start_offset = input.offset;
320    let mut new_geom: Vec<ClipLineResult<M>> = vec![];
321    let mut slice: VectorLineString<M> = vec![];
322    if geom.is_empty() {
323        return new_geom;
324    }
325    let mut last = geom.len() - 1;
326    let intersect = if axis == Axis::X { intersect_x } else { intersect_y };
327
328    let mut cur_offset = start_offset;
329    let mut acc_offset = start_offset;
330    let mut prev_p = &geom[0];
331    let mut first_enter = false;
332
333    let mut i = 0;
334    while i < last {
335        let VectorPoint::<M> { x: ax, y: ay, t: az, m: am, .. } = &geom[i];
336        let VectorPoint::<M> { x: bx, y: by, m: bm, .. } = &geom[i + 1];
337        let a: f64 = if axis == Axis::X { *ax } else { *ay };
338        let b: f64 = if axis == Axis::X { *bx } else { *by };
339        let mut entered = false;
340        let mut exited = false;
341        let mut int_p: Option<VectorPoint<M>> = None;
342
343        // ENTER OR CONTINUE CASES
344        if a < k1 {
345            // ---|-->  | (line enters the clip region from the left)
346            if b > k1 {
347                int_p = Some(intersect(*ax, *ay, *bx, *by, k1, bm));
348                slice.push(int_p.clone().unwrap());
349                entered = true;
350            }
351        } else if a > k2 {
352            // |  <--|--- (line enters the clip region from the right)
353            if b < k2 {
354                int_p = Some(intersect(*ax, *ay, *bx, *by, k2, bm));
355                slice.push(int_p.clone().unwrap());
356                entered = true;
357            }
358        } else {
359            int_p = Some(VectorPoint { x: *ax, y: *ay, z: None, m: am.clone(), t: *az });
360            slice.push(int_p.clone().unwrap());
361        }
362
363        // Update the intersection point and offset if the int_p exists
364        if let Some(int_p) = int_p.as_ref() {
365            // our first enter will change the offset for the line
366            if entered && !first_enter {
367                cur_offset = acc_offset + prev_p.distance(int_p);
368                first_enter = true;
369            }
370        }
371
372        // EXIT CASES
373        if b < k1 && a >= k1 {
374            // <--|---  | or <--|-----|--- (line exits the clip region on the left)
375            int_p = Some(intersect(*ax, *ay, *bx, *by, k1, if bm.is_some() { bm } else { am }));
376            slice.push(int_p.unwrap());
377            exited = true;
378        }
379        if b > k2 && a <= k2 {
380            // |  ---|--> or ---|-----|--> (line exits the clip region on the right)
381            int_p = Some(intersect(*ax, *ay, *bx, *by, k2, if bm.is_some() { bm } else { am }));
382            slice.push(int_p.unwrap());
383            exited = true;
384        }
385
386        // update the offset
387        acc_offset += prev_p.distance(&geom[i + 1]);
388        prev_p = &geom[i + 1];
389
390        // If not a polygon, we can cut it into parts, otherwise we just keep tracking the edges
391        if !is_polygon && exited {
392            new_geom.push(ClipLineResult { line: slice, offset: cur_offset });
393            slice = vec![];
394            first_enter = false;
395        }
396
397        i += 1;
398    }
399
400    // add the last point if inside the clip
401    let last_point = geom[last].clone();
402    let a = if axis == Axis::X { last_point.x } else { last_point.y };
403    if a >= k1 && a <= k2 {
404        slice.push(last_point.clone());
405    }
406
407    // close the polygon if its endpoints are not the same after clipping
408    if !slice.is_empty() && is_polygon {
409        last = slice.len() - 1;
410        let first_p = &slice[0];
411        if last >= 1 && (slice[last].x != first_p.x || slice[last].y != first_p.y) {
412            slice.push(first_p.clone());
413        }
414    }
415
416    // add the final slice
417    if !slice.is_empty() {
418        new_geom.push(ClipLineResult { line: slice, offset: cur_offset });
419    }
420
421    new_geom
422}
423
424/// Find the intersection of two points on the X axis
425fn intersect_x<M: Clone + Default>(
426    ax: f64,
427    ay: f64,
428    bx: f64,
429    by: f64,
430    x: f64,
431    m: &Option<M>,
432) -> VectorPoint<M> {
433    let t = (x - ax) / (bx - ax);
434    VectorPoint::<M> { x, y: ay + (by - ay) * t, z: None, m: m.clone(), t: Some(1.) }
435}
436
437/// Find the intersection of two points on the Y axis
438fn intersect_y<M: Clone + Default>(
439    ax: f64,
440    ay: f64,
441    bx: f64,
442    by: f64,
443    y: f64,
444    m: &Option<M>,
445) -> VectorPoint<M> {
446    let t = (y - ay) / (by - ay);
447    VectorPoint::<M> { x: ax + (bx - ax) * t, y, z: None, m: m.clone(), t: Some(1.) }
448}
gistools/geometry/tools/clip.rs

gistools/geometry/tools/
clip.rs
