Next Previous Contents

3. 2. 2D Polygon Computations

3.1 : How do I find the area of a polygon?

The signed area can be computed in linear time by a simple sum. The key formula is this:

If the coordinates of vertex v_i are x_i and y_i, twice the signed area of a polygon is given by

2 A( P ) = sum_{i=0}^{n-1} (x_i y_{i+1} - y_i x_{i+1}).

Here n is the number of vertices of the polygon. A rearrangement of terms in this equation can save multiplications and operate on coordinate differences, and so may be both faster and more accurate:

2 A( P ) = sum_{i=0}^{n-1} ((x_i + x_{i+1}) (y_{i+1} - y_i)).

References: [O' Rourke (C)] Thm. 1.3.3, p. 21; [Gems II] pp. 5-6: "The Area of a Simple Polygon", Jon Rokne.

To find the area of a planar polygon not in the x-y plane, use:

2 A(P) = abs(N . (sum_{i=0}^{n-1} (v_i x v_{i+1})))

where N is a unit vector normal to the plane. The `.' represents the dot product operator, the `x' represents the cross product operator, and abs() is the absolute value function.

[Gems II] pp. 170-171: "Area of Planar Polygons and Volume of Polyhedra", Ronald N. Goldman.

3.2 : How can the centroid of a polygon be computed?

The centroid (a.k.a. the center of mass, or center of gravity) of a polygon can be computed as the weighted sum of the centroids of a partition of the polygon into triangles. The centroid of a triangle is simply the average of its three vertices, i.e., it has coordinates (x1 + x2 + x3)/3 and (y1 + y2 + y3)/3. This suggests first triangulating the polygon, then forming a sum of the centroids of each triangle, weighted by the area of each triangle, the whole sum normalized by the total polygon area. This indeed works, but there is a simpler method: the triangulation need not be a partition, but rather can use positively and negatively oriented triangles (with positive and negative areas), as is used when computing the area of a polygon. This leads to a very simple algorithm for computing the centroid, based on a sum of triangle centroids weighted with their signed area. The triangles can be taken to be those formed by any fixed point, e.g., the vertex v0 of the polygon, and the two endpoints of consecutive edges of the polygon: (v1,v2), (v2,v3), etc. The area of a triangle with vertices a, b, c is half of this expression:


                (b[X] - a[X]) * (c[Y] - a[Y]) -
                (c[X] - a[X]) * (b[Y] - a[Y]);
   

Code available at ftp://cs.smith.edu/pub/code/centroid.c (3K). Reference: [Gems IV] pp.3-6; also includes code.

3.3 : How do I find if a point lies within a polygon?

The definitive reference is "Point in Polyon Strategies" by Eric Haines [Gems IV] pp. 24-46. The code in the Sedgewick book Algorithms (2nd Edition, p.354) is incorrect.

The essence of the ray-crossing method is as follows. Think of standing inside a field with a fence representing the polygon. Then walk north. If you have to jump the fence you know you are now outside the poly. If you have to cross again you know you are now inside again; i.e., if you were inside the field to start with, the total number of fence jumps you would make will be odd, whereas if you were ouside the jumps will be even.

The code below is from Wm. Randolph Franklin wrf@ecse.rpi.edu with some minor modifications for speed. It returns 1 for strictly interior points, 0 for strictly exterior, and 0 or 1 for points on the boundary. The boundary behavior is complex but determined; in particular, for a partition of a region into polygons, each point is "in" exactly one polygon. See the references below for more detail.


    int pnpoly(int npol, float *xp, float *yp, float x, float y)
    {
      int i, j, c = 0;
      for (i = 0, j = npol-1; i < npol; j = i++) {
        if ((((yp[i]<=y) && (y<yp[j])) ||
             ((yp[j]<=y) && (y<yp[i]))) &&
            (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))

c = !c; } return c; }

The code may be further accelerated, at some loss in clarity, by avoiding the central computation when the inequality can be deduced, and by replacing the division by a multiplication for those processors with slow divides. For code that distinguishes strictly interior points from those on the boundary, see [O'Rourke (C)] pp. 239-245.

References: [Gems IV] pp. 24-46 [O'Rourke (C)] [Glassner:RayTracing]

3.4 : How do I find the intersection of two convex polygons?

Unlike intersections of general polygons, which might produce a quadratic number of pieces, intersection of convex polygons of n and m vertices always produces a polygon of at most (n+m) vertices. There are a variety of algorithms whose time complexity is proportional to this size, i.e., linear. The first, due to Shamos and Hoey, is perhaps the easiest to understand. Let the two polygons be P and Q. Draw a horizontal line through every vertex of each. This partitions each into trapezoids (with at most two triangles, one at the top and one at the bottom). Now scan down the two polygons in concert, one trapezoid at a time, and intersect the trapezoids if they overlap vertically. A more difficult-to-describe algorithm is in [O'Rourke (C)], pp. 242-252. This walks around the boundaries of P and Q in concert, intersecting edges. An implementation is included in [O'Rourke (C)].

3.5 : How do I do a hidden surface test (backface culling) with 2d points?

c = (x1-x2)*(y3-y2)-(y1-y2)*(x3-x2)

x1,y1, x2,y2, x3,y3 = the first three points of the polygon.

If c is positive the polygon is visible. If c is negative the polygon is invisible (or the other way).

3.6 : How do I find a single point inside a simple polygon?

Given a simple polygon, find some point inside it. Here is a method based on the proof that there exists an internal diagonal, in [O'Rourke (C), 13-14]. The idea is that the midpoint of a diagonal is interior to the polygon.

1. Identify a convex vertex v; let its adjacent vertices be a and b. 2. For each other vertex q do: 2a. If q is inside avb, compute distance to v (orthogonal to ab). 2b. Save point q if distance is a new min. 3. If no point is inside, return midpoint of ab, or centroid of avb. 4. Else if some point inside, qv is internal: return its midpoint.

Code for finding a diagonal is in [O'Rourke (C), 35-39], and is part of many other software packages. See Subject 0.07: Where is all the source?

3.7 : How do I find the orientation of a simple polygon?

Compute the signed area (Subject 2.01). The orientation is counter-clockwise if this area is positive.

A slightly faster method is based on the observation that it isn't necessary to compute the area. Find the lowest vertex (or, if there is more than one vertex with the same lowest coordinate, the rightmost of those vertices) and then take the cross product of the edges fore and aft of it. Both methods are O(n) for n vertices, but it does seem a waste to add up the total area when a single cross product (of just the right edges) suffices. Code for this is available at ftp://cs.smith.edu/pub/code/polyorient.C (2K).

The reason that the lowest, rightmost (or any other such extreme) point works is that the internal angle at this vertex is necessarily convex, strictly less than pi (even if there are several equally-lowest points).


Next Previous Contents