Convex Hulls: Incremental Algorithm

Math 282 Computational Geometry

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

Try it out:

leftOf[{1,0},{2,2},{0,5}]

leftOf[{1,1},{2,2},{4,0}]

Incremental 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.

Here is the incremental algorithm for finding a convex hull.
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

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

(* start with the first three points, in counterclockwise order *)
hull=sorted[[;;3]];
If[!leftOf[hull[[1]],hull[[2]],hull[[3]]],hull=Reverse[hull]];

(* loop through the remaining points *)
For[k=4,k≤Length[sorted],k++,
(* find which verticies of the current hull are visible from the kth point *)
{end,start}=findVisible[hull, sorted[[k]]];
(* insert the kth point into the hull *)
hull=If[start < end,
Join[ {sorted[[k]]},hull[[start;;end]] ],
Join[{sorted[[k]]},hull[[start;;]],hull[[;;end]] ] 
]
];

hull (* return the hull *)
]

subroutine: findVisible

To complete the incremental algorithm, we need a subroutine findVisible that finds the part of an existing convex hull that is “visible” from a new point outside of the hull. 
Input: hull given as a list of coordinate pairs, in counterclockwise order; point p given as a coordinate pair outside of the hull
Output: pair {a, b} of indexes of vertices in hull, such that vertices a, b, and everything in between are visible from point p

findVisible[hull_, pt_]:=Module[{end,start,prev,next,i,j},
(* first consider the edge that connects the last vertex of the hull with the first vertex *)
prev=leftOf[ hull[[-1]],hull[[1]],pt ];

(* now loop over all indexes of points in the hull *)
For[i=1,i≤Length[hull],i++,
(* get the next index after i, wrapping around if necessary *)
j=If[i<Length[hull],i+1,1];

(* do the leftOf query *)
next=leftOf[ hull[[i]],hull[[j]],pt ];

(* determine if we hae a transition at vertex i *)
If[prev==True && next==False,end=i];
If[prev==False && next==True, start=i];

prev=next; (* shift next to prev *)
];
{end,start} (* return pair of indexes *)
]

Here is a test case.

testhull = {{1,0},{3,0},{4,3},{2,5},{0,3}};
newpt={6,2};
Graphics[{{LightYellow,EdgeForm[{Black,Thin}],Polygon[testhull]},{PointSize->Medium,Blue,Point[testhull]},{PointSize->Medium,Red,Point[newpt]}},Frame->True,AspectRatio->1]

If findVisible is working correctly, the following line will return {2, 4}.

findVisible[testhull,newpt]

Testing the Incremental 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=incrementalAlg[pts]

Draw it:

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

Investigate Runtime

How quickly does the incremental algorithm run? We can measure the runtime of the incremental algorithm on 1000 points:

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

The output above is the runtime in seconds.

Make a plot of runtime over various numbers of points

First we will make a big list of points:

numPoints = 25000;
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

incRuntime[n_,pts_]:=Timing[ incrementalAlg[ pointsList[[;;n]] ] ][[1]]

Try it out for 500 points:

incRuntime[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,incRuntime[n,pts]},{n,1000,numPoints,1000}]

Next we will plot the values:

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

We can also fit a quadratic function to the runtimes:

f=Fit[runtimes,{1,x,x^2},x]