geo/algorithm/

simplify_vw.rs

use crate::prelude::*;
use crate::{
    Coord, CoordFloat, HasKernel, Line, LineString, MultiLineString, MultiPolygon, Point, Polygon,
    Triangle,
};
use std::cmp::Ordering;
use std::collections::BinaryHeap;

use rstar::{RTree, RTreeNum};

/// Store triangle information
// current is the candidate point for removal
#[derive(Debug)]
struct VScore<T, I>
where
    T: CoordFloat,
{
    left: usize,
    current: usize,
    right: usize,
    area: T,
    // `visvalingam_preserve` uses `intersector`, `visvalingam` does not
    intersector: I,
}

// These impls give us a min-heap
impl<T, I> Ord for VScore<T, I>
where
    T: CoordFloat,
{
    fn cmp(&self, other: &VScore<T, I>) -> Ordering {
        other.area.partial_cmp(&self.area).unwrap()
    }
}

impl<T, I> PartialOrd for VScore<T, I>
where
    T: CoordFloat,
{
    fn partial_cmp(&self, other: &VScore<T, I>) -> Option<Ordering> {
        Some(self.cmp(other))
    }
}

impl<T, I> Eq for VScore<T, I> where T: CoordFloat {}

impl<T, I> PartialEq for VScore<T, I>
where
    T: CoordFloat,
{
    fn eq(&self, other: &VScore<T, I>) -> bool
    where
        T: CoordFloat,
    {
        self.area == other.area
    }
}

/// Simplify a line using the [Visvalingam-Whyatt](http://www.tandfonline.com/doi/abs/10.1179/000870493786962263) algorithm
//
// This method returns the **indices** of the simplified line
// epsilon is the minimum triangle area
// The paper states that:
// If [the new triangle's] calculated area is less than that of the last point to be
// eliminated, use the latter's area instead.
// (This ensures that the current point cannot be eliminated
// without eliminating previously eliminated points)
// (Visvalingam and Whyatt 2013, p47)
// However, this does *not* apply if you're using a user-defined epsilon;
// It's OK to remove triangles with areas below the epsilon,
// then recalculate the new triangle area and push it onto the heap
// based on Huon Wilson's original implementation:
// https://github.com/huonw/isrustfastyet/blob/25e7a68ff26673a8556b170d3c9af52e1c818288/mem/line_simplify.rs
fn visvalingam_indices<T>(orig: &LineString<T>, epsilon: &T) -> Vec<usize>
where
    T: CoordFloat,
{
    // No need to continue without at least three points
    if orig.0.len() < 3 {
        return orig.0.iter().enumerate().map(|(idx, _)| idx).collect();
    }

    let max = orig.0.len();

    // Adjacent retained points. Simulating the points in a
    // linked list with indices into `orig`. Big number (larger than or equal to
    // `max`) means no next element, and (0, 0) means deleted element.
    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();

    // Store all the triangles in a minimum priority queue, based on their area.
    //
    // Invalid triangles are *not* removed if / when points are removed; they're
    // handled by skipping them as necessary in the main loop by checking the
    // corresponding entry in adjacent for (0, 0) values

    // Compute the initial triangles
    let mut pq = orig
        .triangles()
        .enumerate()
        .map(|(i, triangle)| VScore {
            area: triangle.unsigned_area(),
            current: i + 1,
            left: i,
            right: i + 2,
            intersector: (),
        })
        .collect::<BinaryHeap<VScore<T, ()>>>();
    // While there are still points for which the associated triangle
    // has an area below the epsilon
    while let Some(smallest) = pq.pop() {
        // This triangle's area is above epsilon, so skip it
        if smallest.area > *epsilon {
            continue;
        }
        //  This triangle's area is below epsilon: eliminate the associated point
        let (left, right) = adjacent[smallest.current];
        // A point in this triangle has been removed since this VScore
        // was created, so skip it
        if left != smallest.left as i32 || right != smallest.right as i32 {
            continue;
        }
        // We've got a valid triangle, and its area is smaller than epsilon, so
        // remove it from the simulated "linked list"
        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] = (0, 0);

        // Now recompute the adjacent triangle(s), using left and right adjacent points
        let choices = [(ll, left, right), (left, right, rr)];
        for &(ai, current_point, bi) in &choices {
            if ai as usize >= max || bi as usize >= max {
                // Out of bounds, i.e. we're on one edge
                continue;
            }
            let area = Triangle::new(
                orig.0[ai as usize],
                orig.0[current_point as usize],
                orig.0[bi as usize],
            )
            .unsigned_area();
            pq.push(VScore {
                area,
                current: current_point as usize,
                left: ai as usize,
                right: bi as usize,
                intersector: (),
            });
        }
    }
    // Filter out the points that have been deleted, returning remaining point indices
    orig.0
        .iter()
        .enumerate()
        .zip(adjacent.iter())
        .filter_map(|(tup, adj)| if *adj != (0, 0) { Some(tup.0) } else { None })
        .collect::<Vec<usize>>()
}

