2D Convex hull exercise

Question

I'm reading Real World Haskell, and this is my first try with the language.

This is the result of the exercises of chapter 3:

9. Consider three two-dimensional points a, b, and c. If we look at the angle formed by the line segment from a to b and the line segment from b to c, it either turns left, turns right, or forms a straight line. Define a Direction data type that lets you represent these possibilities.
10. Write a function that calculates the turn made by three 2D points and returns a Direction.
11. Define a function that takes a list of 2D points and computes the direction of each successive triple. Given a list of points [a,b,c,d,e], it should begin by computing the turn made by [a,b,c], then the turn made by [b,c,d], then [c,d,e]. Your function should return a list of Direction.
12. Using the code from the preceding three exercises, implement Graham's scan algorithm for the convex hull of a set of 2D points. You can find good description of what a convex hull is, and how the Graham scan algorithm should work.

For short, the Graham scan algorithm is:

1. find the point \$P\$ with the smallest \$y\$
2. sort the points \$A\$ according to the angle \$(Ox, PA)\$ and add P at the end
3. in the sorted list, form the turns \$(A_{i-1}, A_i, A_{i+1})\$ and remove the points \$A_i\$ for which it is a "turn right".

Now my code is the following:

import Data.List
-- 9
data Direction = TurnLeft
               | TurnRight
               | StraightLine
                 deriving (Show, Eq)

-- 10
data Point = Point {
      x :: Float
    , y :: Float
    } deriving (Show, Eq)

turn :: Point -> Point -> Point -> Direction
turn a b c
      | z > 0     = TurnLeft
      | z < 0     = TurnRight
      | otherwise = StraightLine
    where
      bax = x a - x b
      bay = y a - y b -- coordinates of b->a
      bcx = x c - x b
      bcy = y c - y b -- coordinates of b->c
      z   = bax * bcy - bcx * bay

-- 11
turns :: [Point] -> [Direction]
turns []         = []
turns [_]        = []
turns [_,_]      = []
turns (x:y:z:zs) = (turn x y z) : (turns (y:z:zs))

-- 12
angle :: Point -> Point -> Point -> Float
angle a b c = - (bax * bcx + bay * bcy) /
              (sqrt ((bax * bax + bay * bay) * (bcx * bcx + bcy * bcy)))
    -- more efficient than acos for the algorithm and is the same for sorting
    where
      bax = x a - x b
      bay = y a - y b -- coordinates of b->a
      bcx = x c - x b
      bcy = y c - y b -- coordinates of b->c

argmin :: Ord b => (a -> b) -> [a] -> a
argmin _ []  = error "empty list"
argmin f (x:xs)
    | null xs = x
    | (f x) < (f argminrest) = x
    | otherwise              = argminrest
    where argminrest = argmin f xs

hull :: [Point] -> [Point]
hull pts =
    origin : ((map snd (filter keepPt results)) ++ [origin])
     where origin  = argmin y pts
           oxaxis  = Point (1 + x origin) (y origin)
           oxangle = angle oxaxis origin
           compare_angle pt1 pt2 = compare (oxangle pt1) (oxangle pt2)
           list    = reverse (sortBy compare_angle pts)
           -- origin is the first element of the list
           results = zip (turns (list ++ [origin])) (drop 1 list)
           keepPt (dir, _) = dir /= TurnRight

It works, at least according to the following unit test:

hull [(Point (-3) 1),(Point (-4) 1),(Point (-1) 4),(Point 0 0),(Point 2 2),(Point (-1) 3),(Point (-1) 2),(Point 1 0),(Point 3 (-1)),(Point (-1) (-1))]

which gives

[Point {x = -1.0, y = -1.0},Point {x = -3.0, y = 1.0},Point {x = -1.0, y = 2.0},Point {x = -1.0, y = 3.0},Point {x = -1.0, y = 4.0},Point {x = 2.0, y = 2.0},Point {x = 3.0, y = -1.0},Point {x = -1.0, y = -1.0}]

Is it an idiomatic code? How would you improve it?

Answer 1

Terminology

In your hull function, you use the name origin to indicate the chosen point with the minimum y coordinate. That's confusing, as the term usually means the point (0, 0) in the Cartesian plane. For the rest of this review, I'm going to refer to the chosen point as the pivot instead.

Wrong result

The test result that you quoted in your own question shows that the algorithm is incorrect!

Furthermore, when I ran your code with the suggested test case, I got a totally different result than what you claimed, which makes it hard for me to explain exactly what went wrong. (The code produces results in clockwise order, contrary to your claimed answer which lists them in counterclockwise order. I seriously doubt that you posted the question correctly.) I can, however, point to four issues:

1. (drop 1 list) was responsible for removing (-4, 1) from consideration.
2. A correct traversal involves backtracking whenever a right turn occurs. Therefore, you must filter the candidate points while you traverse the list. You cannot simply traverse the list and keep just the vertices where a left turn occurred, since a left turn that occurs at an interior point doesn't mean that it is a vertex of the convex hull.
3. Reporting points that lie collinearly on the convex hull (but not at a vertex) is a difficult task with this algorithm, since it is hard to determine in advance the order in which to consider such points such that they will be guaranteed to be part of a left-turn-only tour.
4. Degenerate cases are poorly handled. hull of a single point lists that point twice. hull of a triangle can also list the pivot twice.

