Skip to main content
LESSON

Convex hull

The essence is in order and maintaining that nature. The first step in online algorithm codes is often to find corner points and their adjacent points, because corner points must be concentrated in convex hull points. But in fact, it is superfluous. Due to the need for order, we first arrange all the points according to coordinates. The first element after the arrangement must also be the corner point.

From the properties of the convex hull point set, we can know the idea of ​​the convex hull algorithm, which constructs the convex hull point set in order and maintains it until the point traversal is completed.

The essence is in order and maintaining that nature. The first step in online algorithm codes is often to find corner points and their adjacent points, because corner points must be concentrated in convex hull points. But in fact, it is superfluous. Due to the need for order, we first arrange all the points according to coordinates. The first element after the arrangement must also be the corner point.

First, we add the first two sorted points to the convex hull point set, then we take out the last two points from the convex hull point set, and process the next points with the edge formed by these two points. Each time the process is processed, the vector is used to determine whether it is on the inside of the edge. If it is on the inside, the point is skipped and deemed to have been processed. If it is on the outside, it means that the edge we are using is not the edge required by the convex hull, so the last point of the current convex hull point set is discarded, and then the last two points are taken and iteratively judged until the point is added. After traversing, a convex hull point set is generated. (Fog) It should be a closed convex hull after traversing in forward order and traversing in reverse order once and returning to the starting point. Otherwise, it will only be half

There are three main algorithms for convex hull type questions:JarvisMarch algorithm,Graham algorithm and Andrew Algorithms, these three algorithms increase in time performance.

1. JarvisMarch algorithm

1.1 Thoughts

The point with the smallest ordinate and then the smallest abscissa must be a point on the convex hull. We mark it as {p_0}, from {p_0} Start by looking for points on the convex hull one by one in a counterclockwise direction, and find a point every step forward, so it is called the step method.

  • How to select the next point:
    Assuming it has been found {p_0}{p_1}, then use {p_0p_1} The point with the smallest angle between vectors is {p_2}。({p_1} Then use {p_0p_1}angle between the vector and the horizontal line)

1.2 Code