// Wrapper for visvalingam_indices, mapping indices back to points
fn visvalingam<T>(orig: &LineString<T>, epsilon: &T) -> Vec<Coord<T>>
where
    T: CoordFloat,
{
    // Epsilon must be greater than zero for any meaningful simplification to happen
    if *epsilon <= T::zero() {
        return orig.0.to_vec();
    }
    let subset = visvalingam_indices(orig, epsilon);
    // filter orig using the indices
    // using get would be more robust here, but the input subset is guaranteed to be valid in this case
    orig.0
        .iter()
        .zip(subset.iter())
        .map(|(_, s)| orig[*s])
        .collect()
}

// Wrap the actual VW function so the R* Tree can be shared.
// this ensures that shell and rings have access to all segments, so
// intersections between outer and inner rings are detected
//
// Constants:
//
// * `INITIAL_MIN`
//   * If we ever have fewer than these, stop immediately
// * `MIN_POINTS`
//   * If we detect a self-intersection before point removal, and we only have `MIN_POINTS` left,
//     stop: since a self-intersection causes removal of the spatially previous point, THAT could
//     lead to a further self-intersection without the possibility of removing more points,
//     potentially leaving the geometry in an invalid state.
fn vwp_wrapper<T, const INITIAL_MIN: usize, const MIN_POINTS: usize>(
    exterior: &LineString<T>,
    interiors: Option<&[LineString<T>]>,
    epsilon: &T,
) -> Vec<Vec<Coord<T>>>
where
    T: CoordFloat + RTreeNum + HasKernel,
{
    let mut rings = vec![];
    // Populate R* tree with exterior and interior samples, if any
    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<_>>(),
    );

    // Simplify shell
    rings.push(visvalingam_preserve::<T, INITIAL_MIN, MIN_POINTS>(
        exterior, epsilon, &mut tree,
    ));
    // Simplify interior rings, if any
    if let Some(interior_rings) = interiors {
        for ring in interior_rings {
            rings.push(visvalingam_preserve::<T, INITIAL_MIN, MIN_POINTS>(
                ring, epsilon, &mut tree,
            ))
        }
    }
    rings
}

