Why doesn't this implementation of Jarvis' March ("Gift wrapping algorithm") work?

Question

I'm trying to implement Jarvis' algorithm for finding the convex hull of a set of points, but for some reason it doesn't work. This is my implementation:

procedure TPointList.ConvexHull(aHull : TPointList); //Return the convex hull of a set of 2D points
var
  vPointOnHull  : TPoint2D;
  vEndpoint     : TPoint2D;
  I             : integer;
begin
  aHull.Clear;
  if Count < 3 then exit;

  vPointOnHull := Self.LeftMostPoint;
  repeat
    aHull.Add(vPointOnHull);
    vEndpoint := Self.Point[0];

    for I := 1 to Self.Count-1 do
      if Orientation(vPointOnHull,vEndpoint,Self.Point[I]) = LeftHandSide then
        vEndpoint := Self.Point[I];

    vPointOnHull := vEndpoint;
  until vEndpoint = aHull.Point[0];
end;

• TPointList is a simple list of points.
• Orientation is a function from Arash Partow's library "FastGEO"
• The implementation is lifted more or less directly from the Wikipedia article on the algorithm

What happens is that the method starts adding the same point to aHull over and over. In one test case I send in the points (200;200) (300;100) (200;50) and (100;100), and the algorithm starts by adding (100;100) to aHull which is correct, but then it starts adding (200;200) over and over again.

Obviously I've done something wrong in my implementation, but for the life of me I can't see what.

UPDATE:

Jonathan Dursi put me on the right track. This line

if Orientation(vPointOnHull,vEndpoint,Self.Point[I]) = LeftHandSide then    

should be replaced with this

if (vPointOnHull = vEndpoint) or (Orientation(vPointOnHull,vEndpoint,Self.Point[I]) = LeftHandSide) then

Works like a charm :-)

Answer 1

The loop adds using this line of code:

aHull.Add(vPointOnHull);

vPointOnHull is only assigned in these lines:

vPointOnHull := Self.LeftMostPoint;
vPointOnHull := vEndpoint;

You already explained that LeftMostPoint is added correctly, so the repeat must come from vEndPoint, which is assigned in these lines:

vEndpoint := Self.Point[0];
vEndpoint := Self.Point[I];

So I guess the last assignment (which is in the below if statement), is never reached.

  if Orientation(vPointOnHull,vEndpoint,Self.Point[I]) = LeftHandSide then
    vEndpoint := Self.Point[I];

--jeroen

Answer 2

It's probably not a conicidence that (200;200) is point 0.

It looks like you're not excluding the current point (vPointOnHull) from being the end point (vEndPoint), and your implementation of Orientation doesn't reject that case; presumably it returns LHS if the cross-product is positive, and if vPointOnHull == vEndPoint, the cross product is zero, so never LHS. So nothing ever replaces Point 0 once Point 0 is selected, et voila.

You could modify Orientation to return "Degenerate" or something in that case, and also reject the point, or you could exclude the current point from ever being the end point. Note that you don't want to do the obvious thing, filter out current CH points from the point set while marching through, because you need to find that the end point is the first point to close the loop.

Update: Looking around a bit at the FastGEO stuff, probably updating Orientation isn't the way to go (although a bit more thought should go into the colinear points case in this algorithm; if there are collinear points on the hull, you really want the closest one first, so you'd like an else if Orientation = Collinear then.. update vEndpoint if new point is closer clause after that if statement).

Easiest might just be to add a couple lines keeping track of the current indicies so you can easily test for equality: something a bit like

iPointOnHull := Self.IndexOfLeftMostPoint;
vPointOnHull := Self.LeftMostPoint
...
vEndpoint := Self.Point[0];
iEndPoint := 0;
if (iPointOnHull = 0) then 
begin
    vEndPoint := Self.Point[1];
    iEndPoint := 1;
end
...
vPointOnHull := vEndPoint;
iPointOnHull := iEndPoint;