Massive Algorithms: Convex Hull | Set 2 (Graham Scan)

Convex Hull | Set 2 (Graham Scan) - GeeksforGeeks

Convex Hull | Set 2 (Graham Scan) - GeeksforGeeks
Given a set of points in the plane. the convex hull of the set is the smallest convex polygon that contains all the points of it.
We have discussed Jarvis’s Algorithm for Convex Hull. Worst case time complexity of Jarvis’s Algorithm is O(n^2). Using Graham’s scan algorithm, we can find Convex Hull in O(nLogn) time.

Let points[0..n-1] be the input array.
1) Find the bottom-most point by comparing y coordinate of all points. If there are two points with same y value, then the point with smaller x coordinate value is considered. Put the bottom-most point at first position.
2) Consider the remaining n-1 points and sort them by polor angle in counterclockwise order around points[0]. If polor angle of two points is same, then put the nearest point first.
3) Create an empty stack ‘S’ and push points[0], points[1] and points[2] to S.
4) Process remaining n-3 points one by one. Do following for every point ‘points[i]‘
        4.1) Keep removing points from stack while orientation of following 3 points is not counterclockwise or they don’t make a left turn.
            a) Point next to top in stack
            b) Point at the top of stack
            c) points[i]
         4.2) Push points[i] to S
5) Print contents of S

Phase 1 (Sort points): We first find the bottom-most point. The idea is to pre-process points be sorting them with respect to the bottom-most point. Once the points are sorted, they form a simple closed path.

What should be the sorting criteria? computation of actual angles would be inefficient since trigonometric functions are not simple to evaluate. The idea is to use the orientation to compare angles without actually computing them

Phase 2 (Accept or Reject Points): Once we have the closed path, the next step is to traverse the path and remove concave points on this path. How to decide which point to remove and which to keep? Again, orientation helps here. The first two points in sorted array are always part of Convex Hull. For remaining points, we keep track of recent three points, and find the angle formed by them. Let the three points be prev(p), curr(c) and next(n). If orientation of these points (considering them in same order) is not counterclockwise, we discard c, otherwise we keep it. Following diagram shows step by step process of this phase (Source of these diagrams is Ref 2).

// A utility function to find next to top in a stack

Point nextToTop(stack<Point> &S)

{

    Point p = S.top();

    S.pop();

    Point res = S.top();

    S.push(p);

    return res;

}

// A utility function to return square of distance between p1 and p2

int dist(Point p1, Point p2)

{

    return (p1.x - p2.x)*(p1.x - p2.x) + (p1.y - p2.y)*(p1.y - p2.y);

}

// To find orientation of ordered triplet (p, q, r).

// The function returns following values

// 0 --> p, q and r are colinear

// 1 --> Clockwise

// 2 --> Counterclockwise

int orientation(Point p, Point q, Point r)

{

    int val = (q.y - p.y) * (r.x - q.x) -

              (q.x - p.x) * (r.y - q.y);

    if (val == 0) return 0;  // colinear

    return (val > 0)? 1: 2; // clock or counterclock wise

}

// A function used by library function qsort() to sort an array of

// points with respect to the first point

int compare(const void *vp1, const void *vp2)

{

   Point *p1 = (Point *)vp1;

   Point *p2 = (Point *)vp2;

   // Find orientation

   int o = orientation(p0, *p1, *p2);

   if (o == 0)

     return (dist(p0, *p2) >= dist(p0, *p1))? -1 : 1;

   return (o == 2)? -1: 1;

}

// Prints convex hull of a set of n points.

void convexHull(Point points[], int n)

{

   // Find the bottommost point

   int ymin = points[0].y, min = 0;

   for (int i = 1; i < n; i++)

   {

     int y = points[i].y;

     // Pick the bottom-most or chose the left most point in case of tie

     if ((y < ymin) || (ymin == y && points[i].x < points[min].x))

        ymin = points[i].y, min = i;

   }

   // Place the bottom-most point at first position

   swap(points[0], points[min]);

   // Sort n-1 points with respect to the first point.  A point p1 comes

   // before p2 in sorted ouput if p2 has larger polar angle (in 

   // counterclockwise direction) than p1

   p0 = points[0];

   qsort(&points[1], n-1, sizeof(Point), compare);

   // Create an empty stack and push first three points to it.

   stack<Point> S;

   S.push(points[0]);

   S.push(points[1]);

   S.push(points[2]);

   // Process remaining n-3 points

   for (int i = 3; i < n; i++)

   {

      // Keep removing top while the angle formed by points next-to-top, 

      // top, and points[i] makes a non-left turn

      while (orientation(nextToTop(S), S.top(), points[i]) != 2)

         S.pop();

      S.push(points[i]);

   }

   // Now stack has the output points, print contents of stack

   while (!S.empty())

   {

       Point p = S.top();

       cout << "(" << p.x << ", " << p.y <<")" << endl;

       S.pop();

   }

}

