From 5d0e037b494683bdd628c284250cbe160387ebc6 Mon Sep 17 00:00:00 2001 From: thehappycheese Date: Wed, 28 Jun 2023 00:37:07 +0800 Subject: [PATCH 01/11] add trait VectorOps & impl for Coord --- geo/src/algorithm/mod.rs | 5 + geo/src/algorithm/vector_ops.rs | 343 ++++++++++++++++++++++++++++++++ 2 files changed, 348 insertions(+) create mode 100644 geo/src/algorithm/vector_ops.rs diff --git a/geo/src/algorithm/mod.rs b/geo/src/algorithm/mod.rs index bba8a95d2c..98f68566a9 100644 --- a/geo/src/algorithm/mod.rs +++ b/geo/src/algorithm/mod.rs @@ -232,6 +232,11 @@ pub mod triangulate_earcut; #[cfg(feature = "earcutr")] pub use triangulate_earcut::TriangulateEarcut; +/// # Vector Operations for 2D coordinates +/// +mod vector_ops; +pub(crate) use vector_ops::Vector2DOps; + /// Calculate the Vincenty distance between two `Point`s. pub mod vincenty_distance; pub use vincenty_distance::VincentyDistance; diff --git a/geo/src/algorithm/vector_ops.rs b/geo/src/algorithm/vector_ops.rs new file mode 100644 index 0000000000..678e040378 --- /dev/null +++ b/geo/src/algorithm/vector_ops.rs @@ -0,0 +1,343 @@ +//! This module defines the [Vector2DOps] trait and implements it for the +//! [Coord] struct. +//! +//! In the future [Vector2DOps] might be implemented on other types: +//! +//! - Based on community discussions it seems like the existing struct +//! [crate::Coord] is just one of many future data structures which may hold +//! coordinate information. For example future data structures may also hold Z +//! and/or M ordinates, or other arbitrary data. +//! - [geotraits::CoordTrait] is a future trait defining accessor methods +//! `.x()` and `.y()` which will facilitate these generic data structures. +//! +//! > Note: [Vector2DOps] is not implemented for [crate::Point] because users +//! > probably expect to see more familiar geospatial functions like +//! > `a.euclidean_distance(b)` at that level and generally not linear algebra +//! > like `a.dot(b)`. +//! +//! For now, it is assumed that users of [Vector2DOps] are satisfied to upcast +//! everything to [CoordFloat]. In future it might make sense to split this +//! trait into 3 flavors supporting progressively more functions: +//! +//! - `trait Vector2DOpsCoordNum` for [CoordNum] +//! - Supports magnitude_squared, dot +//! - `trait Vector2DOpsCoordNumSigned:Vector2DOpsCoordNum` for [CoordNum] + [Signed] +//! - Supports left, right, wedge_product +//! - `trait Vector2DOpsCoordFloat:Vector2DOpsCoordNumSigned` for [CoordFloat] +//! - Supports magnitude, normalize, etc +//! +//! Maybe if these traits were restricted to the future [geotraits::CoordTrait] +//! then they could provide default implementations using the accessors `.x()` +//! and `.y()`? +//! + +use crate::{Coord, CoordFloat, CoordNum}; + +/// Defines vector operations for 2D coordinate types which implement CoordFloat +/// +/// This trait is intended for internal use within the geo crate as a way to +/// bring together the various hand-crafted linear algebra operations used +/// throughout other algorithms and attached to various structs. +/// +/// + +pub trait Vector2DOps +where + Self: Sized, +{ + type NumericType: CoordNum + Send + Sync; + + /// The euclidean distance between this coordinate and the origin + /// + /// `sqrt(x² + y²)` + /// + fn magnitude(self) -> Self::NumericType; + + /// The squared distance between this coordinate and the origin. + /// (Avoids the square root calculation when it is not needed) + /// + /// `x² + y²` + /// + fn magnitude_squared(self) -> Self::NumericType; + + /// Rotate this coordinate around the origin by 90 degrees clockwise. + /// + /// `a.left() => (-a.y, a.x)` + /// + /// Assumes a coordinate system where positive `y` is up and positive `x` is + /// to the right. The described rotation direction is consistent with the + /// documentation for [crate::algorithm::rotate::Rotate]. + fn left(self) -> Self; + + /// Rotate this coordinate around the origin by 90 degrees anti-clockwise. + /// + /// `a.right() => (a.y, -a.x)` + /// + /// Assumes a coordinate system where positive `y` is up and positive `x` is + /// to the right. The described rotation direction is consistent with the + /// documentation for [crate::algorithm::rotate::Rotate]. + fn right(self) -> Self; + + /// The inner product of the coordinate components + /// + /// `a · b = a.x * b.x + a.y * b.y` + /// + fn dot_product(self, other: Rhs) -> Self::NumericType; + + /// The calculates the `wedge product` between two vectors. + /// + /// `a ∧ b = a.x * b.y - a.y * b.x` + /// + /// Also known as: + /// + /// - `exterior product` + /// - because the wedge product comes from 'Exterior Algebra' + /// - `perpendicular product` + /// - because it is equivalent to `a.dot(b.right())` + /// - `2D cross product` + /// - because it is equivalent to the signed magnitude of the + /// conventional 3D cross product assuming `z` ordinates are zero + /// - `determinant` + /// - because it is equivalent to the `determinant` of the 2x2 matrix + /// formed by the column-vector inputs. + /// + /// ## Examples + /// + /// The following list highlights some examples in geo which might be + /// brought together to use this function: + /// + /// 1. [geo_types::Point::cross_prod()] is already defined on + /// [geo_types::Point]... but that it seems to be some other + /// operation on 3 points?? + /// 2. [geo_types::Line] struct also has a [geo_types::Line::determinant()] + /// function which is the same as `line.start.wedge_product(line.end)` + /// 3. The [crate::algorithm::Kernel::orient2d()] trait default + /// implementation uses cross product to compute orientation. It returns + /// an enum, not the numeric value which is needed for line segment + /// intersection. + /// + /// ## Properties + /// + /// - The absolute value of the cross product is the area of the + /// parallelogram formed by the operands + /// - Anti-commutative: The sign of the output is reversed if the operands + /// are reversed + /// - If the operands are colinear with the origin, the value is zero + /// - The sign can be used to check if the operands are clockwise with + /// respect to the origin, or phrased differently: + /// "is a to the left of the line between the origin and b"? + /// - If this is what you are using it for, then please use + /// [crate::algorithm::Kernel::orient2d()] instead as this is more + /// explicit and has a `RobustKernel` option for extra precision. + fn wedge_product(self, other: Rhs) -> Self::NumericType; + + /// Try to find a vector of unit length in the same direction as this + /// vector. + /// + /// Returns `None` if the magnitude of this vector is less than + /// `minimum_magnitude` or the magnitude is not finite. + /// - For f32 the minimum_magnitude can be set to about 1e-30f32. + /// - For F64 the minimum_magnitude can be set to about 2e-301f64. + /// + /// These values should avoid overflowing to Infinity for coordinate values + /// in the range typically relevant for spatial data (+-40e6) which is the + /// approximate length of the earth's great circle in metres. + /// + /// > Note to Reviewer: it may be annoying to have to provide a value for + /// > minimum_magnitude, but that seems to be what `nalgebra` does + /// > (See https://docs.rs/nalgebra/latest/src/nalgebra/base/norm.rs.html#301-307). + /// > Some other parts of the api do not require the user to specify a + /// > value, but I haven't yet figured out how those work because it is + /// > wrapped up in the simba SIMD crate in complicated macros. + /// > + /// > Open to suggestions about how this can be better handled, or the + /// > try_normalize function can just be removed for now. + fn try_normalize(self, minimum_magnitude:Self::NumericType) -> Option; +} + +impl Vector2DOps for Coord +where + T: CoordFloat + Send + Sync, +{ + type NumericType = T; + + fn wedge_product(self, right: Coord) -> Self::NumericType { + self.x * right.y - self.y * right.x + } + + fn dot_product(self, other: Self) -> Self::NumericType { + self.x * other.x + self.y * other.y + } + + fn magnitude(self) -> Self::NumericType { + (self.x * self.x + self.y * self.y).sqrt() + } + + fn magnitude_squared(self) -> Self::NumericType { + self.x * self.x + self.y * self.y + } + + fn left(self) -> Self { + Self { + x: -self.y, + y: self.x, + } + } + + fn right(self) -> Self { + Self { + x: self.y, + y: -self.x, + } + } + + fn try_normalize(self, minimum_magnitude: Self::NumericType) -> Option { + let magnitude = self.magnitude(); + if magnitude.is_finite() && magnitude.abs() > minimum_magnitude { + Some(self / magnitude) + }else{ + None + } + } +} + +#[cfg(test)] +mod test { + // crate dependencies + use crate::coord; + + // private imports + use super::Vector2DOps; + + #[test] + fn test_cross_product() { + // perpendicular unit length + let a = coord! { x: 1f64, y: 0f64 }; + let b = coord! { x: 0f64, y: 1f64 }; + + // expect the area of parallelogram + assert_eq!(a.wedge_product(b), 1f64); + // expect swapping will result in negative + assert_eq!(b.wedge_product(a), -1f64); + + // Add skew; expect results should be the same + let a = coord! { x: 1f64, y: 0f64 }; + let b = coord! { x: 1f64, y: 1f64 }; + + // expect the area of parallelogram + assert_eq!(a.wedge_product(b), 1f64); + // expect swapping will result in negative + assert_eq!(b.wedge_product(a), -1f64); + + // Make Colinear; expect zero + let a = coord! { x: 2f64, y: 2f64 }; + let b = coord! { x: 1f64, y: 1f64 }; + assert_eq!(a.wedge_product(b), 0f64); + } + + #[test] + fn test_dot_product() { + // perpendicular unit length + let a = coord! { x: 1f64, y: 0f64 }; + let b = coord! { x: 0f64, y: 1f64 }; + // expect zero for perpendicular + assert_eq!(a.dot_product(b), 0f64); + + // Parallel, same direction + let a = coord! { x: 1f64, y: 0f64 }; + let b = coord! { x: 2f64, y: 0f64 }; + // expect +ive product of magnitudes + assert_eq!(a.dot_product(b), 2f64); + // expect swapping will have same result + assert_eq!(b.dot_product(a), 2f64); + + // Parallel, opposite direction + let a = coord! { x: 3f64, y: 4f64 }; + let b = coord! { x: -3f64, y: -4f64 }; + // expect -ive product of magnitudes + assert_eq!(a.dot_product(b), -25f64); + // expect swapping will have same result + assert_eq!(b.dot_product(a), -25f64); + } + + #[test] + fn test_magnitude() { + let a = coord! { x: 1f64, y: 0f64 }; + assert_eq!(a.magnitude(), 1f64); + + let a = coord! { x: 0f64, y: 0f64 }; + assert_eq!(a.magnitude(), 0f64); + + let a = coord! { x: -3f64, y: 4f64 }; + assert_eq!(a.magnitude(), 5f64); + } + + #[test] + fn test_magnitude_squared() { + let a = coord! { x: 1f64, y: 0f64 }; + assert_eq!(a.magnitude_squared(), 1f64); + + let a = coord! { x: 0f64, y: 0f64 }; + assert_eq!(a.magnitude_squared(), 0f64); + + let a = coord! { x: -3f64, y: 4f64 }; + assert_eq!(a.magnitude_squared(), 25f64); + } + + #[test] + fn test_left_right() { + let a = coord! { x: 1f64, y: 0f64 }; + let a_left = coord! { x: 0f64, y: 1f64 }; + let a_right = coord! { x: 0f64, y: -1f64 }; + + assert_eq!(a.left(), a_left); + assert_eq!(a.right(), a_right); + assert_eq!(a.left(), -a.right()); + } + + #[test] + fn test_left_right_match_rotate() { + use crate::algorithm::rotate::Rotate; + use crate::Point; + // The aim of this test is to confirm that wording in documentation is + // consistent. + + // when the user is in a coordinate system where the y axis is flipped + // (eg screen coordinates in a HTML canvas), then rotation directions + // will be different to those described in the documentation. + + // The documentation for the Rotate trait says: 'Positive angles are + // counter-clockwise, and negative angles are clockwise rotations' + + let counter_clockwise_rotation_degrees = 90.0; + let clockwise_rotation_degrees = -counter_clockwise_rotation_degrees; + + let a: Point = coord! { x: 1.0, y: 0.0 }.into(); + let origin: Point = coord! { x: 0.0, y: 0.0 }.into(); + + // left is anti-clockwise + assert_relative_eq!( + Point::from(a.0.left()), + a.rotate_around_point(counter_clockwise_rotation_degrees, origin), + ); + // right is clockwise + assert_relative_eq!( + Point::from(a.0.right()), + a.rotate_around_point(clockwise_rotation_degrees, origin), + ); + } + + #[test] + fn test_normalize() { + let a = coord! { x: 1f64, y: 0f64 }; + assert_eq!(a.try_normalize(2e-301f64), Some(a)); + + let a = coord! { x: 1.0 / f64::sqrt(2.0), y: -1.0 / f64::sqrt(2.0) }; + assert_eq!(a.try_normalize(2e-301f64), Some(a)); + + let a = coord! { x: -10.0, y: 8.0 }; + assert_eq!(a.try_normalize(2e-301f64), Some(a/f64::sqrt(10.0*10.0+8.0*8.0))); + } + + +} From 7a9a6d97f5636380a926b5c8c58117478891f434 Mon Sep 17 00:00:00 2001 From: thehappycheese Date: Wed, 28 Jun 2023 00:43:00 +0800 Subject: [PATCH 02/11] fix test for try_normalize --- geo/src/algorithm/vector_ops.rs | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/geo/src/algorithm/vector_ops.rs b/geo/src/algorithm/vector_ops.rs index 678e040378..509022c185 100644 --- a/geo/src/algorithm/vector_ops.rs +++ b/geo/src/algorithm/vector_ops.rs @@ -329,14 +329,17 @@ mod test { #[test] fn test_normalize() { - let a = coord! { x: 1f64, y: 0f64 }; - assert_eq!(a.try_normalize(2e-301f64), Some(a)); + let a = coord! { x: 1.0, y: 0.0 }; + assert_relative_eq!(a.try_normalize(2e-301f64).unwrap(), a); let a = coord! { x: 1.0 / f64::sqrt(2.0), y: -1.0 / f64::sqrt(2.0) }; - assert_eq!(a.try_normalize(2e-301f64), Some(a)); + assert_relative_eq!(a.try_normalize(2e-301f64).unwrap(), a); let a = coord! { x: -10.0, y: 8.0 }; - assert_eq!(a.try_normalize(2e-301f64), Some(a/f64::sqrt(10.0*10.0+8.0*8.0))); + assert_relative_eq!(a.try_normalize(2e-301f64).unwrap(), a/f64::sqrt(10.0*10.0+8.0*8.0)); + + let a = coord! { x: 0.0, y: 1e-301f64 }; + assert_eq!(a.try_normalize(2e-301f64), None); } From b6a2227de2e32f05337b813f8754bf6c1ae2dcc2 Mon Sep 17 00:00:00 2001 From: thehappycheese Date: Wed, 28 Jun 2023 00:46:24 +0800 Subject: [PATCH 03/11] fix format --- geo/src/algorithm/vector_ops.rs | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/geo/src/algorithm/vector_ops.rs b/geo/src/algorithm/vector_ops.rs index 509022c185..8eb9057cde 100644 --- a/geo/src/algorithm/vector_ops.rs +++ b/geo/src/algorithm/vector_ops.rs @@ -138,21 +138,21 @@ where /// `minimum_magnitude` or the magnitude is not finite. /// - For f32 the minimum_magnitude can be set to about 1e-30f32. /// - For F64 the minimum_magnitude can be set to about 2e-301f64. - /// + /// /// These values should avoid overflowing to Infinity for coordinate values /// in the range typically relevant for spatial data (+-40e6) which is the /// approximate length of the earth's great circle in metres. /// - /// > Note to Reviewer: it may be annoying to have to provide a value for + /// > Note to Reviewer: it may be annoying to have to provide a value for /// > minimum_magnitude, but that seems to be what `nalgebra` does /// > (See https://docs.rs/nalgebra/latest/src/nalgebra/base/norm.rs.html#301-307). - /// > Some other parts of the api do not require the user to specify a + /// > Some other parts of the api do not require the user to specify a /// > value, but I haven't yet figured out how those work because it is /// > wrapped up in the simba SIMD crate in complicated macros. - /// > + /// > /// > Open to suggestions about how this can be better handled, or the /// > try_normalize function can just be removed for now. - fn try_normalize(self, minimum_magnitude:Self::NumericType) -> Option; + fn try_normalize(self, minimum_magnitude: Self::NumericType) -> Option; } impl Vector2DOps for Coord @@ -195,7 +195,7 @@ where let magnitude = self.magnitude(); if magnitude.is_finite() && magnitude.abs() > minimum_magnitude { Some(self / magnitude) - }else{ + } else { None } } @@ -336,11 +336,12 @@ mod test { assert_relative_eq!(a.try_normalize(2e-301f64).unwrap(), a); let a = coord! { x: -10.0, y: 8.0 }; - assert_relative_eq!(a.try_normalize(2e-301f64).unwrap(), a/f64::sqrt(10.0*10.0+8.0*8.0)); + assert_relative_eq!( + a.try_normalize(2e-301f64).unwrap(), + a / f64::sqrt(10.0 * 10.0 + 8.0 * 8.0) + ); let a = coord! { x: 0.0, y: 1e-301f64 }; assert_eq!(a.try_normalize(2e-301f64), None); } - - } From 923b8f6a2372ce493893319a4ba750410d083002 Mon Sep 17 00:00:00 2001 From: thehappycheese Date: Thu, 13 Jul 2023 20:17:32 +0800 Subject: [PATCH 04/11] Rename associated type NumericType to Scalar --- geo/src/algorithm/vector_ops.rs | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/geo/src/algorithm/vector_ops.rs b/geo/src/algorithm/vector_ops.rs index 8eb9057cde..6933aeea08 100644 --- a/geo/src/algorithm/vector_ops.rs +++ b/geo/src/algorithm/vector_ops.rs @@ -45,20 +45,20 @@ pub trait Vector2DOps where Self: Sized, { - type NumericType: CoordNum + Send + Sync; + type Scalar: CoordNum + Send + Sync; /// The euclidean distance between this coordinate and the origin /// /// `sqrt(x² + y²)` /// - fn magnitude(self) -> Self::NumericType; + fn magnitude(self) -> Self::Scalar; /// The squared distance between this coordinate and the origin. /// (Avoids the square root calculation when it is not needed) /// /// `x² + y²` /// - fn magnitude_squared(self) -> Self::NumericType; + fn magnitude_squared(self) -> Self::Scalar; /// Rotate this coordinate around the origin by 90 degrees clockwise. /// @@ -82,7 +82,7 @@ where /// /// `a · b = a.x * b.x + a.y * b.y` /// - fn dot_product(self, other: Rhs) -> Self::NumericType; + fn dot_product(self, other: Rhs) -> Self::Scalar; /// The calculates the `wedge product` between two vectors. /// @@ -129,7 +129,7 @@ where /// - If this is what you are using it for, then please use /// [crate::algorithm::Kernel::orient2d()] instead as this is more /// explicit and has a `RobustKernel` option for extra precision. - fn wedge_product(self, other: Rhs) -> Self::NumericType; + fn wedge_product(self, other: Rhs) -> Self::Scalar; /// Try to find a vector of unit length in the same direction as this /// vector. @@ -152,28 +152,28 @@ where /// > /// > Open to suggestions about how this can be better handled, or the /// > try_normalize function can just be removed for now. - fn try_normalize(self, minimum_magnitude: Self::NumericType) -> Option; + fn try_normalize(self, minimum_magnitude: Self::Scalar) -> Option; } impl Vector2DOps for Coord where T: CoordFloat + Send + Sync, { - type NumericType = T; + type Scalar = T; - fn wedge_product(self, right: Coord) -> Self::NumericType { + fn wedge_product(self, right: Coord) -> Self::Scalar { self.x * right.y - self.y * right.x } - fn dot_product(self, other: Self) -> Self::NumericType { + fn dot_product(self, other: Self) -> Self::Scalar { self.x * other.x + self.y * other.y } - fn magnitude(self) -> Self::NumericType { + fn magnitude(self) -> Self::Scalar { (self.x * self.x + self.y * self.y).sqrt() } - fn magnitude_squared(self) -> Self::NumericType { + fn magnitude_squared(self) -> Self::Scalar { self.x * self.x + self.y * self.y } @@ -191,7 +191,7 @@ where } } - fn try_normalize(self, minimum_magnitude: Self::NumericType) -> Option { + fn try_normalize(self, minimum_magnitude: Self::Scalar) -> Option { let magnitude = self.magnitude(); if magnitude.is_finite() && magnitude.abs() > minimum_magnitude { Some(self / magnitude) From 0ba00f0098d665843257039a2b166be2408c4760 Mon Sep 17 00:00:00 2001 From: thehappycheese Date: Fri, 14 Jul 2023 00:17:59 +0800 Subject: [PATCH 05/11] make Vector2DOps public --- geo/src/algorithm/mod.rs | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/geo/src/algorithm/mod.rs b/geo/src/algorithm/mod.rs index 98f68566a9..5aeab5979b 100644 --- a/geo/src/algorithm/mod.rs +++ b/geo/src/algorithm/mod.rs @@ -232,10 +232,9 @@ pub mod triangulate_earcut; #[cfg(feature = "earcutr")] pub use triangulate_earcut::TriangulateEarcut; -/// # Vector Operations for 2D coordinates -/// +/// Vector Operations for 2D coordinates mod vector_ops; -pub(crate) use vector_ops::Vector2DOps; +pub use vector_ops::Vector2DOps; /// Calculate the Vincenty distance between two `Point`s. pub mod vincenty_distance; From fe9df2aafca42d9d293c2825e31b46ba0479b060 Mon Sep 17 00:00:00 2001 From: thehappycheese Date: Fri, 14 Jul 2023 00:37:11 +0800 Subject: [PATCH 06/11] remove commentary intended for reviewers --- geo/src/algorithm/vector_ops.rs | 40 --------------------------------- 1 file changed, 40 deletions(-) diff --git a/geo/src/algorithm/vector_ops.rs b/geo/src/algorithm/vector_ops.rs index 6933aeea08..8189f63012 100644 --- a/geo/src/algorithm/vector_ops.rs +++ b/geo/src/algorithm/vector_ops.rs @@ -1,35 +1,5 @@ //! This module defines the [Vector2DOps] trait and implements it for the //! [Coord] struct. -//! -//! In the future [Vector2DOps] might be implemented on other types: -//! -//! - Based on community discussions it seems like the existing struct -//! [crate::Coord] is just one of many future data structures which may hold -//! coordinate information. For example future data structures may also hold Z -//! and/or M ordinates, or other arbitrary data. -//! - [geotraits::CoordTrait] is a future trait defining accessor methods -//! `.x()` and `.y()` which will facilitate these generic data structures. -//! -//! > Note: [Vector2DOps] is not implemented for [crate::Point] because users -//! > probably expect to see more familiar geospatial functions like -//! > `a.euclidean_distance(b)` at that level and generally not linear algebra -//! > like `a.dot(b)`. -//! -//! For now, it is assumed that users of [Vector2DOps] are satisfied to upcast -//! everything to [CoordFloat]. In future it might make sense to split this -//! trait into 3 flavors supporting progressively more functions: -//! -//! - `trait Vector2DOpsCoordNum` for [CoordNum] -//! - Supports magnitude_squared, dot -//! - `trait Vector2DOpsCoordNumSigned:Vector2DOpsCoordNum` for [CoordNum] + [Signed] -//! - Supports left, right, wedge_product -//! - `trait Vector2DOpsCoordFloat:Vector2DOpsCoordNumSigned` for [CoordFloat] -//! - Supports magnitude, normalize, etc -//! -//! Maybe if these traits were restricted to the future [geotraits::CoordTrait] -//! then they could provide default implementations using the accessors `.x()` -//! and `.y()`? -//! use crate::{Coord, CoordFloat, CoordNum}; @@ -142,16 +112,6 @@ where /// These values should avoid overflowing to Infinity for coordinate values /// in the range typically relevant for spatial data (+-40e6) which is the /// approximate length of the earth's great circle in metres. - /// - /// > Note to Reviewer: it may be annoying to have to provide a value for - /// > minimum_magnitude, but that seems to be what `nalgebra` does - /// > (See https://docs.rs/nalgebra/latest/src/nalgebra/base/norm.rs.html#301-307). - /// > Some other parts of the api do not require the user to specify a - /// > value, but I haven't yet figured out how those work because it is - /// > wrapped up in the simba SIMD crate in complicated macros. - /// > - /// > Open to suggestions about how this can be better handled, or the - /// > try_normalize function can just be removed for now. fn try_normalize(self, minimum_magnitude: Self::Scalar) -> Option; } From 8763ef1ada3c7ce4ded8fc2e024840d14c7ccc27 Mon Sep 17 00:00:00 2001 From: thehappycheese Date: Fri, 14 Jul 2023 00:38:24 +0800 Subject: [PATCH 07/11] improve tests for try_normalize --- geo/src/algorithm/vector_ops.rs | 60 +++++++++++++++++++++++++++------ 1 file changed, 50 insertions(+), 10 deletions(-) diff --git a/geo/src/algorithm/vector_ops.rs b/geo/src/algorithm/vector_ops.rs index 8189f63012..ce52749a39 100644 --- a/geo/src/algorithm/vector_ops.rs +++ b/geo/src/algorithm/vector_ops.rs @@ -289,19 +289,59 @@ mod test { #[test] fn test_normalize() { - let a = coord! { x: 1.0, y: 0.0 }; - assert_relative_eq!(a.try_normalize(2e-301f64).unwrap(), a); - - let a = coord! { x: 1.0 / f64::sqrt(2.0), y: -1.0 / f64::sqrt(2.0) }; - assert_relative_eq!(a.try_normalize(2e-301f64).unwrap(), a); - + let f64_minimum_threshold = 2e-301f64; + // Already Normalized + let a = coord! { + x: 1.0, + y: 0.0 + }; + assert_relative_eq!(a.try_normalize(f64_minimum_threshold).unwrap(), a); + + // Already Normalized + let a = coord! { + x: 1.0 / f64::sqrt(2.0), + y: -1.0 / f64::sqrt(2.0) + }; + assert_relative_eq!(a.try_normalize(f64_minimum_threshold).unwrap(), a); + + // Non trivial example let a = coord! { x: -10.0, y: 8.0 }; assert_relative_eq!( - a.try_normalize(2e-301f64).unwrap(), - a / f64::sqrt(10.0 * 10.0 + 8.0 * 8.0) + a.try_normalize(f64_minimum_threshold).unwrap(), + coord! { x: -10.0, y: 8.0 } / f64::sqrt(10.0 * 10.0 + 8.0 * 8.0) ); - let a = coord! { x: 0.0, y: 1e-301f64 }; - assert_eq!(a.try_normalize(2e-301f64), None); + // Very Small Input - Fails because below threshold + let a = coord! { + x: 0.0, + y: f64_minimum_threshold / 2.0 + }; + assert_eq!(a.try_normalize(f64_minimum_threshold), None); + + // Zero vector - Fails because below threshold + let a = coord! { x: 0.0, y: 0.0 }; + assert_eq!(a.try_normalize(f64_minimum_threshold), None); + + // Large vector Pass; + // Note: this fails because the magnitude calculation overflows to infinity + // this could be avoided my modifying the implementation to use scaling be avoided using scaling + let a = coord! { + x: f64::sqrt(f64::MAX/2.0), + y: f64::sqrt(f64::MAX/2.0) + }; + assert!(a.try_normalize(f64_minimum_threshold).is_some()); + + // Large Vector Fail; + // Note: This test uses the unstable float_next_up_down feature + // let a = coord! { x: f64::sqrt(f64::MAX/2.0), y: f64::next_up(f64::sqrt(f64::MAX/2.0)) }; + // assert!(a.try_normalize(1e-10f64).is_none()); + + // NaN + let a = coord! { x: f64::NAN, y: 0.0 }; + assert_eq!(a.try_normalize(f64_minimum_threshold), None); + + // Infinite + let a = coord! { x: f64::INFINITY, y: 0.0 }; + assert_eq!(a.try_normalize(f64_minimum_threshold), None); } } From 1c132ae7277e579cbc4d46dc694523e600a3e886 Mon Sep 17 00:00:00 2001 From: thehappycheese Date: Fri, 14 Jul 2023 00:38:44 +0800 Subject: [PATCH 08/11] tidy doc comment for try_normalize --- geo/src/algorithm/vector_ops.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/geo/src/algorithm/vector_ops.rs b/geo/src/algorithm/vector_ops.rs index ce52749a39..1d4a3bb811 100644 --- a/geo/src/algorithm/vector_ops.rs +++ b/geo/src/algorithm/vector_ops.rs @@ -106,8 +106,8 @@ where /// /// Returns `None` if the magnitude of this vector is less than /// `minimum_magnitude` or the magnitude is not finite. - /// - For f32 the minimum_magnitude can be set to about 1e-30f32. - /// - For F64 the minimum_magnitude can be set to about 2e-301f64. + /// - For f32 the minimum_magnitude can be set to about `1e-30_f32` + /// - For F64 the minimum_magnitude can be set to about `2e-301_f64` /// /// These values should avoid overflowing to Infinity for coordinate values /// in the range typically relevant for spatial data (+-40e6) which is the From 5706059a977431bcad79fd0c13bcd71f95eec462 Mon Sep 17 00:00:00 2001 From: thehappycheese Date: Fri, 14 Jul 2023 00:42:14 +0800 Subject: [PATCH 09/11] add Vector2DOps to CHANGES.md --- geo/CHANGES.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/geo/CHANGES.md b/geo/CHANGES.md index fa38c4a575..8de905bc70 100644 --- a/geo/CHANGES.md +++ b/geo/CHANGES.md @@ -4,6 +4,8 @@ - Add `TriangulateEarcut` algorithm trait to triangulate polygons with the earcut algorithm. - +- Add `Vector2DOps` trait to algorithims module and implemented it for `Coord` + - ## 0.25.0 From 43d2fefe759398b45e8d3009e6087a2815d9700b Mon Sep 17 00:00:00 2001 From: thehappycheese Date: Wed, 26 Jul 2023 23:40:52 +0800 Subject: [PATCH 10/11] removed threshold from Vector2DOps::try_normalize --- geo/src/algorithm/vector_ops.rs | 100 +++++++++++++++++++------------- 1 file changed, 60 insertions(+), 40 deletions(-) diff --git a/geo/src/algorithm/vector_ops.rs b/geo/src/algorithm/vector_ops.rs index 1d4a3bb811..72feb9ad21 100644 --- a/geo/src/algorithm/vector_ops.rs +++ b/geo/src/algorithm/vector_ops.rs @@ -104,15 +104,17 @@ where /// Try to find a vector of unit length in the same direction as this /// vector. /// - /// Returns `None` if the magnitude of this vector is less than - /// `minimum_magnitude` or the magnitude is not finite. - /// - For f32 the minimum_magnitude can be set to about `1e-30_f32` - /// - For F64 the minimum_magnitude can be set to about `2e-301_f64` - /// - /// These values should avoid overflowing to Infinity for coordinate values - /// in the range typically relevant for spatial data (+-40e6) which is the - /// approximate length of the earth's great circle in metres. - fn try_normalize(self, minimum_magnitude: Self::Scalar) -> Option; + /// Returns `None` if the result is not finite. This can happen when + /// + /// - the vector is really small (or zero length) and the `.magnitude()` + /// calculation has rounded-down to `0.0` + /// - the vector is really large and the `.magnitude()` has rounded-up + /// or 'overflowed' to `f64::INFINITY` + /// - Either x or y are `f64::NAN` or `f64::INFINITY` + fn try_normalize(self) -> Option; + + /// Returns true if both the x and y components are finite + fn is_finite(self) -> bool; } impl Vector2DOps for Coord @@ -151,23 +153,30 @@ where } } - fn try_normalize(self, minimum_magnitude: Self::Scalar) -> Option { + fn try_normalize(self) -> Option { let magnitude = self.magnitude(); - if magnitude.is_finite() && magnitude.abs() > minimum_magnitude { - Some(self / magnitude) + let result = self / magnitude; + // Both the result AND the magnitude must be finite they are finite + // Otherwise very large vectors overflow magnitude to Infinity, + // and the after the division the result would be coord!{x:0.0,y:0.0} + // Note we don't need to check if magnitude is zero, because after the division + // that would have made result non-finite or NaN anyway. + if result.is_finite() && magnitude.is_finite() { + Some(result) } else { None } } + + fn is_finite(self) -> bool { + self.x.is_finite() && self.y.is_finite() + } } #[cfg(test)] mod test { - // crate dependencies - use crate::coord; - - // private imports use super::Vector2DOps; + use crate::coord; #[test] fn test_cross_product() { @@ -288,60 +297,71 @@ mod test { } #[test] - fn test_normalize() { - let f64_minimum_threshold = 2e-301f64; + fn test_try_normalize() { + // Already Normalized let a = coord! { x: 1.0, y: 0.0 }; - assert_relative_eq!(a.try_normalize(f64_minimum_threshold).unwrap(), a); + assert_relative_eq!(a.try_normalize().unwrap(), a); // Already Normalized let a = coord! { x: 1.0 / f64::sqrt(2.0), y: -1.0 / f64::sqrt(2.0) }; - assert_relative_eq!(a.try_normalize(f64_minimum_threshold).unwrap(), a); + assert_relative_eq!(a.try_normalize().unwrap(), a); // Non trivial example let a = coord! { x: -10.0, y: 8.0 }; assert_relative_eq!( - a.try_normalize(f64_minimum_threshold).unwrap(), + a.try_normalize().unwrap(), coord! { x: -10.0, y: 8.0 } / f64::sqrt(10.0 * 10.0 + 8.0 * 8.0) ); + } + + #[test] + fn test_try_normalize_edge_cases() { + use float_next_after::NextAfter; + + // The following tests demonstrate some of the floating point + // edge cases that can cause try_normalize to return None. + + // Zero vector - Normalize returns None + let a = coord! { x: 0.0, y: 0.0 }; + assert_eq!(a.try_normalize(), None); - // Very Small Input - Fails because below threshold + // Very Small Input - Normalize returns None because of + // rounding-down to zero in the .magnitude() calculation let a = coord! { x: 0.0, - y: f64_minimum_threshold / 2.0 + y: 1e-301_f64 }; - assert_eq!(a.try_normalize(f64_minimum_threshold), None); + assert_eq!(a.try_normalize(), None); - // Zero vector - Fails because below threshold - let a = coord! { x: 0.0, y: 0.0 }; - assert_eq!(a.try_normalize(f64_minimum_threshold), None); - - // Large vector Pass; - // Note: this fails because the magnitude calculation overflows to infinity - // this could be avoided my modifying the implementation to use scaling be avoided using scaling + // A large vector where try_normalize returns Some + // Because the magnitude is f64::MAX (Just before overflow to f64::INFINITY) let a = coord! { x: f64::sqrt(f64::MAX/2.0), y: f64::sqrt(f64::MAX/2.0) }; - assert!(a.try_normalize(f64_minimum_threshold).is_some()); + assert!(a.try_normalize().is_some()); - // Large Vector Fail; - // Note: This test uses the unstable float_next_up_down feature - // let a = coord! { x: f64::sqrt(f64::MAX/2.0), y: f64::next_up(f64::sqrt(f64::MAX/2.0)) }; - // assert!(a.try_normalize(1e-10f64).is_none()); + // A large vector where try_normalize returns None + // because the magnitude is just above f64::MAX + let a = coord! { + x: f64::sqrt(f64::MAX / 2.0), + y: f64::sqrt(f64::MAX / 2.0).next_after(f64::INFINITY) + }; + assert_eq!(a.try_normalize(), None); - // NaN + // Where one of the components is NaN try_normalize returns None let a = coord! { x: f64::NAN, y: 0.0 }; - assert_eq!(a.try_normalize(f64_minimum_threshold), None); + assert_eq!(a.try_normalize(), None); - // Infinite + // Where one of the components is Infinite try_normalize returns None let a = coord! { x: f64::INFINITY, y: 0.0 }; - assert_eq!(a.try_normalize(f64_minimum_threshold), None); + assert_eq!(a.try_normalize(), None); } } From 12155aaea512b5ea91a00059942a17c1d63f928f Mon Sep 17 00:00:00 2001 From: thehappycheese Date: Sun, 30 Jul 2023 13:00:29 +0800 Subject: [PATCH 11/11] cargo fmt & check line endings are LF --- geo/src/algorithm/vector_ops.rs | 730 ++++++++++++++++---------------- 1 file changed, 363 insertions(+), 367 deletions(-) diff --git a/geo/src/algorithm/vector_ops.rs b/geo/src/algorithm/vector_ops.rs index 72feb9ad21..90ed6df772 100644 --- a/geo/src/algorithm/vector_ops.rs +++ b/geo/src/algorithm/vector_ops.rs @@ -1,367 +1,363 @@ -//! This module defines the [Vector2DOps] trait and implements it for the -//! [Coord] struct. - -use crate::{Coord, CoordFloat, CoordNum}; - -/// Defines vector operations for 2D coordinate types which implement CoordFloat -/// -/// This trait is intended for internal use within the geo crate as a way to -/// bring together the various hand-crafted linear algebra operations used -/// throughout other algorithms and attached to various structs. -/// -/// - -pub trait Vector2DOps -where - Self: Sized, -{ - type Scalar: CoordNum + Send + Sync; - - /// The euclidean distance between this coordinate and the origin - /// - /// `sqrt(x² + y²)` - /// - fn magnitude(self) -> Self::Scalar; - - /// The squared distance between this coordinate and the origin. - /// (Avoids the square root calculation when it is not needed) - /// - /// `x² + y²` - /// - fn magnitude_squared(self) -> Self::Scalar; - - /// Rotate this coordinate around the origin by 90 degrees clockwise. - /// - /// `a.left() => (-a.y, a.x)` - /// - /// Assumes a coordinate system where positive `y` is up and positive `x` is - /// to the right. The described rotation direction is consistent with the - /// documentation for [crate::algorithm::rotate::Rotate]. - fn left(self) -> Self; - - /// Rotate this coordinate around the origin by 90 degrees anti-clockwise. - /// - /// `a.right() => (a.y, -a.x)` - /// - /// Assumes a coordinate system where positive `y` is up and positive `x` is - /// to the right. The described rotation direction is consistent with the - /// documentation for [crate::algorithm::rotate::Rotate]. - fn right(self) -> Self; - - /// The inner product of the coordinate components - /// - /// `a · b = a.x * b.x + a.y * b.y` - /// - fn dot_product(self, other: Rhs) -> Self::Scalar; - - /// The calculates the `wedge product` between two vectors. - /// - /// `a ∧ b = a.x * b.y - a.y * b.x` - /// - /// Also known as: - /// - /// - `exterior product` - /// - because the wedge product comes from 'Exterior Algebra' - /// - `perpendicular product` - /// - because it is equivalent to `a.dot(b.right())` - /// - `2D cross product` - /// - because it is equivalent to the signed magnitude of the - /// conventional 3D cross product assuming `z` ordinates are zero - /// - `determinant` - /// - because it is equivalent to the `determinant` of the 2x2 matrix - /// formed by the column-vector inputs. - /// - /// ## Examples - /// - /// The following list highlights some examples in geo which might be - /// brought together to use this function: - /// - /// 1. [geo_types::Point::cross_prod()] is already defined on - /// [geo_types::Point]... but that it seems to be some other - /// operation on 3 points?? - /// 2. [geo_types::Line] struct also has a [geo_types::Line::determinant()] - /// function which is the same as `line.start.wedge_product(line.end)` - /// 3. The [crate::algorithm::Kernel::orient2d()] trait default - /// implementation uses cross product to compute orientation. It returns - /// an enum, not the numeric value which is needed for line segment - /// intersection. - /// - /// ## Properties - /// - /// - The absolute value of the cross product is the area of the - /// parallelogram formed by the operands - /// - Anti-commutative: The sign of the output is reversed if the operands - /// are reversed - /// - If the operands are colinear with the origin, the value is zero - /// - The sign can be used to check if the operands are clockwise with - /// respect to the origin, or phrased differently: - /// "is a to the left of the line between the origin and b"? - /// - If this is what you are using it for, then please use - /// [crate::algorithm::Kernel::orient2d()] instead as this is more - /// explicit and has a `RobustKernel` option for extra precision. - fn wedge_product(self, other: Rhs) -> Self::Scalar; - - /// Try to find a vector of unit length in the same direction as this - /// vector. - /// - /// Returns `None` if the result is not finite. This can happen when - /// - /// - the vector is really small (or zero length) and the `.magnitude()` - /// calculation has rounded-down to `0.0` - /// - the vector is really large and the `.magnitude()` has rounded-up - /// or 'overflowed' to `f64::INFINITY` - /// - Either x or y are `f64::NAN` or `f64::INFINITY` - fn try_normalize(self) -> Option; - - /// Returns true if both the x and y components are finite - fn is_finite(self) -> bool; -} - -impl Vector2DOps for Coord -where - T: CoordFloat + Send + Sync, -{ - type Scalar = T; - - fn wedge_product(self, right: Coord) -> Self::Scalar { - self.x * right.y - self.y * right.x - } - - fn dot_product(self, other: Self) -> Self::Scalar { - self.x * other.x + self.y * other.y - } - - fn magnitude(self) -> Self::Scalar { - (self.x * self.x + self.y * self.y).sqrt() - } - - fn magnitude_squared(self) -> Self::Scalar { - self.x * self.x + self.y * self.y - } - - fn left(self) -> Self { - Self { - x: -self.y, - y: self.x, - } - } - - fn right(self) -> Self { - Self { - x: self.y, - y: -self.x, - } - } - - fn try_normalize(self) -> Option { - let magnitude = self.magnitude(); - let result = self / magnitude; - // Both the result AND the magnitude must be finite they are finite - // Otherwise very large vectors overflow magnitude to Infinity, - // and the after the division the result would be coord!{x:0.0,y:0.0} - // Note we don't need to check if magnitude is zero, because after the division - // that would have made result non-finite or NaN anyway. - if result.is_finite() && magnitude.is_finite() { - Some(result) - } else { - None - } - } - - fn is_finite(self) -> bool { - self.x.is_finite() && self.y.is_finite() - } -} - -#[cfg(test)] -mod test { - use super::Vector2DOps; - use crate::coord; - - #[test] - fn test_cross_product() { - // perpendicular unit length - let a = coord! { x: 1f64, y: 0f64 }; - let b = coord! { x: 0f64, y: 1f64 }; - - // expect the area of parallelogram - assert_eq!(a.wedge_product(b), 1f64); - // expect swapping will result in negative - assert_eq!(b.wedge_product(a), -1f64); - - // Add skew; expect results should be the same - let a = coord! { x: 1f64, y: 0f64 }; - let b = coord! { x: 1f64, y: 1f64 }; - - // expect the area of parallelogram - assert_eq!(a.wedge_product(b), 1f64); - // expect swapping will result in negative - assert_eq!(b.wedge_product(a), -1f64); - - // Make Colinear; expect zero - let a = coord! { x: 2f64, y: 2f64 }; - let b = coord! { x: 1f64, y: 1f64 }; - assert_eq!(a.wedge_product(b), 0f64); - } - - #[test] - fn test_dot_product() { - // perpendicular unit length - let a = coord! { x: 1f64, y: 0f64 }; - let b = coord! { x: 0f64, y: 1f64 }; - // expect zero for perpendicular - assert_eq!(a.dot_product(b), 0f64); - - // Parallel, same direction - let a = coord! { x: 1f64, y: 0f64 }; - let b = coord! { x: 2f64, y: 0f64 }; - // expect +ive product of magnitudes - assert_eq!(a.dot_product(b), 2f64); - // expect swapping will have same result - assert_eq!(b.dot_product(a), 2f64); - - // Parallel, opposite direction - let a = coord! { x: 3f64, y: 4f64 }; - let b = coord! { x: -3f64, y: -4f64 }; - // expect -ive product of magnitudes - assert_eq!(a.dot_product(b), -25f64); - // expect swapping will have same result - assert_eq!(b.dot_product(a), -25f64); - } - - #[test] - fn test_magnitude() { - let a = coord! { x: 1f64, y: 0f64 }; - assert_eq!(a.magnitude(), 1f64); - - let a = coord! { x: 0f64, y: 0f64 }; - assert_eq!(a.magnitude(), 0f64); - - let a = coord! { x: -3f64, y: 4f64 }; - assert_eq!(a.magnitude(), 5f64); - } - - #[test] - fn test_magnitude_squared() { - let a = coord! { x: 1f64, y: 0f64 }; - assert_eq!(a.magnitude_squared(), 1f64); - - let a = coord! { x: 0f64, y: 0f64 }; - assert_eq!(a.magnitude_squared(), 0f64); - - let a = coord! { x: -3f64, y: 4f64 }; - assert_eq!(a.magnitude_squared(), 25f64); - } - - #[test] - fn test_left_right() { - let a = coord! { x: 1f64, y: 0f64 }; - let a_left = coord! { x: 0f64, y: 1f64 }; - let a_right = coord! { x: 0f64, y: -1f64 }; - - assert_eq!(a.left(), a_left); - assert_eq!(a.right(), a_right); - assert_eq!(a.left(), -a.right()); - } - - #[test] - fn test_left_right_match_rotate() { - use crate::algorithm::rotate::Rotate; - use crate::Point; - // The aim of this test is to confirm that wording in documentation is - // consistent. - - // when the user is in a coordinate system where the y axis is flipped - // (eg screen coordinates in a HTML canvas), then rotation directions - // will be different to those described in the documentation. - - // The documentation for the Rotate trait says: 'Positive angles are - // counter-clockwise, and negative angles are clockwise rotations' - - let counter_clockwise_rotation_degrees = 90.0; - let clockwise_rotation_degrees = -counter_clockwise_rotation_degrees; - - let a: Point = coord! { x: 1.0, y: 0.0 }.into(); - let origin: Point = coord! { x: 0.0, y: 0.0 }.into(); - - // left is anti-clockwise - assert_relative_eq!( - Point::from(a.0.left()), - a.rotate_around_point(counter_clockwise_rotation_degrees, origin), - ); - // right is clockwise - assert_relative_eq!( - Point::from(a.0.right()), - a.rotate_around_point(clockwise_rotation_degrees, origin), - ); - } - - #[test] - fn test_try_normalize() { - - // Already Normalized - let a = coord! { - x: 1.0, - y: 0.0 - }; - assert_relative_eq!(a.try_normalize().unwrap(), a); - - // Already Normalized - let a = coord! { - x: 1.0 / f64::sqrt(2.0), - y: -1.0 / f64::sqrt(2.0) - }; - assert_relative_eq!(a.try_normalize().unwrap(), a); - - // Non trivial example - let a = coord! { x: -10.0, y: 8.0 }; - assert_relative_eq!( - a.try_normalize().unwrap(), - coord! { x: -10.0, y: 8.0 } / f64::sqrt(10.0 * 10.0 + 8.0 * 8.0) - ); - } - - #[test] - fn test_try_normalize_edge_cases() { - use float_next_after::NextAfter; - - // The following tests demonstrate some of the floating point - // edge cases that can cause try_normalize to return None. - - // Zero vector - Normalize returns None - let a = coord! { x: 0.0, y: 0.0 }; - assert_eq!(a.try_normalize(), None); - - // Very Small Input - Normalize returns None because of - // rounding-down to zero in the .magnitude() calculation - let a = coord! { - x: 0.0, - y: 1e-301_f64 - }; - assert_eq!(a.try_normalize(), None); - - // A large vector where try_normalize returns Some - // Because the magnitude is f64::MAX (Just before overflow to f64::INFINITY) - let a = coord! { - x: f64::sqrt(f64::MAX/2.0), - y: f64::sqrt(f64::MAX/2.0) - }; - assert!(a.try_normalize().is_some()); - - // A large vector where try_normalize returns None - // because the magnitude is just above f64::MAX - let a = coord! { - x: f64::sqrt(f64::MAX / 2.0), - y: f64::sqrt(f64::MAX / 2.0).next_after(f64::INFINITY) - }; - assert_eq!(a.try_normalize(), None); - - // Where one of the components is NaN try_normalize returns None - let a = coord! { x: f64::NAN, y: 0.0 }; - assert_eq!(a.try_normalize(), None); - - // Where one of the components is Infinite try_normalize returns None - let a = coord! { x: f64::INFINITY, y: 0.0 }; - assert_eq!(a.try_normalize(), None); - } -} +//! This module defines the [Vector2DOps] trait and implements it for the +//! [Coord] struct. + +use crate::{Coord, CoordFloat, CoordNum}; + +/// Defines vector operations for 2D coordinate types which implement CoordFloat +/// +/// This trait is intended for internal use within the geo crate as a way to +/// bring together the various hand-crafted linear algebra operations used +/// throughout other algorithms and attached to various structs. +pub trait Vector2DOps +where + Self: Sized, +{ + type Scalar: CoordNum + Send + Sync; + + /// The euclidean distance between this coordinate and the origin + /// + /// `sqrt(x² + y²)` + /// + fn magnitude(self) -> Self::Scalar; + + /// The squared distance between this coordinate and the origin. + /// (Avoids the square root calculation when it is not needed) + /// + /// `x² + y²` + /// + fn magnitude_squared(self) -> Self::Scalar; + + /// Rotate this coordinate around the origin by 90 degrees clockwise. + /// + /// `a.left() => (-a.y, a.x)` + /// + /// Assumes a coordinate system where positive `y` is up and positive `x` is + /// to the right. The described rotation direction is consistent with the + /// documentation for [crate::algorithm::rotate::Rotate]. + fn left(self) -> Self; + + /// Rotate this coordinate around the origin by 90 degrees anti-clockwise. + /// + /// `a.right() => (a.y, -a.x)` + /// + /// Assumes a coordinate system where positive `y` is up and positive `x` is + /// to the right. The described rotation direction is consistent with the + /// documentation for [crate::algorithm::rotate::Rotate]. + fn right(self) -> Self; + + /// The inner product of the coordinate components + /// + /// `a · b = a.x * b.x + a.y * b.y` + /// + fn dot_product(self, other: Rhs) -> Self::Scalar; + + /// The calculates the `wedge product` between two vectors. + /// + /// `a ∧ b = a.x * b.y - a.y * b.x` + /// + /// Also known as: + /// + /// - `exterior product` + /// - because the wedge product comes from 'Exterior Algebra' + /// - `perpendicular product` + /// - because it is equivalent to `a.dot(b.right())` + /// - `2D cross product` + /// - because it is equivalent to the signed magnitude of the + /// conventional 3D cross product assuming `z` ordinates are zero + /// - `determinant` + /// - because it is equivalent to the `determinant` of the 2x2 matrix + /// formed by the column-vector inputs. + /// + /// ## Examples + /// + /// The following list highlights some examples in geo which might be + /// brought together to use this function: + /// + /// 1. [geo_types::Point::cross_prod()] is already defined on + /// [geo_types::Point]... but that it seems to be some other + /// operation on 3 points?? + /// 2. [geo_types::Line] struct also has a [geo_types::Line::determinant()] + /// function which is the same as `line.start.wedge_product(line.end)` + /// 3. The [crate::algorithm::Kernel::orient2d()] trait default + /// implementation uses cross product to compute orientation. It returns + /// an enum, not the numeric value which is needed for line segment + /// intersection. + /// + /// ## Properties + /// + /// - The absolute value of the cross product is the area of the + /// parallelogram formed by the operands + /// - Anti-commutative: The sign of the output is reversed if the operands + /// are reversed + /// - If the operands are colinear with the origin, the value is zero + /// - The sign can be used to check if the operands are clockwise with + /// respect to the origin, or phrased differently: + /// "is a to the left of the line between the origin and b"? + /// - If this is what you are using it for, then please use + /// [crate::algorithm::Kernel::orient2d()] instead as this is more + /// explicit and has a `RobustKernel` option for extra precision. + fn wedge_product(self, other: Rhs) -> Self::Scalar; + + /// Try to find a vector of unit length in the same direction as this + /// vector. + /// + /// Returns `None` if the result is not finite. This can happen when + /// + /// - the vector is really small (or zero length) and the `.magnitude()` + /// calculation has rounded-down to `0.0` + /// - the vector is really large and the `.magnitude()` has rounded-up + /// or 'overflowed' to `f64::INFINITY` + /// - Either x or y are `f64::NAN` or `f64::INFINITY` + fn try_normalize(self) -> Option; + + /// Returns true if both the x and y components are finite + fn is_finite(self) -> bool; +} + +impl Vector2DOps for Coord +where + T: CoordFloat + Send + Sync, +{ + type Scalar = T; + + fn wedge_product(self, right: Coord) -> Self::Scalar { + self.x * right.y - self.y * right.x + } + + fn dot_product(self, other: Self) -> Self::Scalar { + self.x * other.x + self.y * other.y + } + + fn magnitude(self) -> Self::Scalar { + (self.x * self.x + self.y * self.y).sqrt() + } + + fn magnitude_squared(self) -> Self::Scalar { + self.x * self.x + self.y * self.y + } + + fn left(self) -> Self { + Self { + x: -self.y, + y: self.x, + } + } + + fn right(self) -> Self { + Self { + x: self.y, + y: -self.x, + } + } + + fn try_normalize(self) -> Option { + let magnitude = self.magnitude(); + let result = self / magnitude; + // Both the result AND the magnitude must be finite they are finite + // Otherwise very large vectors overflow magnitude to Infinity, + // and the after the division the result would be coord!{x:0.0,y:0.0} + // Note we don't need to check if magnitude is zero, because after the division + // that would have made result non-finite or NaN anyway. + if result.is_finite() && magnitude.is_finite() { + Some(result) + } else { + None + } + } + + fn is_finite(self) -> bool { + self.x.is_finite() && self.y.is_finite() + } +} + +#[cfg(test)] +mod test { + use super::Vector2DOps; + use crate::coord; + + #[test] + fn test_cross_product() { + // perpendicular unit length + let a = coord! { x: 1f64, y: 0f64 }; + let b = coord! { x: 0f64, y: 1f64 }; + + // expect the area of parallelogram + assert_eq!(a.wedge_product(b), 1f64); + // expect swapping will result in negative + assert_eq!(b.wedge_product(a), -1f64); + + // Add skew; expect results should be the same + let a = coord! { x: 1f64, y: 0f64 }; + let b = coord! { x: 1f64, y: 1f64 }; + + // expect the area of parallelogram + assert_eq!(a.wedge_product(b), 1f64); + // expect swapping will result in negative + assert_eq!(b.wedge_product(a), -1f64); + + // Make Colinear; expect zero + let a = coord! { x: 2f64, y: 2f64 }; + let b = coord! { x: 1f64, y: 1f64 }; + assert_eq!(a.wedge_product(b), 0f64); + } + + #[test] + fn test_dot_product() { + // perpendicular unit length + let a = coord! { x: 1f64, y: 0f64 }; + let b = coord! { x: 0f64, y: 1f64 }; + // expect zero for perpendicular + assert_eq!(a.dot_product(b), 0f64); + + // Parallel, same direction + let a = coord! { x: 1f64, y: 0f64 }; + let b = coord! { x: 2f64, y: 0f64 }; + // expect +ive product of magnitudes + assert_eq!(a.dot_product(b), 2f64); + // expect swapping will have same result + assert_eq!(b.dot_product(a), 2f64); + + // Parallel, opposite direction + let a = coord! { x: 3f64, y: 4f64 }; + let b = coord! { x: -3f64, y: -4f64 }; + // expect -ive product of magnitudes + assert_eq!(a.dot_product(b), -25f64); + // expect swapping will have same result + assert_eq!(b.dot_product(a), -25f64); + } + + #[test] + fn test_magnitude() { + let a = coord! { x: 1f64, y: 0f64 }; + assert_eq!(a.magnitude(), 1f64); + + let a = coord! { x: 0f64, y: 0f64 }; + assert_eq!(a.magnitude(), 0f64); + + let a = coord! { x: -3f64, y: 4f64 }; + assert_eq!(a.magnitude(), 5f64); + } + + #[test] + fn test_magnitude_squared() { + let a = coord! { x: 1f64, y: 0f64 }; + assert_eq!(a.magnitude_squared(), 1f64); + + let a = coord! { x: 0f64, y: 0f64 }; + assert_eq!(a.magnitude_squared(), 0f64); + + let a = coord! { x: -3f64, y: 4f64 }; + assert_eq!(a.magnitude_squared(), 25f64); + } + + #[test] + fn test_left_right() { + let a = coord! { x: 1f64, y: 0f64 }; + let a_left = coord! { x: 0f64, y: 1f64 }; + let a_right = coord! { x: 0f64, y: -1f64 }; + + assert_eq!(a.left(), a_left); + assert_eq!(a.right(), a_right); + assert_eq!(a.left(), -a.right()); + } + + #[test] + fn test_left_right_match_rotate() { + use crate::algorithm::rotate::Rotate; + use crate::Point; + // The aim of this test is to confirm that wording in documentation is + // consistent. + + // when the user is in a coordinate system where the y axis is flipped + // (eg screen coordinates in a HTML canvas), then rotation directions + // will be different to those described in the documentation. + + // The documentation for the Rotate trait says: 'Positive angles are + // counter-clockwise, and negative angles are clockwise rotations' + + let counter_clockwise_rotation_degrees = 90.0; + let clockwise_rotation_degrees = -counter_clockwise_rotation_degrees; + + let a: Point = coord! { x: 1.0, y: 0.0 }.into(); + let origin: Point = coord! { x: 0.0, y: 0.0 }.into(); + + // left is anti-clockwise + assert_relative_eq!( + Point::from(a.0.left()), + a.rotate_around_point(counter_clockwise_rotation_degrees, origin), + ); + // right is clockwise + assert_relative_eq!( + Point::from(a.0.right()), + a.rotate_around_point(clockwise_rotation_degrees, origin), + ); + } + + #[test] + fn test_try_normalize() { + // Already Normalized + let a = coord! { + x: 1.0, + y: 0.0 + }; + assert_relative_eq!(a.try_normalize().unwrap(), a); + + // Already Normalized + let a = coord! { + x: 1.0 / f64::sqrt(2.0), + y: -1.0 / f64::sqrt(2.0) + }; + assert_relative_eq!(a.try_normalize().unwrap(), a); + + // Non trivial example + let a = coord! { x: -10.0, y: 8.0 }; + assert_relative_eq!( + a.try_normalize().unwrap(), + coord! { x: -10.0, y: 8.0 } / f64::sqrt(10.0 * 10.0 + 8.0 * 8.0) + ); + } + + #[test] + fn test_try_normalize_edge_cases() { + use float_next_after::NextAfter; + + // The following tests demonstrate some of the floating point + // edge cases that can cause try_normalize to return None. + + // Zero vector - Normalize returns None + let a = coord! { x: 0.0, y: 0.0 }; + assert_eq!(a.try_normalize(), None); + + // Very Small Input - Normalize returns None because of + // rounding-down to zero in the .magnitude() calculation + let a = coord! { + x: 0.0, + y: 1e-301_f64 + }; + assert_eq!(a.try_normalize(), None); + + // A large vector where try_normalize returns Some + // Because the magnitude is f64::MAX (Just before overflow to f64::INFINITY) + let a = coord! { + x: f64::sqrt(f64::MAX/2.0), + y: f64::sqrt(f64::MAX/2.0) + }; + assert!(a.try_normalize().is_some()); + + // A large vector where try_normalize returns None + // because the magnitude is just above f64::MAX + let a = coord! { + x: f64::sqrt(f64::MAX / 2.0), + y: f64::sqrt(f64::MAX / 2.0).next_after(f64::INFINITY) + }; + assert_eq!(a.try_normalize(), None); + + // Where one of the components is NaN try_normalize returns None + let a = coord! { x: f64::NAN, y: 0.0 }; + assert_eq!(a.try_normalize(), None); + + // Where one of the components is Infinite try_normalize returns None + let a = coord! { x: f64::INFINITY, y: 0.0 }; + assert_eq!(a.try_normalize(), None); + } +}