/// Visvalingam-Whyatt with self-intersection detection to preserve topologies
/// this is a port of the technique at https://www.jasondavies.com/simplify/
//
// Constants:
//
// * `INITIAL_MIN`
//   * If we ever have fewer than these, stop immediately
// * `MIN_POINTS`
//   * If we detect a self-intersection before point removal, and we only have `MIN_POINTS` left,
//     stop: since a self-intersection causes removal of the spatially previous point, THAT could
//     lead to a further self-intersection without the possibility of removing more points,
//     potentially leaving the geometry in an invalid state.
fn visvalingam_preserve<T, const INITIAL_MIN: usize, const MIN_POINTS: usize>(
    orig: &LineString<T>,
    epsilon: &T,
    tree: &mut RTree<Line<T>>,
) -> Vec<Coord<T>>
where
    T: CoordFloat + RTreeNum + HasKernel,
{
    if orig.0.len() < 3 || *epsilon <= T::zero() {
        return orig.0.to_vec();
    }
    let max = orig.0.len();
    let mut counter = orig.0.len();

    // Adjacent retained points. Simulating the points in a
    // linked list with indices into `orig`. Big number (larger than or equal to
    // `max`) means no next element, and (0, 0) means deleted element.
    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();
    // Store all the triangles in a minimum priority queue, based on their area.
    //
    // Invalid triangles are *not* removed if / when points are removed; they're
    // handled by skipping them as necessary in the main loop by checking the
    // corresponding entry in adjacent for (0, 0) values

    // Compute the initial triangles
    let mut pq = orig
        .triangles()
        .enumerate()
        .map(|(i, triangle)| VScore {
            area: triangle.unsigned_area(),
            current: i + 1,
            left: i,
            right: i + 2,
            intersector: false,
        })
        .collect::<BinaryHeap<VScore<T, bool>>>();

    // While there are still points for which the associated triangle
    // has an area below the epsilon
    while let Some(mut smallest) = pq.pop() {
        if smallest.area > *epsilon {
            continue;
        }
        if counter <= INITIAL_MIN {
            // we can't remove any more points no matter what
            break;
        }
        let (left, right) = adjacent[smallest.current];
        // A point in this triangle has been removed since this VScore
        // was created, so skip it
        if left != smallest.left as i32 || right != smallest.right as i32 {
            continue;
        }
        // if removal of this point causes a self-intersection, we also remove the previous point
        // that removal alters the geometry, removing the self-intersection
        // HOWEVER if we're within 2 points of the absolute minimum, we can't remove this point or the next
        // because we could then no longer form a valid geometry if removal of next also caused an intersection.
        // The simplification process is thus over.
        smallest.intersector = tree_intersect(tree, &smallest, &orig.0);
        if smallest.intersector && counter <= MIN_POINTS {
            break;
        }
        // We've got a valid triangle, and its area is smaller than epsilon, so
        // remove it from the simulated "linked list"
        adjacent[smallest.current] = (0, 0);
        counter -= 1;
        // Remove stale segments from R* tree
        let left_point = Point::from(orig.0[left as usize]);
        let middle_point = Point::from(orig.0[smallest.current]);
        let right_point = Point::from(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());

        // Restore continuous line segment
        tree.insert(Line::new(left_point, right_point));

        // Now recompute the adjacent triangle(s), using left and right adjacent points
        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 {
                // Out of bounds, i.e. we're on one edge
                continue;
            }
            let new = Triangle::new(
                orig.0[ai as usize],
                orig.0[current_point as usize],
                orig.0[bi as usize],
            );
            // The current point causes a self-intersection, and this point precedes it
            // we ensure it gets removed next by demoting its area to negative epsilon
            let temp_area = if smallest.intersector && (current_point as usize) < smallest.current {
                -*epsilon
            } else {
                new.unsigned_area()
            };
            let new_triangle = VScore {
                area: temp_area,
                current: current_point as usize,
                left: ai as usize,
                right: bi as usize,
                intersector: false,
            };

            // push re-computed triangle onto heap
            pq.push(new_triangle);
        }
    }
    // Filter out the points that have been deleted, returning remaining points
    orig.0
        .iter()
        .zip(adjacent.iter())
        .filter_map(|(tup, adj)| if *adj != (0, 0) { Some(*tup) } else { None })
        .collect()
}

/// is p1 -> p2 -> p3 wound counterclockwise?
#[inline]
fn ccw<T>(p1: Point<T>, p2: Point<T>, p3: Point<T>) -> bool
where
    T: CoordFloat + HasKernel,
{
    let o = <T as HasKernel>::Ker::orient2d(p1.into(), p2.into(), p3.into());
    o == Orientation::CounterClockwise
}

/// checks whether line segments with p1-p4 as their start and endpoints touch or cross
fn cartesian_intersect<T>(p1: Point<T>, p2: Point<T>, p3: Point<T>, p4: Point<T>) -> bool
where
    T: CoordFloat + HasKernel,
{
    (ccw(p1, p3, p4) ^ ccw(p2, p3, p4)) & (ccw(p1, p2, p3) ^ ccw(p1, p2, p4))
}

/// check whether a triangle's edges intersect with any other edges of the LineString
fn tree_intersect<T>(tree: &RTree<Line<T>>, triangle: &VScore<T, bool>, orig: &[Coord<T>]) -> bool
where
    T: CoordFloat + RTreeNum + HasKernel,
{
    let point_a = orig[triangle.left];
    let point_c = orig[triangle.right];
    let bounding_rect = Triangle::new(
        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| {
            // triangle start point, end point
            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::from(point_a), Point::from(point_c))
        })
}

