polysplit 0.1.0

Algorithm that allows to split polylines into segments by the defined list of points not necessary belonging to the polyline
Documentation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
use std::fmt;
use std::error;
use std::cmp::{Ord, Ordering, PartialEq, Eq, PartialOrd};
use std::collections::BinaryHeap;
use std::fmt::Debug;
use std::ops::Add;

/// CutRatioResult presents the closest projection of the point to the segment.
pub enum CutRatioResult {
    /// The closest projection is the start of the segment.
    Begin,
    /// The closest projection is the point on the segment that splits it in
    /// the defined proportion (usually `0.0 < ratio < 1.0` where `0.0` is the start
    /// and `1.0` is the end of the segment).
    Medium(f64),
    /// The closest projection is the end of the segment.
    End,
}

impl PartialEq for CutRatioResult {
    fn eq(&self, other: &Self) -> bool {
        match (self, other) {
            (Self::Begin, Self::Begin) => true,
            (Self::Medium(s), Self::Medium(o)) => s.eq(o),
            (Self::End, Self::End) => true,
            _ => false,
        }
    }
}

impl Eq for CutRatioResult {}

impl PartialOrd for CutRatioResult {
    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
        Some(self.cmp(other))
    }
}

impl Ord for CutRatioResult {
    fn cmp(&self, other: &Self) -> Ordering {
        match (self, other) {
            (Self::Begin, Self::Begin) => Ordering::Equal,
            (Self::Begin, _) => Ordering::Less,
            (Self::Medium(_), Self::Begin) => Ordering::Greater,
            (Self::Medium(s), Self::Medium(o)) => s.total_cmp(o),
            (Self::Medium(_), Self::End) => Ordering::Less,
            (Self::End, Self::End) => Ordering::Equal,
            (Self::End, _) => Ordering::Greater,
        }
    }
}

/// DistanceToSegmentResult presents the projection results of the point to the segment.
pub struct DistanceToSegmentResult<P, D>
where
    P: Copy,
    D: Copy + PartialOrd + Add<Output = D>,
{
    pub cut_ratio: CutRatioResult,
    pub cut_point: P,
    pub distance: D,
}

/// PolySplit defines methods for types that can be used in **polyline_split** method.
pub trait PolySplit<D>
where
    Self: Copy,
    D: Copy + PartialOrd + Add<Output = D>,
{
    /// Returns distance to another point.
    ///
    /// # Arguments
    ///
    /// * `point` - A point distance to should be calculated
    fn distance_to_point(&self, point: &Self) -> D;
    /// Returns projection [results](DistanceToSegmentResult) to the segment.
    ///
    /// # Arguments
    ///
    /// * `segment` - A segment presented by a tuple of points
    fn distance_to_segment(&self, segment: (&Self, &Self)) -> DistanceToSegmentResult<Self, D>;
}

#[derive(Debug, Clone, PartialEq, Eq)]
#[non_exhaustive]
pub enum PolySplitErrorKind {
    InvalidPolyline,
    InvalidPoints,
    PointFarAway,
    CannotSplit,
}

#[derive(Debug)]
pub struct PolySplitError {
    pub(super) kind: PolySplitErrorKind,
    pub message: String,
}

impl PolySplitError {
    pub fn kind(&self) -> &PolySplitErrorKind {
        &self.kind
    }
}

impl fmt::Display for PolySplitError {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
        write!(f, "{}", self.message)
    }
}

impl error::Error for PolySplitError {}

pub type Result<T> = std::result::Result<T, PolySplitError>;

struct CutPoint<P>
where P: std::fmt::Debug {
    segment_index: usize,
    cut_ratio: CutRatioResult,
    cut_point: P,
}

struct Vertex<D> {
    point_index: usize,
    cut_point_index: usize,
    distance_to: D,
}

#[derive(Copy, Clone, PartialEq)]
struct State<D> {
    distance_total: D,
    position: usize,
}

impl<D: PartialOrd> Eq for State<D> {
}