计算几何学(Computational Geometry) - 算法导论第33章

点集Q的凸包是一个最小的凸多边形P，可以把Q中的每个点都想象成是露在一块板外的铁钉，那么凸包就是包围了所有这些铁钉的一条拉紧了的橡皮绳所构成的形状。如图点集Q和凸包CH（Q）。

寻找凸包，可以使用Graham扫描法（Graham’s scan）。借助一个堆栈S，输入集合Q中的每个点都被压入栈一次，非 CH(Q) 中顶点的点最终将被弹出堆栈，算法结束时，堆栈S中仅包含CH（Q）中的顶点，其顺序为各点在边界上出现的逆时针方向排列的顺序。

伪代码：

GRAHAM-SCAN(Q)

1 let p0 be the point in Q with the minimum y-coordinate,

or the leftmost such point in case of a tie

2 let 〈p1, p2, …, pm〉 be the remaining points in Q,

sorted by polar angle in counterclockwise order around p0

(if more than one point has the same angle, remove all but

the one that is farthest from p0)

3 PUSH(p0, S)

4 PUSH(p1, S)

5 PUSH(p2, S)

6 for i ← 3 to m

7 do while the angle formed by points NEXT-TO-TOP(S), TOP(S),

and pi makes a nonleft turn

8 do POP(S)

9 PUSH(pi, S)

10 return S

第一行找到 y 坐标最小的 p₀，第二行根据相对于p0的极角排序（利用叉积），3~5 行堆栈初始化，第 6~9行的迭代对一个点执行一次，意图是：在对点pi进行处理后，在栈S中，按由底到顶得顺序，包含CH（{p₀，p₁, ……，p_i}）中按逆时针方向排列的各个顶点。第 7~8 行由 while 循环把发现不是凸包中的顶点从堆栈中移去。沿逆时针方向通过凸包时，在每个点处应该左转。因此while循环每次发现在一个顶点处没有向左转时，就把该顶点从堆栈中弹出。图示：

如图 a)→b) 过程中，原本p₂在栈S中，到 b) 时，根据p₁p₃，p₁p₂的叉积计算，得知 p₁p₃ 在 p₁p₂ 的顺时针方向行，所以 p₂不在CH（Q）中，所以 p₂ 出栈，p₃入栈，整个算法都循环这个过程，后续完整图示：

typedef struct {

    double x;

    double y;

}Point;

Point p[110], stack[110];

int N, top;

double Direction(Point p1, Point p2, Point p3) {

    return (p2.x - p1.x) * (p3.y - p1.y) - (p2.y - p1.y) * (p3.x - p1.x);

}

double Dis(Point a, Point b) {

    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));

}

int cmp(const void *a, const void *b) {

    Point c = *(Point *)a;

    Point d = *(Point *)b;

    double k = Direction(p[0], c, d);

    if(k < 0 || (!k && Dis(c, p[0]) > Dis(d, p[0]))) {

        return 1;

}

    return -1;

}

void Convex() {

    for(int i = 1; i < N; i++) {

        Point temp;

        if(p[i].y < p[0].y || (p[i].y == p[0].y && p[i].x < p[0].x)) {

            temp = p[i];

            p[i] = p[0];

            p[0] = temp;

}

}

    qsort(p + 1, N - 1, sizeof(p[0]), cmp);

    stack[0] = p[0];        

    stack[1] = p[1];

    stack[2] = p[2];

    top =2;

    for(int i = 3; i < N; i++) {

        while(top >= 2 && Direction(stack[top - 1], stack[top], p[i]) <= 0) {

            top--;

}

        top++;

        stack[top] = p[i];

}

}

