use crate::prelude::*;
use crate::{
Coordinate, Line, LineString, MultiLineString, MultiPolygon, Point, Polygon, Triangle,
};
use num_traits::Float;
use std::cmp::Ordering;
use std::collections::BinaryHeap;
use rstar::{RTree, RTreeNum};
#[derive(Debug)]
struct VScore<T>
where
T: Float,
{
left: usize,
current: usize,
right: usize,
area: T,
intersector: bool,
}
impl<T> Ord for VScore<T>
where
T: Float,
{
fn cmp(&self, other: &VScore<T>) -> Ordering {
other.area.partial_cmp(&self.area).unwrap()
}
}
impl<T> PartialOrd for VScore<T>
where
T: Float,
{
fn partial_cmp(&self, other: &VScore<T>) -> Option<Ordering> {
Some(self.cmp(other))
}
}
impl<T> Eq for VScore<T> where T: Float {}
impl<T> PartialEq for VScore<T>
where
T: Float,
{
fn eq(&self, other: &VScore<T>) -> bool
where
T: Float,
{
self.area == other.area
}
}
#[derive(Debug, Clone, Copy)]
enum GeomType {
Line,
Ring,
}
#[derive(Debug, Clone, Copy)]
struct GeomSettings {
initial_min: usize,
min_points: usize,
geomtype: GeomType,
}
fn visvalingam<T>(orig: &LineString<T>, epsilon: &T) -> Vec<Coordinate<T>>
where
T: Float,
{
if orig.0.len() < 3 {
return orig.0.to_vec();
}
let max = orig.0.len();
let mut adjacent: Vec<(_)> = (0..orig.0.len())
.map(|i| {
if i == 0 {
(-1_i32, 1_i32)
} else {
((i - 1) as i32, (i + 1) as i32)
}
})
.collect();
let mut pq = BinaryHeap::new();
for (i, triangle) in orig.triangles().enumerate() {
pq.push(VScore {
area: triangle.area().abs(),
current: i + 1,
left: i,
right: i + 2,
intersector: false,
});
}
while let Some(smallest) = pq.pop() {
if smallest.area > *epsilon {
continue;
}
let (left, right) = adjacent[smallest.current];
if left as i32 != smallest.left as i32 || right as i32 != smallest.right as i32 {
continue;
}
let (ll, _) = adjacent[left as usize];
let (_, rr) = adjacent[right as usize];
adjacent[left as usize] = (ll, right);
adjacent[right as usize] = (left, rr);
adjacent[smallest.current as usize] = (0, 0);
let choices = [(ll, left, right), (left, right, rr)];
for &(ai, current_point, bi) in &choices {
if ai as usize >= max || bi as usize >= max {
continue;
}
let area = Triangle(
orig.0[ai as usize],
orig.0[current_point as usize],
orig.0[bi as usize],
)
.area()
.abs();
pq.push(VScore {
area: area,
current: current_point as usize,
left: ai as usize,
right: bi as usize,
intersector: false,
});
}
}
orig.0
.iter()
.zip(adjacent.iter())
.filter_map(|(tup, adj)| if *adj != (0, 0) { Some(*tup) } else { None })
.collect::<Vec<Coordinate<T>>>()
}
fn vwp_wrapper<T>(
geomtype: &GeomSettings,
exterior: &LineString<T>,
interiors: Option<&[LineString<T>]>,
epsilon: &T,
) -> Vec<Vec<Coordinate<T>>>
where
T: Float + RTreeNum,
{
let mut rings = vec![];
let mut tree: RTree<Line<_>> = RTree::bulk_load(
exterior
.lines()
.chain(
interiors
.iter()
.flat_map(|ring| *ring)
.flat_map(|line_string| line_string.lines()),
)
.collect::<Vec<_>>(),
);
rings.push(visvalingam_preserve(
geomtype, &exterior, epsilon, &mut tree,
));
if let Some(interior_rings) = interiors {
for ring in interior_rings {
rings.push(visvalingam_preserve(geomtype, &ring, epsilon, &mut tree))
}
}
rings
}
fn visvalingam_preserve<T>(
geomtype: &GeomSettings,
orig: &LineString<T>,
epsilon: &T,
tree: &mut RTree<Line<T>>,
) -> Vec<Coordinate<T>>
where
T: Float + RTreeNum,
{
if orig.0.len() < 3 {
return orig.0.to_vec();
}
let max = orig.0.len();
let mut counter = orig.0.len();
let mut adjacent: Vec<(_)> = (0..orig.0.len())
.map(|i| {
if i == 0 {
(-1_i32, 1_i32)
} else {
((i - 1) as i32, (i + 1) as i32)
}
})
.collect();
let mut pq = BinaryHeap::new();
for (i, triangle) in orig.triangles().enumerate() {
let v = VScore {
area: triangle.area().abs(),
current: i + 1,
left: i,
right: i + 2,
intersector: false,
};
pq.push(v);
}
while let Some(mut smallest) = pq.pop() {
if smallest.area > *epsilon {
continue;
}
if counter <= geomtype.initial_min {
break;
}
let (left, right) = adjacent[smallest.current];
if left as i32 != smallest.left as i32 || right as i32 != smallest.right as i32 {
continue;
}
smallest.intersector = tree_intersect(tree, &smallest, &orig.0);
if smallest.intersector && counter <= geomtype.min_points {
break;
}
adjacent[smallest.current as usize] = (0, 0);
counter -= 1;
let left_point = Point(orig.0[left as usize]);
let middle_point = Point(orig.0[smallest.current]);
let right_point = Point(orig.0[right as usize]);
let line_1 = Line::new(left_point, middle_point);
let line_2 = Line::new(middle_point, right_point);
assert!(tree.remove(&line_1).is_some());
assert!(tree.remove(&line_2).is_some());
tree.insert(Line::new(left_point, right_point));
let (ll, _) = adjacent[left as usize];
let (_, rr) = adjacent[right as usize];
adjacent[left as usize] = (ll, right);
adjacent[right as usize] = (left, rr);
let choices = [(ll, left, right), (left, right, rr)];
for &(ai, current_point, bi) in &choices {
if ai as usize >= max || bi as usize >= max {
continue;
}
let new = Triangle(
orig.0[ai as usize],
orig.0[current_point as usize],
orig.0[bi as usize],
);
let temp_area = if smallest.intersector && (current_point as usize) < smallest.current {
-*epsilon
} else {
new.area().abs()
};
let new_triangle = VScore {
area: temp_area,
current: current_point as usize,
left: ai as usize,
right: bi as usize,
intersector: false,
};
pq.push(new_triangle);
}
}
orig.0
.iter()
.zip(adjacent.iter())
.filter_map(|(tup, adj)| if *adj != (0, 0) { Some(*tup) } else { None })
.collect()
}
fn ccw<T>(p1: Point<T>, p2: Point<T>, p3: Point<T>) -> bool
where
T: Float,
{
(p3.y() - p1.y()) * (p2.x() - p1.x()) > (p2.y() - p1.y()) * (p3.x() - p1.x())
}
fn cartesian_intersect<T>(p1: Point<T>, p2: Point<T>, p3: Point<T>, p4: Point<T>) -> bool
where
T: Float,
{
(ccw(p1, p3, p4) ^ ccw(p2, p3, p4)) & (ccw(p1, p2, p3) ^ ccw(p1, p2, p4))
}
fn tree_intersect<T>(tree: &RTree<Line<T>>, triangle: &VScore<T>, orig: &[Coordinate<T>]) -> bool
where
T: Float + RTreeNum,
{
let point_a = orig[triangle.left];
let point_c = orig[triangle.right];
let bounding_rect = Triangle(
orig[triangle.left],
orig[triangle.current],
orig[triangle.right],
)
.bounding_rect();
let br = Point::new(bounding_rect.min.x, bounding_rect.min.y);
let tl = Point::new(bounding_rect.max.x, bounding_rect.max.y);
tree.locate_in_envelope_intersecting(&rstar::AABB::from_corners(br, tl))
.any(|c| {
let (ca, cb) = c.points();
ca.0 != point_a
&& ca.0 != point_c
&& cb.0 != point_a
&& cb.0 != point_c
&& cartesian_intersect(ca, cb, Point(point_a), Point(point_c))
})
}
pub trait SimplifyVW<T, Epsilon = T> {
fn simplifyvw(&self, epsilon: &T) -> Self
where
T: Float;
}
pub trait SimplifyVWPreserve<T, Epsilon = T> {
fn simplifyvw_preserve(&self, epsilon: &T) -> Self
where
T: Float + RTreeNum;
}
impl<T> SimplifyVWPreserve<T> for LineString<T>
where
T: Float + RTreeNum,
{
fn simplifyvw_preserve(&self, epsilon: &T) -> LineString<T> {
let gt = GeomSettings {
initial_min: 2,
min_points: 4,
geomtype: GeomType::Line,
};
let mut simplified = vwp_wrapper(>, self, None, epsilon);
LineString::from(simplified.pop().unwrap())
}
}
impl<T> SimplifyVWPreserve<T> for MultiLineString<T>
where
T: Float + RTreeNum,
{
fn simplifyvw_preserve(&self, epsilon: &T) -> MultiLineString<T> {
MultiLineString(
self.0
.iter()
.map(|l| l.simplifyvw_preserve(epsilon))
.collect(),
)
}
}
impl<T> SimplifyVWPreserve<T> for Polygon<T>
where
T: Float + RTreeNum,
{
fn simplifyvw_preserve(&self, epsilon: &T) -> Polygon<T> {
let gt = GeomSettings {
initial_min: 4,
min_points: 6,
geomtype: GeomType::Ring,
};
let mut simplified = vwp_wrapper(>, self.exterior(), Some(self.interiors()), epsilon);
let exterior = LineString::from(simplified.remove(0));
let interiors = simplified.into_iter().map(LineString::from).collect();
Polygon::new(exterior, interiors)
}
}
impl<T> SimplifyVWPreserve<T> for MultiPolygon<T>
where
T: Float + RTreeNum,
{
fn simplifyvw_preserve(&self, epsilon: &T) -> MultiPolygon<T> {
MultiPolygon(
self.0
.iter()
.map(|p| p.simplifyvw_preserve(epsilon))
.collect(),
)
}
}
impl<T> SimplifyVW<T> for LineString<T>
where
T: Float,
{
fn simplifyvw(&self, epsilon: &T) -> LineString<T> {
LineString::from(visvalingam(self, epsilon))
}
}
impl<T> SimplifyVW<T> for MultiLineString<T>
where
T: Float,
{
fn simplifyvw(&self, epsilon: &T) -> MultiLineString<T> {
MultiLineString(self.0.iter().map(|l| l.simplifyvw(epsilon)).collect())
}
}
impl<T> SimplifyVW<T> for Polygon<T>
where
T: Float,
{
fn simplifyvw(&self, epsilon: &T) -> Polygon<T> {
Polygon::new(
self.exterior().simplifyvw(epsilon),
self.interiors()
.iter()
.map(|l| l.simplifyvw(epsilon))
.collect(),
)
}
}
impl<T> SimplifyVW<T> for MultiPolygon<T>
where
T: Float,
{
fn simplifyvw(&self, epsilon: &T) -> MultiPolygon<T> {
MultiPolygon(self.0.iter().map(|p| p.simplifyvw(epsilon)).collect())
}
}
#[cfg(test)]
mod test {
use super::{
cartesian_intersect, visvalingam, vwp_wrapper, GeomSettings, GeomType, SimplifyVW,
SimplifyVWPreserve,
};
use crate::{Coordinate, LineString, MultiLineString, MultiPolygon, Point, Polygon};
#[test]
fn visvalingam_test() {
let points = vec![
(5.0, 2.0),
(3.0, 8.0),
(6.0, 20.0),
(7.0, 25.0),
(10.0, 10.0),
];
let points_ls: LineString<_> = points.iter().map(|e| Point::new(e.0, e.1)).collect();
let correct = vec![(5.0, 2.0), (7.0, 25.0), (10.0, 10.0)];
let correct_ls: Vec<_> = correct
.iter()
.map(|e| Coordinate::from((e.0, e.1)))
.collect();
let simplified = visvalingam(&points_ls, &30.);
assert_eq!(simplified, correct_ls);
}
#[test]
fn vwp_intersection_test() {
let a = Point::new(1., 3.);
let b = Point::new(3., 1.);
let c = Point::new(3., 3.);
let d = Point::new(1., 1.);
assert_eq!(cartesian_intersect(a, b, c, d), true);
assert_eq!(cartesian_intersect(b, a, c, d), true);
assert_eq!(cartesian_intersect(a, b, d, c), true);
assert_eq!(cartesian_intersect(b, a, d, c), true);
}
#[test]
fn simple_vwp_test() {
let points = vec![
(10., 60.),
(135., 68.),
(94., 48.),
(126., 31.),
(280., 19.),
(117., 48.),
(300., 40.),
(301., 10.),
];
let points_ls: Vec<_> = points.iter().map(|e| Point::new(e.0, e.1)).collect();
let gt = &GeomSettings {
initial_min: 2,
min_points: 4,
geomtype: GeomType::Line,
};
let simplified = vwp_wrapper(>, &points_ls.into(), None, &668.6);
let correct = vec![
(10., 60.),
(126., 31.),
(280., 19.),
(117., 48.),
(300., 40.),
(301., 10.),
];
let correct_ls: Vec<_> = correct
.iter()
.map(|e| Coordinate::from((e.0, e.1)))
.collect();
assert_eq!(simplified[0], correct_ls);
}
#[test]
fn retained_vwp_test() {
let outer = LineString::from(vec![
(-54.4921875, 21.289374355860424),
(-33.5, 56.9449741808516),
(-22.5, 44.08758502824516),
(-19.5, 23.241346102386135),
(-54.4921875, 21.289374355860424),
]);
let inner = LineString::from(vec![
(-24.451171875, 35.266685523707665),
(-29.513671875, 47.32027765985069),
(-22.869140625, 43.80817468459856),
(-24.451171875, 35.266685523707665),
]);
let poly = Polygon::new(outer.clone(), vec![inner]);
let simplified = poly.simplifyvw_preserve(&95.4);
assert_eq!(simplified.exterior(), &outer);
}
#[test]
fn remove_inner_point_vwp_test() {
let outer = LineString::from(vec![
(-54.4921875, 21.289374355860424),
(-33.5, 56.9449741808516),
(-22.5, 44.08758502824516),
(-19.5, 23.241346102386135),
(-54.4921875, 21.289374355860424),
]);
let inner = LineString::from(vec![
(-24.451171875, 35.266685523707665),
(-40.0, 45.),
(-29.513671875, 47.32027765985069),
(-22.869140625, 43.80817468459856),
(-24.451171875, 35.266685523707665),
]);
let correct_inner = LineString::from(vec![
(-24.451171875, 35.266685523707665),
(-40.0, 45.0),
(-22.869140625, 43.80817468459856),
(-24.451171875, 35.266685523707665),
]);
let poly = Polygon::new(outer.clone(), vec![inner]);
let simplified = poly.simplifyvw_preserve(&95.4);
assert_eq!(simplified.exterior(), &outer);
assert_eq!(simplified.interiors()[0], correct_inner);
}
#[test]
fn very_long_vwp_test() {
let points = include!("test_fixtures/norway_main.rs");
let points_ls: Vec<_> = points.iter().map(|e| Point::new(e[0], e[1])).collect();
let gt = &GeomSettings {
initial_min: 2,
min_points: 4,
geomtype: GeomType::Line,
};
let simplified = vwp_wrapper(>, &points_ls.into(), None, &0.0005);
assert_eq!(simplified[0].len(), 3277);
}
#[test]
fn visvalingam_test_long() {
let points = include!("test_fixtures/vw_orig.rs");
let points_ls: LineString<_> = points.iter().map(|e| Point::new(e[0], e[1])).collect();
let correct = include!("test_fixtures/vw_simplified.rs");
let correct_ls: Vec<_> = correct
.iter()
.map(|e| Coordinate::from((e[0], e[1])))
.collect();
let simplified = visvalingam(&points_ls, &0.0005);
assert_eq!(simplified, correct_ls);
}
#[test]
fn visvalingam_preserve_test_long() {
let points = include!("test_fixtures/vw_orig.rs");
let points_ls: LineString<_> = points.iter().map(|e| Point::new(e[0], e[1])).collect();
let correct = include!("test_fixtures/vw_simplified.rs");
let correct_ls: Vec<_> = correct.iter().map(|e| Point::new(e[0], e[1])).collect();
let simplified = LineString::from(points_ls).simplifyvw_preserve(&0.0005);
assert_eq!(simplified, LineString::from(correct_ls));
}
#[test]
fn visvalingam_test_empty_linestring() {
let vec: Vec<[f32; 2]> = Vec::new();
let compare = Vec::new();
let simplified = visvalingam(&LineString::from(vec), &1.0);
assert_eq!(simplified, compare);
}
#[test]
fn visvalingam_test_two_point_linestring() {
let mut vec = Vec::new();
vec.push(Point::new(0.0, 0.0));
vec.push(Point::new(27.8, 0.1));
let mut compare = Vec::new();
compare.push(Coordinate::from((0.0, 0.0)));
compare.push(Coordinate::from((27.8, 0.1)));
let simplified = visvalingam(&LineString::from(vec), &1.0);
assert_eq!(simplified, compare);
}
#[test]
fn multilinestring() {
let points = vec![
(5.0, 2.0),
(3.0, 8.0),
(6.0, 20.0),
(7.0, 25.0),
(10.0, 10.0),
];
let points_ls: Vec<_> = points.iter().map(|e| Point::new(e.0, e.1)).collect();
let correct = vec![(5.0, 2.0), (7.0, 25.0), (10.0, 10.0)];
let correct_ls: Vec<_> = correct.iter().map(|e| Point::new(e.0, e.1)).collect();
let mline = MultiLineString(vec![LineString::from(points_ls)]);
assert_eq!(
mline.simplifyvw(&30.),
MultiLineString(vec![LineString::from(correct_ls)])
);
}
#[test]
fn polygon() {
let poly = Polygon::new(
LineString::from(vec![
(0., 0.),
(0., 10.),
(5., 11.),
(10., 10.),
(10., 0.),
(0., 0.),
]),
vec![],
);
let poly2 = poly.simplifyvw(&10.);
assert_eq!(
poly2,
Polygon::new(
LineString::from(vec![
Point::new(0., 0.),
Point::new(0., 10.),
Point::new(10., 10.),
Point::new(10., 0.),
Point::new(0., 0.),
]),
vec![],
)
);
}
#[test]
fn multipolygon() {
let mpoly = MultiPolygon(vec![Polygon::new(
LineString::from(vec![
(0., 0.),
(0., 10.),
(5., 11.),
(10., 10.),
(10., 0.),
(0., 0.),
]),
vec![],
)]);
let mpoly2 = mpoly.simplifyvw(&10.);
assert_eq!(
mpoly2,
MultiPolygon(vec![Polygon::new(
LineString::from(vec![(0., 0.), (0., 10.), (10., 10.), (10., 0.), (0., 0.)]),
vec![],
)])
);
}
}