Implementing QuickHull in 2D

QuickHull is a well known algorithm for generating a convex hull given a set of points. The expected complexity is $O(n\cdot \log(n))$, with a worst case complexity of $O(n^2)$. The algorithm generalizes nicely into 3D, but for now let's just focus on the simpler 2D case.

Algorithm

The algorithm is fairly straightforward:

1. Create an edge by finding the min and max points along some axis (in my case I used the x-axis.)
2. For each edge in the hull (initially only one) we assign a set of external points that can "see" the edge.
3. For each set of points, find the point farthest away from its corresponding edge.
4. The point farthest away (out of all the sets) is added to the convex hull by splitting its associated edge.
5. The remaining points in the edge's set are now re-evaluated. They may either be inside the hull or belong to one of the two newly created edges
6. The algorithm terminates when there are no more external points

Implementation

This algorithm has an expected complexity of $O(n\cdot\log(n))$. To actually get this there are some implementation details worth mentioning. First thing to note is that the convex hull is stored as a linked list. This allows for cheap insertion as well as easily traversing the hull.

	QuickHull.Node = function(v)
	{
	this.v = v;
	this.conflicts = [];
	this.next = undefined;
	};
	QuickHull.Node.prototype.Insert = function(v)
	{
	var n = this.next;
	this.next = new QuickHull.Node(v);
	this.next.next = n;
	return this.next;
	};

view raw convex_hull_as_list.js hosted with ❤ by GitHub

To get the $O(n\cdot\log(n))$ complexity we can't keep reiterating over the full set of points. Because of this we divide the set of points into conflict sets. Each edge has its own conflict set of points that are external to the hull and are closest to that edge. Doing this allows us to only iterate over the points we are interested in. To add a new point to the convex hull we iterate over all conflict sets and compute the distance between the edge and the points in its set. The edge with the farthest point will be chosen to split. Here we add caching so that an edge's farthest point is never calculated twice.

	//find farthest point from all segments
	var startNode = undefined;
	var farthestIndex = -1;
	var maxDistance = 0;
	convexNode = convexHull;
	do
	{
	if(convexNode.cachedInfo == undefined) //we only want to compute this if something in the set has changed
	{
	convexNode.cachedInfo = this.FarthestPointToSegment(convexNode.v, convexNode.next.v, convexNode.conflicts);
	}
	var info = convexNode.cachedInfo;
	var d2 = info[1];
	if(d2 >= maxDistance)
	{
	maxDistance = d2;
	farthestIndex = info[0];
	startNode = convexNode;
	}
	convexNode = convexNode.next;
	}while(convexNode != convexHull);
	var newHullPt = startNode.conflicts[farthestIndex];
	startNode.cachedInfo = undefined; //we're splitting the edge so the points info is no longer valid
	startNode.Insert(newHullPt);
	var conflicts = startNode.conflicts;
	conflicts.splice(farthestIndex, 1);
	pointsLeft-= 1;

view raw find_hull_pt.js hosted with ❤ by GitHub

Splitting an edge will (hopefully) put most of its conflict points into the convex hull. At this point we need to re-evaluate the points in the edge's set and determine if they are within the hull, or belong to the set of one of the newly created edges. It's important to remember to invalidate the cache of the edge here since we're technically reusing the same vertex node for a new edge.

	//now fixup conflicts in startNode based on the added triangle
	for(var i=0; i<conflicts.length; i++)
	{
	var pt = conflicts[i];
	var uvw = Barycentric(pt, startNode.v, startNode.next.v, startNode.next.next.v);
	var u = uvw[0];
	var v = uvw[1];
	var w = uvw[2];
	//use barycentric coordinates to determine point location
	// startNode (u =1)
	// /\
	// / \
	// / \
	// /______\ startNode.next (v=1)
	// next.next (w=1)
	if(u >= 0 && v >= 0 && w >= 0) //pt is in new triangle
	{
	//pt is finally covered by hull so remove it from pts
	conflicts.splice(i, 1);
	i--;
	pointsLeft -= 1;
	}else if(v < 0) //edge between next next and start
	{
	//this can only happen in first triangle
	startNode.next.next.conflicts.push(pt);
	conflicts.splice(i,1);
	startNode.next.next.cachedInfo = undefined; //modifying conflicts list so cache is invalid
	i--;
	}else if(u < 0) //edge between next and next next
	{
	startNode.next.conflicts.push(pt);
	conflicts.splice(i,1);
	startNode.next.cachedInfo = undefined; //modifying conflicts list so cache is invalid
	i--;
	}
	}

view raw update_verts.js hosted with ❤ by GitHub

We loop over the above code as long as the number of external points is non-zero. The logarithmic behavior in this algorithm works very well. You can especially see this by generating 100s of random points and hitting the step button. You will see that the number of points get dramatically reduced with each step.

That's it! Hope this has been useful for someone.