The problem:
Given a set of points in a 2D plane determine the smallest encapsulating convex polygon which includes all the points. Conceptually you can think of this as wrapping a rubber band on the set of points and the line segments which make this rubber band is the convex hull. See the figure below:
Applications:
Normally the convex hull is used as a pre-processing step for many computational geometry problems. For example to determine the minimum area convex polygon required to include a given set of points, we first proceed to calculate the convex hull and use the hull as input to a polygon area computation algorithm. We'll see an application of this later.
Algorithm:
There are many algorithms which solve the convex hull problem. Perhaps the simplest is to do a naive brute force. The only way a point belongs to the convex hull is if the point isn't contained in any triangle based off another set of points. Therefore for each point in the set we can iterate all possible 3 points and check, this takes linear time for the set and O(n^3) to check for all possible triangles - hence this approach takes O(n^4).
What we will look at is an efficient yet simple approach which takes O(n log n) time. The algorithm is called Graham's Scan and it makes use of a fundamental concept in Computational Geometry: cross products. We can use cross products to determine whether a point is clockwise or counter-clockwise relative to the origin (0,0). Thus using this we can transform all our points to the left and bottom-most point by doing vector subtraction.
The basis of Graham's scan is to consider one point at a time using a sorted ordering of the points. The sorted criteria is the polar angle from our fixed point P (the left/bottom-most point in the set) to each of the other points in the set. Then we iterate through each of the points in this set and see if the next point makes a counter-clockwise or clockwise turn. If it's clockwise, it means our last point was non-optimal and we backtrack back until the next point we add is counter-clockwise.
A visualisation is shown below (from Wikipedia):
Note that in Step 1, we have our convex hull as P-A-B. We proceed to try point C, as it results in a counter-clockwise turn (otherwise known as a left turn from the relative position of line segment AB) - we add this to our hull. It now becomes P-A-B-C and we proceed to try point D. From the perspective of line segment BC, point D actually is clockwise (we have to turn right to reach it - you can think of it as an ant traversing through the hull and the direction we need to turn in order to get to the next point). Therefore, we pop point C off our hull and re-compute whether or not point D is counter-clockwise to line segment AB which it turns out to be. Hence we add point D to the hull - after the removal of point C we are left with a hull of P-A-B-D (the image in Step 3). One more step is to link D back to the original point P - implementation wise, we usually append the starting point to act as a sentinel. Hence our hull would be P-A-B-D-P.
To calculate the cross-product in the 2D space:
which results in:
Given two line segments P1P2 and P1P3 we can calculate the coefficients of the above expression:
P1P2 = (x2 - x1)i + (y2 - y1)j => a1 = x2 - x1 and b1 = y2 - y1
P1P3 = (x3 - x1)i + (y3 - y1)j => a2 = x3 - x1 and b2 = y3 - y1
As we are only dealing with the 2D plane, a3 = 0 and b3 = 0. Given the above expression we are left with:
P = ((x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1)) k
Using the sign of the expression P, we can determine if the points are collinear (P evaluates to 0), counter-clockwise (P evaluates to a positive number) or clockwise (P evaluates to a negative number). Note that the actual value doesn't matter we are only interested in the sign - this way we can avoid floating point precision problems which is a big bonus in geometry problems.
If you are curious to why the sign of the expression P dictates the counter-clockwise/clockwise nature of the two vectors, simply look at the definition of the cross product:
As a x b is the thing we just calculated and ab is the magnitude of the vectors a and b respectively. Hence ab will always be positive. Therefore if P is negative, then theta will also be negative and hence will be clockwise. Otherwise if P is positive, then theta will also be positive and hence will be counter-clockwise. Note that the absolute value of theta is between 0 and 180 degrees inclusive.
Now back to the actual Graham Scan algorithm: we need to firstly sort the points in accordance to our fixed P (which will always be in the convex hull as it's the lowest and left-most point). We simply override the comparison operator for the Point class which will automatically sort based on angle without actually calculating. To see this, note that the relationship between the counter-clockwise/clockwise is transitive. As in, if point C is CCW to point B and point B is CCW to point A - this also implies that point C is CCW to point A. Hence this is a sufficient sorting criterion.
We process points in this order and calculate whether or not the point is clockwise or counter-clockwise compared to the last line segment in our convex hull. We stop when we reach our starting point (the sentinel) in which case we end up with the hull.
Implementation of the convex hull is shown below:
#define EPS 1e-09
class Point {
public:
double x;
double y;
Point() { }
Point(double a, double b) : x(a), y(b) { }
};
Point first_point;
double cross(const Point& a, const Point& b, const Point& c) {
double res = (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
return res;
}
bool collinear(const Point& a, const Point& b, const Point& c) {
return fabs(cross(a,b,c)) <= EPS;
}
bool ccw(const Point& a, const Point& b, const Point& c) {
return cross(a,b,c) > EPS;
}
bool cw(const Point& a, const Point& b, const Point& c) {
return cross(a,b,c) < -EPS;
}
double distance(const Point& a, const Point& b) {
return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
// ensure first_point is set before calling it for the sort
bool operator<(const Point& p1, const Point& p2) {
if (collinear(first_point, p1, p2)) {
if (distance(first_point, p1) <= distance(first_point, p2)) return true;
return false;
}
if (ccw(first_point,p1,p2)) return true;
return false;
}
using namespace std;
vector<Point> convexHull(vector<Point> ps) {
vector<Point> res;
if (ps.size() <= 3) return res;
int N = ps.size(), k = 0;
first_point.x = 1e50; // sentinel value
for (int i = 0; i < N; i++) {
if (ps[i].x < first_point.x || (ps[i].x == first_point.x
&& ps[i].y < first_point.y)) {
first_point = ps[i];
k = i;
}
}
ps.erase(ps.begin()+k);
res.push_back(first_point);
sort(ps.begin(),ps.end());
res.push_back(ps[0]);
ps.push_back(first_point);
for (int i = 1; i < N; i++) {
while (res.size() >= 2 && !ccw(res[res.size()-2],res[res.size()-1],ps[i]))
res.pop_back();
res.push_back(ps[i]);
}
return res;
}
Applying the Algorithm (SCUD Busters):
To apply and test the convex hull problem, we take a look at SCUD Busters. This problem can be accessed via:
http://acm.pku.edu.cn/JudgeOnline/problem?id=1264
The problem has a set of kingdoms which have a single power station and a set of houses. A set of missiles are then fired which may or may not hit a kingdom, it hits a kingdom if the missile is inside the polygon defined from the set of points belonging to the kingdom. If the missile hits a kingdom - it loses all its power source. We are asked to calculate the total area in which power is lost for all the kingdoms.
First, we need to calculate the convex hull for each kingdom - this is defined by the statement "The wall surrounding a kingdom is the minimal-perimeter wall that completely surrounds all the houses and the power station that comprise a kingdom". Let's say we have calculated the convex hull for each kingdom K, then given a set of missile locations we need to focus on calculating which (if any) kingdom gets hit by the missile. This is determined by iterating through the list of kingdoms and determining if the missile location is inside the convex polygon defined by the convex hull. If it's inside then the kingdom loses power - if it isn't the kingdom is safe. If it's hit we simply need to calculate the area of the convex polygon defined for the kingdom and we need to sum these up to derive the answer. Note that if a kingdom is hit by 2 or more missiles - it doesn't "lose" more power so a checking mechanism is required to make sure we don't double count the areas. A simple visited array suffices.
The implementation is shown below:
#define EPS 1e-09
class Point {
public:
double x;
double y;
Point() { }
Point(double a, double b) : x(a), y(b) { }
};
Point first_point;
double cross(const Point& a, const Point& b, const Point& c) {
double res = (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
return res;
}
bool collinear(const Point& a, const Point& b, const Point& c) {
return fabs(cross(a,b,c)) <= EPS;
}
bool ccw(const Point& a, const Point& b, const Point& c) {
return cross(a,b,c) > EPS;
}
bool cw(const Point& a, const Point& b, const Point& c) {
return cross(a,b,c) < -EPS;
}
double distance(const Point& a, const Point& b) {
return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
// ensure first_point is set before calling it for the sort
bool operator<(const Point& p1, const Point& p2) {
if (collinear(first_point, p1, p2)) {
if (distance(first_point, p1) <= distance(first_point, p2)) return true;
return false;
}
if (ccw(first_point,p1,p2)) return true;
return false;
}
using namespace std;
vector<Point> convexHull(vector<Point> ps) {
vector<Point> res;
if (ps.size() <= 3) return res;
int N = ps.size(), k = 0;
first_point.x = 1e50; // sentinel value
for (int i = 0; i < N; i++) {
if (ps[i].x < first_point.x || (ps[i].x == first_point.x
&& ps[i].y < first_point.y)) {
first_point = ps[i];
k = i;
}
}
ps.erase(ps.begin()+k);
res.push_back(first_point);
sort(ps.begin(),ps.end());
res.push_back(ps[0]);
ps.push_back(first_point);
for (int i = 1; i < N; i++) {
while (res.size() >= 2 && !ccw(res[res.size()-2],res[res.size()-1],ps[i]))
res.pop_back();
res.push_back(ps[i]);
}
return res;
}
// assumes that the polygon ends with the starting point
double polygonArea(vector<Point> p) {
double res = 0.0;
for (int i = 0; i+1 < p.size(); i++) res += (p[i].x*p[i+1].y - p[i+1].x*p[i].y);
return res / 2.0;
}
bool insideConvex(const vector<Point>& hull, Point p) {
bool oddNodes = false;
for (int i = 0; i+1 < hull.size(); i++) {
if ((hull[i].y < p.y && hull[i+1].y >= p.y) ||
(hull[i+1].y < p.y && hull[i].y >= p.y)) {
if (hull[i].x + (p.y - hull[i].y) / (hull[i+1].y - hull[i].y) *
(hull[i+1].x - hull[i].x) < p.x) {
oddNodes = !oddNodes;
}
}
}
return oddNodes;
}
int destroyed[111];
int main() {
vector<vector<Point> > points;
int K;
// process the kingdoms and add the convex hull to points
while (cin >> K && K != -1) {
int x, y;
vector<Point> temp;
for (int i = 0; i < K; i++) {
cin >> x >> y;
temp.push_back(Point(x,y));
}
points.push_back(convexHull(temp));
}
// process the missiles
int mx, my;
double res = 0.0;
while (cin >> mx >> my) {
Point p(mx,my);
for (int j = 0; j < points.size(); j++) {
// if its not destroyed already and inside the kingdom - destroy it
if (!destroyed[j] && insideConvex(points[j], p)) {
res += polygonArea(points[j]);
destroyed[j] = 1;
break;
}
}
}
printf("%.2f\n",res);
return 0;
}
