From d03e58b0cf39d60c91641452810f3b9aee775291 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Fri, 24 Mar 2017 13:57:08 +0000 Subject: [PATCH 01/16] First pass at extreme-finding algorithms --- src/algorithm/extremes.rs | 265 ++++++++++++++++++++++++++++++++++++++ src/algorithm/mod.rs | 2 + src/types.rs | 25 ++++ 3 files changed, 292 insertions(+) create mode 100644 src/algorithm/extremes.rs diff --git a/src/algorithm/extremes.rs b/src/algorithm/extremes.rs new file mode 100644 index 000000000..19d12e279 --- /dev/null +++ b/src/algorithm/extremes.rs @@ -0,0 +1,265 @@ +use num_traits::Float; +use types::{Point, LineString, Polygon}; +use algorithm::convexhull::ConvexHull; +use types::Extremes; + +// Useful direction vectors: +// 1., 0. = largest x +// 0., 1. = largest y +// 0., -1. = smallest y +// -1, 0. = smallest x + +// various tests for vector orientation relative to a direction vector u +fn up(u: &Point, v: &Point) -> bool + where T: Float +{ + u.dot(v) > T::zero() +} + +fn direction_sign(u: &Point, vi: &Point, vj: &Point) -> T + where T: Float +{ + u.dot(&(*vi - *vj)) +} + +// true if Vi is above Vj +fn above(u: &Point, vi: &Point, vj: &Point) -> bool + where T: Float +{ + direction_sign(u, vi, vj) > T::zero() +} + +// true if Vi is below Vj +fn below(u: &Point, vi: &Point, vj: &Point) -> bool + where T: Float +{ + direction_sign(u, vi, vj) < T::zero() +} + +// wrapper for extreme-finding function +fn find_extremes(func: F, polygon: &Polygon, convex: bool, oriented: bool) -> Extremes + where T: Float, + F: Fn(&Point, &Polygon) -> Result, ()> +{ + // TODO: we can't use this until Orient lands + // let mut processed = false; + // if !convex { + // let mut poly = polygon.convex_hull(); + // let processed = true; + // } + // if !oriented and processed { + // poly = poly.orient() + // } else if !oriented and !processed { + // poly = polygon.orient(); + // } else { + // poly = polygon; + // } + let directions = vec![Point::new(T::zero(), -T::one()), + Point::new(T::one(), T::zero()), + Point::new(T::zero(), T::one()), + Point::new(-T::one(), T::zero())]; + directions + .iter() + .map(|p| func(&p, &polygon).unwrap()) + .collect::>>() + .into() +} + +// find a convex, counter-clockwise oriented polygon's maximum vertex in a specified direction +// u: a direction vector. We're using a point to represent this, which is a hack tbh +// this is O(n), because polymax() can't yet calculate minimum x +fn polymax_naive(u: &Point, poly: &Polygon) -> Result, ()> + where T: Float +{ + let vertices = &poly.exterior.0; + let mut max: usize = 0; + for (i, _) in vertices.iter().enumerate() { + // if vertices[i] is above prior vertices[max] + if above(u, &vertices[i], &vertices[max]) { + max = i; + } + } + return Ok(vertices[max]); +} + +// ported from a c++ implementation at +// http://geomalgorithms.com/a14-_extreme_pts.html#polyMax_2D() +// original copyright notice: +// Copyright 2002 softSurfer, 2012-2013 Dan Sunday +// This code may be freely used, distributed and modified for any +// purpose providing that this copyright notice is included with it. +// SoftSurfer makes no warranty for this code, and cannot be held +// liable for any real or imagined damage resulting from its use. +// Users of this code must verify correctness for their application. + +// Original implementation: +// Joseph O'Rourke, Computational Geometry in C (2nd Edition), +// Sect 7.9 "Extreme Point of Convex Polygon" (1998) + +// find a convex, counter-clockwise oriented polygon's maximum vertex in a specified direction +// u: a direction vector. We're using a point to represent this, which is a hack tbh +// this should run in O(log n) time. +// This implementation can't calculate minimum x +// see Section 5.1 at http://codeforces.com/blog/entry/48868 for a discussion +fn polymax(u: &Point, poly: &Polygon) -> Result, ()> + where T: Float +{ + let vertices = &poly.exterior.0; + let n = vertices.len(); + // these are used to divide the vertices slice + // start chain = [0, n] with vertices[n] = vertices[0] + let mut start: usize = 0; + let mut end: usize = n; + let mut mid: usize; + + // edge vectors at vertices[a] and vertices[c] + let mut vec_c: Point; + let vec_a = vertices[1] - vertices[0]; + // test for "up" direction of vec_a + let mut up_a = up(u, &vec_a); + + // test if vertices[0] is a local maximum + if !up_a && !above(u, &vertices[n - 1], &vertices[0]) { + return Ok(vertices[0]); + } + loop { + mid = (start + end) / 2; + vec_c = vertices[mid + 1] - vertices[mid]; + let up_c = up(u, &vec_c); + if !up_c && !above(u, &vertices[mid - 1], &vertices[mid]) { + // vertices[mid] is a local maximum, thus it is a maximum + return Ok(vertices[mid]); + } + // no max yet, so continue with the binary search + // pick one of the two subchains [start, mid] or [mid, end] + if up_a { + // vec_a points up + if !up_c { + // vec_c points down + end = mid; // select [start, mid] + } else { + // vec_c points up + if above(u, &vertices[start], &vertices[mid]) { + // vertices[start] is above vertices[mid] + end = mid; // select [start, mid] + } else { + // vertices[start] is below vertices[mid] + start = mid; // select [mid, end] + up_a = up_c; + } + } + } else { + // vec_a points down + if up_c { + // vec_c points up + start = mid; // select [mid, end] + up_a = up_c; + } else { + // vec_c points down + if below(u, &vertices[start], &vertices[mid]) { + // vertices[start] is below vertices[mid] + end = mid; // select [start, mid] + } else { + // vertices[start] is above vertices[mid] + start = mid; // select [mid, end] + up_a = up_c; + } + } + } + // something went really badly wrong + if end <= start + 1 { + return Err(()); + } + } +} + +pub trait ExtremePoints { + /// Find the extreme `x` and `y` points of a Polygon + /// + /// The polygon must be convex and properly oriented; if you're unsure whether + /// the polygon has these properties: + /// + /// - If you aren't sure whether the polygon is convex, choose `convex=false, oriented=false` + /// - If the polygon is convex but oriented clockwise, choose `convex=true, oriented=false` + /// + /// Convex-hull processing is `O(n log(n))` on average + /// and is thus an upper bound on point-finding, which is otherwise `O(n)` for a `Polygon`. + /// For a `MultiPolygon`, its convex hull must always be calculated first. + /// + /// ``` + /// use geo::{Point, LineString, Polygon}; + /// use geo::extremes::ExtremePoints; + /// // a diamond shape + /// let points_raw = vec![(1.0, 0.0), (2.0, 1.0), (1.0, 2.0), (0.0, 1.0), (1.0, 0.0)]; + /// let points = points_raw.iter().map(|e| Point::new(e.0, e.1)).collect::>(); + /// let poly = Polygon::new(LineString(points), vec![]); + /// // Polygon is both convex and oriented counter-clockwise + /// let extremes = poly.extreme_points(true, true); + /// assert_eq!(extremes.ymin, Point::new(1.0, 0.0)); + /// assert_eq!(extremes.xmax, Point::new(2.0, 1.0)); + /// assert_eq!(extremes.ymax, Point::new(1.0, 2.0)); + /// assert_eq!(extremes.xmin, Point::new(0.0, 1.0)); + /// ``` + fn extreme_points(&self, convex: bool, oriented: bool) -> Extremes; +} + +impl ExtremePoints for Polygon + where T: Float +{ + fn extreme_points(&self, convex: bool, oriented: bool) -> Extremes { + find_extremes(polymax_naive, self, convex, oriented) + } +} + +#[cfg(test)] +mod test { + use types::Point; + use super::*; + #[test] + fn test_polygon_extreme_x() { + // a diamond shape + let points_raw = vec![(1.0, 0.0), (2.0, 1.0), (1.0, 2.0), (0.0, 1.0), (1.0, 0.0)]; + let points = points_raw + .iter() + .map(|e| Point::new(e.0, e.1)) + .collect::>(); + let poly1 = Polygon::new(LineString(points), vec![]); + let min_x = polymax_naive(&Point::new(-1., 0.), &poly1).unwrap(); + let correct = Point::new(0., 1.); + assert_eq!(min_x, correct); + } + #[test] + #[should_panic] + // this test should panic, because the algorithm can't find minimum x + fn test_polygon_extreme_x_fast() { + // a diamond shape + let points_raw = vec![(1.0, 0.0), (2.0, 1.0), (1.0, 2.0), (0.0, 1.0), (1.0, 0.0)]; + let points = points_raw + .iter() + .map(|e| Point::new(e.0, e.1)) + .collect::>(); + let poly1 = Polygon::new(LineString(points), vec![]); + let min_x = polymax(&Point::new(-1., 0.), &poly1).unwrap(); + let correct = Point::new(0., 1.); + assert_eq!(min_x, correct); + } + #[test] + fn test_polygon_extreme_wrapper() { + // a diamond shape with a bump on the top-right edge + let points_raw = + vec![(1.0, 0.0), (2.0, 1.0), (1.75, 1.75), (1.0, 2.0), (0.0, 1.0), (1.0, 0.0)]; + let points = points_raw + .iter() + .map(|e| Point::new(e.0, e.1)) + .collect::>(); + let poly1 = Polygon::new(LineString(points), vec![]); + let extremes = find_extremes(polymax_naive, &poly1, true, true); + let correct = Extremes { + ymin: Point::new(1.0, 0.0), + xmax: Point::new(2.0, 1.0), + ymax: Point::new(1.0, 2.0), + xmin: Point::new(0.0, 1.0), + }; + assert_eq!(extremes, correct); + } +} diff --git a/src/algorithm/mod.rs b/src/algorithm/mod.rs index 23d32c5e4..060c63394 100644 --- a/src/algorithm/mod.rs +++ b/src/algorithm/mod.rs @@ -22,3 +22,5 @@ pub mod simplifyvw; pub mod convexhull; /// Orients a Polygon's exterior and interior rings pub mod orient; +/// calculate Polygon extremes +pub mod extremes; diff --git a/src/types.rs b/src/types.rs index 2a13573e2..7a32393b5 100644 --- a/src/types.rs +++ b/src/types.rs @@ -25,6 +25,31 @@ pub struct Bbox pub ymax: T, } +/// Returns the extreme points (minimum `y`, maximum `x`, maximum `y`, minimum `x`) +/// of a `Polygon` or `MultiPolygon`. +#[derive(PartialEq, Clone, Copy, Debug)] +pub struct Extremes + where T: Float +{ + pub ymin: Point, + pub xmax: Point, + pub ymax: Point, + pub xmin: Point, +} + +impl From>> for Extremes + where T: Float +{ + fn from(original: Vec>) -> Extremes { + Extremes { + ymin: original[0], + xmax: original[1], + ymax: original[2], + xmin: original[3], + } + } +} + #[derive(PartialEq, Clone, Copy, Debug)] pub struct Point (pub Coordinate) where T: Float; From d50468af6538f9eb884e7642c1b372d3f4089d19 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Sun, 9 Apr 2017 22:02:19 +0100 Subject: [PATCH 02/16] Remove binary search implementation for now --- src/algorithm/extremes.rs | 110 +------------------------------------- 1 file changed, 2 insertions(+), 108 deletions(-) diff --git a/src/algorithm/extremes.rs b/src/algorithm/extremes.rs index 19d12e279..b05fabab9 100644 --- a/src/algorithm/extremes.rs +++ b/src/algorithm/extremes.rs @@ -47,9 +47,9 @@ fn find_extremes(func: F, polygon: &Polygon, convex: bool, oriented: bo // let mut poly = polygon.convex_hull(); // let processed = true; // } - // if !oriented and processed { + // if !oriented && processed { // poly = poly.orient() - // } else if !oriented and !processed { + // } else if !oriented && !processed { // poly = polygon.orient(); // } else { // poly = polygon; @@ -82,97 +82,6 @@ fn polymax_naive(u: &Point, poly: &Polygon) -> Result, ()> return Ok(vertices[max]); } -// ported from a c++ implementation at -// http://geomalgorithms.com/a14-_extreme_pts.html#polyMax_2D() -// original copyright notice: -// Copyright 2002 softSurfer, 2012-2013 Dan Sunday -// This code may be freely used, distributed and modified for any -// purpose providing that this copyright notice is included with it. -// SoftSurfer makes no warranty for this code, and cannot be held -// liable for any real or imagined damage resulting from its use. -// Users of this code must verify correctness for their application. - -// Original implementation: -// Joseph O'Rourke, Computational Geometry in C (2nd Edition), -// Sect 7.9 "Extreme Point of Convex Polygon" (1998) - -// find a convex, counter-clockwise oriented polygon's maximum vertex in a specified direction -// u: a direction vector. We're using a point to represent this, which is a hack tbh -// this should run in O(log n) time. -// This implementation can't calculate minimum x -// see Section 5.1 at http://codeforces.com/blog/entry/48868 for a discussion -fn polymax(u: &Point, poly: &Polygon) -> Result, ()> - where T: Float -{ - let vertices = &poly.exterior.0; - let n = vertices.len(); - // these are used to divide the vertices slice - // start chain = [0, n] with vertices[n] = vertices[0] - let mut start: usize = 0; - let mut end: usize = n; - let mut mid: usize; - - // edge vectors at vertices[a] and vertices[c] - let mut vec_c: Point; - let vec_a = vertices[1] - vertices[0]; - // test for "up" direction of vec_a - let mut up_a = up(u, &vec_a); - - // test if vertices[0] is a local maximum - if !up_a && !above(u, &vertices[n - 1], &vertices[0]) { - return Ok(vertices[0]); - } - loop { - mid = (start + end) / 2; - vec_c = vertices[mid + 1] - vertices[mid]; - let up_c = up(u, &vec_c); - if !up_c && !above(u, &vertices[mid - 1], &vertices[mid]) { - // vertices[mid] is a local maximum, thus it is a maximum - return Ok(vertices[mid]); - } - // no max yet, so continue with the binary search - // pick one of the two subchains [start, mid] or [mid, end] - if up_a { - // vec_a points up - if !up_c { - // vec_c points down - end = mid; // select [start, mid] - } else { - // vec_c points up - if above(u, &vertices[start], &vertices[mid]) { - // vertices[start] is above vertices[mid] - end = mid; // select [start, mid] - } else { - // vertices[start] is below vertices[mid] - start = mid; // select [mid, end] - up_a = up_c; - } - } - } else { - // vec_a points down - if up_c { - // vec_c points up - start = mid; // select [mid, end] - up_a = up_c; - } else { - // vec_c points down - if below(u, &vertices[start], &vertices[mid]) { - // vertices[start] is below vertices[mid] - end = mid; // select [start, mid] - } else { - // vertices[start] is above vertices[mid] - start = mid; // select [mid, end] - up_a = up_c; - } - } - } - // something went really badly wrong - if end <= start + 1 { - return Err(()); - } - } -} - pub trait ExtremePoints { /// Find the extreme `x` and `y` points of a Polygon /// @@ -229,21 +138,6 @@ mod test { assert_eq!(min_x, correct); } #[test] - #[should_panic] - // this test should panic, because the algorithm can't find minimum x - fn test_polygon_extreme_x_fast() { - // a diamond shape - let points_raw = vec![(1.0, 0.0), (2.0, 1.0), (1.0, 2.0), (0.0, 1.0), (1.0, 0.0)]; - let points = points_raw - .iter() - .map(|e| Point::new(e.0, e.1)) - .collect::>(); - let poly1 = Polygon::new(LineString(points), vec![]); - let min_x = polymax(&Point::new(-1., 0.), &poly1).unwrap(); - let correct = Point::new(0., 1.); - assert_eq!(min_x, correct); - } - #[test] fn test_polygon_extreme_wrapper() { // a diamond shape with a bump on the top-right edge let points_raw = From d89d1441da91298a594efe706e863d85b11e9e66 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Mon, 10 Apr 2017 00:48:20 +0100 Subject: [PATCH 03/16] Return indices, not points --- src/algorithm/extremes.rs | 34 +++++++++++++++++----------------- src/types.rs | 20 ++++++++------------ 2 files changed, 25 insertions(+), 29 deletions(-) diff --git a/src/algorithm/extremes.rs b/src/algorithm/extremes.rs index b05fabab9..2d70f56d0 100644 --- a/src/algorithm/extremes.rs +++ b/src/algorithm/extremes.rs @@ -37,9 +37,9 @@ fn below(u: &Point, vi: &Point, vj: &Point) -> bool } // wrapper for extreme-finding function -fn find_extremes(func: F, polygon: &Polygon, convex: bool, oriented: bool) -> Extremes +fn find_extremes(func: F, polygon: &Polygon, convex: bool, oriented: bool) -> Extremes where T: Float, - F: Fn(&Point, &Polygon) -> Result, ()> + F: Fn(&Point, &Polygon) -> Result { // TODO: we can't use this until Orient lands // let mut processed = false; @@ -61,14 +61,14 @@ fn find_extremes(func: F, polygon: &Polygon, convex: bool, oriented: bo directions .iter() .map(|p| func(&p, &polygon).unwrap()) - .collect::>>() + .collect::>() .into() } // find a convex, counter-clockwise oriented polygon's maximum vertex in a specified direction // u: a direction vector. We're using a point to represent this, which is a hack tbh // this is O(n), because polymax() can't yet calculate minimum x -fn polymax_naive(u: &Point, poly: &Polygon) -> Result, ()> +fn polymax_naive(u: &Point, poly: &Polygon) -> Result where T: Float { let vertices = &poly.exterior.0; @@ -79,11 +79,11 @@ fn polymax_naive(u: &Point, poly: &Polygon) -> Result, ()> max = i; } } - return Ok(vertices[max]); + return Ok(max); } pub trait ExtremePoints { - /// Find the extreme `x` and `y` points of a Polygon + /// Find the extreme `x` and `y` indices of a Polygon /// /// The polygon must be convex and properly oriented; if you're unsure whether /// the polygon has these properties: @@ -104,18 +104,18 @@ pub trait ExtremePoints { /// let poly = Polygon::new(LineString(points), vec![]); /// // Polygon is both convex and oriented counter-clockwise /// let extremes = poly.extreme_points(true, true); - /// assert_eq!(extremes.ymin, Point::new(1.0, 0.0)); - /// assert_eq!(extremes.xmax, Point::new(2.0, 1.0)); - /// assert_eq!(extremes.ymax, Point::new(1.0, 2.0)); - /// assert_eq!(extremes.xmin, Point::new(0.0, 1.0)); + /// assert_eq!(extremes.ymin, 0); + /// assert_eq!(extremes.xmax, 1); + /// assert_eq!(extremes.ymax, 2); + /// assert_eq!(extremes.xmin, 3); /// ``` - fn extreme_points(&self, convex: bool, oriented: bool) -> Extremes; + fn extreme_points(&self, convex: bool, oriented: bool) -> Extremes; } impl ExtremePoints for Polygon where T: Float { - fn extreme_points(&self, convex: bool, oriented: bool) -> Extremes { + fn extreme_points(&self, convex: bool, oriented: bool) -> Extremes { find_extremes(polymax_naive, self, convex, oriented) } } @@ -134,7 +134,7 @@ mod test { .collect::>(); let poly1 = Polygon::new(LineString(points), vec![]); let min_x = polymax_naive(&Point::new(-1., 0.), &poly1).unwrap(); - let correct = Point::new(0., 1.); + let correct = 3_usize; assert_eq!(min_x, correct); } #[test] @@ -149,10 +149,10 @@ mod test { let poly1 = Polygon::new(LineString(points), vec![]); let extremes = find_extremes(polymax_naive, &poly1, true, true); let correct = Extremes { - ymin: Point::new(1.0, 0.0), - xmax: Point::new(2.0, 1.0), - ymax: Point::new(1.0, 2.0), - xmin: Point::new(0.0, 1.0), + ymin: 0, + xmax: 1, + ymax: 3, + xmin: 4, }; assert_eq!(extremes, correct); } diff --git a/src/types.rs b/src/types.rs index 7a32393b5..ad37fee10 100644 --- a/src/types.rs +++ b/src/types.rs @@ -25,22 +25,18 @@ pub struct Bbox pub ymax: T, } -/// Returns the extreme points (minimum `y`, maximum `x`, maximum `y`, minimum `x`) +/// Returns the extreme vertex indices (minimum `y`, maximum `x`, maximum `y`, minimum `x`) /// of a `Polygon` or `MultiPolygon`. #[derive(PartialEq, Clone, Copy, Debug)] -pub struct Extremes - where T: Float -{ - pub ymin: Point, - pub xmax: Point, - pub ymax: Point, - pub xmin: Point, +pub struct Extremes { + pub ymin: usize, + pub xmax: usize, + pub ymax: usize, + pub xmin: usize, } -impl From>> for Extremes - where T: Float -{ - fn from(original: Vec>) -> Extremes { +impl From> for Extremes { + fn from(original: Vec) -> Extremes { Extremes { ymin: original[0], xmax: original[1], From 4f5df74f80edc42d3d3a759fcc91d3ffe07c8593 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Mon, 10 Apr 2017 11:35:08 +0100 Subject: [PATCH 04/16] Simplify matching, remove allocations, add another test --- src/algorithm/extremes.rs | 77 ++++++++++++++++++++++++++++----------- 1 file changed, 55 insertions(+), 22 deletions(-) diff --git a/src/algorithm/extremes.rs b/src/algorithm/extremes.rs index 2d70f56d0..d8dea8ae0 100644 --- a/src/algorithm/extremes.rs +++ b/src/algorithm/extremes.rs @@ -1,6 +1,7 @@ use num_traits::Float; use types::{Point, LineString, Polygon}; use algorithm::convexhull::ConvexHull; +use algorithm::orient::{Orient, Direction}; use types::Extremes; // Useful direction vectors: @@ -41,32 +42,37 @@ fn find_extremes(func: F, polygon: &Polygon, convex: bool, oriented: bo where T: Float, F: Fn(&Point, &Polygon) -> Result { - // TODO: we can't use this until Orient lands - // let mut processed = false; - // if !convex { - // let mut poly = polygon.convex_hull(); - // let processed = true; - // } - // if !oriented && processed { - // poly = poly.orient() - // } else if !oriented && !processed { - // poly = polygon.orient(); - // } else { - // poly = polygon; - // } let directions = vec![Point::new(T::zero(), -T::one()), Point::new(T::one(), T::zero()), Point::new(T::zero(), T::one()), Point::new(-T::one(), T::zero())]; - directions - .iter() - .map(|p| func(&p, &polygon).unwrap()) - .collect::>() - .into() + match (convex, oriented) { + (false, _) => { + directions + .iter() + .map(|p| func(&p, &polygon.convex_hull()).unwrap()) + .collect::>() + .into() + } + (true, false) => { + directions + .iter() + .map(|p| func(&p, &polygon.orient(Direction::Default)).unwrap()) + .collect::>() + .into() + } + _ => { + directions + .iter() + .map(|p| func(&p, &polygon).unwrap()) + .collect::>() + .into() + } + } } // find a convex, counter-clockwise oriented polygon's maximum vertex in a specified direction -// u: a direction vector. We're using a point to represent this, which is a hack tbh +// u: a direction vector. We're using a point to represent this, which is a hack but works fine // this is O(n), because polymax() can't yet calculate minimum x fn polymax_naive(u: &Point, poly: &Polygon) -> Result where T: Float @@ -139,15 +145,42 @@ mod test { } #[test] fn test_polygon_extreme_wrapper() { - // a diamond shape with a bump on the top-right edge - let points_raw = + // non-convex, with a bump on the top-right edge + let points_raw = vec![(1.0, 0.0), + (1.3, 1.), + (2.0, 1.0), + (1.75, 1.75), + (1.0, 2.0), + (0.0, 1.0), + (1.0, 0.0)]; + let points = points_raw + .iter() + .map(|e| Point::new(e.0, e.1)) + .collect::>(); + let poly1 = Polygon::new(LineString(points), vec![]); + let extremes = find_extremes(polymax_naive, &poly1, false, false); + let correct = Extremes { + ymin: 0, + xmax: 1, + ymax: 3, + xmin: 4, + }; + assert_eq!(extremes, correct); + } + #[test] + fn test_polygon_extreme_wrapper_convex() { + // convex, with a bump on the top-right edge + let mut points_raw = vec![(1.0, 0.0), (2.0, 1.0), (1.75, 1.75), (1.0, 2.0), (0.0, 1.0), (1.0, 0.0)]; + // orient the vector clockwise + points_raw.reverse(); let points = points_raw .iter() .map(|e| Point::new(e.0, e.1)) .collect::>(); let poly1 = Polygon::new(LineString(points), vec![]); - let extremes = find_extremes(polymax_naive, &poly1, true, true); + // specify convexity, wrong orientation + let extremes = find_extremes(polymax_naive, &poly1, true, false); let correct = Extremes { ymin: 0, xmax: 1, From 392f38aecdb22016600c5892c82e820d168cddc6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Mon, 10 Apr 2017 11:57:29 +0100 Subject: [PATCH 05/16] Add MultiPolygon and MultiPoint impls --- src/algorithm/extremes.rs | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/src/algorithm/extremes.rs b/src/algorithm/extremes.rs index d8dea8ae0..d2d471f9a 100644 --- a/src/algorithm/extremes.rs +++ b/src/algorithm/extremes.rs @@ -1,5 +1,5 @@ use num_traits::Float; -use types::{Point, LineString, Polygon}; +use types::{Point, LineString, Polygon, MultiPoint, MultiPolygon}; use algorithm::convexhull::ConvexHull; use algorithm::orient::{Orient, Direction}; use types::Extremes; @@ -73,7 +73,7 @@ fn find_extremes(func: F, polygon: &Polygon, convex: bool, oriented: bo // find a convex, counter-clockwise oriented polygon's maximum vertex in a specified direction // u: a direction vector. We're using a point to represent this, which is a hack but works fine -// this is O(n), because polymax() can't yet calculate minimum x +// this is O(n) if the polygon is convex. fn polymax_naive(u: &Point, poly: &Polygon) -> Result where T: Float { @@ -126,6 +126,24 @@ impl ExtremePoints for Polygon } } +impl ExtremePoints for MultiPolygon + where T: Float +{ + fn extreme_points(&self, convex: bool, oriented: bool) -> Extremes { + // we can disregard the input because convex-hull processing always orients + find_extremes(polymax_naive, &self.convex_hull(), true, true) + } +} + +impl ExtremePoints for MultiPoint + where T: Float +{ + fn extreme_points(&self, convex: bool, oriented: bool) -> Extremes { + // we can disregard the input because convex-hull processing always orients + find_extremes(polymax_naive, &self.convex_hull(), true, true) + } +} + #[cfg(test)] mod test { use types::Point; From 32d789ddaaee47867099912086d99aae46f02e48 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Mon, 10 Apr 2017 12:08:19 +0100 Subject: [PATCH 06/16] Better docs --- src/algorithm/mod.rs | 4 ++-- src/types.rs | 2 -- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/algorithm/mod.rs b/src/algorithm/mod.rs index 060c63394..da1483d79 100644 --- a/src/algorithm/mod.rs +++ b/src/algorithm/mod.rs @@ -20,7 +20,7 @@ pub mod simplify; pub mod simplifyvw; /// Calculates the convex hull of a geometry. pub mod convexhull; -/// Orients a Polygon's exterior and interior rings +/// Orients a Polygon's exterior and interior rings. pub mod orient; -/// calculate Polygon extremes +/// Returns the extreme indices of a `Polygon`, `MultiPolygon`, or `MultiPoint`. pub mod extremes; diff --git a/src/types.rs b/src/types.rs index ad37fee10..85c8cc345 100644 --- a/src/types.rs +++ b/src/types.rs @@ -25,8 +25,6 @@ pub struct Bbox pub ymax: T, } -/// Returns the extreme vertex indices (minimum `y`, maximum `x`, maximum `y`, minimum `x`) -/// of a `Polygon` or `MultiPolygon`. #[derive(PartialEq, Clone, Copy, Debug)] pub struct Extremes { pub ymin: usize, From f2f2218ae25de7908c32c2f965a0f56ef9b5908e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Sat, 15 Apr 2017 00:26:30 +0100 Subject: [PATCH 07/16] We're returning indices, so make that clear --- src/algorithm/extremes.rs | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/algorithm/extremes.rs b/src/algorithm/extremes.rs index d2d471f9a..59e523241 100644 --- a/src/algorithm/extremes.rs +++ b/src/algorithm/extremes.rs @@ -109,19 +109,19 @@ pub trait ExtremePoints { /// let points = points_raw.iter().map(|e| Point::new(e.0, e.1)).collect::>(); /// let poly = Polygon::new(LineString(points), vec![]); /// // Polygon is both convex and oriented counter-clockwise - /// let extremes = poly.extreme_points(true, true); + /// let extremes = poly.extreme_indices(true, true); /// assert_eq!(extremes.ymin, 0); /// assert_eq!(extremes.xmax, 1); /// assert_eq!(extremes.ymax, 2); /// assert_eq!(extremes.xmin, 3); /// ``` - fn extreme_points(&self, convex: bool, oriented: bool) -> Extremes; + fn extreme_indices(&self, convex: bool, oriented: bool) -> Extremes; } impl ExtremePoints for Polygon where T: Float { - fn extreme_points(&self, convex: bool, oriented: bool) -> Extremes { + fn extreme_indices(&self, convex: bool, oriented: bool) -> Extremes { find_extremes(polymax_naive, self, convex, oriented) } } @@ -129,7 +129,7 @@ impl ExtremePoints for Polygon impl ExtremePoints for MultiPolygon where T: Float { - fn extreme_points(&self, convex: bool, oriented: bool) -> Extremes { + fn extreme_indices(&self, convex: bool, oriented: bool) -> Extremes { // we can disregard the input because convex-hull processing always orients find_extremes(polymax_naive, &self.convex_hull(), true, true) } @@ -138,7 +138,7 @@ impl ExtremePoints for MultiPolygon impl ExtremePoints for MultiPoint where T: Float { - fn extreme_points(&self, convex: bool, oriented: bool) -> Extremes { + fn extreme_indices(&self, convex: bool, oriented: bool) -> Extremes { // we can disregard the input because convex-hull processing always orients find_extremes(polymax_naive, &self.convex_hull(), true, true) } From 9e6cf9bc0e643c8033d4a7a6a0c4d1c10a26c2a6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Sat, 15 Apr 2017 00:29:28 +0100 Subject: [PATCH 08/16] Rename trait --- src/algorithm/extremes.rs | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/algorithm/extremes.rs b/src/algorithm/extremes.rs index 59e523241..c14ef472a 100644 --- a/src/algorithm/extremes.rs +++ b/src/algorithm/extremes.rs @@ -88,7 +88,7 @@ fn polymax_naive(u: &Point, poly: &Polygon) -> Result return Ok(max); } -pub trait ExtremePoints { +pub trait ExtremeIndices { /// Find the extreme `x` and `y` indices of a Polygon /// /// The polygon must be convex and properly oriented; if you're unsure whether @@ -103,7 +103,7 @@ pub trait ExtremePoints { /// /// ``` /// use geo::{Point, LineString, Polygon}; - /// use geo::extremes::ExtremePoints; + /// use geo::extremes::ExtremeIndices; /// // a diamond shape /// let points_raw = vec![(1.0, 0.0), (2.0, 1.0), (1.0, 2.0), (0.0, 1.0), (1.0, 0.0)]; /// let points = points_raw.iter().map(|e| Point::new(e.0, e.1)).collect::>(); @@ -118,7 +118,7 @@ pub trait ExtremePoints { fn extreme_indices(&self, convex: bool, oriented: bool) -> Extremes; } -impl ExtremePoints for Polygon +impl ExtremeIndices for Polygon where T: Float { fn extreme_indices(&self, convex: bool, oriented: bool) -> Extremes { @@ -126,7 +126,7 @@ impl ExtremePoints for Polygon } } -impl ExtremePoints for MultiPolygon +impl ExtremeIndices for MultiPolygon where T: Float { fn extreme_indices(&self, convex: bool, oriented: bool) -> Extremes { @@ -135,7 +135,7 @@ impl ExtremePoints for MultiPolygon } } -impl ExtremePoints for MultiPoint +impl ExtremeIndices for MultiPoint where T: Float { fn extreme_indices(&self, convex: bool, oriented: bool) -> Extremes { From 8c959bda8814a9734e41360cc949f9f0bbb6a9df Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Sat, 15 Apr 2017 23:55:58 +0100 Subject: [PATCH 09/16] Fix comment --- src/algorithm/extremes.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/src/algorithm/extremes.rs b/src/algorithm/extremes.rs index c14ef472a..099cc3517 100644 --- a/src/algorithm/extremes.rs +++ b/src/algorithm/extremes.rs @@ -73,7 +73,6 @@ fn find_extremes(func: F, polygon: &Polygon, convex: bool, oriented: bo // find a convex, counter-clockwise oriented polygon's maximum vertex in a specified direction // u: a direction vector. We're using a point to represent this, which is a hack but works fine -// this is O(n) if the polygon is convex. fn polymax_naive(u: &Point, poly: &Polygon) -> Result where T: Float { From c29a01ed0338c8146585229162c2b5094a607cb8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Sun, 16 Apr 2017 00:15:08 +0100 Subject: [PATCH 10/16] Better function names --- src/algorithm/extremes.rs | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/algorithm/extremes.rs b/src/algorithm/extremes.rs index 099cc3517..3250640d8 100644 --- a/src/algorithm/extremes.rs +++ b/src/algorithm/extremes.rs @@ -38,7 +38,7 @@ fn below(u: &Point, vi: &Point, vj: &Point) -> bool } // wrapper for extreme-finding function -fn find_extremes(func: F, polygon: &Polygon, convex: bool, oriented: bool) -> Extremes +fn find_extreme_indices(func: F, polygon: &Polygon, convex: bool, oriented: bool) -> Extremes where T: Float, F: Fn(&Point, &Polygon) -> Result { @@ -73,7 +73,7 @@ fn find_extremes(func: F, polygon: &Polygon, convex: bool, oriented: bo // find a convex, counter-clockwise oriented polygon's maximum vertex in a specified direction // u: a direction vector. We're using a point to represent this, which is a hack but works fine -fn polymax_naive(u: &Point, poly: &Polygon) -> Result +fn polymax_naive_indices(u: &Point, poly: &Polygon) -> Result where T: Float { let vertices = &poly.exterior.0; @@ -121,7 +121,7 @@ impl ExtremeIndices for Polygon where T: Float { fn extreme_indices(&self, convex: bool, oriented: bool) -> Extremes { - find_extremes(polymax_naive, self, convex, oriented) + find_extreme_indices(polymax_naive_indices, self, convex, oriented) } } @@ -130,7 +130,7 @@ impl ExtremeIndices for MultiPolygon { fn extreme_indices(&self, convex: bool, oriented: bool) -> Extremes { // we can disregard the input because convex-hull processing always orients - find_extremes(polymax_naive, &self.convex_hull(), true, true) + find_extreme_indices(polymax_naive_indices, &self.convex_hull(), true, true) } } @@ -139,7 +139,7 @@ impl ExtremeIndices for MultiPoint { fn extreme_indices(&self, convex: bool, oriented: bool) -> Extremes { // we can disregard the input because convex-hull processing always orients - find_extremes(polymax_naive, &self.convex_hull(), true, true) + find_extreme_indices(polymax_naive_indices, &self.convex_hull(), true, true) } } @@ -156,7 +156,7 @@ mod test { .map(|e| Point::new(e.0, e.1)) .collect::>(); let poly1 = Polygon::new(LineString(points), vec![]); - let min_x = polymax_naive(&Point::new(-1., 0.), &poly1).unwrap(); + let min_x = polymax_naive_indices(&Point::new(-1., 0.), &poly1).unwrap(); let correct = 3_usize; assert_eq!(min_x, correct); } @@ -175,7 +175,7 @@ mod test { .map(|e| Point::new(e.0, e.1)) .collect::>(); let poly1 = Polygon::new(LineString(points), vec![]); - let extremes = find_extremes(polymax_naive, &poly1, false, false); + let extremes = find_extreme_indices(polymax_naive_indices, &poly1, false, false); let correct = Extremes { ymin: 0, xmax: 1, @@ -197,7 +197,7 @@ mod test { .collect::>(); let poly1 = Polygon::new(LineString(points), vec![]); // specify convexity, wrong orientation - let extremes = find_extremes(polymax_naive, &poly1, true, false); + let extremes = find_extreme_indices(polymax_naive_indices, &poly1, true, false); let correct = Extremes { ymin: 0, xmax: 1, From 1ec7e715fc05c1cd220dcbae2554dae204b3a75a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Sun, 16 Apr 2017 00:37:37 +0100 Subject: [PATCH 11/16] Simplify ExtremeIndices trait since this trait can't return the indices of a non-convex polygon, there's no need to specify what to do with an input polygon This makes the trait publicly useless, since there's no simple way of guaranteeing the convex property of the input geometry. However, the trait is still useful as a helper for returning extreme points, as an internal helper for returning extreme indices for polygon-to-polygon distance there's no simple way to restrict a trait's visibility until 1.18, so we're currently back to the drawing board --- src/algorithm/extremes.rs | 67 +++++++++++---------------------------- 1 file changed, 18 insertions(+), 49 deletions(-) diff --git a/src/algorithm/extremes.rs b/src/algorithm/extremes.rs index 3250640d8..d7edeab0a 100644 --- a/src/algorithm/extremes.rs +++ b/src/algorithm/extremes.rs @@ -38,7 +38,7 @@ fn below(u: &Point, vi: &Point, vj: &Point) -> bool } // wrapper for extreme-finding function -fn find_extreme_indices(func: F, polygon: &Polygon, convex: bool, oriented: bool) -> Extremes +fn find_extreme_indices(func: F, polygon: &Polygon) -> Extremes where T: Float, F: Fn(&Point, &Polygon) -> Result { @@ -46,29 +46,11 @@ fn find_extreme_indices(func: F, polygon: &Polygon, convex: bool, orien Point::new(T::one(), T::zero()), Point::new(T::zero(), T::one()), Point::new(-T::one(), T::zero())]; - match (convex, oriented) { - (false, _) => { - directions - .iter() - .map(|p| func(&p, &polygon.convex_hull()).unwrap()) - .collect::>() - .into() - } - (true, false) => { - directions - .iter() - .map(|p| func(&p, &polygon.orient(Direction::Default)).unwrap()) - .collect::>() - .into() - } - _ => { - directions - .iter() - .map(|p| func(&p, &polygon).unwrap()) - .collect::>() - .into() - } - } + directions + .iter() + .map(|p| func(&p, &polygon).unwrap()) + .collect::>() + .into() } // find a convex, counter-clockwise oriented polygon's maximum vertex in a specified direction @@ -88,17 +70,9 @@ fn polymax_naive_indices(u: &Point, poly: &Polygon) -> Result { - /// Find the extreme `x` and `y` indices of a Polygon - /// - /// The polygon must be convex and properly oriented; if you're unsure whether - /// the polygon has these properties: - /// - /// - If you aren't sure whether the polygon is convex, choose `convex=false, oriented=false` - /// - If the polygon is convex but oriented clockwise, choose `convex=true, oriented=false` + /// Find the extreme `x` and `y` indices of a convex Polygon /// - /// Convex-hull processing is `O(n log(n))` on average - /// and is thus an upper bound on point-finding, which is otherwise `O(n)` for a `Polygon`. - /// For a `MultiPolygon`, its convex hull must always be calculated first. + /// The polygon **must be convex and properly (ccw) oriented**. /// /// ``` /// use geo::{Point, LineString, Polygon}; @@ -108,38 +82,36 @@ pub trait ExtremeIndices { /// let points = points_raw.iter().map(|e| Point::new(e.0, e.1)).collect::>(); /// let poly = Polygon::new(LineString(points), vec![]); /// // Polygon is both convex and oriented counter-clockwise - /// let extremes = poly.extreme_indices(true, true); + /// let extremes = poly.extreme_indices(); /// assert_eq!(extremes.ymin, 0); /// assert_eq!(extremes.xmax, 1); /// assert_eq!(extremes.ymax, 2); /// assert_eq!(extremes.xmin, 3); /// ``` - fn extreme_indices(&self, convex: bool, oriented: bool) -> Extremes; + fn extreme_indices(&self) -> Extremes; } impl ExtremeIndices for Polygon where T: Float { - fn extreme_indices(&self, convex: bool, oriented: bool) -> Extremes { - find_extreme_indices(polymax_naive_indices, self, convex, oriented) + fn extreme_indices(&self) -> Extremes { + find_extreme_indices(polymax_naive_indices, self) } } impl ExtremeIndices for MultiPolygon where T: Float { - fn extreme_indices(&self, convex: bool, oriented: bool) -> Extremes { - // we can disregard the input because convex-hull processing always orients - find_extreme_indices(polymax_naive_indices, &self.convex_hull(), true, true) + fn extreme_indices(&self) -> Extremes { + find_extreme_indices(polymax_naive_indices, &self.convex_hull()) } } impl ExtremeIndices for MultiPoint where T: Float { - fn extreme_indices(&self, convex: bool, oriented: bool) -> Extremes { - // we can disregard the input because convex-hull processing always orients - find_extreme_indices(polymax_naive_indices, &self.convex_hull(), true, true) + fn extreme_indices(&self) -> Extremes { + find_extreme_indices(polymax_naive_indices, &self.convex_hull()) } } @@ -175,7 +147,7 @@ mod test { .map(|e| Point::new(e.0, e.1)) .collect::>(); let poly1 = Polygon::new(LineString(points), vec![]); - let extremes = find_extreme_indices(polymax_naive_indices, &poly1, false, false); + let extremes = find_extreme_indices(polymax_naive_indices, &poly1.convex_hull()); let correct = Extremes { ymin: 0, xmax: 1, @@ -189,15 +161,12 @@ mod test { // convex, with a bump on the top-right edge let mut points_raw = vec![(1.0, 0.0), (2.0, 1.0), (1.75, 1.75), (1.0, 2.0), (0.0, 1.0), (1.0, 0.0)]; - // orient the vector clockwise - points_raw.reverse(); let points = points_raw .iter() .map(|e| Point::new(e.0, e.1)) .collect::>(); let poly1 = Polygon::new(LineString(points), vec![]); - // specify convexity, wrong orientation - let extremes = find_extreme_indices(polymax_naive_indices, &poly1, true, false); + let extremes = find_extreme_indices(polymax_naive_indices, &poly1.convex_hull()); let correct = Extremes { ymin: 0, xmax: 1, From 867c67903a8bd0c088cea71a1c19c291b6f442bd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Sun, 16 Apr 2017 01:30:02 +0100 Subject: [PATCH 12/16] Implement ExtremePoints trait Like ExtremeIndices, but for points --- src/algorithm/extremes.rs | 54 ++++++++++++++++++++++++++++++++++++++- src/types.rs | 10 ++++++++ 2 files changed, 63 insertions(+), 1 deletion(-) diff --git a/src/algorithm/extremes.rs b/src/algorithm/extremes.rs index d7edeab0a..afacf754d 100644 --- a/src/algorithm/extremes.rs +++ b/src/algorithm/extremes.rs @@ -2,7 +2,7 @@ use num_traits::Float; use types::{Point, LineString, Polygon, MultiPoint, MultiPolygon}; use algorithm::convexhull::ConvexHull; use algorithm::orient::{Orient, Direction}; -use types::Extremes; +use types::{Extremes, ExtremePoint}; // Useful direction vectors: // 1., 0. = largest x @@ -42,6 +42,7 @@ fn find_extreme_indices(func: F, polygon: &Polygon) -> Extremes where T: Float, F: Fn(&Point, &Polygon) -> Result { + // TODO: test for non-convexity, and return a Result let directions = vec![Point::new(T::zero(), -T::one()), Point::new(T::one(), T::zero()), Point::new(T::zero(), T::one()), @@ -115,6 +116,44 @@ impl ExtremeIndices for MultiPoint } } +pub trait ExtremePoints { + /// Find the extreme `x` and `y` points of a Geometry + /// + /// This trait is available to any struct implementing both `ConvexHull` amd `ExtremeIndices` + /// + /// ``` + /// use geo::{Point, LineString, Polygon}; + /// use geo::extremes::ExtremePoints; + /// let points_raw = vec![(1.0, 0.0), (2.0, 1.0), (1.0, 2.0), (0.0, 1.0), (1.0, 0.0)]; + /// let points = points_raw + /// .iter() + /// .map(|e| Point::new(e.0, e.1)) + /// .collect::>(); + /// let poly1 = Polygon::new(LineString(points), vec![]); + /// let extremes = poly1.extreme_points(); + /// let correct = Point::new(0.0, 1.0); + /// assert_eq!(extremes.xmin, correct); + /// ``` + fn extreme_points(&self) -> ExtremePoint; +} + +impl ExtremePoints for G + where T: Float, + G: ConvexHull + ExtremeIndices +{ + // Any Geometry implementing `ConvexHull` and `ExtremeIndices` gets this automatically + fn extreme_points(&self) -> ExtremePoint { + let ch = self.convex_hull(); + let indices = ch.extreme_indices(); + ExtremePoint { + ymin: ch.exterior.0[indices.ymin], + xmax: ch.exterior.0[indices.xmax], + ymax: ch.exterior.0[indices.ymax], + xmin: ch.exterior.0[indices.xmin], + } + } +} + #[cfg(test)] mod test { use types::Point; @@ -175,4 +214,17 @@ mod test { }; assert_eq!(extremes, correct); } + #[test] + fn test_polygon_extreme_point_x() { + // a diamond shape + let points_raw = vec![(1.0, 0.0), (2.0, 1.0), (1.0, 2.0), (0.0, 1.0), (1.0, 0.0)]; + let points = points_raw + .iter() + .map(|e| Point::new(e.0, e.1)) + .collect::>(); + let poly1 = Polygon::new(LineString(points), vec![]); + let extremes = poly1.extreme_points(); + let correct = Point::new(0.0, 1.0); + assert_eq!(extremes.xmin, correct); + } } diff --git a/src/types.rs b/src/types.rs index 85c8cc345..9e54037b3 100644 --- a/src/types.rs +++ b/src/types.rs @@ -44,6 +44,16 @@ impl From> for Extremes { } } +#[derive(PartialEq, Clone, Copy, Debug)] +pub struct ExtremePoint + where T: Float + { + pub ymin: Point, + pub xmax: Point, + pub ymax: Point, + pub xmin: Point, +} + #[derive(PartialEq, Clone, Copy, Debug)] pub struct Point (pub Coordinate) where T: Float; From 7d9b73a29559f4f53fce294d372c5bab262484e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Mon, 17 Apr 2017 00:23:02 +0100 Subject: [PATCH 13/16] Add convexity check to ExtremeIndices trait The input Polygon is now checked for convexity. This functionality adds considerable complexity to the trait, but I feel like it's worth it. --- src/algorithm/extremes.rs | 128 +++++++++++++++++++++++++++++++------- 1 file changed, 106 insertions(+), 22 deletions(-) diff --git a/src/algorithm/extremes.rs b/src/algorithm/extremes.rs index afacf754d..99ed85aef 100644 --- a/src/algorithm/extremes.rs +++ b/src/algorithm/extremes.rs @@ -1,10 +1,10 @@ -use num_traits::Float; +use num_traits::{Float, Signed}; use types::{Point, LineString, Polygon, MultiPoint, MultiPolygon}; use algorithm::convexhull::ConvexHull; use algorithm::orient::{Orient, Direction}; use types::{Extremes, ExtremePoint}; -// Useful direction vectors: +// Useful direction vectors, aligned with x and y axes: // 1., 0. = largest x // 0., 1. = largest y // 0., -1. = smallest y @@ -37,21 +37,79 @@ fn below(u: &Point, vi: &Point, vj: &Point) -> bool direction_sign(u, vi, vj) < T::zero() } +// used to check the sign of a vec of floats +#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)] +enum ListSign { + Empty, + Positive, + Negative, + Mixed, +} + +// Wrap-around previous-vertex +impl Polygon + where T: Float +{ + fn previous_vertex(&self, current_vertex: &usize) -> usize + where T: Float + { + (current_vertex + (self.exterior.0.len() - 1) - 1) % (self.exterior.0.len() - 1) + } +} + +// positive implies a -> b -> c is counter-clockwise, negative implies clockwise +pub fn cross_prod(p_a: &Point, p_b: &Point, p_c: &Point) -> T + where T: Float +{ + (p_b.x() - p_a.x()) * (p_c.y() - p_a.y()) - (p_b.y() - p_a.y()) * (p_c.x() - p_a.x()) +} + // wrapper for extreme-finding function -fn find_extreme_indices(func: F, polygon: &Polygon) -> Extremes - where T: Float, +fn find_extreme_indices(func: F, polygon: &Polygon) -> Result + where T: Float + Signed, F: Fn(&Point, &Polygon) -> Result { - // TODO: test for non-convexity, and return a Result + // For each consecutive pair of edges of the polygon (each triplet of points), + // compute the z-component of the cross product of the vectors defined by the + // edges pointing towards the points in increasing order. + // Take the cross product of these vectors + // The polygon is convex if the z-components of the cross products are either + // all positive or all negative. Otherwise, the polygon is non-convex. + // see: http://stackoverflow.com/a/1881201/416626 + let convex = polygon + .exterior + .0 + .iter() + .enumerate() + .map(|(idx, _)| { + // get previous and previous - 1 index offsets + let prev_1 = polygon.previous_vertex(&idx); + let prev_2 = polygon.previous_vertex(&prev_1); + // now make sure they're all positive or all negative, which implies convexity + cross_prod(&polygon.exterior.0[prev_2], + &polygon.exterior.0[prev_1], + &polygon.exterior.0[idx]) + }) + // positive implies ccw convexity, negative implies cw convexity + .fold(ListSign::Empty, |acc, n| { + match (acc, n.is_positive()) { + (ListSign::Empty, true) | (ListSign::Positive, true) => ListSign::Positive, + (ListSign::Empty, false) | (ListSign::Negative, false) => ListSign::Negative, + _ => ListSign::Mixed + } + }); + if convex == ListSign::Mixed { + return Err(()); + } let directions = vec![Point::new(T::zero(), -T::one()), Point::new(T::one(), T::zero()), Point::new(T::zero(), T::one()), Point::new(-T::one(), T::zero())]; - directions + Ok(directions .iter() .map(|p| func(&p, &polygon).unwrap()) .collect::>() - .into() + .into()) } // find a convex, counter-clockwise oriented polygon's maximum vertex in a specified direction @@ -70,7 +128,7 @@ fn polymax_naive_indices(u: &Point, poly: &Polygon) -> Result { +pub trait ExtremeIndices { /// Find the extreme `x` and `y` indices of a convex Polygon /// /// The polygon **must be convex and properly (ccw) oriented**. @@ -83,35 +141,35 @@ pub trait ExtremeIndices { /// let points = points_raw.iter().map(|e| Point::new(e.0, e.1)).collect::>(); /// let poly = Polygon::new(LineString(points), vec![]); /// // Polygon is both convex and oriented counter-clockwise - /// let extremes = poly.extreme_indices(); + /// let extremes = poly.extreme_indices().unwrap(); /// assert_eq!(extremes.ymin, 0); /// assert_eq!(extremes.xmax, 1); /// assert_eq!(extremes.ymax, 2); /// assert_eq!(extremes.xmin, 3); /// ``` - fn extreme_indices(&self) -> Extremes; + fn extreme_indices(&self) -> Result; } impl ExtremeIndices for Polygon - where T: Float + where T: Float + Signed { - fn extreme_indices(&self) -> Extremes { + fn extreme_indices(&self) -> Result { find_extreme_indices(polymax_naive_indices, self) } } impl ExtremeIndices for MultiPolygon - where T: Float + where T: Float + Signed { - fn extreme_indices(&self) -> Extremes { + fn extreme_indices(&self) -> Result { find_extreme_indices(polymax_naive_indices, &self.convex_hull()) } } impl ExtremeIndices for MultiPoint - where T: Float + where T: Float + Signed { - fn extreme_indices(&self) -> Extremes { + fn extreme_indices(&self) -> Result { find_extreme_indices(polymax_naive_indices, &self.convex_hull()) } } @@ -138,13 +196,14 @@ pub trait ExtremePoints { } impl ExtremePoints for G - where T: Float, + where T: Float + Signed, G: ConvexHull + ExtremeIndices { // Any Geometry implementing `ConvexHull` and `ExtremeIndices` gets this automatically fn extreme_points(&self) -> ExtremePoint { let ch = self.convex_hull(); - let indices = ch.extreme_indices(); + // safe to unwrap, since we're guaranteeing the polygon's convexity + let indices = ch.extreme_indices().unwrap(); ExtremePoint { ymin: ch.exterior.0[indices.ymin], xmax: ch.exterior.0[indices.xmax], @@ -172,7 +231,32 @@ mod test { assert_eq!(min_x, correct); } #[test] - fn test_polygon_extreme_wrapper() { + #[should_panic] + fn test_extreme_indices_bad_polygon() { + // non-convex, with a bump on the top-right edge + let points_raw = vec![(1.0, 0.0), + (1.3, 1.), + (2.0, 1.0), + (1.75, 1.75), + (1.0, 2.0), + (0.0, 1.0), + (1.0, 0.0)]; + let points = points_raw + .iter() + .map(|e| Point::new(e.0, e.1)) + .collect::>(); + let poly1 = Polygon::new(LineString(points), vec![]); + let extremes = find_extreme_indices(polymax_naive_indices, &poly1).unwrap(); + let correct = Extremes { + ymin: 0, + xmax: 1, + ymax: 3, + xmin: 4, + }; + assert_eq!(extremes, correct); + } + #[test] + fn test_extreme_indices_good_polygon() { // non-convex, with a bump on the top-right edge let points_raw = vec![(1.0, 0.0), (1.3, 1.), @@ -186,7 +270,7 @@ mod test { .map(|e| Point::new(e.0, e.1)) .collect::>(); let poly1 = Polygon::new(LineString(points), vec![]); - let extremes = find_extreme_indices(polymax_naive_indices, &poly1.convex_hull()); + let extremes = find_extreme_indices(polymax_naive_indices, &poly1.convex_hull()).unwrap(); let correct = Extremes { ymin: 0, xmax: 1, @@ -198,14 +282,14 @@ mod test { #[test] fn test_polygon_extreme_wrapper_convex() { // convex, with a bump on the top-right edge - let mut points_raw = + let points_raw = vec![(1.0, 0.0), (2.0, 1.0), (1.75, 1.75), (1.0, 2.0), (0.0, 1.0), (1.0, 0.0)]; let points = points_raw .iter() .map(|e| Point::new(e.0, e.1)) .collect::>(); let poly1 = Polygon::new(LineString(points), vec![]); - let extremes = find_extreme_indices(polymax_naive_indices, &poly1.convex_hull()); + let extremes = find_extreme_indices(polymax_naive_indices, &poly1.convex_hull()).unwrap(); let correct = Extremes { ymin: 0, xmax: 1, From 81ef980b301447fc0165c13786188615e7484daa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Mon, 17 Apr 2017 01:02:53 +0100 Subject: [PATCH 14/16] Better comments [ci skip] --- src/algorithm/extremes.rs | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/algorithm/extremes.rs b/src/algorithm/extremes.rs index 99ed85aef..d4f652283 100644 --- a/src/algorithm/extremes.rs +++ b/src/algorithm/extremes.rs @@ -82,15 +82,15 @@ fn find_extreme_indices(func: F, polygon: &Polygon) -> Result ListSign::Positive, @@ -106,10 +106,10 @@ fn find_extreme_indices(func: F, polygon: &Polygon) -> Result>() - .into()) + .iter() + .map(|p| func(&p, &polygon).unwrap()) + .collect::>() + .into()) } // find a convex, counter-clockwise oriented polygon's maximum vertex in a specified direction From 6a81d835ecc24ff44ef64c32830dd2e1fd80afd4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Mon, 17 Apr 2017 10:41:34 +0100 Subject: [PATCH 15/16] cross product doesn't need to be pub --- src/algorithm/extremes.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithm/extremes.rs b/src/algorithm/extremes.rs index d4f652283..501d1e093 100644 --- a/src/algorithm/extremes.rs +++ b/src/algorithm/extremes.rs @@ -58,7 +58,7 @@ impl Polygon } // positive implies a -> b -> c is counter-clockwise, negative implies clockwise -pub fn cross_prod(p_a: &Point, p_b: &Point, p_c: &Point) -> T +fn cross_prod(p_a: &Point, p_b: &Point, p_c: &Point) -> T where T: Float { (p_b.x() - p_a.x()) * (p_c.y() - p_a.y()) - (p_b.y() - p_a.y()) * (p_c.x() - p_a.x()) From a969dabc6ed4e20f7f2a90fa304aa48291dff060 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Mon, 17 Apr 2017 11:36:09 +0100 Subject: [PATCH 16/16] Make points raw --- src/algorithm/extremes.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithm/extremes.rs b/src/algorithm/extremes.rs index 501d1e093..3f81f0533 100644 --- a/src/algorithm/extremes.rs +++ b/src/algorithm/extremes.rs @@ -282,7 +282,7 @@ mod test { #[test] fn test_polygon_extreme_wrapper_convex() { // convex, with a bump on the top-right edge - let points_raw = + let mut points_raw = vec![(1.0, 0.0), (2.0, 1.0), (1.75, 1.75), (1.0, 2.0), (0.0, 1.0), (1.0, 0.0)]; let points = points_raw .iter()