Skip to content

Commit

Permalink
Merge #1025
Browse files Browse the repository at this point in the history 
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 <[email protected]>
  • Loading branch information
@bors @thehappycheese
bors[bot] and thehappycheese authored Jul 31, 2023
2 parents 2ab2b27 + 10b7f08 commit 707ffa0
Show file tree
Hide file tree
Showing 3 changed files with 372 additions and 0 deletions.
2 changes: 2 additions & 0 deletions geo/CHANGES.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
4 changes: 4 additions & 0 deletions geo/src/algorithm/mod.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
366 changes: 366 additions & 0 deletions geo/src/algorithm/vector_ops.rs
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,366 @@
//! 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
// Annotation to disable bad clippy lint; It is not good to use
// `&self` as clippy suggests since Coord is Copy
#[allow(clippy::wrong_self_convention)]
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);
}
}

0 comments on commit 707ffa0

Please sign in to comment.