1
\$\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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.