nsys-math-utils 0.1.0

//! Convex hulls

use std;

use crate::*;

/// 2D convex hull
#[derive(Clone, Debug, PartialEq)]
pub struct Hull2 <S> {
  /// Set of unique points sorted in counter-clockwise order for efficient
  /// minimum bounding box algorithm
  points : Vec <Point2 <S>>
}

/// 3D convex hull
#[derive(Clone, Debug, PartialEq)]
pub struct Hull3 <S> {
  points : Vec <Point3 <S>>
}

impl <S> Hull2 <S> {
  /// Create a new 2D convex hull from a bag of points.
  ///
  /// Uses Graham scan algorithm.
  ///
  /// Input must contain at least 1 point:
  /// ```
  /// # use math_utils::geometry::Hull2;
  /// assert!(Hull2::<f32>::from_points (&[]).is_none());
  /// ```
  pub fn from_points (points : &[Point2 <S>]) -> Option <Self> where S : Real {
    // code adapted from scirs2-spatial:
    // <https://github.com/cool-japan/scirs/blob/a176b462aca55e73bd4b25eea83aad10a9f4a2b0/scirs2-spatial/src/convex_hull.rs>
    use std::cmp::Ordering;
    match points.len() {
      0 => return None,
      1 => return Some (Hull2 { points: points.to_vec() }),
      2 => {
        let mut points = points.to_vec();
        points.dedup();
        return Some (Hull2 { points})
      }
      _ => {} // continue
    }
    let mut indexed_points = points.iter().cloned().enumerate()
      .collect::<Vec <(usize, Point2 <S>)>>();
    // find bottom-most (lowest y-coordinate), then left-most x
    let start_index = indexed_points.iter().min_by (|a, b|{
      let cmp = a.1.0.y.partial_cmp (&b.1.0.y).unwrap();
      if cmp == Ordering::Equal {
        a.1.0.x.partial_cmp (&b.1.0.x).unwrap()
      } else {
        cmp
      }
    }).unwrap().0;
    let start_point = indexed_points[start_index].1;
    // sort points by polar angle with respect to start point
    indexed_points.sort_by (|a, b| {
      if a.0 == start_index {
        return Ordering::Less
      }
      if b.0 == start_index {
        return Ordering::Greater
      }
      let angle_a = (a.1.0.y - start_point.0.y)
        .atan2 (a.1.0.x - start_point.0.x);
      let angle_b = (b.1.0.y - start_point.0.y)
        .atan2 (b.1.0.x - start_point.0.x);
      let angle_cmp = angle_a.partial_cmp (&angle_b).unwrap();
      if angle_cmp == Ordering::Equal {
        // if angles are equal sort by distance
        let dist_a = (a.1.0.x - start_point.0.x).powi (2) +
          (a.1.0.y - start_point.0.y).powi (2);
        let dist_b = (b.1.0.x - start_point.0.x).powi (2) +
          (b.1.0.y - start_point.0.y).powi (2);
        dist_a.partial_cmp (&dist_b).unwrap()
      } else {
        angle_cmp
      }
    });
    // graham scan
    let mut stack : Vec <usize> = Vec::new();
    for (point_index, point) in indexed_points {
      while stack.len() >= 2 {
        let top = stack[stack.len() - 1];
        let second = stack[stack.len() -2];
        let p1 = points[second];
        let p2 = points[top];
        let p3 = point;
        let det = Matrix2 {
          cols: vector2 (p2 - p1, p3 - p1)
        }.determinant();
        if det <= S::zero() {
          stack.pop();
        } else {
          break
        }
      }
      stack.push (point_index)
    }
    let p = points;
    let points = stack.iter().map (|i| p[*i]).collect();
    Some (Hull2 { points })
  }

  /// Points sorted in counter-clockwise order
  pub fn points (&self) -> &[Point2 <S>] {
    &self.points
  }
}

impl <S> Hull3 <S> {
  pub fn from_points (_points : &[Point3 <S>]) -> Option <Self> where S : Real {
    // code adapted from GeometricTools:
    // <https://github.com/davideberly/GeometricTools/blob/f78dd0b65bc3a0a4723586de6991dd2c339b08ad/GTE/Mathematics/ConvexHull3.h>
    unimplemented!("TODO")
  }

  pub fn points (&self) -> &[Point3 <S>] {
    &self.points
  }
}

#[cfg(test)]
mod tests {
  use super::*;
  #[test]
  fn hull2() {
    use crate::*;
    // dedup unique points
    let points : Vec <Point2 <f32>> = [
      [ 1.0,  1.0],
      [ 1.0,  1.0]
    ].map (Point2::from).into_iter().collect();
    let hull = Hull2::from_points (points.as_slice()).unwrap();
    assert_eq!(hull.points(), &[[1.0, 1.0].into()]);
    // interior point at origin is excluded
    let points : Vec <Point2 <f32>> = [
      [ 0.0,  0.0],
      [ 1.0,  3.0],
      [ 2.0, -3.0],
      [-3.0, -1.0]
    ].map (Point2::from).into_iter().collect();
    let hull = Hull2::from_points (points.as_slice()).unwrap();
    // points are in counter-clockwise order
    let points : Vec <Point2 <f32>> = [
      [ 2.0, -3.0],
      [ 1.0,  3.0],
      [-3.0, -1.0]
    ].map (Point2::from).into_iter().collect();
    assert_eq!(hull.points(), points);
    // colinear point on edge is excluded
    let points : Vec <Point2 <f32>> = [
      [ 0.0,  2.0],
      [-2.0, -2.0],
      [ 0.0, -2.0],
      [ 2.0, -2.0]
    ].map (Point2::from).into_iter().collect();
    let hull = Hull2::from_points (points.as_slice()).unwrap();
    // points are in counter-clockwise order
    let points : Vec <Point2 <f32>> = [
      [-2.0, -2.0],
      [ 2.0, -2.0],
      [ 0.0,  2.0]
    ].map (Point2::from).into_iter().collect();
    assert_eq!(hull.points(), points);
    // multiple edge and interior points are excluded
    let points : Vec <Point2 <f32>> = [
      // perimeter points
      [ 0.0,  6.0],
      [ 1.0,  5.0],
      [ 2.0,  4.0],
      [ 3.0,  3.0],
      [ 3.0,  1.0],
      [ 3.0, -1.0],
      [ 3.0, -3.0],
      [ 1.0, -3.0],
      [-1.0, -3.0],
      [-3.0, -3.0],
      [-3.0, -1.0],
      [-3.0,  1.0],
      [-3.0,  3.0],
      [-2.0,  4.0],
      [-1.0,  5.0],
      // interior points
      [-1.0,  2.0],
      [ 2.0, -1.0],
      [ 0.0,  3.0],
      [-2.0, -2.0]
    ].map (Point2::from).into_iter().collect();
    let hull = Hull2::from_points (points.as_slice()).unwrap();
    // points are in counter-clockwise order
    let points : Vec <Point2 <f32>> = [
      [-3.0, -3.0],
      [ 3.0, -3.0],
      [ 3.0,  3.0],
      [ 0.0,  6.0],
      [-3.0,  3.0]
    ].map (Point2::from).into_iter().collect();
    assert_eq!(hull.points(), points);
  }
}