Convex Hulls: Incremental Algorithm

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

Try it out:

signedArea[{0,0},{1,0},{1,1}]

signedArea[{0,0},{1,0},{1,-1}]

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)

Try them out:

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

leftOn[{1,1},{2,2},{3,3}]

collinear[{1,2},{2,0},{3,-2}]

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

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

How could we avoid this?

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.

findVisible

First, we will write 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_, p_]:=Module[{ltor,rtol,prev,cur,i,j},
prev=leftOf[ Last[hull],First[hull],p ];
For[i=1,i≤Length[hull],i++,
j=Mod[i,Length[hull]]+1;
cur=leftOf[ hull[[i]],hull[[j]],p ];
If[prev==True && cur==False,ltor=i];
If[prev==False && cur==True, rtol=i];
(*Print[{prev,cur}];*)
prev=cur;
];
{ltor,rtol}
]

Try it out:

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

findVisible[testhull,newpt]

The full algorithm

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

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? Here is one way to find out:

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

Your Turn: Time the incremental algorithm for various input sizes. Can you verify that this algorithm has quadratic time complexity?