Convex Hull Algorithms

Math 282 Computational Geometry

“Left Of” Queries

For points a, b, c, in the plane, the signed area of triangle abc is half the cross product (b-a)×(c-a).

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

We can use this to define Boolean functions (i.e., functions that return true or false) to determine whether point c is to the left of directed segment ab or collinear to ab.

leftOf[a_,b_,c_]:=(signedArea[a,b,c] > 0)
leftOn[a_,b_,c_]:=(signedArea[a,b,c] >= 0)
collinear[a_,b_,c_]:=(signedArea[a,b,c]==0)

Note that we might sometimes get the wrong answer due to round-off error.

collinear[{1.45,2},{2,0},{2.55,-2}]

Gift Wrapping Algorithm

First, a module to find the point with the smallest y-coordinate. If there is a tie, return the rightmost point. This module returns the index of the point, not the point itself.

findLowest[pts_]:=Module[{n=Length[pts],i,low},
low=1;
Do[
If[(pts[[i,2]]<pts[[low,2]] )|| (pts[[i,2]]==pts[[low,2]] && pts[[i,1]]>pts[[low,1]]),
low=i],
{i,2,n}];
low
]

Now a module to implement the gift wrapping algorithm.

giftWrapping[pts_]:=Module[{startIndex,hull,curIndex,nextIndex,i},
(* start at the "lowest" point *)
startIndex = findLowest[pts];
hull={pts[[startIndex]]};
curIndex=startIndex;

While[True,
(* search for the point that makes the largest angle *)
nextIndex=1;
For[i=2,i≤Length[pts],i++,
If[ leftOf[pts[[curIndex]],pts[[i]],pts[[nextIndex]] ],
nextIndex=i];
];
curIndex=nextIndex;

(* quit if we are back to the starting point *)
If[curIndex==startIndex, Break[]];

(* otherwise, append the next point to the hull *)
AppendTo[hull,pts[[nextIndex]]]

(* TESTING: print the partial hull 
Print[Graphics[{{Purple,Line[hull]},{Blue,PointSize[0.02],Point[pts]}},Frame->True,ImageSize->Small]]*)
];
hull
]

Try it out:

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

hull=giftWrapping[pts]

Draw the points and their convex hull:

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

Divide and Conquer Algorithm

First a module to determine whether a polygon has a left turn at vertex index i. Assume vertices are given in counterclockwise order around the polygon. This module uses modular arithmetic to wrap around the polygon.

leftTurn[poly_,i_]:=Module[{n=Length[poly]},
(*Print[{"leftTurn",poly,i}];*)
leftOf[poly[[Mod[i-2,n]+1]],poly[[Mod[i-1,n]+1]],poly[[Mod[i,n]+1]]]
]

Next a module to remove right turns from a polygon formed from the merge of two hulls. Starts checking for right turns at indexes i and i+1. This module requires i>0.

removeRightTurns[poly_,index_]:=Module[{p=poly,i=index},
While[!leftTurn[p,i] || !leftTurn[p,i+1],
If[!leftTurn[p,i+1],p=Delete[p,i+1]];
If[!leftTurn[p,i],p=Delete[p,i];i--]
];
p
]

Now a module to implement the recursive hull algorithm.

recursiveHull[pts_]:=Module[{sorted,mid,left,right,a,b,c,hull},
(* base case *)
If[Length[pts]≤2,!! complete this code  !!  ];
If[Length[pts]==3, !! complete this code  !!  ];

(* 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 value *)
]

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]