Graham’s Scan Convex Hull Algorithm
Was spending my free time working through Real World Haskell by O’Sullivan, Stewart, and Goerzen. I can appreciate how many of the exercises end up being very mathy. For instance, the Chapter 3 exercises culminate in an implementation of Graham’s scan algorithm for finding the Convex Hull of a finite set of points in the plane. Here’s my implementation.
Input: A finite list of points \( (x,y) \) in the plane.
Sample Input: (0,0), (0.5, 1), (1,0), (1,1), (0,1), (10,10)
Output: The (ordered) list of vertices of the smallest convex polygon containing all of the input points.
Sample Output: (10,10), (1,0), (1,1), (0,1)
By ordered we mean that traversing the list of vertices should be analogous to taking a counterclockwise walk around the outside edge of the polygon, without cutting through the interior.
Strategy: By extremal point we mean a point in the input set that is a vertex of the final convex polygon. Our task is to identify all of the extremal points and order them.
We first find an obviouslyextremal point, then we use a method not unlike ray tracing to order the remaining points. Once the remaining points are ordered, we can determine which points are extremal by taking a walk along the points, making sure we only ever take left turns (by throwing away points when we don’t make a left).

Remove any duplicates from the input list.

Find the leftmost point on the input list (using lowest \(y\)coordinate to break ties). Call it
b
. The description ofb
makes it clear thatb
is extremal. 
For each of the remaining points
p
, find the slope and distance betweenb
andp
. If multiplep
s form the same slope withb
, then obviously the ones closer tob
can’t be extremal, so only keep the furthest.At this step, we worry about what will happen to points directly above
b
, since the slope will be undefined. Haskell’sDouble
implementation will treat these slopes asInfinity
. 
Order the remaining points by increasing slope. Then append
b
to the front of the ordered list. We now have a listps
of points that starts withb
, is ordered by increasing slope, and contains at most one point for each slope.A quick check in GHCi shows
(1.0 / 0.0) > 1.0
evaluates toTrue
, so if there were any points directly aboveb
, the furthest such point will be the last entry ofps
. 
Now, we take a walk along our list, looking at three consecutive points at a time. For each triple of points
p1
,p2
,p3
we imagine that we travel fromp1
throughp2
top3
, and we calculate the direction we turn atp2
.Since our list is ordered in increasing slope, anything other than a left turn at
p2
indicates thatp2
can’t be extremal, so in that case we throw it away. If we do make a left turn atp2
, then because of how our list is orderedp2
is extremal, and we can then move on to considering the trip fromp2
throughp3
top4
, and so on.
Let’s take a look at the code. The whole file is 98 lines, and you can see it in its entirety on GitHub, but we dissect it below.
First, imports and a general convenience function.
1
2
3
4
5
6
7
8
9
import Data.List (groupBy, sortBy)
removeDuplicates :: Eq a => [a] > [a]
 ^ Why isn't this in Prelude?
removeDuplicates = foldr skipIfElem []
where
skipIfElem x ys = if x `elem` ys
then ys
else x : ys
removeDuplicates
does exactly what it sounds like it does.
Edit: The function I’m looking for is
nub
fromData.List
.
Next we introduce a datatype for points in the plane and a few functions that we’ll need later.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
  Type for encoding points in the Cartesian plane.
data Point = Point { xProj :: Double, yProj :: Double }
deriving (Eq, Read)
  Show instance displays points in usual mathematical notation.
instance Show Point where
show (Point x y) = "(" ++ show x ++ "," ++ show y ++ ")"
p :: (Double, Double) > Point
 ^ Function for creating points in usual mathematical notation.
p (x, y) = Point x y
slope :: Point > Point > Double
 ^ Calculates the slope between two points.
slope (Point x1 y1) (Point x2 y2) = (y2  y1) / (x2  x1)
norm :: Point > Point > Double
 ^ Calculates the taxicab distance between two points.
norm (Point x1 y1) (Point x2 y2) = abs (x2  x1) + abs (y2  y1)
coordinateSort :: [Point] > [Point]
 ^ Sorts by lowest xcoord, then by lowest ycoord.
