Convex Hulls: Gift Wrapping Algorithm

MATH 261 Computational Geometry

Preliminaries

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). This is computed as follows:

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

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]).

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

Gift Wrapping Algorithm

For now, assume that points are in general position (i.e., all x-coordinates are unique and no three points are collinear), and that numerical errors don’t happen.

findLowest

First, we need a subroutine that returns a point guaranteed to be on the convex hull. We choose to find the point with the smallest y-coordinate. If multiple points have the smallest y-coordinate, then we choose the rightmost of them. 
Input: a list of points, given as ordered pairs
Output: the index of the (rightmost) lowest point in the list

findLowest[pts_]:=Module[{n=Length[pts],low},
low=1; (* start with the first point *) 
Do[ (* loop over the other points; if we find a new lowest point, then update low *)
If[(pts[[i,2]]<pts[[low,2]] )||(pts[[i,2]]==pts[[low,2]] && pts[[i,1]]>pts[[low,1]]),
low=i],
{i,2,n}];
(* return the index of the lowest point *) 
low
]

Try it out:

findLowest[{{1,0},{2,3},{4,-1},{5,-1},{3,3}}]

The gift wrapping algorithm

The following Module implements the gift wrapping algorithm.
Input: a list of points, given as ordered pairs
Output: convex hull of the points, given as a list of ordered pairs, in counterclockwise order

Your Task: Complete the code below.

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 is most to the right, from the perspective of pts[[curIndex]]. 
Let nextIndex be the index of this point. *)
nextIndex=If[curIndex==1,2,1];

For[i=2,i<=Length[pts],i++,
If[i==nextIndex||i==curIndex,Continue[]];
If[leftOf[pts[[curIndex]],pts[[i]],pts[[nextIndex]]],nextIndex=i]
];

(* quit if we are back to the starting point *)
If[nextIndex==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]];*)

(* advance to the next index *)
curIndex=nextIndex
];
hull
]

Try it out:

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

giftWrapping[pts]

Compute convex hulls using the gift wrapping algorithm

Create a set of random points:

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

Draw the points, just to see what we are dealing with:

Graphics[{PointSize->Medium, Blue, Point[pts]},Frame->True,AspectRatio->1]

Now find the convex hull of the points:

hull=giftWrapping[pts]

Draw it:

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

Investigate Runtime

Sample code to generate some random points and time the gift wrapping algorithm.

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

Your Task: Time the gift wrapping algorithm for various input sizes.

First we will make a big list of points:

numPoints = 10000;
pointsList=RandomReal[{0,10},{numPoints,2}];

Now we will define a little function that takes the first n points and measures the time necessary to compute their convex hull

giftwrapRuntime[n_,pts_]:=Timing[ giftWrapping[ pointsList[[;;n]] ] ][[1]]

Try it out for 500 points:

giftwrapRuntime[500,pts]

Now we will make a list of lots of pairs (n, t), where t is the runtime required to compute the convex hull of n points. We will use values of n that are multiples of 1000.

runtimes=Table[{n,giftwrapRuntime[n,pts]},{n,1000,numPoints,1000}]

Now we will plot the runtime values:

ListPlot[runtimes,AxesLabel->{"number of points","runtime (seconds)"}]

Points on a circle

numPoints=1000;
circPts=RandomPoint[Circle[],{1,numPoints}][[1]];
Graphics[{PointSize->Medium, Blue, Point[circPts]},Frame->True,AspectRatio->1]

giftwrapCircleRuntime[n_,pts_]:=Timing[ giftWrapping[ circPts[[;;n]] ] ][[1]]
circRuntimes=Table[{n,giftwrapCircleRuntime[n,pts]},{n,100,numPoints,100}]

ListPlot[circRuntimes,AxesLabel->{"number of points","runtime (seconds)"}]