Implementation

Stuff like this is hard to follow:

angle :: Point -> Point -> Point -> Float
angle a b c = - (bax * bcx + bay * bcy) /
              (sqrt ((bax * bax + bay * bay) * (bcx * bcx + bcy * bcy)))
    -- more efficient than acos for the algorithm and is the same for sorting
    where
      bax = x a - x b
      bay = y a - y b -- coordinates of b->a
      bcx = x c - x b
      bcy = y c - y b -- coordinates of b->c

A remedy is to use a vector package to allow you to think at a more abstract level.

angle a b c :: Vector2 -> Vector2 -> Vector2 -> Double
angle a b c = (ab `vdot` bc) / ((vmag ab) * (vmag bc))
  where
    ab = b - a
    bc = c - b

In the implementation below, I've rewritten the function further, restricting it to only compare a single vector against \$\mathbf{\hat{x}}\$, and renaming it to cosAngleX to be more honest.

You don't really need to define the argmin function. argmin y pts could just be written as minimumBy (comparing y) pts, using support from Data.List.

Suggested solution

import Data.List (minimumBy, sortBy)
import Data.Ord (comparing)

-- https://hackage.haskell.org/package/AC-Vector
import Data.Vector.Class (vdot, vmag)
import Data.Vector.V2 (Vector2(..), v2x, v2y)
import Data.Vector.V3 (Vector3(..), v3z, vcross)

data Direction = TurnLeft
               | TurnRight
               | StraightLine
                 deriving (Show, Eq)

-- Direction of the turn when travelling from a to b to c
turn :: Vector2 -> Vector2 -> Vector2 -> Direction
turn a b c
  | z > 0     = TurnLeft
  | z < 0     = TurnRight
  | otherwise = StraightLine
  where
    ba = a - b
    bc = c - b
    promote23 (Vector2 x y) = Vector3 x y 0
    z = v3z ((promote23 bc) `vcross` (promote23 ba))

-- Cosine of the angle that a vector forms with the X axis.
-- We can sort by this instead of sorting by the angle as long as all vectors
-- have a nonnegative y coordinate (which is the case, since we picked a point
-- with the minimum y coordinate as the pivot).
-- For (Vector2 0 0), the result is Infinity
cosAngleX :: Vector2 -> Double
cosAngleX v
  | norm == 0  = 1/0
  | otherwise  = (v `vdot` xHat) / norm
  where
    xHat = Vector2 1 0
    norm = vmag v

giftWrap :: [Vector2] -> [Vector2]
giftWrap vs = result
  where
    -- Decides whether to keep or reject vertex b.
    -- Vertex b is provisionally kept if a left turn occurs there; the result
    -- is finalized only when a closed left-turn-only tour is completed.
    giftWrap' (a:b:c:vs)
      | vs == [] && dir == TurnLeft  = (True, [b, c])
      | vs == []                     = (False, [c])
      | dir == TurnLeft && final     = (True, (b:gw))
      | dir == TurnLeft              = giftWrap' (a:gw)
      | otherwise                    = (False, (a : c : vs))
      where
        dir = turn a b c
        (final, gw) = giftWrap' (b:c:vs)
    (_, result) = giftWrap' (vs ++ [head vs])

-- Chooses the point with the smallest y coordinate (breaking ties by smallest
-- x coordinate), then presents the points by the angle between the x axis,
-- the pivot, and each point.  The pivot will be the last element of the
-- resulting list.
chain :: [Vector2] -> [Vector2]
chain pts = sortBy (comparing byAngleX) pts
  where
    comparing2 primaryKey secondaryKey a b
      | firstCmp == EQ  = comparing secondaryKey a b
      | otherwise       = firstCmp
      where firstCmp = comparing primaryKey a b
    pivot      = minimumBy (comparing2 v2y v2x) pts
    byAngleX v = -cosAngleX (v - pivot)

-- Convex hull in some counter-clockwise order
hull :: [Vector2] -> [Vector2]
hull (a:[])    = [a]
hull (a:b:[])  = [a, b]
hull vs        = giftWrap $ chain vs

The test case:

*Main> hull [(Vector2 (-3)   1),  (Vector2 (-4) 1),   (Vector2 (-1)  4),
             (Vector2   0    0),  (Vector2   2  2),   (Vector2 (-1)  3),
             (Vector2 (-1)   2),  (Vector2   1  0),   (Vector2   3 (-1)),
             (Vector2 (-1) (-1))]
[Vector2 {v2x = 3.0, v2y = -1.0}, Vector2 {v2x = 2.0, v2y = 2.0},
 Vector2 {v2x = -1.0, v2y = 4.0}, Vector2 {v2x = -4.0, v2y = 1.0},
 Vector2 {v2x = -1.0, v2y = -1.0}]