Merge #1025

1025: Vector2DOps Trait - Proposal r=michaelkirk a=thehappycheese

- [x] I agree to follow the project's [code of conduct](https://github.com/georust/geo/blob/main/CODE_OF_CONDUCT.md).
- [x] I added an entry to `CHANGES.md` if knowledge of this change could be valuable to users.
---

Proposal to define `Vector2DOps` trait and implement it for the `Coord` struct.
Intention is to facilitate general purpose linear algebra operations such as 'dot product' etc all in one place.

Most of these functions are are already implemented in various places throughout the crate or sometimes hard-coded in place. This would bring them all together.

For example

- `wedge_product`( aka cross/perp/exterior product) `geotypes::Line::determinant`,
  - `geotypes::Point::cross_prod` (which takes three arguments for some reason) etc.
  - Also used in the default implementation of Kernel orient2d, although this serves a different purpose as it returns an enum
- `dot_product` `geotypes::Point::dot`,
- `magnitude` in `Line::line_euclidean_length`
- `magnitude_squared` hard coded ( [here](https://github.com/thehappycheese/geo/blob/9b477fc7b24fcb8731bd5e01d559b613e79ada0c/geo/src/algorithm/line_locate_point.rs#L61C9-L61C29))
- `normalize` ... interestingly it is hard to find an example of this operation in the existing crate. Possibly because it is generally best avoided due to floating point instability.

For more discussion about motivation for this proposal, see [this comment](#935 (comment)) on PR #935 . 

Possibly the timing of this proposal is weird, I am sure that the `geotraits::CoordTrait` will affect how this should be implemented. I suspect it may allow the `Vector2DOps` to define default implementations?

Note:

 - Currently `Vector2DOps` isn't used anywhere so it will trigger dead code warnings.
 - For an example that illustrates how I imagine it would be used in the OffsetCurve trait, see [this snippet](https://github.com/thehappycheese/geo/blob/d95420f1f8d6b0d2dd1e6d803fe3d3e2a2b3dd13/geo/src/algorithm/offset_curve/offset_line_raw.rs#L46C1-L51C1) from this PR #935 
   - Note that PR #935 is a separate, older branch, in that branch I was calling it `VectorExtensions`. In this PR I renamed it to `Vector2DOps` so that it is more similar to the existing `AffineOps`

Many thanks for your time and any feedback you can provide


Co-authored-by: thehappycheese <the.nicholas.archer@gmail.com>

geo/CHANGES.md

@@ -11,6 +11,8 @@
 
 - Add `TriangulateEarcut` algorithm trait to triangulate polygons with the earcut algorithm.
   - <https://github.com/georust/geo/pull/1007>
+- Add `Vector2DOps` trait to algorithims module and implemented it for `Coord<T::CoordFloat>`
+  - <https://github.com/georust/geo/pull/1025>
 
 ## 0.25.0
 

geo/src/algorithm/mod.rs

@@ -236,6 +236,10 @@ pub mod triangulate_earcut;
 #[cfg(feature = "earcutr")]
 pub use triangulate_earcut::TriangulateEarcut;
 
+/// Vector Operations for 2D coordinates
+mod vector_ops;
+pub use vector_ops::Vector2DOps;
+
 /// Calculate the Vincenty distance between two `Point`s.
 pub mod vincenty_distance;
 pub use vincenty_distance::VincentyDistance;

geo/src/algorithm/vector_ops.rs

@@ -0,0 +1,367 @@
+//! 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<Rhs = Self>
+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<Self>;
+
+    /// Returns true if both the x and y components are finite
+    fn is_finite(self) -> bool;
+}
+
+impl<T> Vector2DOps for Coord<T>
+where
+    T: CoordFloat + Send + Sync,
+{
+    type Scalar = T;
+
+    fn wedge_product(self, right: Coord<T>) -> 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<Self> {
+        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);
+    }
+}