Determining Convex Hulls with Quickhull Algorithm

Taran Besant and Spencer Eanes

The basic idea behind this project is to find the set of points in a 2d plane that contains all other points in a convex shape in the least number of points.
The algorithm we choose for this is the quickhull algorithm. More information about this algorithm can be found on wikipedia, and in our presentation slides.
Motivation for this project can be found in submitted presentation slides.

XtremeP

Calculate the farthest left and farthest right most points in a given data set. in the cast that two points share the same extreme x value, the function returns the one closer to the midpoint of the data set, for simplicity.
   @arg hull - a set of points to create a convex hull around

xtremeP[hull_]:=Module[{x,y},
x=hull[[1]];
For[i=1,i< Length[hull],i++,
If[hull[[i+1]][[1]]>x[[1]], x=hull[[i+1]],If[hull[[i+1]][[1]]==x[[1]],
If[Abs[hull[[i+1]][[2]]-50]>Abs[x[[1]]-50], ,x=hull[[i+1]]
];
];
];
];

y=hull[[1]];
For[i=1,i< Length[hull],i++,
If[hull[[i+1]][[1]]<y[[1]], y=hull[[i+1]],If[hull[[i+1]][[1]]==y[[1]],
If[Abs[hull[[i+1]][[2]]-50]<Abs[y[[2]]-50],y=hull[[i+1]],
];
];
];
];
Return[{x,y}];
];

HullPartition

A module that partitions the remaining points not in the convex hull into two groups : those above a line segment and those below based on y value. This runs into a problem with vertical lines, because the slope is returned as a complex infinity which we can' t work with in the future (See Hull Partition I).
   @arg allpts - the set of points to be partitioned.
   @arg x - the first point used to create the partitioning line.
   @arg y - the second point used to create the partitioning line.

hullPartition[allpts_,x_,y_]:=Module[{m,b,Coordy,tempList,listAbove,listBelow},
m=(y[[2]]-x[[2]])/(y[[1]]-x[[1]]);
b=x[[2]]-m*x[[1]];
listAbove={};
listBelow={};
For[i=0,i≤ 100,i++,
Coordy=m*i+b;
tempList={};
For[l=1,l≤ Length[allpts],l++,
If[allpts[[l]][[1]]==i,AppendTo[tempList,allpts[[l]]],
];
];
For[k=1,k≤ Length[tempList],k++,
If[tempList[[k]][[2]]>Coordy, AppendTo[listAbove,tempList[[k]]],
If[tempList[[k]][[2]]==Coordy,,AppendTo[listBelow,tempList[[k]]]
];
];
];
];
Return[{listAbove,listBelow}];
];

PtTrans

This takes a 2d point, pt, and and two more points, P and Q, which are endpoints of a line.  
It finds the equation of the line through P and Q using the fit function, which returns the RHS of the equation of a line in the form y = mx + b. 
It then find a rotation matrix using Mathematica’s built in function by taking the ArcTangent of the slope of the line between P and Q.
	Note the adjustment here, and taking the inverse outside. This is because a standard rotation matrix rotates counter-clockwise, and we want to rotate clockwise.
We then adjust the point down by the y-intercept of the line, and apply the rotation matrix to the point.

This function also has two special cases, if the line determined by P and Q is horizontal or vertical.
In these cases the result of the Fit function does not play nicely with the rest of the code, so we create these special cases.
If the line is vertical, we rotate 3π/2 radians. If they are horizontal, we only adjust down by the y intercept.

PtTrans[pt_,{P_,Q_}]:=Module[{point = pt,f1,yint,A,x},
If[P[[1]]-Q[[1]]==0,
(*Print["Line is vertical"];*)
A=Inverse@RotationMatrix[π/2];
Return[A.point];
];
If[P[[2]]-Q[[2]]==0,
(*Print["Line is horizontal"];*)
yint = Q[[2]];
point-={0,yint};
Return[point];
];
f1=Fit[{P,Q},{1,x},x];
yint=f1/.x->0;
point-={0,yint};
(*Print[point];*)
A=Inverse@RotationMatrix[(ArcTan[1,f1[[2]]/.x->1])];

(*Print[A.point];*)

Return[A.point];
]

HullPartition1

This is a more general version of HullPartition that works in all cases.
This function implements the PtTrans method on all the points in the passed to it in the first argument, pts.
It uses the map function to apply PtTrans to all  the points and stores their output, where the line determined by P and Q is treated as the x-axis in transPts.
It constructs a list of two sublists, the first of which is the indices of points above, those with transformed y-values greater than 0, and those below, with y-value less than 0.
It returns the this sorted set of two lists.

Note that if a point is exactly on the line, it is ignored. This is intended behavior, because if a point is on the line we don’t want to include it.

