## Graham scan algorithm in Haskell

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)) )
where
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)))
where
pivot = minimum points
compareAnglesWith p a b = compare (sineAngle a p b) 0.0

convexHull points = hullify (extend (sortPoints points))
where
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.

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.

```*Main> convexHull points