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 Monotone Chain algorithm -- also known as Andrew's algorithm.
This algorithm was first introduced by A.M. Andrew in 1979 in the paper "Another Efficient Algorithm for Convex Hulls in Two Dimensions." This algorithm also has $O(n \log n)$ complexity which again is mostly due to sorting. The algorithm will compute the convex hull by computing two pieces: an "upper" and "lower" hull. The upper hull is the piece of the hull you would see looking down from above and the lower hull is the piece of the hull you would see looking up from below. One benefit this algorithm has over the Graham scan is that it will only require us to sort according to x-value instead of by angle -- This is a computationally simpler problem.
TODO: Reread some of this and rephrase some steps.
We now need to split the points into our two categories: $x_L$ and $x_U$. We can do this by taking each x-coordinate and computing the y-coordinate that corresponded to being on the line between $x_L$ and $x_U$. If the y-coordinate of the point is above the line then it belongs to $x_U$ and if it is below then it belongs to $x_L$.
We write a function which splits the points into $x_L$ and $x_U$ below.
In [2]:
function split_xL_xU{T<:Real}(psort::Vector{Point{2, T}})
# Allocate space to identify where points belong
# We make endpoints true because we want them in both sets
npts = length(psort)
xL_bool = Array(Bool, npts)
xL_bool[1] = true; xL_bool[end] = true
xU_bool = Array(Bool, npts)
xU_bool[1] = true; xU_bool[end] = true
# Get the endpoints and slope
pmin = psort[1]; pmax = psort[end]
ls = LineSegment(pmin, pmax)
# We check whether it is above or below the line between
# the points, but could also check whether it was clockwise
# or counter-clockwise (above/below seems like less operations)
# Iterate through all points except endpoints
@inbounds for i=2:npts-1
# Pull out x and y values for ith point
p_i = psort[i]
xi, yi = p_i[1], p_i[2]
# Compute where y on line would be
yi_line = evaluatey(ls, xi)
if yi_line <= yi
xL_bool[i] = false
xU_bool[i] = true
else
xL_bool[i] = true
xU_bool[i] = false
end
end
xL = psort[xL_bool]
xU = psort[xU_bool]
return xL, xU
end
Out[2]:
In this step, we will use an approach that is similar to the Graham scan algorithm. If we are finding the upper (lower) hull then we will begin with $x_{\text{min}}$ ($x_{\text{max}}$) and work our way towards $x_{\text{max}}$ ($x_{\text{min}}$).
Just as in the Graham scan, we will need to identify the direction of the turn generated by three points. Imagine we have three points $p_{i-1}, p_{i}, p_{i+1}$.
We follow until we reach the end point which will be $x_{\text{max}}$ when we are finding the lower hull and $x_{\text{min}}$ when we are finding the upper hull.
Notice that this is the same function from the convex hull notebook that we used in the Graham scan algorithm.
In [3]:
function _monotonechain{T<:Real}(p::Vector{Point{2, T}})
# Sort points
psort = sort!(p)
# Split into upper and lower
xl, xu = split_xL_xU(psort)
lh = wrappoints(xl)
uh = wrappoints(reverse!(xu))
return [lh[1:end-1]; uh[1:end-1]]
end
Out[3]:
In [4]:
srand(1112016)
p = [Point(randn(), randn()) for i=1:250]
ep = _monotonechain(p)
# Add the first point to end to make plot pretty
push!(ep, ep[1]);
In [5]:
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[5]: