Convex hull algorithms fall into the category of “computational geometry”. Computational geometry has a wide range of applications in many different fields, ranging from computer graphics, to ecology, city planning, and much, much more.
In this post I will cover one of the simplest ways to compute the convex hull of a polygon: the gift wrapping method, also known as the Jarvis March for the programmer who developed it. But what is a convex hull anyway?
The Convex Hull
To put it simply, the convex hull is the set points that form a polygon containing all of the given points in a data set.
Gift wrapping is a line sweeping algorithm, that is, starting from the point with the smallest y-axis value, we “sweep” a horizontal line until in whichever direction is positive until we find the next point on the hull, and restart the process. This repeats until we arrive back at our starting point with the set of points encountered on our journey forming the convex hull of our data set.
Data Structures and Sample Data
Since we are dealing with geometry and points on a plane, it only makes sense that we would have a data structure to represent a point. This is a simply plain old data structure consisting of an x value, a y value, and a label:
struct point {
int x, y;
char c;
};
Our sample set is going to consist of 16 points:
A B C D E F G H I J K L M N O P
x: 3 11 6 4 5 8 1 7 9 14 10 16 15 13 3 12
y: 9 1 8 3 15 11 6 4 7 5 13 14 2 16 12 10
Determining which direction is positive
As i mentioned above, to perform the “line sweep” we need to know which direction to sweep in, as it must be positive. To do this we will be determining the angle between two sets of points using a function, theta().
float theta(point p1, point p2) {
int dx, dy, ax, ay;
float t;
dx = p2.x - p1.x; dy = p2.y - p1.y;
ax = abs(dx); ay = abs(dy);
t = (ax + ay == 0) ? 0 : (float)dy/(ax+ay);
if (dx < 0) t = 2-t;
else if (dy < 0) t = 4+t;
return t*90.0;
}
The Gift Wrapping Algorithm (The Jarvis March)
Jarvis derived the gift wrapping method in the early 1970’s to compute the convex hull of a set of points, inspired by the way one would wrap a gift in wrapping paper.
The algorithm makes use of a sentinel node, consisting of the point with smallest y value, placed in the N+1 position of the array of points.
vector<point> wrap(vector<point>& p, int N) {
int i, min, M;
float th, sweepAngle;
vector<point> hull;
for (min = 0, i = 1; i < N; i++)
if (p[i].y < p[min].y) min = i;
p[N] = p[min];
We’re now ready to find the M points beling to the convex hull, remember the first point on the convex hull is the point with the lowest y-value, thus on the first go through of the loop, were doing little more than marking our sentinel node as being the first point on the convex hull.
for (M = 0; M < N; M++) {
swap(p[M], p[min]);
hull.push_back(p[M]);
Then, we want to sweep our line up, looking for the next point on the convex hull. This accomplised by computing the “sweep angle”, which is the angle formed from the horizontal of the line between p[M-1] and p[M]. We are looking for the point that has the smallest angle in the set of points with angles greater than our sweepAngle value.
min = N; sweepAngle = th; th = 360.0;
for (i = M + 1; i <= N; i++) {
if (theta(p[M], p[i]) > sweepAngle) {
if (theta(p[M], p[i]) < th) {
min = i; th = theta(p[M], p[min]);
}
}
}
if (min == N) return hull;
}
}
And thats our whole algorithm. When it encounters the point we first started at, it ends and returns the set of points in the convex hull. So lets put it all together and see what we come up with:
#include <iostream>
#include <iomanip>
#include <vector>
using namespace std;
const int maxN = 16;
struct point {
int x, y;
char c;
};
float theta(point p1, point p2) {
int dx, dy, ax, ay;
float t;
dx = p2.x - p1.x; dy = p2.y - p1.y;
ax = abs(dx); ay = abs(dy);
t = (ax + ay == 0) ? 0 : (float)dy/(ax+ay);
if (dx < 0) t = 2-t;
else if (dy < 0) t = 4+t;
return t*90.0;
}
vector<point> wrap(vector<point>& p, int N) {
int i, min, M;
float th, sweepAngle;
vector<point> hull;
for (min = 0, i = 1; i < N; i++)
if (p[i].y < p[min].y) min = i;
p[N] = p[min]; th = 0.0;
for (M = 0; M < N; M++) {
swap(p[M], p[min]);
hull.push_back(p[M]);
min = N; sweepAngle = th; th = 360.0;
for (i = M + 1; i <= N; i++) {
if (theta(p[M], p[i]) > sweepAngle) {
if (theta(p[M], p[i]) < th) {
min = i; th = theta(p[M], p[min]);
}
}
}
if (min == N) return hull;
}
}
void printPoints(const std::vector<point>& data, int n) {
cout<<" ";
for (int i = 0; i < n; i++)
cout<<setw(2)<<data[i].c<<" ";
cout<<endl<<"x: ";
for (int i = 0; i < n; i++)
cout<<setw(2)<<data[i].x<<" ";
cout<<endl<<"y: ";
for (int i = 0; i < n; i++)
cout<<setw(2)<<data[i].y<<" ";
cout<<endl;
}
int main() {
// maxN + 1 so there is room for sentinel node;
vector<point> data = { {3, 9, 'A'}, {11, 1, 'B'}, {6, 8, 'C'}, {4, 3, 'D'}, {5, 15, 'E'},
{8, 11, 'F'}, {1, 6, 'G'}, {7, 4, 'H'}, {9, 7, 'I'}, {14, 5, 'J'},
{10, 13, 'K'},{16,14,'L'}, {15, 2, 'M'},{13, 16, 'N'},{3, 12, 'O'},
{12, 10, 'P'}};
printPoints(data, maxN);
cout<<"Convex Hull: \n";
vector<point> hull = wrap(data, maxN);
printPoints(hull, hull.size());
return 0;
}A B C D E F G H I J K L M N O P
x: 3 11 6 4 5 8 1 7 9 14 10 16 15 13 3 12
y: 9 1 8 3 15 11 6 4 7 5 13 14 2 16 12 10
Convex Hull:
B M L N E O G D
x: 11 15 16 13 5 3 1 4
y: 1 2 14 16 15 12 6 3
One can see by the arrangement of the loops, that each point on the hull (h) must check against each point in the set of points (n), thus this algorithm runs in O(nh). In the worst case, such as an empty circle of points, the algorithm devolves to exponential time instead of linear. There are several other convex hull algorithms with O(n log n) performance, but as i mentioned in the beginning, this is the simplest of the bunch, both to implement and understand.
Learning More…
The algorithms and helper functions discussed in this article are covered comprehensively in the 2nd edition (1992) of “Algorithms in C++” By Robert Sedgewick. This book contains several chapters on computational geometry amongst its many other data structure and algorithm topics and is well worth the read.