// Heap order depends on `Ord`, so implementing trait to make min-heap instead of max-heap
impl<D: Copy + PartialOrd> Ord for State<D> {
    fn cmp(&self, other: &Self) -> Ordering {
        // Notice that the we flip the ordering on costs.
        // In case of a tie we compare positions - this step is necessary
        // to make implementations of `PartialEq` and `Ord` consistent.
        other.distance_total
            .partial_cmp(&self.distance_total)
            .unwrap()
            .then_with(|| self.position.cmp(&other.position))
    }
}

// It is required `PartialOrd` to be also implemented
impl<D: Copy + PartialOrd> PartialOrd for State<D> {
    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
        Some(self.cmp(other))
    }
}

/// Splits polyline into segments by the defined list of points.
///
/// # Examples
///
/// ```
/// use polysplit::euclidean::Point;
/// use polysplit::polyline_split;
///
/// let polyline = vec![
///     Point(0.0, 0.0),
///     Point(10.0, 0.0),
///     Point(20.0, 0.0),
/// ];
/// let points = vec![
///     Point(1.0, 1.0),
///     Point(19.0, 1.0),
/// ];
///
/// let segments = polyline_split(&polyline, &points, None).unwrap();
///
/// assert_eq!(segments.len(), 1);
/// println!("{:?}", segments[0]);
/// ```
pub fn polyline_split<P, D>(
    polyline: &[P],
    points: &[P],
    distance_threshold: Option<D>,
) -> Result<Vec<Vec<P>>>
where
    P: PolySplit<D> + std::fmt::Debug,
    D: Copy + PartialOrd + Add<Output = D>,
{
    if polyline.len() <= 1 {
        return Err(PolySplitError{
            kind: PolySplitErrorKind::InvalidPolyline,
            message: "polyline has not enough points".to_string(),
        });
    }

    if points.len() <= 1 {
        return Err(PolySplitError{
            kind: PolySplitErrorKind::InvalidPoints,
            message: "number of points are not enough".to_string(),
        });
    }

    let segments_len = polyline.len() - 1;
    let points_len = points.len();

    // Collecting all possible cut points
    let mut cut_points: Vec<CutPoint<P>> = Vec::new();

    for segment_index in 0..segments_len {
        let segment_a = &polyline[segment_index];
        let segment_b = &polyline[segment_index + 1];

        let mut is_start_added = false;
        let mut is_end_added = false;

        for point in points.iter() {
            let psd: DistanceToSegmentResult<P, D> = point.distance_to_segment((segment_a, segment_b));
            if let Some(dt) = distance_threshold {
                if psd.distance > dt {
                    continue;
                }
            }

            match psd.cut_ratio {
                CutRatioResult::Begin => {
                    if segment_index == 0 && !is_start_added {
                        cut_points.push(CutPoint {
                            segment_index,
                            cut_ratio: psd.cut_ratio,
                            cut_point: *segment_a,
                        });

                        is_start_added = true;
                    }
                }

                CutRatioResult::End => {
                    if !is_end_added {
                        cut_points.push(CutPoint {
                            segment_index,
                            cut_ratio: psd.cut_ratio,
                            cut_point: *segment_b,
                        });

                        is_end_added = true;
                    }
                },

                _ => {
                    cut_points.push(CutPoint {
                        segment_index,
                        cut_ratio: psd.cut_ratio,
                        cut_point: psd.cut_point,
                    });
                }
            }
        }
    }

    cut_points.sort_unstable_by(|a, b| {
        match a.segment_index.cmp(&b.segment_index) {
            Ordering::Equal => a.cut_ratio.partial_cmp(&b.cut_ratio).unwrap(),
            v => v,
        }
    });

    // Building graph
    let mut vertexes: Vec<Vertex<D>> = Vec::new();
    let mut edges: Vec<(usize, usize)> = Vec::new();

    let mut last_reachable_cut_point_index = 0;

    for (point_index, point) in points.iter().enumerate() {
        let start_position = vertexes.len();
        let mut first_match_cut_point_index = None;

        for (cut_point_index, cut_point) in cut_points.iter().enumerate().skip(last_reachable_cut_point_index) {
            let distance_to = point.distance_to_point(&cut_point.cut_point);
            if let Some(dt) = distance_threshold {
                if distance_to > dt {
                    continue;
                }
            }

            if first_match_cut_point_index.is_none() {
                first_match_cut_point_index = Some(cut_point_index);
            }

            vertexes.push(Vertex {
                point_index,
                cut_point_index,
                distance_to,
            });
        }

        let end_position = vertexes.len();
        if start_position == end_position {
            return Err(PolySplitError{
                kind: PolySplitErrorKind::PointFarAway,
                message: "point has no closest segments".to_string(),
            });
        }

        edges.push((start_position, end_position));
        last_reachable_cut_point_index = first_match_cut_point_index.unwrap_or_default();
    }

    // Initializing start points
    let vertexes_len = vertexes.len();
    let mut dist: Vec<Option<D>> = (0..vertexes_len).map(|_| None).collect();
    let mut prev: Vec<Option<usize>> = (0..vertexes_len).map(|_| None).collect();
    let mut priority_queue = BinaryHeap::new();

    for idx in edges[0].0..edges[0].1 {
        let vertex = &vertexes[idx];

        dist[idx] = Some(vertex.distance_to);
        prev[idx] = None;
        priority_queue.push(State {
            distance_total: vertex.distance_to,
            position: idx,
        });
    }

    // Searching for shortest path using Dijkstra's algorithm
    let mut destination = None;
    while let Some(State { distance_total, position }) = priority_queue.pop() {
        let current_vertex = &vertexes[position];

        // Goal is reached
        if current_vertex.point_index + 1 == points_len {
            destination = Some(position);
            break;
        }

        // Useless state because there is better one
        if let Some(d) = dist[position] {
            if distance_total > d {
                continue;
            }
        }

        // Iterating connected vertexes
        let (from_idx, to_idx) = edges[current_vertex.point_index + 1];
        for idx in from_idx..to_idx {
            let neighbour_vertex = &vertexes[idx];

            if current_vertex.cut_point_index > neighbour_vertex.cut_point_index {
                continue;
            }

            let relaxed_distance_total = distance_total + neighbour_vertex.distance_to;
            if dist[idx].map_or(true, |d| d > relaxed_distance_total) {
                dist[idx] = Some(relaxed_distance_total);
                prev[idx] = Some(position);
                priority_queue.push(State {
                    distance_total: relaxed_distance_total,
                    position: idx,
                });
            }
        }
    }

    // Restoring path
    let mut path = Vec::new();
    while let Some(idx) = destination {
        path.push(idx);
        destination = prev[idx];
    }

    if path.is_empty() {
        return Err(PolySplitError{
            kind: PolySplitErrorKind::CannotSplit,
            message: "cannot split polyline".to_string(),
        });
    }

    path.reverse();

    // Building sub-segments
    let mut segments: Vec<_> = Vec::with_capacity(segments_len);
    let mut current_vertex = &vertexes[path[0]];
    let mut current_cut_point = &cut_points[current_vertex.cut_point_index];

    for next_idx in path[1..].iter() {
        let next_vertex: &Vertex<D> = &vertexes[*next_idx];
        let next_cut_point = &cut_points[next_vertex.cut_point_index];
        let mut segment: Vec<_> = Vec::new();

        if !matches!(current_cut_point.cut_ratio, CutRatioResult::End) {
            segment.push(current_cut_point.cut_point);
        }

        for segment_idx in current_cut_point.segment_index..next_cut_point.segment_index {
            segment.push(polyline[segment_idx + 1]);
        }

        if !matches!(next_cut_point.cut_ratio, CutRatioResult::Begin) {
            segment.push(next_cut_point.cut_point);
        }

        // Two points are matched to same cut point
        // So adding same point to be valid segment
        if segment.len() == 1 {
            segment.push(segment[0]);
        }

        segments.push(segment);
        current_vertex = next_vertex;
        current_cut_point = &cut_points[current_vertex.cut_point_index];
    }

    Ok(segments)
}