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).
Point nextToTop(stack<Point> &S)
{
Point p = S.top();
S.pop();
Point res = S.top();
S.push(p);
return res;
}
int dist(Point p1, Point p2)
{
return (p1.x - p2.x)*(p1.x - p2.x) + (p1.y - p2.y)*(p1.y - p2.y);
}
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;
return (val > 0)? 1: 2;
}
int compare(const void *vp1, const void *vp2)
{
Point *p1 = (Point *)vp1;
Point *p2 = (Point *)vp2;
int o = orientation(p0, *p1, *p2);
if (o == 0)
return (dist(p0, *p2) >= dist(p0, *p1))? -1 : 1;
return (o == 2)? -1: 1;
}
void convexHull(Point points[], int n)
{
int ymin = points[0].y, min = 0;
for (int i = 1; i < n; i++)
{
int y = points[i].y;
if ((y < ymin) || (ymin == y && points[i].x < points[min].x))
ymin = points[i].y, min = i;
}
swap(points[0], points[min]);
p0 = points[0];
qsort(&points[1], n-1, sizeof(Point), compare);
stack<Point> S;
S.push(points[0]);
S.push(points[1]);
S.push(points[2]);
for (int i = 3; i < n; i++)
{
while (orientation(nextToTop(S), S.top(), points[i]) != 2)
S.pop();
S.push(points[i]);
}
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 坐标最小的 p0，第二行根据相对于p0的极角排序（利用叉积），3~5 行堆栈初始化，第 6~9行的迭代对一个点执行一次，意图是：在对点pi进行处理后，在栈S中，按由底到顶得顺序，包含CH（{p0，p1, ……，pi}）中按逆时针方向排列的各个顶点。第 7~8 行由 while 循环把发现不是凸包中的顶点从堆栈中移去。沿逆时针方向通过凸包时，在每个点处应该左转。因此while循环每次发现在一个顶点处没有向左转时，就把该顶点从堆栈中弹出。图示：
如图 a)→b) 过程中，原本p2在栈S中，到 b) 时，根据p1p3，p1p2的叉积计算，得知 p1p3 在 p1p2 的顺时针方向行，所以 p2 不在CH（Q）中，所以 p2 出栈，p3入栈，整个算法都循环这个过程，后续完整图示：
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) {
int N = pts.length;
Point2D[] points = new Point2D[N];
for (int i = 0; i < N; i++)
points[i] = pts[i];
Arrays.sort(points);
Arrays.sort(points, 1, N, points[0].polarOrder());
hull.push(points[0]);
int k1;
for (k1 = 1; k1 < N; k1++)
if (!points[0].equals(points[k1])) break;
if (k1 == N) return;
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]);
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();
}
public Iterable<Point2D> hull() {
Stack<Point2D> s = new Stack<Point2D>();
for (Point2D p : hull) s.push(p);
return s;
}
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;
}
public final class Point2D implements Comparable<Point2D> {
public static final Comparator<Point2D> X_ORDER = new XOrder();
public static final Comparator<Point2D> Y_ORDER = new YOrder();
public static final Comparator<Point2D> R_ORDER = new ROrder();
private final double x;
private final double y;
public double r() {
return Math.sqrt(x*x + y*y);
}
public double theta() {
return Math.atan2(y, x);
}
private double angleTo(Point2D that) {
double dx = that.x - this.x;
double dy = that.y - this.y;
return Math.atan2(dy, dx);
}
@param
@param
@param
@return
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;
}
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);
}
public double distanceTo(Point2D that) {
double dx = this.x - that.x;
double dy = this.y - that.y;
return Math.sqrt(dx*dx + dy*dy);
}
public double distanceSquaredTo(Point2D that) {
double dx = this.x - that.x;
double dy = this.y - that.y;
return dx*dx + dy*dy;
}
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;
}
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;
}
}
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;
}
}
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;
else if (dy2 >= 0 && dy1 < 0) return +1;
else if (dy1 == 0 && dy2 == 0) {
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);
}
}
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.
-
- int oIndex = 0;
- for (int i = 1; i < points.size(); i++) {
- if (points.get(i).y < points.get(oIndex).y)
- oIndex = i;
- }
-
- final Point origin = points.get(oIndex);
- points.remove(oIndex);
-
-
- Stack<Point> hullStack = new Stack<Point>();
- hullStack.push(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;
- }
- });
-
- points.add(origin);
- for (Point p : points) {
- Point pre = hullStack.pop();
- while(!hullStack.isEmpty() && !isLeftTurn(hullStack.peek(), pre, p)) {
- pre = hullStack.pop();
- }
- hullStack.push(pre);
- hullStack.push(p);
- }
-
- hullStack.pop();
- while(!hullStack.isEmpty()) {
- result.add(hullStack.pop());
- }
- private static float getPolar(Point start, Point end) {
-
- 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.
Related:
Find Simple Closed Path for a given set of points
Read full article from
Convex Hull | Set 2 (Graham Scan) - GeeksforGeeks
Even you only have trees in a line, you need to use rope to enclose them.
