I am reading the excellent book Real World Haskell by Bryan O’Sullivan, Don Stewart, and John Goerzen, boh the paper version and the on-line version. I did a lot of functional programming at university (until 1993), but more recently, XSLT is the only functional language that I have used for real projects. I hope to use Haskell or Erlang in one of my (near) future projects.

One of the exercises in chapter 3 of Real World Haskell is programming the Graham scan algorithm for computing a two-dimensional complex hull. Some comments in the on-line versions of the book give good solutions, like the ones by Jarkko Sakkinen, Greg Burri, and Hynek Schlawack.

I added my solution to the comments, but unfortunately, the indenting was not preserved. There were also some slight improvements I made later. So here it is, with indents.

Beware: Some commenters have found a bug in this code. Thank you, enoksrd, and Nathan. I’ll revise my code as soon as I find the time for it.

import Data.List

-- The cross-product equals |a|*|b|*sin(apb), so sin(apb) = (a x b) / (|a| * |b|)
sineAngle (ax, ay) (px, py) (bx, by) =
  (cross (ax - px, ay - py) (bx - px, by - py)) /
  ( (length (ax - px, ay - py)) * (length (bx - px, by - py)) )
    cross (x1, y1) (x2, y2) = x1*y2 - x2*y1
    length (x, y) = sqrt ((x^2) + (y^2))

sortPoints points =
  pivot : (sortBy (compareAnglesWith pivot) (delete pivot (nub points)))
    pivot = minimum points
    compareAnglesWith p a b = compare (sineAngle a p b) 0.0

convexHull points = hullify (extend (sortPoints points))
    extend ps = [last ps] ++ ps ++ [head ps]
    hullify (p1:p2:p3:ps)
      -- If sineAngle is greater than zero, we turn right.
      | (sineAngle p1 p2 p3) > 0 = p2:(hullify (p2:p3:ps))
      | otherwise = (hullify (p2:p3:ps))
    hullify _ = []

I have not used the Direction datatype, as suggested in exercise 10, because to my taste it just added unnecessary fluff.

The solution is made of three functions:

SineAngle computes the sine of the angle between two segments AP and PB, using the cross product formula. Note that we take the cross-product of two 2-dimensional vectors, which results in a (perpendicular) vector in the third dimension.

SortPoints sorts the points according to the Graham scan algorithm. The pivot point (minimum points) is the one with lowest x-coordinate (and the lowest y-coordinate if there is a tie). The other points are sorted according to their angle w.r.t. the pivot. The pivot must be the first point in the list, and we nub the list to remove duplicate points.

To compute the convex hull, we sort the points. Then we append the first and prepend the last points, because to determine the turn (left, right, straight) at a point, we need its predecessor and successor, which would otherwise not be available for the first and the last point. Then we determine the turn for each point, and keep every point that is at a right turn (hullify).

My implementation fails for an empty set of points, and for a set of points on a line.

To illustrate how it works, I used

points = [(1.0,1.0), (0.0,4.0), (0.0, 6.0), (3.0,5.0), (4.0,4.0), (4.0,1.0), (3.0,3.0), (2.0,2.0), (5.0, 5.0)]
computing the convex hull of a set of points, using the Graham scan.

computing the convex hull of a set of points, using the Graham scan.

This set of points is shown in the figure. The red filled point is the pivot. The red circled points are in the convex hull. The green points are not part of the convex hull.

Indeed, Haskell tells me

*Main> convexHull points

I travel clockwise,so I only include points at a right turn in the hull.

Nathan asked what I used to make the figure. I made it using Inkscape, a great open source SVG editor. I use Inkscape a lot, because it is simple yet powerful. It runs on Windows, MacOS and Linuxes. Definitely recommended!