diff --git a/geo/CHANGES.md b/geo/CHANGES.md index 00a9f1a41..b56b1662c 100644 --- a/geo/CHANGES.md +++ b/geo/CHANGES.md @@ -13,6 +13,9 @@ shouldn't break for any common numeric types, but if you are using something exotic you'll need to manually implement `GeoNum` for your numeric type. * +* POSSIBLY BREAKING: `minimum_rotated_rect` is about 33%, but might return + slightly different results. + * * POSSIBLY BREAKING: `SimplifyVwPreserve` trait implementation moved from `geo_types::CoordNum` to `geo::GeoNum` as a consequence of introducing the `GeoNum::total_cmp`. This shouldn't break anything for common numeric diff --git a/geo/Cargo.toml b/geo/Cargo.toml index f38fb5cbc..ec2e31711 100644 --- a/geo/Cargo.toml +++ b/geo/Cargo.toml @@ -114,3 +114,7 @@ harness = false [[bench]] name = "triangulate" harness = false + +[[bench]] +name = "minimum_rotated_rect" +harness = false diff --git a/geo/src/algorithm/minimum_rotated_rect.rs b/geo/src/algorithm/minimum_rotated_rect.rs index 73ecc3aba..213443c9c 100644 --- a/geo/src/algorithm/minimum_rotated_rect.rs +++ b/geo/src/algorithm/minimum_rotated_rect.rs @@ -1,9 +1,7 @@ -use num_traits::Float; +use geo_types::LineString; +use num_traits::Bounded; -use crate::{ - algorithm::{centroid::Centroid, rotate::Rotate, BoundingRect, CoordsIter}, - Area, ConvexHull, CoordFloat, GeoFloat, GeoNum, LinesIter, Polygon, -}; +use crate::{algorithm::CoordsIter, ConvexHull, CoordFloat, GeoFloat, GeoNum, Polygon}; /// Return the minimum bounding rectangle(MBR) of geometry /// reference: /// minimum rotated rect is the rectangle that can enclose all points given @@ -15,16 +13,16 @@ use crate::{ /// ``` /// use geo_types::{line_string, polygon, LineString, Polygon}; /// use geo::MinimumRotatedRect; -/// let poly: Polygon = polygon![(x: 3.3, y: 30.4), (x: 1.7, y: 24.6), (x: 13.4, y: 25.1), (x: 14.4, y: 31.0),(x:3.3,y:30.4)]; +/// let poly: Polygon = polygon![(x: 3.3, y: 30.4), (x: 1.7, y: 24.6), (x: 13.4, y: 25.1), (x: 14.4, y: 31.0), (x:3.3, y:30.4)]; /// let mbr = MinimumRotatedRect::minimum_rotated_rect(&poly).unwrap(); /// assert_eq!( /// mbr.exterior(), /// &LineString::from(vec![ -/// (1.7000000000000006, 24.6), -/// (1.4501458363715918, 30.446587428904767), -/// (14.4, 31.0), -/// (14.649854163628408, 25.153412571095235), -/// (1.7000000000000006, 24.6), +/// (1.7000000000000004, 24.600000000000005), +/// (14.64985416362841, 25.153412571095235), +/// (14.400000000000002, 31.000000000000004), +/// (1.4501458363715916, 30.446587428904774), +/// (1.7000000000000004, 24.600000000000005), /// ]) /// ); /// ``` @@ -41,24 +39,74 @@ where type Scalar = T; fn minimum_rotated_rect(&self) -> Option> { - let convex_poly = ConvexHull::convex_hull(self); - let mut min_area: T = Float::max_value(); - let mut min_angle: T = T::zero(); - let mut rect_poly: Option> = None; - let rotate_point = convex_poly.centroid(); - for line in convex_poly.exterior().lines_iter() { - let (ci, cii) = line.points(); - let angle = (cii.y() - ci.y()).atan2(cii.x() - ci.x()).to_degrees(); - let rotated_poly = Rotate::rotate_around_point(&convex_poly, -angle, rotate_point?); - let tmp_poly = rotated_poly.bounding_rect()?.to_polygon(); - let area = tmp_poly.unsigned_area(); + let hull = ConvexHull::convex_hull(self); + + // We take unit vectors along and perpendicular to each edge, then use + // them to build a rotation matrix and apply it to the convex hull, + // keeping track of the best AABB found. + // + // See also the discussion at + // https://gis.stackexchange.com/questions/22895/finding-minimum-area-rectangle-for-given-points/22934 + let mut min_area = ::max_value(); + let mut best_edge = None; + let (mut best_min_x, mut best_max_x) = (T::zero(), T::zero()); + let (mut best_min_y, mut best_max_y) = (T::zero(), T::zero()); + for edge in hull.exterior().lines() { + let (dx, dy) = edge.delta().x_y(); + let norm = dx.hypot(dy); + if norm.is_zero() { + continue; + } + let edge = (dx / norm, dy / norm); + + let (mut min_x, mut max_x) = (::max_value(), ::min_value()); + let (mut min_y, mut max_y) = (::max_value(), ::min_value()); + for point in hull.exterior().points() { + let x = point.x() * edge.0 + point.y() * edge.1; + let y = point.x() * -edge.1 + point.y() * edge.0; + + min_x = min_x.min(x); + max_x = max_x.max(x); + min_y = min_y.min(y); + max_y = max_y.max(y); + } + + let area = (max_y - min_y) * (max_x - min_x); if area < min_area { min_area = area; - min_angle = angle; - rect_poly = Some(tmp_poly); + best_edge = Some(edge); + best_min_x = min_x; + best_max_x = max_x; + best_min_y = min_y; + best_max_y = max_y; } } - Some(rect_poly?.rotate_around_point(min_angle, rotate_point?)) + + if let Some((dx, dy)) = best_edge { + let p1 = ( + best_min_x * dx + -best_min_y * dy, + best_min_x * dy + best_min_y * dx, + ); + let p2 = ( + best_max_x * dx + -best_min_y * dy, + best_max_x * dy + best_min_y * dx, + ); + let p3 = ( + best_max_x * dx + -best_max_y * dy, + best_max_x * dy + best_max_y * dx, + ); + let p4 = ( + best_min_x * dx + -best_max_y * dy, + best_min_x * dy + best_max_y * dx, + ); + let rectangle = Polygon::new( + LineString(vec![p1.into(), p2.into(), p3.into(), p4.into(), p1.into()]), + vec![], + ); + Some(rectangle) + } else { + None + } } } @@ -75,11 +123,11 @@ mod test { assert_eq!( mbr.exterior(), &LineString::from(vec![ - (1.7000000000000006, 24.6), - (1.4501458363715918, 30.446587428904767), - (14.4, 31.0), - (14.649854163628408, 25.153412571095235), - (1.7000000000000006, 24.6), + (1.7000000000000004, 24.600000000000005), + (14.64985416362841, 25.153412571095235), + (14.400000000000002, 31.000000000000004), + (1.4501458363715916, 30.446587428904774), + (1.7000000000000004, 24.600000000000005), ]) ); } @@ -90,11 +138,11 @@ mod test { assert_eq!( mbr.exterior(), &LineString::from(vec![ - (1.7000000000000006, 24.6), - (1.4501458363715918, 30.446587428904767), - (14.4, 31.0), - (14.649854163628408, 25.153412571095235), - (1.7000000000000006, 24.6), + (1.7000000000000004, 24.600000000000005), + (14.64985416362841, 25.153412571095235), + (14.400000000000002, 31.000000000000004), + (1.4501458363715916, 30.446587428904774), + (1.7000000000000004, 24.600000000000005), ]) ); }