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
abs() strips the sign out. –
Nike (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 –
Silurid p = np.asarray([[1306, 157], [1330, 178], [1347, 126], [1322, 106]]) this wrongly gives 24526.5 (should give 1600) –
Daisydaitzman def area(polygon): return abs(numpy.cross(polygon, numpy.roll(polygon, -1, 0)).sum() / 2) –
Daisydaitzman 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
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.
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
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.
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.
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;
}
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:
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.
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);
}
cross() function takes 3 parameters, but is given only 2? –
Rubin 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.
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;
}
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;
}
As described here: http://www.wikihow.com/Calculate-the-Area-of-a-Polygon
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
© 2022 - 2024 — McMap. All rights reserved.