An essential topic in computational geometry is convex hulls, which have found use in many diverse applications: statistics, soil sciences, medical imaging, character recognition, parts inspection, and robotics systems, to name a few. Additionally, finding the convex hull is the first step in the solution to a number of other problems in computational geometry.

The Convex Hull

Convexity may be an abstract mathematical concept, but it is fairly easy to recognize if an object has it or not. A shape is convex if the line connecting any two points in the shape is itself completely contained in the shape. Intuitively, we realize that convex shapes are basically ball-like; they have no indentations or holes.

The convex hull of a set of points is the smallest convex shape that contains all the points. Again, intuitively, we realize that the convex hull must be a polygon (a shape with straight edges), and that the vertices of the polygon must be points in the set.

In fact, you can easily visualize the convex hull of a set of points by doing the following thought experiment: Let each point be represented by a nail, stretch a rubber band around all the nails, and then let go. The shape of the rubber band is the convex hull of the original set of points.

Cross Products Defined

The Swiss-army knife of computational geometry is arguably the cross product, a computation borrowed from the study of vectors. Vectors are fundamental to the disciplines of physics, and you will have encountered them in linear algebra and calculus courses. In this article, I will use this versatile tool in two different ways as well as provide examples of its many uses.

A vector has direction as well as magnitude (length). For example, consider the two vectors P₁ and P₂ in Figure 1. P₁ goes from the origin to the point (x₁,y₁) and P₂ goes from the origin to the point (x₂,y₂). The sum P₁+P₂=P₃ is another vector that goes from the origin to the point (x₁+x₂,y₁+y₂).

The cross product P_1XP2 can then be interpreted as the signed area of the parallelogram formed by the origin, P₁, P₂, and P₃ (for more details, see the accompanying text box entitled "Finding the Area of a Triangle"). The value of the cross product is given by the formula in Example 1(a). Finally, if the point common to P₁ and P₂ is P₀(x₀,y₀) instead of the origin, you simply translate the origin to P₀. In this translation, the formula for the cross product becomes as in Example 1(b).

Cross Products Applied

How can you use the cross product to solve problems in computational geometry? Recall that the cross product is the area of a parallelogram, which is twice the area of the triangle formed by the three points P₀, P₁, and P₂. Thus, with just five subtracts and two multiplies you can compute the area of any triangle given the coordinates of its three vertices.

An extremely useful observation is that P₁XP₂ will be positive whenever P₁ is clockwise from P₂, and negative whenever P₁ is counterclockwise from P₂. (Since the cross product is an area, if the three points P₀, P₁, and P₂ are collinear you can anticipate that the cross product will be 0.)

Therefore, using just the sign of the cross product, you can tell if the point P₂ is on the left or right side of the line going from P₀ to P₁.

As an example, using just this orientation information provided by the cross product, you can easily determine if a point P(x,y) is inside or outside of the triangle formed by the three points P₀, P₁, and P₂.

Here's how: As you go clockwise around the triangle, if the point P is inside, it will always be to the right of the edge you are on. If it is not inside, the point P will be on the left of some edge but to the right of another.

Thus, the point P is inside the triangle only if the three cross products (P₁-P₀)X (P-P₀), (P₂-P₁)X(P-P₁), and (P₀-P₂)X(P-P₂) all have the same sign. A moment's reflection will convince you that this result is correct regardless of the direction you go around the triangle. (If the point P is on the perimeter of the triangle, one of the cross products will be 0; it is your call as to whether such a point should be classified as inside, outside, or treated in some other special way.)

Graham's Scan

Because of its fundamental nature, the problem of finding the convex hull of a set of points has received much attention, and as a result several efficient algorithms now exist.

The problem of sorting N numbers is linear-time transformable to finding the convex hull of N points, and therefore finding the convex hull must require at least O(N log N) time to compute the solution.

Graham's Scan, the algorithm presented here, achieves this lower bound. Additionally, it is known to be optimal in the worst-case sense. Jarvis's March is another popular algorithm that is asymptotically faster than Graham's Scan, but its worst-case performance is O(N²). For this reason Graham's Scan may be the better choice for a general-purpose convex hull algorithm.

I will first describe the algorithm for Graham's Scan by giving an overview, ignoring any nasty details and special cases that might get in the way. I will then use a small example to expose the heart of the algorithm, and finally the code and essential details will be presented.

Algorithm Overview

Graham's Scan begins by locating the lowest point; call it the pivot point. Since none of the other points are below the pivot, sweeping in the counterclockwise direction from 0 to 180 degrees you know you will encounter them all.

As each point is encountered, it is added to a candidate convex hull list. The newest point may reveal that one or more points added earlier cannot be part of the convex hull, and so they are removed from the list. This step is the heart of the algorithm.

