This notebook is the beginning of a sequence of notebooks that discuss numerical convex hull algorithms. We will focus on computing convex hulls for 2 dimensions. These are intended to give an idea of what is involved in finding the convex hull. If needed one can find other more technically complete notes (some of which are mentioned in the references) that will be helpful in dealing with more general cases. All related code to these operations can be found on github under the CHull2D.jl package.

Convex Hulls

In this notebook we present the convex hulls and discuss briefly why they matter. Additionally, we will compare the performance of two (eventually we will add a third) algorithms that are discussed in other notebooks of this series.

What is a convex set

Mathematically defined, a set $X$ is convex if the mixture of any two points of the set is also in the set.

More succinctly, $X$ is convex if $\forall \lambda \in (0, 1)$ and $\forall$ $x, y \in X$ then $z \equiv \lambda x + (1-\lambda)y \in X$.

Less formally, we can think of this as the straight line drawn between two points as being completely inside the set $X$.

The extreme points of a convex set is the minimal set of points that can be used to generate the entire set. Mathematically, we would say that an extreme point of a convex set, $S$, is a point in $s \in S$ which does not lie in any open line segment joining two points of $S$. For shapes with finite edges, we can think about these being the vertices.

For example, the extreme points of a square in $\mathcal{R}^2$ would simply be its 4 corners. On the other hand, the extreme points of a circle are its entire circumference.

Why do we care

Convex sets are front and center for many important results within many fields -- This includes our field of economics. Representing a convex set on paper --or in your mind-- is simple, but representing such a set on a computer requires a little more thought. In this notebook, we will show you one way a convex set can be represented on a computer and how to generate the convex hull of a given set of points.

More generally, representing a set of points in a computer is very difficult because in order to keep track of the points the computer needs to represent every point from the set. This becomes impossible as soon as the set of points isn't discrete (we can't keep track of an infinite number of points in a computer). Luckily, convex sets can be represented using only their extreme points -- This greatly reduces the amount of information that we need to keep.

In later notebooks, we will use the concept of a convex hull to represent sets of continuation values for dynamic games.

Numerical Convex Hull Algorithms

Numerical algorithms to generate convex hulls complete the following exercise: Given a set of points, tell me which points are necessary to represent the convex hull of that set. This exercise is simple to do on paper, you draw the least distance line possible that circles the set of points -- i.e. you draw lines around the outside of the set to connect the extreme points. This task becomes slightly harder once we move to a computer because we must find a practical way to do this generally.

In this sequence of notebooks, we will cover three algorithms that are used to numerically compute the convex hull:

• The Graham Scan
• Monotone Chain
• Quickhull (This will be implemented at a later date)

At the end of this notebook we will run some simple comparisons of the three algorithms. Before doing this, we will introduce some prerequisite material that will be useful for developing these algorithms.

Useful Prerequisite Material

Here we introduce the previously mentioned prerequisite material that will simplify the exposition of the algorithms later.

In our exposition below, we will consider a set of points $P$. Let $a = (a_1, a_2)$, $b = (b_1, b_2)$, and $c = (c_1, c_2)$ be three points in $P$. We will think of these points as being ordered by $(a, b, c)$. We will think about two vectors created by these points: $u = \vec{ab}$ and $v = \vec{ac}$. Finally, we will consider the triangle formed by the three points $\triangle abc$.

NOTE: Some of the following material will be repetitive for some readers, feel free to skip as much of this section as you feel comfortable skipping.

Angles between Points

At the base of much of what we do is determine the angle between three points. We will use this to sort points by their angle relative to some initial point and also to determine whether three points make a clockwise or counter-clockwise turn.

Consider our three points $(a, b, c)$ and the two vectors formed by them $\vec{u}$ and $\vec{v}$. One way we could find the angle formed by the three points is to find the angle, $\theta$, formed between our two vectors $u$ and $v$. This can be computed by

$$\theta = \arccos \left( \frac{u \cdot v}{\|u\| \|v\|} \right)$$