hullPartition1[pts_,{P_,Q_}]:=Module[{transPts,topPtInd={},botPtInd={},topPts,botPts},
transPts = Map[PtTrans[#,{P,Q}]&,pts];
(*Print["transPts: ",transPts];*)
Do[
If[transPts[[ix]][[2]]>0,
AppendTo[topPtInd,ix];
];
If[transPts[[ix]][[2]]<0,
AppendTo[botPtInd,ix];
],
{ix,1,Length[pts]}
];
(*indices of all top and bottom points found*)
(*if a point is on the line it is ignored*)
topPts = pts[[topPtInd]];
botPts = pts[[botPtInd]];

Return[{topPts,botPts}];
]

Distance and Sorting

PointLineDistance takes a a point pt, and a tuple of points, s1 and s2. It uses mathematica’s built in Euclidean distance and projection functions to find the shortest distance from the point to the line.

PointLineDistance[pt_,{s1_,s2_}]:=With[{tp=s1-pt},EuclideanDistance[tp,Projection[tp,s2-s1]]];

This function implements PointLineDistance on a set of points, pts, and returns the index of the point furthest from the line.

Furthest[pts_,{s1_,s2_}]:=Module[{dists},
dists=Map[PointLineDistance[#,{s1,s2}]&,pts];
Return[Position[dists,Max[dists]]];
];

This functions sorts the points passed to it by their x-value.

SortPts[pts_]:=Module[{},
Return@SortBy[pts,First];
]

Find Hull

This is a recursive function that implements hullPartition1 to find the complex hull using divide and conquer.
This is the main implementation of the quickhull algorithm.
It’s input are:
	S - the list of all points in the hull.
	P, Q - the two most extreme x and y points.
	cHull - all the points in the hull so far, this is mainly used so that recursive calls can track what has been found so far.
	tb - this is an indicator that should be 1 or 2, which tells if we are finding the hull for the top or bottom half. It is used to index the return from hullPartition1.
If S is empty, we do nothing.
We then find the furthest point in S from the line determined by P and Q and add it to the convexHull.
Then, we find all the points that are in S but not within the polygon determined by the points of the convexHull. 
In order for the previous step to work, the set of points in the convex hull needs to be sorted into the correct order so that the polygon is correct, so we first use the SortPts function.
We only add points if they are not region members of the polygon, using mathematica’s built in function.
Then, if there are points in S that are not contained in the polygon, we partition them into those above and below the two new line segments created with the new furthest point and the two previous line segments.
We then recursively run the FindHull algorithm on the correct half of of the partitioned points. 
This will terminate once there are no points outside the polygon, and return the set of points found.

FindHull[S_,{P_,Q_},cHull_,tb_:1]:=Module[{convexHull = cHull, furthestIndex,maxDist=0,S1,S2,newPts},
If[S≠{},
(*Print[S];*)
furthestIndex=Furthest[S,{P,Q}][[1]][[1]];
(*Print[furthestIndex];*)
(*Print[S[[furthestIndex]]];*)
AppendTo[convexHull,S[[furthestIndex]]];
(*Print[convexHull];*)
(*Print["convexHull: ", convexHull];*)
newPts={};
convexHull=SortPts[convexHull];
Do[
If[!RegionMember[Polygon[convexHull],S[[ix]]],
AppendTo[newPts,S[[ix]]];
],
{ix,1,Length[S]}
];

(*Print["newPts: ",newPts];*)
If[newPts≠{},
S1=hullPartition1[newPts,{S[[furthestIndex]],P}][[tb]];
S2=hullPartition1[newPts,{Q,S[[furthestIndex]]}][[tb]];

convexHull =FindHull[S1,{S[[furthestIndex]],P},convexHull,tb];
convexHull = FindHull[S2,{Q,S[[furthestIndex]]},convexHull,tb];
];
];
Return@convexHull
]

Quickhull Driver

This function acts like a C++ driver program. It calls all constituent parts of the quickhull function:
	find the extreme x points.
	partition the set on the line segment determined by these points.
	Fund the upper and lower hull using the FindHull function, and finally
	Return these two lists.

We include the delete case clauses because of a simple error that can be caused. If a point had exactly the same x-coordinate as one of the extreme points, the algorithm would sometimes sort them in the incorrect order causing the correct hull to display as if it was incorrect. By doing deletecases and then prepending and appending the first and last hull points, we make sure that everything is in the order that we want.

QuickHull[pts_]:=Module[{xpts,part,hullPtsTop,hullPtsBot},
xpts =xtremeP[pts];
xpts=SortPts[xpts];
part = hullPartition1[pts,xpts];

hullPtsTop=FindHull[part[[1]],xpts,xpts,1];
hullPtsBot = FindHull[part[[2]],xpts,xpts,2];

hullPtsTop=DeleteCases[hullPtsTop,xpts[[1]]];
hullPtsTop = DeleteCases[hullPtsTop,xpts[[2]]];
PrependTo[hullPtsTop,xpts[[1]]];
AppendTo[hullPtsTop,xpts[[2]]];

hullPtsBot=DeleteCases[hullPtsBot,xpts[[1]]];
hullPtsBot = DeleteCases[hullPtsBot,xpts[[2]]];
PrependTo[hullPtsBot,xpts[[1]]];
AppendTo[hullPtsBot,xpts[[2]]];

Return[{hullPtsTop,hullPtsBot}];
]

Output

In our output, a random selection of points is generated in a circle of radius 10 with center at the origin. This is because with other large amounts of randomly generated points the hull was very boxlike, so this gives a more interesting visual output. The QuickHull module returns two lists of points : the points making up the top half of the convex polygon and the points making up to bottom half. These are shown individually, and then together.

Interestingly, we notice that as the number of points is increased, the hull more and more precisely approximates the exact boundaries of the region specified in the point selection. Cool!

n=100; (*total number of points*)
pts = RandomPoint[Polygon[{{0,0},{0,1},{1,1},{1,0}}],n];
(*an alternative input that is in a circle*)
(*pts = RandomPoint[Disk[{0,0},1],n];*)
Hull=QuickHull[pts]

a =ListPlot[pts,PlotRange->All];
b=ListLinePlot[Hull[[1]]];
c=ListLinePlot[Hull[[2]]];
Show[a,b]
Show[a,c]
Show[a,b,c]

Note that mathematica also has a built in function that will do this. Mathematica’s function also works in 3d. Though this doesn’t have an option to print specific points in the hull, we can see on close observation that the vertices appear to be in the same locations.

ConvexHullMesh[pts]

As we can see, mathematica’s implementation is several orders of magnitude faster than ours... Oh well!

AbsoluteTiming[QuickHull[pts]][[1]]
AbsoluteTiming[ConvexHullMesh[pts]][[1]]