From e9b3012520ba580b6e798512b2318df48d45c3b7 Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Thu, 28 May 2020 14:42:37 -0400 Subject: euclid: Implement Euclidean distance --- src/euclid.rs | 390 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/lib.rs | 2 + 2 files changed, 392 insertions(+) create mode 100644 src/euclid.rs (limited to 'src') diff --git a/src/euclid.rs b/src/euclid.rs new file mode 100644 index 0000000..fa82e9c --- /dev/null +++ b/src/euclid.rs @@ -0,0 +1,390 @@ +//! Euclidean space. + +use crate::coords::{Coordinates, CoordinateMetric, CoordinateProximity}; +use crate::distance::{Distance, Metric, Proximity, Value}; + +use num_traits::zero; + +use std::cmp::Ordering; +use std::convert::TryFrom; + +/// A point in Euclidean space. +/// +/// This wrapper equips any [coordinate space] with the [Euclidean distance] metric. +/// +/// [coordinate space]: [Coordinates] +/// [Euclidean distance]: https://en.wikipedia.org/wiki/Euclidean_distance +#[derive(Clone, Copy, Debug, Eq, PartialEq)] +pub struct Euclidean(pub T); + +impl Euclidean { + /// Wrap a point. + pub fn new(point: T) -> Self { + Self(point) + } + + /// Unwrap a point. + pub fn inner(&self) -> &T { + &self.0 + } + + /// Unwrap a point. + pub fn into_inner(self) -> T { + self.0 + } +} + +impl Coordinates for Euclidean { + type Value = T::Value; + + fn dims(&self) -> usize { + self.0.dims() + } + + fn coord(&self, i: usize) -> Self::Value { + self.0.coord(i) + } +} + +/// Compute the Euclidean distance between two points. +pub fn euclidean_distance(x: T, y: U) -> EuclideanDistance +where + T: Coordinates, + U: Coordinates, +{ + debug_assert!(x.dims() == y.dims()); + + let mut sum = zero(); + for i in 0..x.dims() { + let diff = x.coord(i) - y.coord(i); + sum += diff * diff; + } + + EuclideanDistance::from_squared(sum) +} + +/// The Euclidean distance function. +impl Proximity for Euclidean +where + T: Coordinates, + EuclideanDistance: Distance, +{ + type Distance = EuclideanDistance; + + fn distance(&self, other: &Self) -> Self::Distance { + euclidean_distance(self, other) + } +} + +impl Proximity for Euclidean +where + T: Coordinates, + EuclideanDistance: Distance, +{ + type Distance = EuclideanDistance; + + fn distance(&self, other: &T) -> Self::Distance { + euclidean_distance(self, other) + } +} + +impl Proximity> for T +where + T: Coordinates, + EuclideanDistance: Distance, +{ + type Distance = EuclideanDistance; + + fn distance(&self, other: &Euclidean) -> Self::Distance { + euclidean_distance(self, other) + } +} + +/// Euclidean distance is a metric. +impl Metric for Euclidean +where + T: Coordinates, + EuclideanDistance: Distance, +{} + +impl Metric for Euclidean +where + T: Coordinates, + EuclideanDistance: Distance, +{} + +impl Metric> for T +where + T: Coordinates, + EuclideanDistance: Distance, +{} + +impl CoordinateProximity for Euclidean +where + T: Coordinates, + EuclideanDistance: Distance, +{ + type Distance = EuclideanDistance; + + fn distance_to_coords(&self, coords: &[T::Value]) -> Self::Distance { + euclidean_distance(self, coords) + } +} + +impl CoordinateMetric for Euclidean +where + T: Coordinates, + EuclideanDistance: Distance, +{} + +/// A [Euclidean distance]. +/// +/// This type stores the squared value of the Euclidean distance, to avoid computing expensive +/// square roots until absolutely necessary. +/// +/// # use acap::distance::Distance; +/// # use acap::euclid::EuclideanDistance; +/// # use std::convert::TryFrom; +/// let a = EuclideanDistance::try_from(3).unwrap(); +/// let b = EuclideanDistance::try_from(4).unwrap(); +/// let c = EuclideanDistance::from_squared(a.squared_value() + b.squared_value()); +/// assert!(a < c && b < c); +/// assert_eq!(c.value(), 5.0f32); +/// +/// [Euclidean distance]: https://en.wikipedia.org/wiki/Euclidean_distance +#[derive(Clone, Copy, Debug, PartialEq, PartialOrd)] +pub struct EuclideanDistance(T); + +impl EuclideanDistance { + /// Createa a `EuclideanDistance` from an already-squared value. + pub fn from_squared(value: T) -> Self { + debug_assert!(!value.is_negative()); + Self(value) + } + + /// Get the squared distance value. + pub fn squared_value(self) -> T { + self.0 + } +} + +/// Error type for failed conversions from negative numbers to [EuclideanDistance]. +#[derive(Debug)] +pub struct NegativeDistanceError; + +/// Implement EuclideanDistance for a floating-point type. +macro_rules! float_distance { + ($f:ty) => { + impl TryFrom<$f> for EuclideanDistance<$f> { + type Error = NegativeDistanceError; + + #[inline] + fn try_from(value: $f) -> Result { + if value >= 0.0 { + Ok(Self(value * value)) + } else { + Err(NegativeDistanceError) + } + } + } + + impl From> for $f { + #[inline] + fn from(value: EuclideanDistance<$f>) -> $f { + value.0.sqrt() + } + } + + impl PartialOrd<$f> for EuclideanDistance<$f> { + #[inline] + fn partial_cmp(&self, other: &$f) -> Option { + if let Ok(rhs) = Self::try_from(*other) { + self.partial_cmp(&rhs) + } else { + Some(Ordering::Greater) + } + } + } + + impl PartialOrd> for $f { + #[inline] + fn partial_cmp(&self, other: &EuclideanDistance<$f>) -> Option { + if let Ok(lhs) = EuclideanDistance::try_from(*self) { + lhs.partial_cmp(other) + } else { + Some(Ordering::Less) + } + } + } + + impl PartialEq<$f> for EuclideanDistance<$f> { + #[inline] + fn eq(&self, other: &$f) -> bool { + self.partial_cmp(other) == Some(Ordering::Equal) + } + } + + impl PartialEq> for $f { + #[inline] + fn eq(&self, other: &EuclideanDistance<$f>) -> bool { + self.partial_cmp(other) == Some(Ordering::Equal) + } + } + + impl Distance for EuclideanDistance<$f> { + type Value = $f; + } + } +} + +float_distance!(f32); +float_distance!(f64); + +/// Implement EuclideanDistance for an integer type. +macro_rules! int_distance { + ($i:ty, $f:ty, $ff:ty) => { + impl TryFrom<$i> for EuclideanDistance<$i> { + type Error = NegativeDistanceError; + + #[inline] + fn try_from(value: $i) -> Result { + if value >= 0 { + Ok(Self(value * value)) + } else { + Err(NegativeDistanceError) + } + } + } + + impl From> for $f { + #[inline] + fn from(value: EuclideanDistance<$i>) -> Self { + (value.0 as $ff).sqrt() as $f + } + } + + impl PartialOrd<$i> for EuclideanDistance<$i> { + #[inline] + fn partial_cmp(&self, other: &$i) -> Option { + if let Ok(rhs) = Self::try_from(*other) { + self.partial_cmp(&rhs) + } else { + Some(Ordering::Greater) + } + } + } + + impl PartialOrd> for $i { + #[inline] + fn partial_cmp(&self, other: &EuclideanDistance<$i>) -> Option { + if let Ok(lhs) = EuclideanDistance::try_from(*self) { + lhs.partial_cmp(other) + } else { + Some(Ordering::Less) + } + } + } + + impl PartialEq<$i> for EuclideanDistance<$i> { + #[inline] + fn eq(&self, other: &$i) -> bool { + self.partial_cmp(other) == Some(Ordering::Equal) + } + } + + impl PartialEq> for $i { + #[inline] + fn eq(&self, other: &EuclideanDistance<$i>) -> bool { + self.partial_cmp(other) == Some(Ordering::Equal) + } + } + + impl PartialOrd<$f> for EuclideanDistance<$i> { + #[inline] + fn partial_cmp(&self, other: &$f) -> Option { + if *other >= 0.0 { + let lhs = self.0 as $ff; + let mut rhs = *other as $ff; + rhs *= rhs; + lhs.partial_cmp(&rhs) + } else { + Some(Ordering::Greater) + } + } + } + + impl PartialOrd> for $f { + #[inline] + fn partial_cmp(&self, other: &EuclideanDistance<$i>) -> Option { + if *other >= 0.0 { + let mut lhs = *self as $ff; + lhs *= lhs; + let rhs = other.0 as $ff; + lhs.partial_cmp(&rhs) + } else { + Some(Ordering::Greater) + } + } + } + + impl PartialEq<$f> for EuclideanDistance<$i> { + #[inline] + fn eq(&self, other: &$f) -> bool { + self.partial_cmp(other) == Some(Ordering::Equal) + } + } + + impl PartialEq> for $f { + #[inline] + fn eq(&self, other: &EuclideanDistance<$i>) -> bool { + self.partial_cmp(other) == Some(Ordering::Equal) + } + } + + impl Distance for EuclideanDistance<$i> { + type Value = $f; + } + } +} + +int_distance!(i16, f32, f32); +int_distance!(i32, f32, f64); +int_distance!(i64, f64, f64); +int_distance!(isize, f64, f64); + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_i32() { + let five = euclidean_distance([0, 0], [3, 4]); + assert_eq!(five, EuclideanDistance::from_squared(25)); + assert_eq!(five, 5.0f32); + + let thirteen = Euclidean([0, 0]).distance(&Euclidean([5, 12])); + assert_eq!(thirteen, EuclideanDistance::from_squared(169)); + assert_eq!(thirteen, 13.0f32); + + assert!(five < thirteen); + assert!(five < 13); + assert!(5 < thirteen); + assert!(-5 < thirteen); + } + + #[test] + fn test_f64() { + let five = euclidean_distance([0.0, 0.0], [3.0, 4.0]); + assert_eq!(five, EuclideanDistance::from_squared(25.0)); + assert_eq!(five, 5.0); + + let thirteen = Euclidean([0.0, 0.0]).distance(&Euclidean([5.0, 12.0])); + assert_eq!(thirteen, EuclideanDistance::from_squared(169.0)); + assert_eq!(thirteen, 13.0); + + assert!(five < thirteen); + assert!(five < 13.0); + assert!(5.0 < thirteen); + assert!(-5.0 < thirteen); + } +} diff --git a/src/lib.rs b/src/lib.rs index 1532211..e18025b 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -4,6 +4,8 @@ pub mod coords; pub mod distance; +pub mod euclid; pub use coords::Coordinates; pub use distance::{Distance, Metric, Proximity}; +pub use euclid::{euclidean_distance, Euclidean, EuclideanDistance}; -- cgit v1.2.3