This would tell us the exact angle. It also tells us clockwise or counter-clockwise by the angle being greater or less than 180 degrees.

While this approach would certainly work, it requires more operations than we need. For most of what we need to do, it will suffice to determnine whether the three points generate a clockwise or counter-clockwise turn. This can be done by using a cross product of the vectors $u$ and $v$.

Wikipedia tells us,

In computational geometry of the plane, the cross product is used to determine the sign of the acute angle defined by three points $p_1=(x_1,y_1)$, $p_2=(x_2,y_2)$ and $p_3=(x_3,y_3)$. It corresponds to the direction of the cross product of the two coplanar vectors defined by the pairs of points $p_1, p_2$ and $p_1, p_3$, i.e., by the sign of the expression $P = (x_2-x_1)(y_3-y_1)-(y_2-y_1)(x_3-x_1)$. In the "right-handed" coordinate system, if the result is 0, the points are collinear; if it is positive, the three points constitute a positive angle of rotation around $p_1$ from $p_2$ to $p_3$, otherwise a negative angle. From another point of view, the sign of $P$ tells whether $p_3$ lies to the left or to the right of line $p_1, p_2$.

We can use the formula they give us to compute the last (and relevant) element of the cross-product.

$$P = (b_1 - a_1)(c_2 - a_2) - (b_2 - a_2)(c_1 - a_1)$$

If $P$ is positive then we moved counter-clockwise and if $P$ is negative then we moved clockwise.

Point Relative to Triangle

Another important piece of what we do will be to determine whether a point lies inside of or outside of a triangle. We will do this by using Barycentric Coordinates.

Consider a point $x$ and our triangle $\triangle abc$. We want to know whether $x \in \triangle abc$. The key insight will be that triangles are convex sets and so we can represent any point within the triangle as a convex combination of the three points that make up the triangle -- In our case $a$, $b$, and $c$.

Thus $x \in \triangle abc$ if and only if $\exists$ $\lambda_1, \lambda_2, \lambda_3$ where $\lambda_i \in [0, 1]$ and $\sum_i \lambda_i = 1$ such that $x = \lambda_1 a + \lambda_2 b + \lambda_3 c$.

If we use $\lambda_3 = 1 - \lambda_1 - \lambda_2$ and rearrange this equation then we get:

$$\begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} a_1 - c_1 & b_1 - c_1 \\ a_2 - c_2 & b_2 - c_2 \end{bmatrix} \begin{bmatrix} \lambda_1 \\ \lambda_2 \end{bmatrix} + \begin{bmatrix} c_1 \\ c_2 \end{bmatrix}$$

Define $T = \begin{bmatrix} a_1 - c_1 & b_1 - c_1 \\ a_2 - c_2 & b_2 - c_2 \end{bmatrix}$, then we can rearrange this again to get an expression for $\lambda_1$ and $\lambda_2$:

$$ \begin{bmatrix} \lambda_1 & \lambda_2 \end{bmatrix} = T^{-1} \begin{bmatrix} x_1 - c_1 \\ x_2 - c_2 \end{bmatrix}$$

This equation pins down $\lambda_1$ and $\lambda_2$ and once we have these then we are able to determine whether the conditions on $\lambda_i$ hold. If

• $\lambda_1 \in [0, 1]$
• $\lambda_2 \in [0, 1]$
• $\lambda_1 + \lambda_2 \leq 1$

then we have found the $(\lambda_1, \lambda_2, \lambda_3)$ that we need for our conditions to hold and we know $x \in \triangle abc$.

Wrapping a Sequence of Points

In both the monotone chain algorithm and the Graham scan, we will be required to take a sequence of sorted points and wrap our convex hull around them in the right way. We can do this by assuring that the direction of the turn between any three points in the convex hull is the same direction (worth drawing a picture to convince yourself of this fact) -- We will use the counter-clockwise direction.

