By Jason Z


2009-01-16 18:18:08 8 Comments

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 comments

@abe312 2017-04-02 22:47:44

I'm going to give a few simple functions for calculating area of 2d polygon. This works for both convex and concave polygons. we simply divide the polygon into many sub-triangles.

//don't forget to include cmath for abs function
struct Point{
  double x;
  double y;
}
// cross_product
double cp(Point a, Point b){ //returns cross product
  return a.x*b.y-a.y*b.x;
}

double area(Point * vertices, int n){  //n is number of sides
  double sum=0.0;
  for(i=0; i<n; i++){
    sum+=cp(vertices[i], vertices[(i+1)%n]); //%n is for last triangle
  }
  return abs(sum)/2.0;
}

@user9315861 2018-07-09 17:03:35

cp takes two arguments, yet you're calling it with one.

@firelynx 2017-09-29 06:29:58

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

@Darius Bacon 2009-01-16 18:39:05

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

@David Lehavi 2009-01-17 06:44:43

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_permieter(-y dx + x dy) = integral_area((-(-dy)/dy+dx/dx)dydyx = 2 Area

@Yakov 2011-07-11 12:42:19

The link in the post is dead.Does anyone have another?

@Andrew Андрей Листочкин 2011-09-07 12:10:33

The linked discussion on [email protected] mailing list is unavalable to me. I copied the message from Google Cache: gist.github.com/1200393

@imsc 2012-02-17 10:09:37

Here is another link for explanation. mathopenref.com/coordpolygonarea.html

@Sibbs Gambling 2013-11-09 08:34:20

does the point order matter? I mean must it be in counterclockwise or clockwise or doesn't matter?

@Darius Bacon 2013-11-19 09:38:24

@perfectionm1ng switching directions would flip the sign in the sum, but abs() strips the sign out.

@OneWorld 2016-10-12 07:01:34

Limitations: This method will produce the wrong answer for self-intersecting polygons, where one side crosses over another, as shown on the right. It will work correctly however for triangles, regular and irregular polygons, convex or concave polygons. (mathopenref.com/coordpolygonarea.html)

@Darius Bacon 2016-10-15 03:35:54

The 'wrong' answer is a matter of definition. For example, if you computed the area of a pentagram this way, it'd leave out the area of the central pentagon. But if you rendered the filled pentagram in e.g. SVG, it'd also leave that pentagon unfilled! I agree that it's worth bringing up that this answer gives the absolute value of the signed area.

@PHPirate 2018-06-25 12:37:02

This is also known at Wikipedia as Shoelace formula.

@Nathan 2019-02-01 14:57:27

I started to write this question, then realized the answer to it was already in Darius' answer. Hopefully my comment can help someone else who has this question; given as few as 4 points, a concave polygon is ambiguous; consider the points (0,0), (-1,-1),(1,0), and (0,1); the area is dependent on which 2 points are connected to (0,0) in the middle. Darius answered this in the "implicitly wrapping around from the last vertex to the first," though, because the order tells you which 2 other points (0,0) is adjacent to

@Joseph 2015-12-02 21:27:54

C way of doing that:

float areaForPoly(const int numVerts, const Point *verts)
{
    Point v2;
    float area = 0.0f;

    for (int i = 0; i<numVerts; i++){
        v2 = verts[(i + 1) % numVerts];
        area += verts[i].x*v2.y - verts[i].y*v2.x;
    }

    return area / 2.0f;
}

@Mahdi 2015-06-20 06:07:23

Implementation of Shoelace formula could be done in Numpy. Assuming these vertices:

import numpy as np
x = np.arange(0,1,0.001)
y = np.sqrt(1-x**2)

We can define the following function to find area:

def PolyArea(x,y):
    return 0.5*np.abs(np.dot(x,np.roll(y,1))-np.dot(y,np.roll(x,1)))

And getting results:

print PolyArea(x,y)
# 0.26353377782163534

Avoiding loop makes this function ~50X faster than PolygonArea:

%timeit PolyArea(x,y)
# 10000 loops, best of 3: 42 µs per loop
%timeit PolygonArea(zip(x,y))
# 100 loops, best of 3: 2.09 ms per loop

Note: I had written this answer for another question, I just mention this here to have a complete list of solutions.

@unutbu 2013-11-09 11:48:55

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

@physicsmichael 2014-06-11 17:03:53

Thanks for the NumPy version.

@Steve 2012-10-04 15:05:38

To calc the area of the polygon

http://community.topcoder.com/tc?module=Static&d1=tutorials&d2=geometry1#polygon_area

int cross(vct a,vct b,vct c)
{
    vct ab,bc;
    ab=b-a;
    bc=c-b;
    return ab.x*bc.y-ab.y*bc.x;
}    
double area(vct p[],int n)
{ 
    int ar=0;
    for(i=1;i+1<n;i++)
    {
        vct a=p[i]-p[0];
        vct b=p[i+1]-p[0];
        area+=cross(a,b);
    }
    return abs(area/2.0);
}    

@Mark Taylor 2012-10-04 21:05:45

This is a 3 year old question with 34 upvotes on the accepted answer. Tell us how your answer is better than any of the other answers already posted.

@underdoeg 2015-08-27 14:11:27

It is an example in c and not python. Not better but nice to have it in different languages

@MSN 2009-01-16 18:40:31

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.

@user836725 2011-08-21 15:06:08

language independent solution:

GIVEN: a polygon can ALWAYS be composed by n-2 triangles that do not overlap (n = number of points OR sides). 1 triangle = 3 sided polygon = 1 triangle; 1 square = 4 sided polygon = 2 triangles; etc ad nauseam QED

therefore, a polygon can be reduced by "chopping off" triangles and the total area will be the sum of the areas of these triangles. try it with a piece of paper and scissors, it is best if you ca visualize the process before following.

if you take any 3 consecutive points in a polygons path and create a triangle with these points, you will have one and only one of three possible scenarios:

  1. resulting triangle is completely inside original polygon
  2. resulting triangle is totally outside original polygon
  3. resulting triangle is partially contained in original polygon

we are interested only in cases that fall in the first option (totally contained).

every time we find one of these, we chop it off, calculate its area (easy peasy, wont explain formula here) and make a new polygon with one less side (equivalent to polygon with this triangle chopped off). until we have only one triangle left.

how to implement this programatically:

create an array of (consecutive) points that represent the path AROUND the polygon. start at point 0. run the array making triangles (one at a time) from points x, x+1 and x+2. transform each triangle from a shape to an area and intersect it with area created from polygon. IF the resulting intersection is identical to the original triangle, then said triangle is totally contained in polygon and can be chopped off. remove x+1 from the array and start again from x=0. otherwise (if triangle is outside [partially or completely] polygon), move to next point x+1 in array.

additionally if you are looking to integrate with mapping and are starting from geopoints, you must convert from geopoints to screenpoints FIRST. this requires deciding a modelling and formula for earths shape (though we tend to think of the earth as a sphere, it is actually an irregular ovoid (eggshape), with dents). there are many models out there, for further info wiki. an important issue is whether or not you will consider the area to be a plane or to be curved. in general, "small" areas, where the points are up to some km apart, will not generate significant error if consider planar and not convex.

@chmike 2009-04-04 16:38:57

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;

@Cygon 2012-04-29 09:19:03

I cannot imagine the second code snippet will work. It's pretty obvious that, the further the polygon is on the X axis, the larger its area would be.

@chmike 2012-05-14 09:12:05

It is a correct mathematical rearrangement of the algorithm described above saving some multiplications. You are right, but the areas defined by other vertexes will subtract. But this may indeed lead to precision degradation.

@Cygon 2012-05-15 14:58:10

Sorry, but I don't think this is right. One just needs to look at the loop to see that it cannot work. Here's an example program that demonstrates that the result of the second algorithm changes completely depending on the translation of the polygon on the X axis: pastebin.com/Mb8uQpz5

@chmike 2012-05-15 15:45:26

I verified with a sample program using your values and don't see a difference in area. I do see a difference between your code and the second algorithm: The for loop condition is i <= N and not i < n. This may explain why you don't get the expected result.

@chmike 2012-05-15 15:52:04

What you overlooked is that the addition has always some negative terms because of the y subtraction. Consider any 2d polygonal shape and compare y values of consecutive vertexes. You'll see that some subtraction will yield a negative value and some positive.

@Cygon 2012-05-15 16:07:26

Indeed, that last paragraph is what I couldn't wrap my mind around! With i <= N it works. Thanks for your patience, I take everything back :)

@chmike 2012-05-16 11:45:58

No problem. I'm glad to be helpful. That is the purpose of StackExchange. Beside, you could have found a bug or a wrong answer, in which case you would have been helpful to me. :)

@NightElfik 2014-11-06 06:57:08

On a side note, the area returned by the algorithm is "signed" (negative or positive based on order of points) so if you want always positive area just use absolute value.

@ploth 2018-09-19 10:34:22

The first algorithm is dangerous. After the first iteration you may access behind the array if there are only 4 points.

@David Hanak 2009-01-16 19:34:35

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;
}

