I recently obtained the "Dutch" computational geometry book to supplement the other computational geometry book that I have. Instead of doing what I normally do which is skip directly to the topic that I'm interested in, I've resolved to work through the entire book systematically implementing each of the algorithms in F# or some other functional language. The purpose is two fold; I know almost nothing about the topic and it's hard to find anything about computational geometry implementations in a functional language. Who knows, someone may actually find my floundering helpful.
Chapter one is all about convex hulls and after doing a little bit of reading I determined that the "Jarvis March" algorithm looked like a good place to start. I've only included core part of the algorithm here in the blog post. The entire source including the visualization can be found here.
The algorithm is pretty simple I start first by sorting the point list by the Y coordinate and use the minimum as the starting point for my hull. I start working through the list trying to find the left most point in relation to the previous point added to the hull; this is the fold_left bit of the _jarvis function. Once I've found the left most point I append it to the head of the hull list and make a recursive call to _jarvis. When I reach the point with the maximum Y value I know that I have successfully built the left side of the hull. I then repeat the entire process but this time instead of searching for the left most point I search for the right most point building up the right side of the hull. When done I append these two lists together to form the complete hull.
type Point = {
x :double;
y :double;
} with
static member ( = ) a b =
a.x = b.x && a.y = b.y;
static member ( <> ) a b =
not(a = b);
end;;
type Edge = {
orig:Point;
dest:Point;
} with
static member ( = ) a b =
a.orig = b.orig && a.dest = b.dest;
static member ( <> ) a b =
not(a = b);
end;;
let sort point_list =
List.sort( fun p1 p2 ->
match Pervasives.compare p1.y p2.y with
| 0 -> Pervasives.compare p1.x p2.x;
| _ as v -> v;
) point_list;;
// Input: three points P0, P1, and P2
// Return: >0 for P2 left of the line through P0 and P1
// =0 for P2 on the line
// <0 for P2 right of the line
let testSide p0 p1 p2 =
(p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);;
let isLeft p0 p1 p2 =
(testSide p0 p1 p2) > 0.;;
let isRight p0 p1 p2 =
(testSide p0 p1 p2) < 0.;;
let rec _jarvis cmp max_point (hull_list:Point list) (input_list:Point list) =
match input_list with
| p::tl ->
let last = List.hd hull_list in
if last = max_point then hull_list
else (
let best = List.fold_left(fun p1 p2 ->
if (cmp last p1 p2) then p2
else p1;
) p input_list in
_jarvis cmp max_point (best::hull_list) tl
)
| [] -> hull_list;;
let jarvis point_list =
if List.length point_list < 3 then failwith "3 points are required for a polygon";
let point_list = sort point_list in
let miny = point_list |> List.hd in
let maxy = point_list |> List.rev |> List.hd in
let left = _jarvis (isLeft) maxy [miny] point_list in
let right = _jarvis (isRight) maxy [miny] point_list in
right @ left |> List.rev;;