• Start the convex hull with the first two elements of the sequence of points.
• For each element $p_i \in P$ take the convex hull $C$ and the number of points currently in $C$, $m = |C|$, as given.
• Check whether $(C_{m-1}, C_{m}, p_i)$ makes up a clockwise or counter-clockwise turn.
  • If this term is clockwise then $C_{m}$ cannot be a memeber of the convex hull, so we discard it from $C$.
    • If after removing $C_{m}$ from the convex hull there is only one element of the convex hull then add $p_i$ to the convex hull and move to the next point, $p_{i+1}$
    • If there are more than one elements left in the convex hull then check the next element of the convex hull as well -- Check whether $(C_{m-2}, C_{m-1}, p_i)$ makes up a clockwise turn and repeat until there is only one element in the convex hull, or the last two points of the convex hull and $p_i$ make a counter-clockwise turn.
  • If the turn is counter-clockwise then add $p_i$ to the convex hull and move to $p_{i+1}$.

As long as the points are ordered as described in the Graham scan algorithm or the monotone chain algorithm, this will determine the extreme points of the set.

Pruning: Akl–Toussaint Heuristic

Since it is applicable to the algorithms we discuss in the other notebooks, we present the Akl-Toussaint Heuristic. This tool will allow us to reduce the number of points that we need to consider when computing the convex hull by removing points that we know are not elements of the convex hull. This is important because the algorithms that we present are all $O(n \log n)$ complexity meaning that the number of points that we are passing to our actual algorithms will have a significant effect on compute time.

This was introduced in work done by Selim Akl and G. T. Toussaint in their 1978 paper, "A fast convex hull algorithm." Their key insight is that if you draw a quadrilateral using four points that you know are elements of the convex hull, then you are able to immediately disregard a large number of points that you know are not in the convex hull. The four points they propose are the maximal and minimal $x$ and $y$ coordinates.

Algorithm Comparisons

These comparisons were run by drawing $np$ points from a standard normal distribution and then computing the convex hull $nr$ times to make sure each algorithm had a chance to perform well. There are obvious short comings of this comparison because it only is comparing them for a specific type of set -- a set with many points but few extreme points. More careful analysis has been done by others, though I may eventually edit this notebook to include analysis on sets with different properties. The code to do this benchmark can be found in the file benchmark.jl.

Comparisons were done using a 2.2 GhZ laptop running Linux.

Algorithm \ Number of points		25	250	2,500	25,000	250,000	2,500,000
Graham scan		1.6e-5 s	1.6e-4 s	2e-3 s	2e-2 s	3.4e-1 s	4e0 s
Graham scan (with Akl Toussaint)		3.1e-5 s	1.1e-4 s	9e-4 s	8e-3 s	8.5e-2 s	8.7e-1 s
Monotone chain		2.7e-6 s	3.7e-5 s	5e-4 ms	6e-3 s	7.7e-2 s	9.4e-1 s
Monotone chain (with Akl Toussaint)		2.0e-5 s	9.8e-5 s	9e-4 ms	8e-3 s	8.4e-2 s	8.7e-1 s
Quickhull		NA s	NA s	NA s	NA s	NA s	NA s

The first thing to notice it that the algorithms scale just slower than linearly -- This is because of the complexity $O(n \log n)$ which is slightly more than linear.

One thing that I found surprising was that the Akl-Toussaint heuristic did not speed things up until the number of points were very large. It was more powerful when applied to the Graham scan because the Graham scan is slightly slower than the monotone chain algorithm for large numbers of points.

TODO: Try running other experiments like taking sets in which every point is an extreme point etc...

References

TODO: Clean up references.

• http://geomalgorithms.com/a10-_hull-1.html
• https://github.com/intdxdt/convexhull.jl
• Wikipedia -- Specific articles are cited within text
• http://keithbriggs.info/documents/Skarakis_MSc.pdf
• http://www.toptal.com/python/computational-geometry-in-python-from-theory-to-implementation