How do I calculate the area of a 2d polygon?

Viewed 83417

Assuming a series of points in 2d space that do not self-intersect, what is an efficient method of determining the area of the resulting polygon?

As a side note, this is not homework and I am not looking for code. I am looking for a description I can use to implement my own method. I have my ideas about pulling a sequence of triangles from the list of points, but I know there are a bunch of edge cases regarding convex and concave polygons that I probably won't catch.

16 Answers

Here is the standard method, AFAIK. Basically sum the cross products around each vertex. Much simpler than triangulation.

Python code, given a polygon represented as a list of (x,y) vertex coordinates, implicitly wrapping around from the last vertex to the first:

def area(p):
    return 0.5 * abs(sum(x0*y1 - x1*y0
                         for ((x0, y0), (x1, y1)) in segments(p)))

def segments(p):
    return zip(p, p[1:] + [p[0]])

David Lehavi comments: It is worth mentioning why this algorithm works: It is an application of Green's theorem for the functions −y and x; exactly in the way a planimeter works. More specifically:

Formula above =
integral_over_perimeter(-y dx + x dy) =
integral_over_area((-(-dy)/dy+dx/dx) dy dx) =
2 Area

The cross product is a classic.

If you have zillion of such computation to do, try the following optimized version that requires half less multiplications:

area = 0;
for( i = 0; i < N; i += 2 )
   area += x[i+1]*(y[i+2]-y[i]) + y[i+1]*(x[i]-x[i+2]);
area /= 2;

I use array subscript for clarity. It is more efficient to use pointers. Though good compilers will do it for you.

The polygon is assumed to be "closed", which means you copy the first point as point with subscript N. It also assume the polygon has an even number of points. Append an additional copy of the first point if N is not even.

The algorithm is obtained by unrolling and combining two successive iterations of the classic cross product algorithm.

I'm not so sure how the two algorithms compare regarding numerical precision. My impression is that the above algorithm is better than the classic one because the multiplication tend to restore the loss of precision of the subtraction. When constrained to use floats, as with GPU, this can make a significant difference.

EDIT: "Area of Triangles and Polygons 2D & 3D" describes an even more efficient method

// "close" polygon
x[N] = x[0];
x[N+1] = x[1];
y[N] = y[0];
y[N+1] = y[1];

// compute area
area = 0;
for( size_t i = 1; i <= N; ++i )
  area += x[i]*( y[i+1] - y[i-1] );
area /= 2;

This page shows that the formula

enter image description here

can be simplified to:

enter image description here

If you write out a few terms and group them according to common factors of xi, the equality is not hard to see.

The final summation is more efficient since it requires only n multiplications instead of 2n.

def area(x, y):
    return abs(sum(x[i] * (y[i + 1] - y[i - 1]) for i in xrange(-1, len(x) - 1))) / 2.0

I learned this simplification from Joe Kington, here.


If you have NumPy, this version is faster (for all but very small arrays):

def area_np(x, y):        
    x = np.asanyarray(x)
    y = np.asanyarray(y)
    n = len(x)
    shift_up = np.arange(-n+1, 1)
    shift_down = np.arange(-1, n-1)    
    return (x * (y.take(shift_up) - y.take(shift_down))).sum() / 2.0

To expand on the triangulate and sum triangle areas, those work if you happen to have a convex polygon OR you happen to pick a point that doesn't generate lines to every other point that intersect the polygon.

For a general non-intersecting polygon, you need to sum the cross product of the vectors (reference point, point a), (reference point, point b) where a and b are "next" to each other.

Assuming you have a list of points that define the polygon in order (order being points i and i+1 form a line of the polygon):

Sum(cross product ((point 0, point i), (point 0, point i + 1)) for i = 1 to n - 1.

Take the magnitude of that cross product and you have the surface area.

This will handle concave polygons without having to worry about picking a good reference point; any three points that generate a triangle that is not inside the polygon will have a cross product that points in the opposite direction of any triangle that is inside the polygon, so the areas get summed correctly.

Or do a contour integral. Stokes' Theorem allows you to express an area integral as a contour integral. A little Gauss quadrature and Bob's your uncle.

  1. Set a base point (the most convex point). This will be your pivot point of the triangles.
  2. Calculate the most-left point (arbitrary), other than your base point.
  3. Calculate the 2nd-most-left point to complete your triangle.
  4. Save this triangulated area.
  5. Shift over one point to the right each iteration.
  6. Sum the triangulated areas

Better than summing triangles is summing trapezoids in the Cartesian space:

area = 0;
for (i = 0; i < n; i++) {
  i1 = (i + 1) % n;
  area += (vertex[i].y + vertex[i1].y) * (vertex[i1].x - vertex[i].x) / 2.0;
}

My inclination would be to simply start slicing off triangles. I don't see how anything else could avoid being awfully hairy.

Take three sequential points that comprise the polygon. Ensure the angle is less than 180. You now have a new triangle which should be no problem to calculate, delete the middle point from the polygon's list of points. Repeat until you have only three points left.

Python code

As described here: http://www.wikihow.com/Calculate-the-Area-of-a-Polygon

With pandas

import pandas as pd

df = pd.DataFrame({'x': [10, 20, 20, 30, 20, 10, 0], 'y': [-10, -10, -10, 0, 10, 30, 20]})
df = df.append(df.loc[0])

first_product = (df['x'].shift(1) * df['y']).fillna(0).sum()
second_product = (df['y'].shift(1) * df['x']).fillna(0).sum()

(first_product - second_product) / 2
600
Related