Convex Hulls: Divide and Conquer Algorithm

MATH 261 Computational Geometry

Preliminaries: “Left Of” Queries

For points a, b, c, in the plane, the signed area of parallelogram determined by (OverVector[ab]) and (OverVector[ac]) is the third component of the cross product (b-a)×(c-a). We can use this to define a Boolean function (i.e., function that returns true or false) to determine whether point c is to the left of directed segment (OverVector[ab]) or collinear to (OverVector[ab]).

signedArea[a_,b_,c_]:=((b[[1]]-a[[1]])(c[[2]]-a[[2]])-(b[[2]]-a[[2]])(c[[1]]-a[[1]]))
leftOf[a_,b_,c_]:=(signedArea[a,b,c] > 0)

Divide and Conquer Algorithm

First, a module to determine whether a polygon has a left turn at vertex index i. Returns True if there is a left turn at vertex i, and False otherwise. This module assumes vertices are given in counterclockwise order around the polygon and uses modular arithmetic to wrap around the polygon.

leftTurn[poly_,i_]:=Module[{n=Length[poly]},
leftOf[poly[[Mod[i-2,n]+1]],poly[[Mod[i-1,n]+1]],poly[[Mod[i,n]+1]]]
]

Here is the main module that implements the divide and conquer algorithm. This module calls itself recursively. It also calls a module removeRightTurns to remove right turns after each merge of two partial hulls.

recursiveHull[pts_]:=Module[{sorted,mid,left,right,a,b,c,hull},
(* base cases *)
(* make sure points are in counterclockwise order *)
If[Length[pts]≤2,Return[pts]  ];
If[Length[pts]==3, If[leftTurn[pts,2],Return[pts],Return[Reverse[pts]]]  ];

(* sort the points by x-coordinate *)
sorted=Sort[pts];

(* partition the points into two halves, and find the convex hull of each *)
mid=Floor[Length[sorted]/2];
left=recursiveHull[ sorted[[;;mid]] ];
right=recursiveHull[ sorted[[mid+1;;]] ];

(* merge the two hulls *)
a=Ordering[left,-1][[1]];  (* index of rightmost vertex in left *)
b=Ordering[right,1][[1]];  (* index of leftmost vertex in right *)
hull=Join[ left[[;;a]],right[[b;;]],right[[;;b]],left[[a;;]] ];
hull=removeRightTurns[hull,a]; (* eliminate right turns on botom *)
c=Length[hull]-(Length[left]-a+1);
hull=removeRightTurns[hull,c]; (* eliminate right turns on top *)

(* TESTING: print the partial hull *)
Print[Graphics[{{LightYellow,EdgeForm[Black],Polygon[hull]},{Blue,PointSize[0.02],Point[pts]}},Frame->True,PlotRange->{{0,10},{0,10}},ImageSize->Small]];

hull  (* return the hull *)
]

Finally, we need a module to remove right turns from a polygon formed from the merge of two hulls. This module should start by checking for right turns at indexes i and i+1.

removeRightTurns[poly_,index_]:=Module[{p=poly,i=index},
While[!leftTurn[p,i]|| !leftTurn[p,i+1],
(* remove a vertex *)
If[!leftTurn[p,i],
p=Delete[p,i];i=i-1,
(* else *)
p=Delete[p,i+1]
]
];

p   (* return the new hull *)
]

Your Task: Fill in the missing code above.

Try it out:

pts=RandomReal[{0,10},{20,2}]

Graphics[{Blue,PointSize[0.02],Point[pts]},Frame->True]

hull=recursiveHull[pts]

Graphics[{{LightYellow,EdgeForm[Black],Polygon[hull]},{Blue,PointSize[0.02],Point[pts]}},Frame->True]

Investigate Runtime

How quickly does the divide and conquer algorithm run? Here is one way to find out:

pts=RandomReal[{0,10},{1000,2}];
Timing[recursiveHull[pts]]

Your Task: Time the divide and conquer algorithm for various input sizes.