Graham's scan

I just realized today that I have not written on that blog for two years. It is definitely time that I take it more seriously.

For the last months, I have started using the functional programming language Haskell. Since I have never been an expert and since it has been a while, I would like to share with you the steps I have been following in order to implement Graham's scan which computes the convex hull of a set of points in the plane.

First, since we handle 2D points, it is useful to define de corresponding data type. I assume that the coordinates are doubles (Double in Haskell):

-- 2D points
data Point2D = Point2D Double Double
               deriving (Show, Eq)

I recall that, thanks to the expression deriving (Show, Eq), it is possible to display (Show) a point when we use the ghci interpreter and also to test for point equality (Eq), which is done coordinate-wise.

A useful function is the square of the distance between two points:

-- Square of distance between points
squareDistance :: Point2D -> Point2D -> Double
squareDistance (Point2D a b) (Point2D c d)
    = (c - a) * (c - a) + (d - b) * (d - b)

We will also need a sample to test our algorithm:

-- Points
points = [Point2D 0 0, Point2D 2 9, Point2D 3 4,
          Point2D 0 5, Point2D 6 2, Point2D 7 3,
          Point2D 4 1, Point2D 5 7]

A graphical representation can be obtained with the Sage Cell Server :

Now, let me list below the 3 main steps of Graham's algorithm:

1. One identifies the lowest point. If two points share the same y-coordinate, then one chooses the leftmost one. The minimal point is called the origin.
2. One sorts the points in increasing order of angle with respect to the origin. Note that, by choice of the origin, the angle is between 0^∘ and 180^∘. If two points have the same angle, the the comparison is done on the distance with respect to the origin.
3. For each successive triple of points, one verifies if they form a left turn, a right turn or a straight line.
   • If it is a left turn, we move on to the next triple;
   • Otherwise, we eliminate the central point from the convex hull. Then, we keep the two remaining points in the next triple.

Step 1 is easy. It suffices to define a lexicographic order on the points with respect to the y-coordinate first, then the x coordinate.

-- Lexicographic order on points (with respect to y and then to x)
comparePoints :: Point2D -> Point2D -> Ordering
comparePoints (Point2D a b) (Point2D c d)
    | b < d  = LT
    | b > d  = GT
    | b == d = compare a c

This function returns the value LT, EQ or GT (of Ordering type) accordingly. The origin is then retrieved with the help of the expression l'expression

List.minimumBy comparePoints points

To complete step 2, it suffices to introduce an order relation based on the angle of each point in the polar coordinates system, where the origin is the lowest-leftmost point computed before. Notice that when the angle is equal, then the distance with respect to the origin is used. This ensures that the origin is always the first point in the resulting list:

-- Angle comparator between two points with respect to some origin
compareAngle :: Point2D -> Point2D -> Point2D -> Ordering
compareAngle origin@(Point2D a b) p@(Point2D c d) q@(Point2D e f)
        | cross > 0  = LT
        | cross == 0 = Ord.comparing (squareDistance origin) p q
        | cross < 0  = GT
    where cross = (c - a) * (f - b) - (d - b) * (e - a)

The first parameter is the origin. The next two parameters are the two points we wish to compare. This test is very standard in algorithmic geometry: it consists in verifying the sign of the cross product of two vectors induced by the three points. According to whether the sign is positive, negative or null, one deduces that the angle of the first point is smaller, greater or equal to the second one. The reader less familiar with linear algebra might be interested in toying with the following interact to verify the property:

Hence, in the main function, we first sort the points, then we call grahamScan to complete the processing:

-- Convex hull computation
convexHull :: [Point2D] -> [Point2D]
convexHull points = grahamScan (List.sortBy (compareAngle origin) points)
    where origin = List.minimumBy comparePoints points

Notice how short is the scan:

-- Graham scan
grahamScan []         = []
grahamScan [p]        = [p]
grahamScan [p,q]      = [p,q]
grahamScan (p:q:r:ps)
    | compareAngle q p r == GT = p:(grahamScan (q:r:ps))
    | otherwise                = grahamScan (p:r:ps)

We may verify that the function returns what is expected when applied to the set of points mentioned above:

Main> convexHull points
[Point2D 0.0 0.0,Point2D 4.0 1.0,Point2D 6.0 2.0,Point2D 7.0 3.0,
 Point2D 5.0 7.0,Point2D 2.0 9.0,Point2D 0.0 5.0]

Here is the complete code:

-- Imports
import qualified Data.List as List
import qualified Data.Ord as Ord

-- 2D points
data Point2D = Point2D Double Double
               deriving (Show, Eq)

-- Square of distance between points
squareDistance :: Point2D -> Point2D -> Double
squareDistance (Point2D a b) (Point2D c d)
    = (c - a) * (c - a) + (d - b) * (d - b)

-- Points
points = [Point2D 0 0, Point2D 2 9, Point2D 3 4,
          Point2D 0 5, Point2D 6 2, Point2D 7 3,
          Point2D 4 1, Point2D 5 7]

-- Angle comparator between two points with respect to some origin
compareAngle :: Point2D -> Point2D -> Point2D -> Ordering
compareAngle origin@(Point2D a b) p@(Point2D c d) q@(Point2D e f)
        | cross > 0  = LT
        | cross == 0 = Ord.comparing (squareDistance origin) p q
        | cross < 0  = GT
    where cross = (c - a) * (f - b) - (d - b) * (e - a)

-- Lexicographic order on points (with respect to y and then to x)
comparePoints :: Point2D -> Point2D -> Ordering
comparePoints (Point2D a b) (Point2D c d)
    | b < d  = LT
    | b > d  = GT
    | b == d = compare a c

-- Convex hull computation
convexHull :: [Point2D] -> [Point2D]
convexHull points = grahamScan (List.sortBy (compareAngle origin) points)
    where origin = List.minimumBy comparePoints points

-- Graham scan
grahamScan []         = []
grahamScan [p]        = [p]
grahamScan [p,q]      = [p,q]
grahamScan (p:q:r:ps)
    | compareAngle q p r == GT = p:(grahamScan (q:r:ps))
    | otherwise                = grahamScan (p:r:ps)

It is interesting to see that the relations using more space are those about order relations. Once they are done, the core of the algorithm does not need more than 10 lines, including commentaries.