/// Simplifies a geometry.
///
/// Polygons are simplified by running the algorithm on all their constituent rings.  This may
/// result in invalid Polygons, and has no guarantee of preserving topology. Multi* objects are
/// simplified by simplifying all their constituent geometries individually.
///
/// An epsilon less than or equal to zero will return an unaltered version of the geometry.
pub trait SimplifyVw<T, Epsilon = T> {
    /// Returns the simplified representation of a geometry, using the [Visvalingam-Whyatt](http://www.tandfonline.com/doi/abs/10.1179/000870493786962263) algorithm
    ///
    /// See [here](https://bost.ocks.org/mike/simplify/) for a graphical explanation
    ///
    /// # Examples
    ///
    /// ```
    /// use geo::SimplifyVw;
    /// use geo::line_string;
    ///
    /// let line_string = line_string![
    ///     (x: 5.0, y: 2.0),
    ///     (x: 3.0, y: 8.0),
    ///     (x: 6.0, y: 20.0),
    ///     (x: 7.0, y: 25.0),
    ///     (x: 10.0, y: 10.0),
    /// ];
    ///
    /// let simplified = line_string.simplify_vw(&30.0);
    ///
    /// let expected = line_string![
    ///     (x: 5.0, y: 2.0),
    ///     (x: 7.0, y: 25.0),
    ///     (x: 10.0, y: 10.0),
    /// ];
    ///
    /// assert_eq!(expected, simplified);
    /// ```
    fn simplify_vw(&self, epsilon: &T) -> Self
    where
        T: CoordFloat;
}

/// Simplifies a geometry, returning the retained _indices_ of the output
///
/// This operation uses the Visvalingam-Whyatt algorithm,
/// and does **not** guarantee that the returned geometry is valid.
///
/// An epsilon less than or equal to zero will return an unaltered version of the geometry.
pub trait SimplifyVwIdx<T, Epsilon = T> {
    /// Returns the simplified representation of a geometry, using the [Visvalingam-Whyatt](http://www.tandfonline.com/doi/abs/10.1179/000870493786962263) algorithm
    ///
    /// See [here](https://bost.ocks.org/mike/simplify/) for a graphical explanation
    ///
    /// # Examples
    ///
    /// ```
    /// use geo::SimplifyVwIdx;
    /// use geo::line_string;
    ///
    /// let line_string = line_string![
    ///     (x: 5.0, y: 2.0),
    ///     (x: 3.0, y: 8.0),
    ///     (x: 6.0, y: 20.0),
    ///     (x: 7.0, y: 25.0),
    ///     (x: 10.0, y: 10.0),
    /// ];
    ///
    /// let simplified = line_string.simplify_vw_idx(&30.0);
    ///
    /// let expected = vec![
    ///     0_usize,
    ///     3_usize,
    ///     4_usize,
    /// ];
    ///
    /// assert_eq!(expected, simplified);
    /// ```
    fn simplify_vw_idx(&self, epsilon: &T) -> Vec<usize>
    where
        T: CoordFloat;
}