Notice that you always turn to the left as you move counterclockwise from one point to the next on the convex hull. Thus, if you ever turn right, the point you turned right from cannot be part of the convex hull; it must be on the inside.

If this happens, remove the offending point from the list, then check to see if you turn left or right to get from the new head of the list to the new point. Repeat this step until the left-turn-only property is restored.

After all points have been visited, the list will contain all the points that lie on the convex hull of the original set and no others.

Algorithm Example

To illustrate the Graham's Scan algorithm, consider the set of points in Figure 3. The pivot point is labeled P₀, and the other points are labeled in increasing order by the angle they make with P₀ and the positive x-axis.

Both P₀ and P₁ must be on the convex hull, and that P₂ cannot be right of the line from P₀ to P₁. This is important: The algorithm begins with two points guaranteed to be on the convex hull and with the left-turn-only property intact.

So, begin with points P₀, P₁, and P₂ in the convex hull candidate list, and consider point P₃. Since you turn left from P₂ to get to P₃ it is added to the list; the same goes for P₄.

However, P₅ is on the right side of the line from P₃ to P₄, so you now know that P₃ cannot be part of the convex hull after all. After removing P₄ from the list you again consider P₅. This time P₅ is on the left side of the top two points on the list (P₂ and P₃), so it is added to the list.

Points P₆, P₇, and P₈ all get added to our list, since they do not cause us to violate the left-turn-only rule.

When you get to point P₉, however, you are led to remove P₈ from the list. After doing so, you find that P₉ is still to the right of the top two points on the list (P₆ and P₇), and so P₇ must be removed also. This restores the left-turn-only property to our list, so P₉ can now be added.

Finally, point P₁₀ is added to the list without incident and you are done. The dashed lines in Figure 3 show the points that were added to the convex hull candidate list that later had to be deleted.

Algorithm Details

Assume you are given as input an array S of N points, and must produce as output the set of points of S that are on the convex hull of S in clockwise order.

With the example fresh in mind, Listing One presents the Graham Scan algorithm proper. Assume that the pivot point has already been located and placed in S[0], and that the rest of the array has been sorted in order of increasing polar angle.

It turns out to be most convenient to implement the convex hull candidate list as a stack. To begin, you initialize the list by simply pushing the first three points onto the stack. The rest of the points P₃ through P_N-1 are then processed in order within the for loop.

The first two items in the parameter list of the call to the CrossProduct function require explanation. StackTop and Stack2nd are macros defined to reference the point most recently pushed onto the stack and the point pushed before that, respectively; that is, they provide access to the top two points on the stack.

For each new point, if the left-turn-only rule is violated, the while loop removes points from the stack until it is restored, and then the new point is pushed (recall that the cross product provides the necessary orientation information).

The rest of the code is fairly straightforward. Listing Two shows the point data type and the function that locates the pivot point. If there is more than one point with the smallest y-coordinate, the leftmost of these is chosen as the tiebreaker.

Listing Three provides the stack defines, data structures, and functions. Since every point in the original set may be on the convex hull, the stack must be able to accommodate N points.

Listing Four contains the function to compute the cross product and the function used in sorting the points by their polar angle they make with the pivot and the positive x-axis. I used the standard library qsort function, and this requires you to provide an item comparison function. I again use the cross product, this time to determine which of two points has the greater polar angle.

Finally, Listing Five shows everything glued together. The pivot is global so that it can be accessed by the point comparison function. Its value is set at the same time that the pivot point is swapped into array location 0.

One Small Detail

If you study Listing One closely you may realize that the "left-turn-only" rule as implemented should really be called the "no-right-turn" rule; a cross product of zero does not cause any points to be removed from the list.

This detail allows some messy special case handling to be sidestepped. The effect is that all points on the convex hull are returned, even if they are collinear, which may be undesirable in some applications.

Happily, it is a simple and straightforward matter to remove any in-between points after the fact by taking a final lap around the hull, deleting the middle point of any triple having a cross product of zero.

A Refinement

The Graham Scan achieves the theoretical lower bound of O(N log N). Locating the pivot point and extracting the hull points from the sorted list take only linear time; it is sorting the points by their polar angles that takes up the bulk of the time.

The sort's comparison function will be executed about N log N times, so it would help if the cross product computation could be removed. An obvious approach is to calculate the actual angles. This, however, requires you to use the arctangent function. Not only is this computationally expensive, but you are suddenly forced into using floating-point numbers.

If your application is such that the point coordinates must be floating-point numbers, this may not be a problem. Nonetheless, there is a beautiful device by which the angle computation can be avoided.

To set the stage, locate the leftmost and rightmost points, call them P and Q, and let L be the line through them. Let S_A be the points above L, and S_B be the points below L. It is clear that the points in the upper part of the convex hull must come from S_A, and the points in lower part must come from S_B.