Java Implementaion
http://algs4.cs.princeton.edu/99hull/GrahamScan.java.html

    private Stack<Point2D> hull = new Stack<Point2D>();

    public GrahamScan(Point2D[] pts) {

        // defensive copy
        int N = pts.length;
        Point2D[] points = new Point2D[N];
        for (int i = 0; i < N; i++)
            points[i] = pts[i];

        // preprocess so that points[0] has lowest y-coordinate; break ties by x-coordinate
        // points[0] is an extreme point of the convex hull
        // (alternatively, could do easily in linear time)
        Arrays.sort(points);

        // sort by polar angle with respect to base point points[0],
        // breaking ties by distance to points[0]
        Arrays.sort(points, 1, N, points[0].polarOrder());

        hull.push(points[0]);       // p[0] is first extreme point

        // find index k1 of first point not equal to points[0]
        int k1;
        for (k1 = 1; k1 < N; k1++)
            if (!points[0].equals(points[k1])) break;
        if (k1 == N) return;        // all points equal

        // find index k2 of first point not collinear with points[0] and points[k1]
        int k2;
        for (k2 = k1 + 1; k2 < N; k2++)
            if (Point2D.ccw(points[0], points[k1], points[k2]) != 0) break;
        hull.push(points[k2-1]);    // points[k2-1] is second extreme point

        // Graham scan; note that points[N-1] is extreme point different from points[0]
        for (int i = k2; i < N; i++) {
            Point2D top = hull.pop();
            while (Point2D.ccw(hull.peek(), top, points[i]) <= 0) {
                top = hull.pop();
            }
            hull.push(top);
            hull.push(points[i]);
        }

        assert isConvex();
    }

    // return extreme points on convex hull in counterclockwise order as an Iterable
    public Iterable<Point2D> hull() {
        Stack<Point2D> s = new Stack<Point2D>();
        for (Point2D p : hull) s.push(p);
        return s;
    }

    // check that boundary of hull is strictly convex
    private boolean isConvex() {
        int N = hull.size();
        if (N <= 2) return true;

        Point2D[] points = new Point2D[N];
        int n = 0;
        for (Point2D p : hull()) {
            points[n++] = p;
        }

        for (int i = 0; i < N; i++) {
            if (Point2D.ccw(points[i], points[(i+1) % N], points[(i+2) % N]) <= 0) {
                return false;
            }
        }
        return true;
    }

http://algs4.cs.princeton.edu/99hull/Point2D.java.html

public final class Point2D implements Comparable<Point2D> {

    /**
     * Compares two points by x-coordinate.
     */
    public static final Comparator<Point2D> X_ORDER = new XOrder();

    /**
     * Compares two points by y-coordinate.
     */
    public static final Comparator<Point2D> Y_ORDER = new YOrder();

    /**
     * Compares two points by polar radius.
     */
    public static final Comparator<Point2D> R_ORDER = new ROrder();

    private final double x;    // x coordinate
    private final double y;    // y coordinate

     * Returns the polar radius of this point.
    public double r() {
        return Math.sqrt(x*x + y*y);
    }
     * Returns the angle of this point in polar coordinates.
    public double theta() {
        return Math.atan2(y, x);
    }

     * Returns the angle between this point and that point.
    private double angleTo(Point2D that) {
        double dx = that.x - this.x;
        double dy = that.y - this.y;
        return Math.atan2(dy, dx);
    }

     * Is a->b->c a counterclockwise turn?

     * @return { -1, 0, +1 } if a->b->c is a { clockwise, collinear; counterclocwise } turn.
    public static int ccw(Point2D a, Point2D b, Point2D c) {
        double area2 = (b.x-a.x)*(c.y-a.y) - (b.y-a.y)*(c.x-a.x);
        if      (area2 < 0) return -1;
        else if (area2 > 0) return +1;
        else                return  0;
    }

     * Is a->b->c a counterclockwise turn?
     * @param a first point
     * @param b second point
     * @param c third point
     * @return { -1, 0, +1 } if a->b->c is a { clockwise, collinear; counterclocwise } turn.
     */
    public static int ccw(Point2D a, Point2D b, Point2D c) {
        double area2 = (b.x-a.x)*(c.y-a.y) - (b.y-a.y)*(c.x-a.x);
        if      (area2 < 0) return -1;
        else if (area2 > 0) return +1;
        else                return  0;
    }
     * Returns twice the signed area of the triangle a-b-c.
    public static double area2(Point2D a, Point2D b, Point2D c) {
        return (b.x-a.x)*(c.y-a.y) - (b.y-a.y)*(c.x-a.x);
    }
     * Returns the Euclidean distance between this point and that point.
    public double distanceTo(Point2D that) {
        double dx = this.x - that.x;
        double dy = this.y - that.y;
        return Math.sqrt(dx*dx + dy*dy);
    }
     * Returns the square of the Euclidean distance between this point and that point.
    public double distanceSquaredTo(Point2D that) {
        double dx = this.x - that.x;
        double dy = this.y - that.y;
        return dx*dx + dy*dy;
    }
     * Compares this point to that point by y-coordinate, breaking ties by x-coordinate.
    public int compareTo(Point2D that) {
        if (this.y < that.y) return -1;
        if (this.y > that.y) return +1;
        if (this.x < that.x) return -1;
        if (this.x > that.x) return +1;
        return 0;
    }

