From 5429bd7b61401c77e0163cf663f111b3fb6f3447 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lauren=C8=9Biu=20Nicola?= Date: Fri, 20 Oct 2023 11:47:36 +0300 Subject: [PATCH 1/3] Avoid atan2 in mrr --- geo/src/algorithm/minimum_rotated_rect.rs | 118 +++++++++++++++------- 1 file changed, 83 insertions(+), 35 deletions(-) diff --git a/geo/src/algorithm/minimum_rotated_rect.rs b/geo/src/algorithm/minimum_rotated_rect.rs index 73ecc3aba..7964561cd 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.649854163628412, 25.153412571095238), +/// (14.400000000000002, 31.000000000000007), +/// (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.y() * edge.1 + point.x() * edge.0; + let y = point.y() * edge.0 - point.x() * edge.1; + + 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(e) = best_edge { + let p1 = ( + best_min_x * e.0 - best_min_y * e.1, + best_min_x * e.1 + best_min_y * e.0, + ); + let p2 = ( + best_max_x * e.0 - best_min_y * e.1, + best_max_x * e.1 + best_min_y * e.0, + ); + let p3 = ( + best_max_x * e.0 - best_max_y * e.1, + best_max_x * e.1 + best_max_y * e.0, + ); + let p4 = ( + best_min_x * e.0 - best_max_y * e.1, + best_min_x * e.1 + best_max_y * e.0, + ); + 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.649854163628412, 25.153412571095238), + (14.400000000000002, 31.000000000000007), + (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.649854163628412, 25.153412571095238), + (14.400000000000002, 31.000000000000007), + (1.4501458363715916, 30.446587428904774), + (1.7000000000000004, 24.600000000000005), ]) ); } From 626dce664d147969ea94338bfb0ab22a915a1ef6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lauren=C8=9Biu=20Nicola?= Date: Fri, 20 Oct 2023 19:11:35 +0300 Subject: [PATCH 2/3] Use FMA --- geo/src/algorithm/minimum_rotated_rect.rs | 32 +++++++++++------------ 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/geo/src/algorithm/minimum_rotated_rect.rs b/geo/src/algorithm/minimum_rotated_rect.rs index 7964561cd..0098eaed5 100644 --- a/geo/src/algorithm/minimum_rotated_rect.rs +++ b/geo/src/algorithm/minimum_rotated_rect.rs @@ -19,8 +19,8 @@ use crate::{algorithm::CoordsIter, ConvexHull, CoordFloat, GeoFloat, GeoNum, Pol /// mbr.exterior(), /// &LineString::from(vec![ /// (1.7000000000000004, 24.600000000000005), -/// (14.649854163628412, 25.153412571095238), -/// (14.400000000000002, 31.000000000000007), +/// (14.64985416362841, 25.153412571095235), +/// (14.400000000000002, 31.000000000000004), /// (1.4501458363715916, 30.446587428904774), /// (1.7000000000000004, 24.600000000000005), /// ]) @@ -62,8 +62,8 @@ where 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.y() * edge.1 + point.x() * edge.0; - let y = point.y() * edge.0 - point.x() * edge.1; + let x = point.x().mul_add(edge.0, point.y() * edge.1); + let y = point.x().mul_add(-edge.1, point.y() * edge.0); min_x = min_x.min(x); max_x = max_x.max(x); @@ -84,20 +84,20 @@ where if let Some(e) = best_edge { let p1 = ( - best_min_x * e.0 - best_min_y * e.1, - best_min_x * e.1 + best_min_y * e.0, + best_min_x.mul_add(e.0, -best_min_y * e.1), + best_min_x.mul_add(e.1, best_min_y * e.0), ); let p2 = ( - best_max_x * e.0 - best_min_y * e.1, - best_max_x * e.1 + best_min_y * e.0, + best_max_x.mul_add(e.0, -best_min_y * e.1), + best_max_x.mul_add(e.1, best_min_y * e.0), ); let p3 = ( - best_max_x * e.0 - best_max_y * e.1, - best_max_x * e.1 + best_max_y * e.0, + best_max_x.mul_add(e.0, -best_max_y * e.1), + best_max_x.mul_add(e.1, best_max_y * e.0), ); let p4 = ( - best_min_x * e.0 - best_max_y * e.1, - best_min_x * e.1 + best_max_y * e.0, + best_min_x.mul_add(e.0, -best_max_y * e.1), + best_min_x.mul_add(e.1, best_max_y * e.0), ); let rectangle = Polygon::new( LineString(vec![p1.into(), p2.into(), p3.into(), p4.into(), p1.into()]), @@ -124,8 +124,8 @@ mod test { mbr.exterior(), &LineString::from(vec![ (1.7000000000000004, 24.600000000000005), - (14.649854163628412, 25.153412571095238), - (14.400000000000002, 31.000000000000007), + (14.64985416362841, 25.153412571095235), + (14.400000000000002, 31.000000000000004), (1.4501458363715916, 30.446587428904774), (1.7000000000000004, 24.600000000000005), ]) @@ -139,8 +139,8 @@ mod test { mbr.exterior(), &LineString::from(vec![ (1.7000000000000004, 24.600000000000005), - (14.649854163628412, 25.153412571095238), - (14.400000000000002, 31.000000000000007), + (14.64985416362841, 25.153412571095235), + (14.400000000000002, 31.000000000000004), (1.4501458363715916, 30.446587428904774), (1.7000000000000004, 24.600000000000005), ]) From 3ed0cfa221964e6204064c28037d7e40082ff1e7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lauren=C8=9Biu=20Nicola?= Date: Sat, 3 Feb 2024 18:45:34 +0200 Subject: [PATCH 3/3] Update changelog --- geo/CHANGES.md | 3 +++ geo/Cargo.toml | 4 ++++ geo/src/algorithm/minimum_rotated_rect.rs | 22 +++++++++++----------- 3 files changed, 18 insertions(+), 11 deletions(-) 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 0098eaed5..213443c9c 100644 --- a/geo/src/algorithm/minimum_rotated_rect.rs +++ b/geo/src/algorithm/minimum_rotated_rect.rs @@ -62,8 +62,8 @@ where 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().mul_add(edge.0, point.y() * edge.1); - let y = point.x().mul_add(-edge.1, point.y() * edge.0); + 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); @@ -82,22 +82,22 @@ where } } - if let Some(e) = best_edge { + if let Some((dx, dy)) = best_edge { let p1 = ( - best_min_x.mul_add(e.0, -best_min_y * e.1), - best_min_x.mul_add(e.1, best_min_y * e.0), + best_min_x * dx + -best_min_y * dy, + best_min_x * dy + best_min_y * dx, ); let p2 = ( - best_max_x.mul_add(e.0, -best_min_y * e.1), - best_max_x.mul_add(e.1, best_min_y * e.0), + best_max_x * dx + -best_min_y * dy, + best_max_x * dy + best_min_y * dx, ); let p3 = ( - best_max_x.mul_add(e.0, -best_max_y * e.1), - best_max_x.mul_add(e.1, best_max_y * e.0), + best_max_x * dx + -best_max_y * dy, + best_max_x * dy + best_max_y * dx, ); let p4 = ( - best_min_x.mul_add(e.0, -best_max_y * e.1), - best_min_x.mul_add(e.1, best_max_y * e.0), + 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()]),