You are now ready. Let the point at y=- be the pivot point. Then the lines connecting the pivot to the points in S_A are vertical, and so the ordering by polar angles is the same as the ordering by the x-coordinates.

Thus, you sort the points in S_A by their x-coordinates in descending order, and apply Graham's Scan starting at Q to determine the upper portion of the convex hull. The same thing is done with the points in S_B, except you sort them in ascending order and apply Graham's Scan starting at P to determine the upper portion of the convex hull.

Conclusion

You now have a very efficient algorithm for computing the convex hull of a set of points to add to your library. Additionally, the cross product was used for two different purposes here, and other uses were hinted at. You will likely find use for it in solving other problems in computational geometry.

References

F.P. Preparata and M.I. Shamos, Computational Geometry, 1985.

T.H. Cormen, C. E. Leiserson, and R.L. Rivest, Introduction to Algorithms, 1994.

DDJ

Listing One

//***********************************************************************
// This is the actual Graham scan algorithm. It is called after the pivot 
// point has been located and the rest of the points have been sorted by 
// their polar angle.

void Graham_Scan(point S[], int N)
{
   int i;
   point p;

   Push(S[0]);
   Push(S[1]);
   Push(S[2]);

   for (i = 3; i < N; i++)
   {
      while (CrossProduct(Stack2nd, StackTop, S[i]) < 0)
      {
         Pop(&p);
      };
      Push(S[i]);
   }
}

Back to Article

Listing Two

typedef struct
{
   int x, y;
} point;

// Find point with smallest y coordinate; if there is more
// than one select the one with the smallest x coordinate.

int LocatePivot(point S[], int N)
{
   int i, ndx, minx, miny;

// initialize minimum coordinates found so far
   ndx  = 0;
   minx = S[ndx].x;
   miny = S[ndx].y;

// locate lowest point in set
   for (i = 1; i < N; i++)
   {
      if (S[i].y > miny) continue;

//    found a candidate for a new pivot
      if (S[i].y < miny)
      {
         ndx = i;
         miny = S[i].y;
      }
      else if (S[i].x < minx)
      {
         ndx = i;
         minx = S[i].x;
      }
   }

// all done - return index of pivot point
   return(ndx);
}

Back to Article

Listing Three

#define sizeof_Stack 100
#define StackTop     stack[TOS_ptr-1]
#define Stack2nd     stack[TOS_ptr-2]

int TOS_ptr = 0;           // Top Of Stack pointer
point stack[sizeof_Stack]; // the actual point stack array

//*********************************************************
// pushes point p onto the stack
// if successful Push() returns 1
// on stack overflow Push() returns 0

int Push(point p)
{
   if (TOS_ptr == sizeof_Stack) return(0);

   stack[TOS_ptr++] = p;
   return(1);
}
//*********************************************************
// pops a point off of the stack
// if successful Pop() returns 1
// on stack underflow Pop() returns 0

int Pop(point *p)
{
   if (TOS_ptr == 0) return(0);
   *p = stack[--TOS_ptr];
   return(1);
}

Back to Article

Listing Four

//*************************************************************************
// Computes the cross product (p1-p0)x(p2-p0). If p2 is to the left of the 
// line from p0 to p1 the result will be positive.

int CrossProduct(point p0, point p1, point p2)
{
   int xprod;
   xprod = (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);
   return(xprod);
}
//*********************************************************
// Support function for qsort()
// a < b if a makes a smaller angle with Pivot than does b

int compare_points(const void *a, const void *b)
{
   int xprod;
   point *point_a, *point_b;

   point_a = (point *)a;
   point_b = (point *)b;

   xprod = CrossProduct(Pivot, *point_a, *point_b);

   if (xprod < 0) return(+1);    // a > b
   else
   if (xprod > 0) return(-1);    // a < b

   return(0);                    // a = b
}

Back to Article

Listing Five

point Pivot, S[100];

void main(void)
{
   int i, N;
   point p;

// make up some points
   N = 12;
   randomize();
   for (i = 0; i < N; i++)
   {
      S[i].x = random(10);
      S[i].y = random(10);
      printf("(%3d, %3d)\n", S[i].x, S[i].y);
   }
// find the convex hull
   i = LocatePivot(S, N);
   Pivot = S[i];
   S[i]  = S[0];
   S[0]  = Pivot;

   qsort((void *)&S[1], N-1, sizeof(S[0]), compare_points);

   Graham_Scan(S, N);

// output points of convex hull in clockwise order
   for (i = 0; Pop(&p) != 0; i++)
   {
      printf("(%3d, %3d)\n", p.x, p.y);
   }
   printf("%d points in convex hull\n", i);
}

Back to Article