# 2D Convex hull exercise

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:

1. 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.
2. Write a function that calculates the turn made by three 2D points and returns a Direction.
3. 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.
4. 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?

## 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)

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

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.
| 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

-- 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}]

• Dec 24 '14 at 23:08