/****************************************************************** Jarvis March’s stepping algorithm
Algorithm complexity: O(nH). (where n is the total number of points and H is the number of points on the convex hull)
******************************************************************/ //#include <bits/stdc++.h>
#include <queue>
#include <cstdio>
#include <cmath>
#include <algorithm> using namespace std; const int MAXN = 10005;
const double MAXD = 1e9;
const double ACCUR = 1e-9; struct node
{ double x, y; //The coordinates of the point bool operator < (const node n) const { if (abs(y-n.y) < ACCUR) return x < n.x; else return y < n.y; } bool operator == (const node n) const { return (abs(x-n.x) < ACCUR) && (abs(y-n.y) < ACCUR); } void operator = (const node n) { x = n.x; y = n.y; }
}; struct vect
{ double x, y; void operator = (const vect v) { x = v.x; y = v.y; } double operator *(const vect v) const { return x*v.x + y*v.y; }
}; bool equal (const double d1, const double d2)
{ return abs(d1-d2) < ACCUR;
}
vect vform(const node n1, const node n2)
{ vect tmpv; tmpv.x = n2.x - n1.x; tmpv.y = n2.y - n1.y; return tmpv;
}
double vlen(const vect v)
{ return sqrt(v.x*v.x+v.y*v.y);
}
double vcos(const vect v1, const vect v2)
{ return (v1*v2)/(vlen(v1)*vlen(v2));
}
double area (const node n1, const node n2, const node n3)
{ double b1, b2, b3; b1 = vlen(vform(n1,n2)); b2 = vlen(vform(n2,n3)); b3 = vlen(vform(n3,n1)); double b = (b1+b2+b3)/2; return sqrt(b*(b-b1)*(b-b2)*(b-b3));
} node p[MAXN]; //point set
queue <node> bq; //convex hull vertex set int main()
{ int n; while(scanf("%d", &n) == 1) { if(n == 0) { break; } /*[Note] The first point is not marked first, as the loop end condition (that is, the first point is finally found)*/ int f[MAXN] = {0}; //Point set mark array; vect v; //v represents the vector formed by the previous two points. node p0, p1; //p0 represents the first point, p1 represents the previous point. p0.x = p0.y = MAXD; //Initialization for (int i = 0; i < n; ++i) { scanf("%lf%lf", &(p[i].x), &(p[i].y)); if (p[i] < p0) { p0 = p[i]; } } p1 = p0; //Initialize the previous point //[Note] The selection of the initialization vector is related to the selection of the first point. //If the first point has the smallest abscissa and then the smallest ordinate, the initial vector is the vertical unit vector. v.x = 1; v.y = 0; //The initial vector is the horizontal unit vector. do { node p2; //Point to be determined. vect v1; //Vector to be determined int j; //The subscript of the point with judgment double minvcos = -1, minvlen = MAXD; //Initial maximum angle and minimum vector length. for (int i = 0; i < n; ++i) { if (!f[i]) //Judge whether the point is already on the convex hull { vect tmpv; tmpv.x = p[i].x-p1.x; tmpv.y = p[i].y-p1.y; if (vcos(v,tmpv) > minvcos) { p2 = p[i]; v1 = tmpv; j = i; minvcos = vcos(v,tmpv); minvlen = vlen(tmpv); } else if (equal(vcos(v,tmpv),minvcos) && vlen(tmpv) < minvlen) { p2 = p[i]; v1 = tmpv; j = i; minvcos = vcos(v,tmpv); minvlen = vlen(tmpv); } } } bq.push(p2); p1 = p2; v = v1; f[j] = 1; //printf("minvcos=%f,minvlen=%f\n", minvcos, minvlen); }while(!(p1==p0)); /* while(!bq.empty()) { node tmpp = bq.front(); printf("(%f,%f)\n", tmpp.x, tmpp.y); bq.pop(); } */ //convex hull perimeter double ans = 0; node fp, ep; fp = p0; while(!bq.empty()) { ep = bq.front(); bq.pop(); ans += vlen(vform(fp, ep)); fp = ep; } printf("%.2f\n", ans); /* //convex hull area double ans = 0; node fp = bq.front(); bq.pop(); node np = bq.front(); bq.pop(); while(!bq.empty()) { node ep = bq.front(); bq.pop(); ans += area(fp,np,ep); np = ep; //printf("(%f,%f)\n", tmpp.x, tmpp.y); } printf("%d\n", (int)ans/50); */ } return 0;
}

2. Graham’s algorithm

2.1 Thoughts

Put all the points in the two-dimensional coordinate system, then the point with the smallest ordinate must be the point on the convex hull, recorded as {p_0} . Calculate the relative values of each point {p_0} Argument of , sort the points in order from small to large. (When the angles are the same, the distance {p_0} The closer ones are ranked first), then the point with the smallest argument angle and the largest point must be on the convex hull. Take the point with the smallest argument angle and record it as {p_1}, will {p_0}{p_1} Push onto the stack. Connect the point at the top of the stack and the point at the top of the next stack to get a straight line l, to see whether the current point is on the right or left of the straight line. If it is on the right, the top element of the stack is not a point on the convex hull, pop it up, and return to continue execution. If it is on the left, the current point is a point on the convex hull. Until the point with the largest argument angle.

  • cross product principle
    cross product of two vectors {P_1 imes P_2 = x_1y_2 - x_2y_1}, where the positive or negative of the result represents the direction of the cross product result. The essence of this formula is two three-dimensional vectors (z axis component is 0) the result of the cross product (the original result is{(x_1y_2 - x_2y_1)dot ec{k}}, among which {ec{k}} Yes z axis unit positive vector).

2.2 Code

/************************************************************************ Graham Scan algorithm
Time complexity: O(nlogn). The Scan process is O(n), and the preprocessing sorting is O(nlogn).
Preprocessing sorting: polar angle sorting.
************************************************************************/ //#include <bits/stdc++.h>
#include <stack>
#include <cstdio>
#include <cmath>
#include <algorithm> using namespace std; const int MAXN = 10005;
const double MAXD = 1e9;
const double ACCUR = 1e-9; struct node
{ double x, y; //The coordinates of the point bool operator < (const node n) const { if (abs(y-n.y) < ACCUR) return x < n.x; else return y < n.y; } bool operator == (const node n) const { return (abs(x-n.x) < ACCUR) && (abs(y-n.y) < ACCUR); } void operator = (const node n) { x = n.x; y = n.y; }
}; struct vect
{ double x, y; void operator = (const vect v) { x = v.x; y = v.y; } //cross product double operator *(const vect v) const { return x*v.y - y*v.x; }
}; node p0; //The point with the smallest ordinate bool equal (const double d1, const double d2)
{ return abs(d1-d2) < ACCUR;
}
vect vform(const node n1, const node n2)
{ vect tmpv; tmpv.x = n2.x - n1.x; tmpv.y = n2.y - n1.y; return tmpv;
}
double vlen(const vect v)
{ return sqrt(v.x*v.x+v.y*v.y);
}
double vcos(const vect v1, const vect v2)
{ return (v1*v2)/(vlen(v1)*vlen(v2));
}
//Polar angle sorting
bool cmpp (const node p1, const node p2)
{ vect v1, v2; v1 = vform(p0, p1); v2 = vform(p0, p2); if (equal(v1*v2,0)) { return vlen(v1) < vlen(v2); } else { return v1*v2 > 0; }
}
//The cross product judgment point (the end point of v2) is to the left or right of v1
bool cmpv (const vect v1, const vect v2)
{ return (v1*v2 > 0) || equal(v1*v2,0);
}
double area (const node n1, const node n2, const node n3)
{ /* //Helen's formula double b1, b2, b3; b1 = vlen(vform(n1,n2)); b2 = vlen(vform(n2,n3)); b3 = vlen(vform(n3,n1)); double b = (b1+b2+b3)/2; return sqrt(b*(b-b1)*(b-b2)*(b-b3)); */ //Cross product formula (cross product is the area of a parallelogram) vect v1, v2; v1 = vform(n1, n2); v2 = vform(n1, n3); return abs(v1*v2)/2;
} node p[MAXN]; //point set
stack <node> bs; //convex hull vertex set int main()
{ int n; p0.x = p0.y = MAXD; //Initialize the first point scanf("%d", &n); for(int i = 0; i < n; ++i) { scanf("%lf%lf", &(p[i].x), &(p[i].y)); if(p[i] < p0) { p0 = p[i]; } } sort(p,p+n,cmpp); bs.push(p[0]); bs.push(p[1]); int j = 2; while(j < n) { //Remove the top of the stack and the top of the sub-stack node p1, p2; p2 = bs.top(); bs.pop(); p1 = bs.top(); //Construct cross product vector vect v1, v2; v1 = vform(p1,p2); v2 = vform(p2,p[j]); if(cmpv(v1,v2)) { bs.push(p2); bs.push(p[j]); ++j; } } /* while(!bs.empty()) { node tmpp = bs.top(); printf("(%f,%f)\n", tmpp.x, tmpp.y); bs.pop(); } */ /* //Calculate perimeter double ans = 0; node fp, ep; fp = p[0]; while(!bs.empty()) { ep = bs.top(); bs.pop(); ans += vlen(vform(fp, ep)); fp = ep; } printf("%.2f\n", ans); */ //Calculate area double ans = 0; node fp, np, ep; fp = bs.top(); bs.pop(); np = bs.top(); bs.pop(); while(!bs.empty()) { ep = bs.top(); bs.pop(); ans += area(fp,np,ep); np = ep; } printf("%d\n", (int)ans/50); return 0;
}

3. Andrew’s algorithm

3.1 Thoughts

The preprocessing sorting is changed to horizontal sorting, sorting according to the abscissa from small to large, and if the abscissas are the same, the ordinates are sorted from small to large. follow graham Algorithmic thinking starts from {p_0}{p_1} Scan all points to get the lower convex hull, and then {p_{n-1}}{p_{n-2}} Scan all points to get the upper convex hull, and the combination of the two is the entire convex hull. (Note: here {p_1} not necessarily in the convex hull)

3.2 Code

/************************************************************************ Andrew algorithm (Graham Scan algorithm variant)
Time complexity: O(nlogn). The Scan process is O(n), and the preprocessing sorting is O(nlogn).
Preprocessing sorting: horizontal sorting sorting.
************************************************************************/ //#include <bits/stdc++.h>
#include <stack>
#include <cstdio>
#include <cmath>
#include <algorithm> using namespace std; const int MAXN = 10005;
const double MAXD = 1e9;
const double ACCUR = 1e-9; struct node
{ double x, y; //The coordinates of the point //Horizontal sorting (different from polar angular sorting, it can only determine that p0 and pn-1 are within the convex hull) bool operator < (const node n) const { if (abs(x-n.x) < ACCUR) return y < n.y; else return x < n.x; } bool operator == (const node n) const { return (abs(x-n.x) < ACCUR) && (abs(y-n.y) < ACCUR); } void operator = (const node n) { x = n.x; y = n.y; }
}; struct vect
{ double x, y; void operator = (const vect v) { x = v.x; y = v.y; } //cross product double operator *(const vect v) const { return x*v.y - y*v.x; }
}; bool equal (const double d1, const double d2)
{ return abs(d1-d2) < ACCUR;
}
vect vform(const node n1, const node n2)
{ vect tmpv; tmpv.x = n2.x - n1.x; tmpv.y = n2.y - n1.y; return tmpv;
}
//Calculate vector length
double vlen(const vect v)
{ return sqrt(v.x*v.x+v.y*v.y);
}
//Calculate the cosine of the angle between vectors
double vcos(const vect v1, const vect v2)
{ return (v1*v2)/(vlen(v1)*vlen(v2));
}
/*
//Polar angle sorting
bool cmpp (const node p1, const node p2)
{ vect v1, v2; v1 = vform(p0, p1); v2 = vform(p0, p2); if (equal(v1*v2,0)) { return vlen(v1) < vlen(v2); } else { return v1*v2 > 0; }
}
*/
//The cross product judgment point (the end point of v2) is to the left or right of v1
bool cmpv (const vect v1, const vect v2)
{ return (v1*v2 > 0) || equal(v1*v2,0);
}
double area (const node n1, const node n2, const node n3)
{ /* //Helen's formula double b1, b2, b3; b1 = vlen(vform(n1,n2)); b2 = vlen(vform(n2,n3)); b3 = vlen(vform(n3,n1)); double b = (b1+b2+b3)/2; return sqrt(b*(b-b1)*(b-b2)*(b-b3)); */ //Cross product formula (cross product is the area of a parallelogram) vect v1, v2; v1 = vform(n1, n2); v2 = vform(n1, n3); return abs(v1*v2)/2;
} node p[MAXN]; //point set
stack <node> bs; //convex hull vertex set int main()
{ int n; scanf("%d", &n); for(int i = 0; i < n; ++i) { scanf("%lf%lf", &(p[i].x), &(p[i].y)); } sort(p,p+n); //Forward scan (upper convex hull) bs.push(p[0]); bs.push(p[1]); int j = 2; while(j < n) { //Remove the top of the stack and the top of the sub-stack node p1, p2; p2 = bs.top(); bs.pop(); //p1 is not necessarily in the convex hull if(bs.empty()) { bs.push(p2); bs.push(p[j]); ++j; } else { p1 = bs.top(); //Construct cross product vector vect v1, v2; v1 = vform(p1,p2); v2 = vform(p2,p[j]); if(cmpv(v1,v2)) { bs.push(p2); bs.push(p[j]); ++j; } } } //Reverse scan (lower convex hull) int k = n-2; while(k >= 0) { //Remove the top of the stack and the top of the sub-stack node p1, p2; p2 = bs.top(); bs.pop(); p1 = bs.top(); //Construct cross product vector vect v1, v2; v1 = vform(p1,p2); v2 = vform(p2,p[k]); if(cmpv(v1,v2)) { bs.push(p2); bs.push(p[k]); --k; } } bs.pop(); //p0 is pushed onto the stack repeatedly /* while(!bs.empty()) { node tmpp = bs.top(); printf("(%f,%f)\n", tmpp.x, tmpp.y); bs.pop(); } */ /* //Calculate perimeter double ans = 0; node fp, ep; fp = p[0]; while(!bs.empty()) { ep = bs.top(); bs.pop(); ans += vlen(vform(fp, ep)); fp = ep; } printf("%.2f\n", ans); */ //Calculate area double ans = 0; node fp, np, ep; fp = bs.top(); bs.pop(); np = bs.top(); bs.pop(); while(!bs.empty()) { ep = bs.top(); bs.pop(); ans += area(fp,np,ep); np = ep; } printf("%d\n", (int)ans/50); return 0;
}