In [1]:
using CHull2D
using FixedSizeArrays
using PyPlot
This notebook is a piece of a sequence of notebooks that discuss numerical convex hull algorithms. The goal of these algorithms is to take a set of points $P$ and output a set of extreme points $EP$ such that every point in $P$ is contained in $\text{co}(EP)$.
In this notebook we present the Graham Scan algorithm.
This algorithm was initially introduced by Ronald Graham in 1972 in "An Efficient Algorithm for Determining the Convex Hull of a Finite Planar Set." The algorithm will have $O(n \log(n))$ complexity mostly due to the required sorting done within the algorithm. As you will observe, the way in which the algorithm traverses the points will lead for it to be a better fit for sets of points which have many extreme points.
In our implementation, we will use the point with the minimum $y$ value as the initial point. Additionally, our rotation direction will be counter-clockwise.
Call our initial point $a$. What we really need to know is whether for any two points $b, c \in P$ whether the angle formed by $\vec{ab}$ and $\vec{bc}$ is greater than or less than 180 degrees. This can be determined by using the cross-product. If the cross-product is negative then we have moved clockwise; if the cross-product is positive then we have moved counter-clockwise (as desired). This leaves us to simply sort by comparing the appropriate cross-product between vectors -- It turns out we can directly write this as:
$$\vec{ab} \times \vec{bc} = (b_1 - a_1) (c_2 - a_2) - (b_2 - a_2) (c_1 - a_1)$$We wrote a function to do this in the ConvexHulls notebook.
As previously mentioned, we will begin with the point which has the smallest $y$ value. We will then work through our ordered points three at a time.
Imagine we have three points $p_{i-1}, p_{i}, p_{i+1}$.
We follow until we reach the final point (which is also our starting point).
This is the wrappoints
function from the ConvexHulls notebook.
To summarize the algorithm we could write:
In [2]:
function _grahamscan{T<:Real}(p::Vector{Point{2, T}})
# Get yminind
xminind = indmin(p)
xmin = p[xminind]
# Create function to sort by angles
lt(a, b) = xmin == a ? true : xmin == b ? false : ccw(xmin, a, b)
psort = sort!(p, lt=lt)
# Add the starting point at end so we go full circle
push!(psort, xmin)
# Call function that works around points
ep = wrappoints(psort)
return ep[1:end-1]
end
Out[2]:
In [3]:
srand(1112016)
p = [Point(randn(), randn()) for i=1:250]
ep = _grahamscan(p)
# Add the first point to end to make plot pretty
push!(ep, ep[1]);
In [7]:
fig, ax = subplots()
ax[:scatter]([el[1] for el in p], [el[2] for el in p], color="r")
ax[:plot]([el[1] for el in ep], [el[2] for el in ep], color="k")
Out[7]: