How do I calculate the area of a 2d polygon?
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.
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
can be simplified to:
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