  // compare points according to their polar radius
    private static class ROrder implements Comparator<Point2D> {
        public int compare(Point2D p, Point2D q) {
            double delta = (p.x*p.x + p.y*p.y) - (q.x*q.x + q.y*q.y);
            if (delta < 0) return -1;
            if (delta > 0) return +1;
            return 0;
        }
    }
 
    // compare other points relative to atan2 angle (bewteen -pi/2 and pi/2) they make with this Point
    private class Atan2Order implements Comparator<Point2D> {
        public int compare(Point2D q1, Point2D q2) {
            double angle1 = angleTo(q1);
            double angle2 = angleTo(q2);
            if      (angle1 < angle2) return -1;
            else if (angle1 > angle2) return +1;
            else                      return  0;
        }
    }

    // compare other points relative to polar angle (between 0 and 2pi) they make with this Point
    private class PolarOrder implements Comparator<Point2D> {
        public int compare(Point2D q1, Point2D q2) {
            double dx1 = q1.x - x;
            double dy1 = q1.y - y;
            double dx2 = q2.x - x;
            double dy2 = q2.y - y;

            if      (dy1 >= 0 && dy2 < 0) return -1;    // q1 above; q2 below
            else if (dy2 >= 0 && dy1 < 0) return +1;    // q1 below; q2 above
            else if (dy1 == 0 && dy2 == 0) {            // 3-collinear and horizontal
                if      (dx1 >= 0 && dx2 < 0) return -1;
                else if (dx2 >= 0 && dx1 < 0) return +1;
                else                          return  0;
            }
            else return -ccw(Point2D.this, q1, q2);     // both above or below

            // Note: ccw() recomputes dx1, dy1, dx2, and dy2
        }
    }

    // compare points according to their distance to this point
    private class DistanceToOrder implements Comparator<Point2D> {
        public int compare(Point2D p, Point2D q) {
            double dist1 = distanceSquaredTo(p);
            double dist2 = distanceSquaredTo(q);
            if      (dist1 < dist2) return -1;
            else if (dist1 > dist2) return +1;
            else                    return  0;
        }
    }

Java Program to Implement Graham Scan Algorithm to Find the Convex Hull
http://www.dcs.gla.ac.uk/~pat/52233/slides/Hull1x1.pdf

http://yuanhsh.iteye.com/blog/2194148
Convex hull: the smallest polygon that covers all the points in the given set.

//find the point with the lowest y
int oIndex = 0;
for (int i = 1; i < points.size(); i++) {
if (points.get(i).y < points.get(oIndex).y)
oIndex = i;
}
//remove the origin points and put it into the stack
final Point origin = points.get(oIndex);
points.remove(oIndex);
//stack that we use to iterate through the
Stack<Point> hullStack = new Stack<Point>();
hullStack.push(origin); //put origin point in stack
//sort the points according to their polarity relative to origin
Collections.sort(points, new Comparator<Point>() {
@Override
public int compare(Point p1, Point p2) {
float polar1 = getPolar(origin, p1);
float polar2 = getPolar(origin, p2);
if (polar1 > polar2)
return -1;
else if (polar1 < polar2)
return 1;
return 0;
}
});
//add origin to as the last element
points.add(origin);
for (Point p : points) {
Point pre = hullStack.pop();
while(!hullStack.isEmpty() && !isLeftTurn(hullStack.peek(), pre, p)) {
pre = hullStack.pop();//discard the pre and take new top of the stack as pre
}
hullStack.push(pre);
hullStack.push(p);
}
//remove the duplicate origin
hullStack.pop();
while(!hullStack.isEmpty()) {
result.add(hullStack.pop());
}
private static float getPolar(Point start, Point end) {
//get the cosin between end and start
float r = ((end.x - start.x)*(end.x - start.x) + (end.y - start.y)*(end.y - start.y));
return (float) ((end.x - start.x)/Math.pow(r, 0.5));
}
private static boolean isLeftTurn(Point a, Point b, Point c) {
if ((c.y - a.y)/(c.x - a.x) > (b.y - a.y)/(b.x - a.x))
return true;
return false;
}

Two things need to pay attention:
1. the original point is the one with the lowest y coordinate. So all the other point was above it. We want the points to be sorted according to the polar angle relative to the original point, in the order of counterclockwise. So we should use cosin which on the uppper part of coordinate system, decrease in the order of counterclock direction.

2. To find out whether we make a left (counter clock wise) turn in a three point a->b->c, there is a simple way. We should check if (c.y – a.y)*(b.x – a.x) – (b.y – a.y)*(c.x – a.x). If it is equals to 0, then a, b and c are on the same line. if positive, it is making a left turn, otherwise it is making a right turn.