@Joe Phillips 2009-01-16 18:25:57

  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

@recursive 2009-01-17 20:24:42

Make sure you negate the triangle area if the next point is moving "backwards".

@mbeckish 2009-01-16 18:34:18

A set of points without any other constraints don't necessarily uniquely define a polygon.

So, first you have to decide what polygon to build from these points - perhaps the convex hull? http://en.wikipedia.org/wiki/Convex_hull

Then triangulate and calculate area. http://www.mathopenref.com/polygonirregulararea.html

@duffymo 2009-01-16 18:33:40

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.

@MattK 2009-01-16 18:24:40

One way to do it would be to decompose the polygon into triangles, compute the area of the triangles, and take the sum as the area of the polygon.

@Loren Pechtel 2009-01-16 18:23:55

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.

@Richard 2012-11-17 19:44:45

The hairy part about this is that if your three consecutive points define a triangle outside or partially outside of the polygon, then you have a problem.

@Loren Pechtel 2012-11-18 20:28:21

@Richard: That's why the qualification about 180 degrees. If you slice off a triangle outside the polygon you'll end up with too many degrees.

@Richard 2012-11-18 22:09:14

you may need to better describe how you are finding the angle. There is no way in plane geometry to have 3 points as part of a triangle and have any angle or combination of angles exceed 180 degrees - the check would seem to be meaningless.

@Loren Pechtel 2012-11-19 19:08:24

@Richard: On your polygon you have the angle of every junction. If the relevant triangle would lie outside the polygon the angle between the two segments will be greater than 180 degrees.

@Richard 2012-11-19 20:36:35

You mean the interior angle of the two adjacent edge segments would be greater than 180 degrees.

@Loren Pechtel 2012-11-20 21:17:47

Yup--I'm looking at the angles as seen from the inside of the polygon. If the prospective triangle actually lies outside the polygon you will have three angles that are greater than 180 as seen from this vantage point. Testing the first is enough, though.

Related Questions

Sponsored Content

54 Answered Questions

20 Answered Questions

[SOLVED] How do CSS triangles work?

4 Answered Questions

5 Answered Questions

[SOLVED] Distance from a point to a polygon

36 Answered Questions

[SOLVED] How can I pair socks from a pile efficiently?

6 Answered Questions

[SOLVED] How to determine a diagonal is in or out of a concave polygon?

  • 2009-03-29 00:09:52
  • Kamran Bigdely
  • 5131 View
  • 4 Score
  • 6 Answer
  • Tags:   geometry polygon

9 Answered Questions

3 Answered Questions

[SOLVED] Intersection of a line and a concave polygon 3D

Sponsored Content