/// Simplifies a geometry, preserving its topology by removing self-intersections
///
/// An epsilon less than or equal to zero will return an unaltered version of the geometry.
pub trait SimplifyVwPreserve<T, Epsilon = T> {
    /// Returns the simplified representation of a geometry, using a topology-preserving variant of the
    /// [Visvalingam-Whyatt](http://www.tandfonline.com/doi/abs/10.1179/000870493786962263) algorithm.
    ///
    /// See [here](https://www.jasondavies.com/simplify/) for a graphical explanation.
    ///
    /// The topology-preserving algorithm uses an [R* tree](../../../rstar/struct.RTree.html) to
    /// efficiently find candidate line segments which are tested for intersection with a given triangle.
    /// If intersections are found, the previous point (i.e. the left component of the current triangle)
    /// is also removed, altering the geometry and removing the intersection.
    ///
    /// In the example below, `(135.0, 68.0)` would be retained by the standard algorithm,
    /// forming triangle `(0, 1, 3),` which intersects with the segments `(280.0, 19.0),
    /// (117.0, 48.0)` and `(117.0, 48.0), (300,0, 40.0)`. By removing it,
    /// a new triangle with indices `(0, 3, 4)` is formed, which does not cause a self-intersection.
    ///
    /// **Note**: it is possible for the simplification algorithm to displace a Polygon's interior ring outside its shell.
    ///
    /// **Note**: if removal of a point causes a self-intersection, but the geometry only has `n + 2`
    /// points remaining (4 for a `LineString`, 6 for a `Polygon`), the point is retained and the
    /// simplification process ends. This is because there is no guarantee that removal of two points will remove
    /// the intersection, but removal of further points would leave too few points to form a valid geometry.
    ///
    /// # Examples
    ///
    /// ```
    /// use approx::assert_relative_eq;
    /// use geo::SimplifyVwPreserve;
    /// use geo::line_string;
    ///
    /// let line_string = line_string![
    ///     (x: 10., y: 60.),
    ///     (x: 135., y: 68.),
    ///     (x: 94., y: 48.),
    ///     (x: 126., y: 31.),
    ///     (x: 280., y: 19.),
    ///     (x: 117., y: 48.),
    ///     (x: 300., y: 40.),
    ///     (x: 301., y: 10.),
    /// ];
    ///
    /// let simplified = line_string.simplify_vw_preserve(&668.6);
    ///
    /// let expected = line_string![
    ///     (x: 10., y: 60.),
    ///     (x: 126., y: 31.),
    ///     (x: 280., y: 19.),
    ///     (x: 117., y: 48.),
    ///     (x: 300., y: 40.),
    ///     (x: 301., y: 10.),
    /// ];
    ///
    /// assert_relative_eq!(expected, simplified, epsilon = 1e-6);
    /// ```
    fn simplify_vw_preserve(&self, epsilon: &T) -> Self
    where
        T: CoordFloat + RTreeNum;
}

impl<T> SimplifyVwPreserve<T> for LineString<T>
where
    T: CoordFloat + RTreeNum + HasKernel,
{
    fn simplify_vw_preserve(&self, epsilon: &T) -> LineString<T> {
        let mut simplified = vwp_wrapper::<_, 2, 4>(self, None, epsilon);
        LineString::from(simplified.pop().unwrap())
    }
}

impl<T> SimplifyVwPreserve<T> for MultiLineString<T>
where
    T: CoordFloat + RTreeNum + HasKernel,
{
    fn simplify_vw_preserve(&self, epsilon: &T) -> MultiLineString<T> {
        MultiLineString::new(
            self.0
                .iter()
                .map(|l| l.simplify_vw_preserve(epsilon))
                .collect(),
        )
    }
}