coordinateSort = sortBy (\(Point x1 _) (Point x2 _) > compare x1 x2)
. sortBy (\(Point _ y1) (Point _ y2) > compare y1 y2)
Edit: Using a type alias, such as
type Point = (Double, Double)
instead of adata
declaration allows us to use the defaultOrd
,Read
, andWrite
instances, greatly simplifying the code. These changes are reflected on the GitHub version ofexB12.hs
.
We could implement this algorithm using angle instead of slope, but slope is easier to calculate and ends up being equivalent.
I.e. if we were to calculate the angle that each point forms with the ray moving in the positive \(x\) direction emanating from a fixed point, we’d end up needing to calculate atan
of the slope, anyway, and since atan
is monotonic, sorting on atan . slope
is the same as sorting on slope
, so why compute atan
?
Likewise, we chose to use the taxicab metric instead of the usual distance formula, since the taxicab metric is easier to calculate (we avoid sqrt
) and results in the same sort.
Next we define a datatype for relative directions and a function that discerns the direction we need to turn at p2
in order to proceed to p3
from p1
.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
data Direction = GoLeft  GoStraight  GoRight  GoBackwards  GoNowhere
deriving (Eq, Read, Show)
direction :: Point > Point > Point > Direction
 ^ Given points p1, p2, and p3, returns the direction to turn at `p2`
 when traveling from `p1` through `p2` to `p3`.
direction p1@(Point x1 y1) p2@(Point x2 y2) p3@(Point x3 y3) = do
 Check a few easy edge cases.
if p1 == p2  p2 == p3  direction is not welldefined
then GoNowhere
else if p1 == p3  direction is backwards
then GoBackwards
else do
 Here, we do a coordinate transformation that moves `p1` to the
 origin and moves `p2` to (1,0). This puts `p3'` somewhere in the
 plane, and then we branch on the trichotomy of `y3'`.
let det = (x2  x1) ** 2 + (y2  y1) ** 2
a = (x2  x1) / det
b = (y2  y1) / det
c =  (y2  y1)
d = x2  x1
x3' = a * (x3  x1) + b * (y3  y1)
y3' = c * (x3  x1) + d * (y3  y1)
if y3' > 0  `p3'` is in the upper half plane
then GoLeft
else if y3' < 0  `p3'` is in the lower half plane
then GoRight
else if x3' > 1  `p3'` is on the xaxis, beyond (1,0)
then GoStraight
else GoBackwards  `p3'` is on the xaxis behind (1, 0)
Lines 17 through 23 perform a coordinate transformation on p3
.
We could pull the coordinate transformation out of direction
and make it a distinct toplevel function, but I don’t see any reason to at this point.
Wikipedia suggests embedding our points in three space and using the \(z\)coordinate of the cross product of the vector \(\vec{p_1 p_2}\) with the vector \(\vec{p_2 p_3}\) to determine direction. That approach works, but it misses a good number of edge cases, and taking care of those edge cases results in a function that’s not any simpler than the one at hand.
We’re ready for the algorithm itself:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
grahamScan :: [Point] > [Point]
 ^ Takes a list of points in the plane and returns the vertices of
 the smallest convex polygon containing the input points, ordered
 counterclockwise around the perimeter.
grahamScan input = walkPerimeter . sort $ input
where
sort ps = (b :)
. map (\(p,_,_) > p)
 ^ forget slopes and norms, and push `b`
. map last
 ^ take the furthest from each slope group
. map (sortBy (\(_,_,n1) (_,_,n2) > compare n1 n2))
 ^ sort each slope group by norm
. groupBy (\(_,s1,_) (_,s2,_) > s1 == s2)
 ^ group by slope
. sortBy (\(_,s1,_) (_,s2,_) > compare s1 s2)
 ^ sort by slope
. map (\p > (p, slope b p, norm b p)) $ bs
 ^ calculate slope and norm for remaning points
where
(b : bs) = coordinateSort . removeDuplicates $ ps
 ^ remove duplicates and find the lowestleftmost point `b`
walkPerimeter (p1 : p2 : p3 : ps) =
 ^ examine three points at a time
 notably, assumes p1 is good
if direction p1 p2 p3 == GoLeft
then p1 : (walkPerimeter (p2 : p3 : ps))  GoLeft > p2 is good
else walkPerimeter (p1 : p3 : ps)  not GoLeft > p2 is bad
walkPerimeter ps = ps  base case, fewer than three points > end
The algorithm flows through two phases: sort
massages our input data into a form that is easy to recurse on, walkPerimeter
process the sorted data, three points at a time.
I’ve annotated sort
and walkPerimeter
, and you’ll see they implement the fivestep strategy described at the start of this post.
One thing to remember while reading the annotations, though: read from bottom to top.
This is functional language, and sort
is a composition of a bunch of functions.
The rightmost (here, lowest) function is the first to touch our input data.
walkPerimeter
assumes that its first argument is a point that we intend to keep (it has no way of deleting the first point it’s fed), and this is okay since we know b
is an extremal point.
Let’s give it a test drive in GHCi.
Let’s noiseup our data a little and see what happens.
1
2
3
4
noiseEm :: [Point] > [Point]
noiseEm = concatMap noiseIt
where
noiseIt (Point x y) = [Point (x + e) (y + e)  e < [0.1, 0, 0.1]]
I’d like to write some robust some testing suites, but I don’t really know how. Anyway, that’s enough for tonight.