gistools/geometry/wm/
clip.rs

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

gistools/geometry/wm/
clip.rs