impl<T> SimplifyVwPreserve<T> for Polygon<T>
where
    T: CoordFloat + RTreeNum + HasKernel,
{
    fn simplify_vw_preserve(&self, epsilon: &T) -> Polygon<T> {
        let mut simplified =
            vwp_wrapper::<_, 4, 6>(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: CoordFloat + RTreeNum + HasKernel,
{
    fn simplify_vw_preserve(&self, epsilon: &T) -> MultiPolygon<T> {
        MultiPolygon::new(
            self.0
                .iter()
                .map(|p| p.simplify_vw_preserve(epsilon))
                .collect(),
        )
    }
}

impl<T> SimplifyVw<T> for LineString<T>
where
    T: CoordFloat,
{
    fn simplify_vw(&self, epsilon: &T) -> LineString<T> {
        LineString::from(visvalingam(self, epsilon))
    }
}

impl<T> SimplifyVwIdx<T> for LineString<T>
where
    T: CoordFloat,
{
    fn simplify_vw_idx(&self, epsilon: &T) -> Vec<usize> {
        visvalingam_indices(self, epsilon)
    }
}

impl<T> SimplifyVw<T> for MultiLineString<T>
where
    T: CoordFloat,
{
    fn simplify_vw(&self, epsilon: &T) -> MultiLineString<T> {
        MultiLineString::new(self.iter().map(|l| l.simplify_vw(epsilon)).collect())
    }
}

impl<T> SimplifyVw<T> for Polygon<T>
where
    T: CoordFloat,
{
    fn simplify_vw(&self, epsilon: &T) -> Polygon<T> {
        Polygon::new(
            self.exterior().simplify_vw(epsilon),
            self.interiors()
                .iter()
                .map(|l| l.simplify_vw(epsilon))
                .collect(),
        )
    }
}

impl<T> SimplifyVw<T> for MultiPolygon<T>
where
    T: CoordFloat,
{
    fn simplify_vw(&self, epsilon: &T) -> MultiPolygon<T> {
        MultiPolygon::new(self.iter().map(|p| p.simplify_vw(epsilon)).collect())
    }
}

#[cfg(test)]
mod test {
    use super::{cartesian_intersect, visvalingam, vwp_wrapper, SimplifyVw, SimplifyVwPreserve};
    use crate::{
        line_string, point, polygon, Coord, LineString, MultiLineString, MultiPolygon, Point,
        Polygon,
    };

    #[test]
    fn visvalingam_test() {
        // this is the PostGIS example
        let ls = line_string![
            (x: 5.0, y: 2.0),
            (x: 3.0, y: 8.0),
            (x: 6.0, y: 20.0),
            (x: 7.0, y: 25.0),
            (x: 10.0, y: 10.0)
        ];

        let correct = vec![(5.0, 2.0), (7.0, 25.0), (10.0, 10.0)];
        let correct_ls: Vec<_> = correct.iter().map(|e| Coord::from((e.0, e.1))).collect();

        let simplified = visvalingam(&ls, &30.);
        assert_eq!(simplified, correct_ls);
    }
    #[test]
    fn vwp_intersection_test() {
        // does the intersection check always work
        let a = point!(x: 1., y: 3.);
        let b = point!(x: 3., y: 1.);
        let c = point!(x: 3., y: 3.);
        let d = point!(x: 1., y: 1.);
        // cw + ccw
        assert!(cartesian_intersect(a, b, c, d));
        // ccw + ccw
        assert!(cartesian_intersect(b, a, c, d));
        // cw + cw
        assert!(cartesian_intersect(a, b, d, c));
        // ccw + cw
        assert!(cartesian_intersect(b, a, d, c));
    }
    #[test]
    fn simple_vwp_test() {
        // this LineString will have a self-intersection if the point with the
        // smallest associated area is removed
        // the associated triangle is (1, 2, 3), and has an area of 668.5
        // the new triangle (0, 1, 3) self-intersects with triangle (3, 4, 5)
        // Point 1 must also be removed giving a final, valid
        // LineString of (0, 3, 4, 5, 6, 7)
        let ls = line_string![
            (x: 10., y:60.),
            (x: 135., y: 68.),
            (x: 94.,  y: 48.),
            (x: 126., y: 31.),
            (x: 280., y: 19.),
            (x: 117., y: 48.),
            (x: 300., y: 40.),
            (x: 301., y: 10.)
        ];
        let simplified = vwp_wrapper::<_, 2, 4>(&ls, None, &668.6);
        // this is the correct, non-intersecting LineString
        let correct = vec![
            (10., 60.),
            (126., 31.),
            (280., 19.),
            (117., 48.),
            (300., 40.),
            (301., 10.),
        ];
        let correct_ls: Vec<_> = correct.iter().map(|e| Coord::from((e.0, e.1))).collect();
        assert_eq!(simplified[0], correct_ls);
    }
    #[test]
    fn retained_vwp_test() {
        // we would expect outer[2] to be removed, as its associated area
        // is below epsilon. However, this causes a self-intersection
        // with the inner ring, which would also trigger removal of outer[1],
        // leaving the geometry below min_points. It is thus retained.
        // Inner should also be reduced, but has points == initial_min for the Polygon type
        let outer = line_string![
            (x: -54.4921875, y: 21.289374355860424),
            (x: -33.5, y: 56.9449741808516),
            (x: -22.5, y: 44.08758502824516),
            (x: -19.5, y: 23.241346102386135),
            (x: -54.4921875, y: 21.289374355860424)
        ];
        let inner = line_string![
            (x: -24.451171875, y: 35.266685523707665),
            (x: -29.513671875, y: 47.32027765985069),
            (x: -22.869140625, y: 43.80817468459856),
            (x: -24.451171875, y: 35.266685523707665)
        ];
        let poly = Polygon::new(outer.clone(), vec![inner]);
        let simplified = poly.simplify_vw_preserve(&95.4);
        assert_relative_eq!(simplified.exterior(), &outer, epsilon = 1e-6);
    }
    #[test]
    fn remove_inner_point_vwp_test() {
        // we would expect outer[2] to be removed, as its associated area
        // is below epsilon. However, this causes a self-intersection
        // with the inner ring, which would also trigger removal of outer[1],
        // leaving the geometry below min_points. It is thus retained.
        // Inner should be reduced to four points by removing inner[2]
        let outer = line_string![
            (x: -54.4921875, y: 21.289374355860424),
            (x: -33.5, y: 56.9449741808516),
            (x: -22.5, y: 44.08758502824516),
            (x: -19.5, y: 23.241346102386135),
            (x: -54.4921875, y: 21.289374355860424)
        ];
        let inner = line_string![
            (x: -24.451171875, y: 35.266685523707665),
            (x: -40.0, y: 45.),
            (x: -29.513671875, y: 47.32027765985069),
            (x: -22.869140625, y: 43.80817468459856),
            (x: -24.451171875, y: 35.266685523707665)
        ];
        let correct_inner = line_string![
            (x: -24.451171875, y: 35.266685523707665),
            (x: -40.0, y: 45.0),
            (x: -22.869140625, y: 43.80817468459856),
            (x: -24.451171875, y: 35.266685523707665)
        ];
        let poly = Polygon::new(outer.clone(), vec![inner]);
        let simplified = poly.simplify_vw_preserve(&95.4);
        assert_eq!(simplified.exterior(), &outer);
        assert_eq!(simplified.interiors()[0], correct_inner);
    }
    #[test]
    fn very_long_vwp_test() {
        // simplify an 8k-point LineString, eliminating self-intersections
        let points_ls = geo_test_fixtures::norway_main::<f64>();
        let simplified = vwp_wrapper::<_, 2, 4>(&points_ls, None, &0.0005);
        assert_eq!(simplified[0].len(), 3278);
    }

    #[test]
    fn visvalingam_test_long() {
        // simplify a longer LineString
        let points_ls = geo_test_fixtures::vw_orig::<f64>();
        let correct_ls = geo_test_fixtures::vw_simplified::<f64>();
        let simplified = visvalingam(&points_ls, &0.0005);
        assert_eq!(simplified, correct_ls.0);
    }
    #[test]
    fn visvalingam_preserve_test_long() {
        // simplify a longer LineString using the preserve variant
        let points_ls = geo_test_fixtures::vw_orig::<f64>();
        let correct_ls = geo_test_fixtures::vw_simplified::<f64>();
        let simplified = points_ls.simplify_vw_preserve(&0.0005);
        assert_relative_eq!(simplified, correct_ls, epsilon = 1e-6);
    }
    #[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 vec = vec![Point::new(0.0, 0.0), Point::new(27.8, 0.1)];
        let compare = vec![Coord::from((0.0, 0.0)), Coord::from((27.8, 0.1))];
        let simplified = visvalingam(&LineString::from(vec), &1.0);
        assert_eq!(simplified, compare);
    }

    #[test]
    fn multilinestring() {
        // this is the PostGIS example
        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::new(vec![LineString::from(points_ls)]);
        assert_relative_eq!(
            mline.simplify_vw(&30.),
            MultiLineString::new(vec![LineString::from(correct_ls)]),
            epsilon = 1e-6
        );
    }

    #[test]
    fn polygon() {
        let poly = polygon![
            (x: 0., y: 0.),
            (x: 0., y: 10.),
            (x: 5., y: 11.),
            (x: 10., y: 10.),
            (x: 10., y: 0.),
            (x: 0., y: 0.),
        ];

        let poly2 = poly.simplify_vw(&10.);

        assert_relative_eq!(
            poly2,
            polygon![
                (x: 0., y: 0.),
                (x: 0., y: 10.),
                (x: 10., y: 10.),
                (x: 10., y: 0.),
                (x: 0., y: 0.),
            ],
            epsilon = 1e-6
        );
    }

    #[test]
    fn multipolygon() {
        let mpoly = MultiPolygon::new(vec![Polygon::new(
            LineString::from(vec![
                (0., 0.),
                (0., 10.),
                (5., 11.),
                (10., 10.),
                (10., 0.),
                (0., 0.),
            ]),
            vec![],
        )]);

        let mpoly2 = mpoly.simplify_vw(&10.);

        assert_relative_eq!(
            mpoly2,
            MultiPolygon::new(vec![Polygon::new(
                LineString::from(vec![(0., 0.), (0., 10.), (10., 10.), (10., 0.), (0., 0.)]),
                vec![],
            )]),
            epsilon = 1e-6
        );
    }
}