simples 0.1.1

N-dimensional linestring smoothing and simplification
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

//! Simplify a linestring using the [Ramer-Douglas-Peucker](https://en.wikipedia.org/wiki/Ramer%E2%80%93Douglas%E2%80%93Peucker_algorithm) algorithm.
use crate::Precision;
use nalgebra::{distance_squared, Point};

fn proj_dist2<const D: usize>(
    start: &Point<Precision, D>,
    end: &Point<Precision, D>,
    p: &Point<Precision, D>,
    length_sq: Precision,
) -> Precision {
    let lensq_along = (p - start).dot(&(end - start));
    let proj;
    if lensq_along <= 0.0 {
        proj = start;
    } else if lensq_along >= length_sq {
        proj = end;
    } else {
        return distance_squared(p, &(start + lensq_along / length_sq * (end - start)));
    }
    distance_squared(p, proj)
}

fn rdp_keep_inner<const D: usize>(
    line: &[Point<Precision, D>],
    epsilon_sq: Precision,
    offset: usize,
) -> Vec<usize> {
    if line.len() <= 2 {
        return vec![];
    }

    let first = line.first().unwrap();
    let last = line.last().unwrap();
    let length_sq = distance_squared(first, last);

    let mut greatest_dist2 = (0, -Precision::INFINITY);
    for (idx, point) in line.iter().enumerate().skip(1).take(line.len() - 2) {
        let d2 = proj_dist2(first, last, point, length_sq);
        if d2 > greatest_dist2.1 {
            greatest_dist2 = (idx, d2)
        }
    }

    let mut to_keep = vec![];
    if greatest_dist2.1 > epsilon_sq {
        if greatest_dist2.0 > 1 {
            to_keep.append(&mut rdp_keep_inner(
                &line[0..=greatest_dist2.0],
                epsilon_sq,
                offset,
            ));
        }
        to_keep.push(greatest_dist2.0);
        if greatest_dist2.0 < line.len() - 2 {
            to_keep.append(&mut rdp_keep_inner(
                &line[greatest_dist2.0..line.len()],
                epsilon_sq,
                offset + greatest_dist2.0,
            ));
        }
    }

    to_keep
}

/// Return the indices of the points in the line which would be kept if simplified using RDP.
pub fn rdp_keep<const D: usize>(line: &[Point<Precision, D>], epsilon: Precision) -> Vec<usize> {
    let epsilon_sq = epsilon * epsilon;
    let mut out = vec![0];
    out.append(&mut rdp_keep_inner(line, epsilon_sq, 0));
    out.push(line.len() - 1);
    out
}

/// Decimate the linestring using RDP.
pub fn rdp_reduce<const D: usize>(
    line: &[Point<Precision, D>],
    epsilon: Precision,
) -> Vec<Point<Precision, D>> {
    let kept = rdp_keep(line, epsilon);
    kept.into_iter().map(|idx| line[idx]).collect()
}