Convex Hull Algorithms

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

Incremental Algorithm

First, a module to find the part of hull visible to a new point. Assume general position.

findVisible[hull_, pt_]:=Module[{ltor,rtol,prev,cur,i,j},
prev=leftOf[ Last[hull],First[hull],pt ];
For[i=1,i≤Length[hull],i++,
j=Mod[i,Length[hull]]+1;
cur=leftOf[ hull[[i]],hull[[j]],pt ];
If[prev==True && cur==False,ltor=i];
If[prev==False && cur==True, rtol=i];
(*Print[{prev,cur}];*)
prev=cur;
];
{ltor,rtol}
]

Now the incremental algorithm:

incrementalAlg[pts_]:=Module[{sorted,hull,k,end,start},
(* sort the points by x-coordinate *)
sorted=Sort[pts];

(* form hull of the first three points *)
hull=sorted[[;;3]];
If[!leftOf[hull[[1]],hull[[2]],hull[[3]]],hull=Reverse[hull]];

(* incrementally add each vertex to the hull *)
For[k=4,k≤Length[sorted],k++,
{end,start}=findVisible[hull, sorted[[k]]];
hull=If[start < end,
Join[ {sorted[[k]]},hull[[start;;end]] ],
Join[{sorted[[k]]},hull[[start;;]],hull[[;;end]] ] 
]
];

hull
]

Try it out:

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

hull=incrementalAlg[pts]

Draw the points and their convex hull:

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

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. Note that the while loop has been updated since last class to avoid i being equal to curIndex or nextIndex.

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=If[curIndex==1,2,1];  (* ensure nextIndex is different than curIndex *)
For[i=2,i≤Length[pts],i++,
If[i==curIndex||i==nextIndex,Continue[]];
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 point to the hull *)
AppendTo[hull,pts[[curIndex]]]

(* 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]},
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,Return[pts]];
If[Length[pts]==3, (* make sure triangle is counerclockwise *)
If[!leftOf[pts[[1]],pts[[2]],pts[[3]]],Return[Reverse[pts]],Return[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 value *)
]

Try it out:

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

hull=recursiveHull[pts]

Draw the points and their convex hull:

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

Runtime Comparison

Random points in the plane

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

incTimes=Table[ {1000n,Timing[incrementalAlg[pts[[;;1000n]]]][[1]]},{n,1,20}]

gwTimes=Table[ {1000n,Timing[giftWrapping[pts[[;;1000n]]]][[1]]},{n,1,20}]

recurTimes=Table[{1000n, Timing[recursiveHull[pts[[;;1000n]]]][[1]]},{n,1,20}]

ListLinePlot[{incTimes,gwTimes,recurTimes},PlotLegends->{"incremental","gift wrapping","recursive"},AxesLabel->{"number of points","time (seconds)"},PlotLabel->"Runtime of Convex Hull Algorithms"]

Random points in a circle

cpts=RandomPoint[Circle[],{1,10000}][[1]];

incCTimes=Table[ {200n,Timing[incrementalAlg[cpts[[;;200n]]]][[1]]},{n,1,5}]

gwCTimes=Table[ {200n,Timing[giftWrapping[cpts[[;;200n]]]][[1]]},{n,1,4}]

recCTimes=Table[{500n, Timing[recursiveHull[cpts[[;;500n]]]][[1]]},{n,1,20}]

ListLinePlot[{incCTimes,gwCTimes,recCTimes},PlotRange->{0,4},PlotLegends->{"incremental","gift wrapping","recursive"},AxesLabel->{"points","time (seconds)"},PlotLabel->"Runtime of Convex Hull Algorithms"]

A 3D Example

The following example is from Joseph O’Rourke, Computational Geometry in C, Figure 4.7.

verticesToSolid[v_List,f_List]:=Polygon[v[[#]]]&/@f

verts1={{-20,10,5},{-20,20,5},{-5,10,5},{-5,20,5},{-5,10,8},{-5,20,8}};
faces1={{1,2,4,3},{1,2,6,5},{1,3,5},{2,4,6},{3,4,6,5}};

verts2={{10,18,20},{20,18,20},{20,12,20},{10,12,20},{10,18,-10},{20,18,-10},{20,12,-10},{10,12,-10}};
faces2={{1,2,3,4},{5,6,7,8},{1,2,6,5},{3,4,8,7},{1,4,8,5},{2,3,7,6}};

Show[Graphics3D[{{Opacity[0.8],verticesToSolid[verts1,faces1]},{Opacity[0.8],verticesToSolid[verts2,faces2]}}],Axes->True,ViewPoint->{-1,-1,1}]