\$\begingroup\$

Continuing the Algorithm in a Nutshell series, here is the code:

use point::{Point, sort_points, Direction}; // see http://i.imgur.com/C2zng5r.png // I have done this from a slightly different perspective, // i.e. intead of using the lowest point as the head, I used the leftmost. fn graham_scan(points: &mut Vec<Point>) -> Vec<Point> { let mut hull: Vec<Point> = Vec::new(); sort_points(points); hull.push(points[0]); hull.push(points[1]); for i in 2..points.len() { loop { println!("{:?}", &hull); let m1 = hull.len() - 1; let m0 = m1 - 1; let direction = hull[m0].direction(&hull[m1], &points[i]); match direction { Direction::Left => { hull.push(points[i]); break; }, Direction::Ahead =>{ hull.pop(); hull.push(points[i]); break; }, _ => { hull.pop(); () } } } } return hull; } #[cfg(test)] mod test { use super::graham_scan; use point::Point; #[test] fn test_graham_scan() { let mut points: Vec<Point> = Vec::new(); // These points form a triangle, so only the 3 vertices should be in the convex hull. for i in 1..10 { points.push(Point::new(i as f64, i as f64)); points.push(Point::new(i as f64, (-i) as f64)); points.push(Point::new(i as f64, 0.0)); } points.push(Point::new(0.0, 0.0)); let hull = graham_scan(&mut points); let hull_should_be = vec![ Point::new(0.0, 0.0), Point::new(9.0, -9.0), Point::new(9.0, 9.0), ]; assert_eq!(hull, hull_should_be); } }

The point module:

use std::ops::{Add, Sub, Mul}; use std::cmp::Ordering; use super::definite_num::DefinitelyANumber; #[derive(Debug, Copy, Clone, PartialEq, Eq, PartialOrd, Ord)] pub struct Point { pub x: DefinitelyANumber, pub y: DefinitelyANumber, } impl Point { pub fn new(x: f64, y: f64) -> Point { Point { x: DefinitelyANumber::new(x).expect("X coordinate cannot be NaN!"), y: DefinitelyANumber::new(y).expect("Y coordinate cannot be NaN!"), } } // Euclidean distance pub fn distance(&self, other: &Point) -> f64 { ((self.x - other.x).to_f64().powi(2) + (self.y - other.y).to_f64().powi(2)).sqrt() } // Draw a horizontal line through this point, connect this point with the other, and measure the angle between these two lines. pub fn angle(&self, other: &Point) -> f64 { if self == other { 0.0 } else { (other.y - self.y).to_f64().atan2((other.x - self.x).to_f64()) } } pub fn magnitude(&self) -> f64 { (self.x.to_f64().powi(2) + self.y.to_f64().powi(2)).sqrt() } pub fn sin_cos(&self) -> (f64, f64) { let mag = self.magnitude(); (self.y.to_f64() / mag, self.x.to_f64() / mag) } pub fn rotate(&self, theta: f64) -> Point { let x = self.x.to_f64(); let y = self.y.to_f64(); let cosine = theta.cos(); let sine = theta.sin(); let x_cos_theta = x * cosine; let x_sin_theta = x * sine; let y_cos_theta = y * cosine; let y_sin_theta = y * sine; let x1 = x_cos_theta - y_sin_theta; let y1 = x_sin_theta + y_cos_theta; Point::new(x1, y1) } pub fn direction(&self, p1: &Point, p2: &Point) -> Direction { let v1 = *p1 - *self; let v2 = *p2 - *self; let x1 = v1.x.to_f64(); let x2 = v2.x.to_f64(); let y1 = v1.y.to_f64(); let y2 = v2.y.to_f64(); let det = x1 * y2 - y1 * x2; if det < 0.0 { Direction::Right } else if det > 0.0 { Direction::Left } else { Direction::Ahead } } } #[derive(Debug, PartialEq)] pub enum Direction { Left, Right, Ahead, } impl Add for Point { type Output = Point; fn add(self, rhs: Point) -> Point { Point { x: self.x + rhs.x, y: self.y + rhs.y, } } } impl Sub for Point { type Output = Point; fn sub(self, rhs: Point) -> Point { Point { x: self.x - rhs.x, y: self.y - rhs.y, } } } // dot product impl Mul for Point { type Output = f64; fn mul(self, rhs: Point) -> f64 { (self.x * rhs.x + self.y * rhs.y).to_f64() } } // sort by angle to head pub fn sort_points(points: &mut Vec<Point>) { // sort by coordinates so that the first point is the leftmost points.sort(); let head = points[0]; // sort by the angle with the first point // when that is equal, sort by x // when that is equal, sort by y points.sort_by(|a, b| { // head always comes first. if a == &head { return Ordering::Less; } if b == &head { return Ordering::Greater } let angle_a = head.angle(&a); let angle_b = head.angle(&b); let angle_cmp = angle_a.partial_cmp(&angle_b).unwrap(); if angle_cmp == Ordering::Equal { a.cmp(&b) } else { angle_cmp } }); } #[cfg(test)] mod test { use point::Point; use super::Direction; #[test] fn test_point() { use std::f64::consts::PI; let p1 = Point::new(0.0, 0.0); let p2 = Point::new(0.0, 1.0); assert_eq!(p1.angle(&p2), PI / 2.0); assert_eq!(p1.distance(&p2), 1.0); let p1 = Point::new(0.0, 0.0); let p2 = Point::new(1.0, 1.0); assert_eq!(p1.angle(&p2), PI / 4.0); assert_eq!(p1.distance(&p2), 2.0f64.sqrt()); let p1 = Point::new(0.0, 0.0); let p2 = Point::new(1.0, -1.0); assert_eq!(p1.angle(&p2), -PI / 4.0); assert_eq!(p1.distance(&p2), 2.0f64.sqrt()); } #[test] fn test_direction() { let p1 = Point::new(1.0, 1.0); let p2 = Point::new(2.0, 2.0); let p3 = Point::new(3.0, 3.0); assert_eq!(p1.direction(&p2, &p3), Direction::Ahead); let p1 = Point::new(1.0, 1.0); let p2 = Point::new(2.0, 2.0); let p3 = Point::new(3.0, 2.5); assert_eq!(p1.direction(&p2, &p3), Direction::Right); let p1 = Point::new(1.0, 1.0); let p2 = Point::new(2.0, 2.0); let p3 = Point::new(3.0, 3.5); assert_eq!(p1.direction(&p2, &p3), Direction::Left); let p1 = Point::new(1.0, -1.0); let p2 = Point::new(2.0, -2.0); let p3 = Point::new(3.0, -3.0); assert_eq!(p1.direction(&p2, &p3), Direction::Ahead); let p1 = Point::new(1.0, -1.0); let p2 = Point::new(2.0, -2.0); let p3 = Point::new(3.0, -2.5); assert_eq!(p1.direction(&p2, &p3), Direction::Left); let p1 = Point::new(1.0, -1.0); let p2 = Point::new(2.0, -2.0); let p3 = Point::new(3.0, -3.5); assert_eq!(p1.direction(&p2, &p3), Direction::Right); let p3 = Point::new(1.0, -1.0); let p2 = Point::new(2.0, -2.0); let p1 = Point::new(3.0, -3.5); assert_eq!(p1.direction(&p2, &p3), Direction::Left); let p3 = Point::new(1.0, -1.0); let p2 = Point::new(2.0, -2.0); let p1 = Point::new(3.0, -2.5); assert_eq!(p1.direction(&p2, &p3), Direction::Right); } }

The DefinitelyANumber trait was provided by @Shepmaster in my previous post.

All suggestions are welcome.