Subscribe
to eAlerts
for Geom Site
Updates

CLICK HERE

by Dan Sunday

Computing the area of a planar polygon is a basic geometry calculation and can be found in many introductory texts. However, there are several different methods for computing planar areas depending on the information available.

Triangles

Ancient Triangles

Before Pythagoras, the area of the parallelogram (including the rectangle and the square) had been known to equal the product of its base times its height. Further, two copies of the same triangle paste together to form a parallelogram, and thus the area of a triangle is half of its base b times its height h. So, for these simple but commonly occurring cases, we have:

Parallelogram:  Triangle:

However, except in special situations, finding the height of a triangle at an arbitrary orientation usually requires also computing the perpendicular distance of the top vertex from the base.

For example, if one knows the lengths of two sides, a and b, of a triangle and also the angle between them, then this is enough to determine the triangle and its area. Using trigonometry, the height of the triangle over the base b is given by , and thus the area is:

Another frequently used computation is derived from the fact that triangles with equal sides are congruent, and thus have the same area. This observation from Euclid (~300 BC) culminated in Heron's formula (~50 AD) for area as a function of the lengths of its three sides [Note: some historians attribute this result to Archimedes (~250 BC)]; namely:

where a,b,c are the lengths of the sides, and s is the semiperimeter. There are interesting algebraic variations of this formula; such as:

which efficiently avoids calculating the 3 square roots to explicitly get the lengths a,b,c from the triangle's vertex coordinates. Other variations on Heron's formula can be found at Eric Weisstein's Triangle page.

The remaining classical triangle congruence is when two angles and one side are known. Knowing two angles gives all three, so we can assume the angles and are both adjacent to the known base b. Then the formula for area is:

Modern Triangles

More recently, starting in the 17-th century with Descartes and Fermat, linear algebra produced new simple formulas for area. In 3 dimensional space (3D), the area of a planar parallelogram or triangle can be expressed by the magnitude of the cross-product of two edge vectors, since where is the angle between the two vectors v and w. Thus for a 3D triangle with vertices putting and , one gets:

A 2D vector (x, y) can be viewed as embedded in 3D by adding a third z component set = 0. This lets one take the cross-product of 2D vectors, and use it to compute area. Given a triangle with vertices for i = 0,1,2, we can compute that:

And the absolute value of the third z-component is twice the absolute area of the triangle. However, it is useful to not take the absolute value here, and instead let the area be a signed quantity.

This formula for area is a very efficient computation with no roots and no trigonometric functions involved - just 2 multiplications and 5 additions, and possibly 1 division by 2 (which can sometimes be avoided).

Note that the signed area will be positive if the vertices V₀V₁V₂ are oriented counterclockwise around the triangle, and will be negative if the triangle is oriented clockwise; and so, this area computation can be used to test for a triangle's orientation. This signed area can also be used to test whether the point V₂ is to the left (positive) or the right (negative) of the directed line segment V₀V₁ . So this value is a very useful primitive, and it's great to have such an efficient formula for it.

Quadrilaterals

The Greeks singled out certain quadrilaterals (also called quadrangles) for special treatment, including the square, the rectangle, the parallelogram, and the trapezium. Then, given an arbitrary quadrilateral, they showed how to construct a parallelogram [Euclid, Book I, Prop 45] or square [Euclid, Book II, Prop 14] with an equal area. And the area of the parallelogram was equal to its base times its height. But there was no general formula for the quadrilateral's area.

An extension of Heron's triangle area formula to quadrilaterals was discovered by the Hindu geometer Brahmagupta (620 AD) [Coxeter,1967, Section 3.2]. However, it only works for cyclic quadrilaterals where all four vertices lie on the same circle. For a cyclic quadrilateral , let the lengths of the four sides be a, b, c, d, and the semiperimeter be s = (a+b+c+d)/2. Then, the area of is given by:

which is an amazing symmetric formula. If one side is zero length, say d = 0, then we have a triangle (which is always cyclic) and this formula reduces to Heron's one.

In modern linear algebra, as already noted, the area of a planar parallelogram is the magnitude of the cross product of two adjacent edge vectors. So, for any 3D planar parallelogram , we have:

In 2D, with vertices for , this becomes:

which is again a signed area, just as we had for triangles. It also indicates orientation.

Next, for an arbitrary quadrilateral, one can compute its area using a parallelogram discovered by Pierre Varignon (first published in 1731). It is amazing that the Greeks missed Varignon's simple result which was discovered 2000 years after Euclid! Given any quadrilateral, one can take the midpoints of its 4 edges to get 4 vertices which form a new quadrilateral. It is then easy to show that this midpoint quadrilateral is always a parallelogram, called the "Varignon parallelogram", and that its area is exactly one-half the area of the original quadrilateral [Coxeter,1967, Section 3.1]. To see this, for any quadrilateral , let the midpoint vertices be as shown in the diagram:

From elementary geometry, we know that in triangle V₀V₁V₂ the midpoint line M₀M₁ is parallel to the base V₀V₂. In triangle V₀V₂V₃, the line M₂M₃ is parallel to that same base V₀V₂. Thus, M₀M₁ and M₂M₃ are parallel to each other. Similarly, M₀M₃ and M₁M₂ are parallel, which shows that is a parallelogram. The area relation is also easy to demonstrate. And we can then compute the area as:

which is one-half the magnitude of the cross-product of the two diagonals of the quadrilateral. This result was noted by [Van Gelder,1995] who used a different proof. This formula holds for any 3D planar quadrilateral. When restricted to 2D with , this becomes a formula for any quadrilateral’s signed area:

This formula for an arbitrary quadrilateral is just as efficient as the one for an arbitrary triangle, using only 2 multiplications and 5 additions. For simple quadrilaterals, the area is positive when the vertices are oriented counterclockwise, and negative when they are clockwise. However, it also works for nonsimple quadrilaterals and is equal to the difference in area of the two regions the quadrilateral bounds. For example, in the following diagram where I is the self-intersection point of a nonsimple quadrilateral , we have:

Polygons

2D Polygons

A 2D polygon can be decomposed into triangles. For computing area, there is a very easy decomposition method for simple polygons (i.e. ones without self intersections). Let a polygon be defined by its vertices for with . Also, let P be any point; and for each edge of , form the triangle . Then, the area of is equal to the sum of the signed areas of all the triangles for ; and we have:

Notice that, for a counterclockwise oriented polygon, when the point P is on the "inside" left side of an edge , then the area of is positive; whereas, when P is on the "outside" right side of an edge , then has a negative area. If instead the polygon is oriented clockwise, then the signs are reversed, and inside triangles become negative.

For example, in the above diagram, the triangles and have positive area, and contribute positively to the total area of polygon . However, as one can see, only part of and are actually inside and there is a part of each triangle that is also exterior. On the other hand, the triangles and have negative area, and this cancels out the exterior excesses of positive area triangles. In the final analysis, the exterior areas all get canceled, and one is left with exactly the area of the polygon .

One can make the formula more explicit by picking a specific point P and expanding the terms. By selecting , the area formula of each triangle reduces to . This yields:

A little algebra shows that the more efficient second and third summations are equal to the first [Sunday, 2002]. For a polygon with n vertices, the first summation uses 2n multiplications and (2n–1) additions; the second uses n multiplications and (3n–1) additions; and the third uses only n multiplications and (2n–1) additions. So, the third is preferred for efficiency. And, to avoid any overhead from computing the index i(mod n), one can either extend the polygon array up to , or simply put the final term outside of the summation loop. We give an efficient implementation below in the routine area2D_Polygon().

This computation gives a signed area for a polygon; and, similar to the signed area of a triangle, is positive when the vertices are oriented counterclockwise around the polygon, and negative when oriented clockwise. So, this computation can be used to test for a polygon's global orientation. However, there are other more efficient algorithms for determining polygon orientation. The easiest is to find the rightmost lowest vertex of the polygon, and then test the orientation of the entering and leaving edges at this vertex. This test can be made by checking if the end vertex of the leaving edge is to the left of the entering edge, which means that the orientation is counterclockwise, otherwise it is clockwise.

3D Planar Polygons

An important generalization is for planar polygons embedded in 3D space [Goldman, 1994]. We have already shown that the area of a 3D triangle is given by half the magnitude of the cross product of two edge vectors; namely, .

The Standard Formula

There is a classic standard formula for the area of a 3D polygon [Goldman, 1994] that extends the cross-product formula for a triangle. It can be derived from Stokes Theorem. However, we show here how to derive it from a 3D triangular decomposition that is geometrically more intuitive.

A general 3D planar polygon has vertices for with , where all the vertices lie on the same 3D plane P  which has a unit normal vector n. Now, as in the 2D case, let P be any 3D point (not generally on the plane P ); and for each edge of , form the 3D triangle . We would like to relate the sum of the areas of all these triangles to the area of the polygon in the plane P . But what we have is a pyramidal cone with P as an apex over the polygon as a base. We are going to project the triangular sides of this cone onto the plane P  of the base polygon, and compute signed areas of the projected triangles. Then the sum of the projected areas will equal the total area of the planar polygon.

To achieve this, start by associating to each triangle an area vector , which is perpendicular to , and whose magnitude we know is equal to that triangle's area. Next, drop a perpendicular from P to a point P₀ on P , and consider the projected triangle . Then drop a perpendicular P₀B_i from P₀ to B_i on the edge . Since PP₀ is also perpendicular to e_i, the three points PP₀B_i define a plane that is perpendicular to e_i, and thus PB_i is a perpendicular from P to e_i. Thus |PB_i| is the height of , and |P₀B_i| is the height of T_i. Further, the angle between these two altitudes = = the angle between n and since a 90° rotation (in the PP₀B_i plane) results in congruence. This gives:

This signed area computation is positive if the vertices of T_i are oriented counterclockwise when we look at the plane P  from the side pointed to by n. As in the 2D case, we can now add together the signed areas of all the triangles T_i to get the area of the polygon . Writing this down, we have:

Finally, by selecting P = (0,0,0), we have PV_i = V_i and this produces the concise formula:

which uses 6n+3 multiplications and 4n+2 additions

Similar to the 2D case, this is a signed area which is positive when the vertices are oriented counterclockwise around the polygon when viewed from the side of P  pointed to by n.

Quadrilateral Decomposition

[Van Gelder, 1995] has shown how to significantly speed up this computation by using a decomposition into quadrilaterals instead of triangles. As we have already shown, the area of a 3D planar quadrilateral can be computed in terms of the cross-product of its diagonals; namely as: , which reduces four expensive cross-product computations to just one!

Then, any polygon (with n > 4 vertices) can be decomposed into quadrilaterals formed by V₀ and three other sequential vertices V_2i–1, V_2i, and V_2i+1 for i = 1,h where h = the greatest integer . When n is odd, the decomposition ends with a triangle. This gives:

where k = 0 for n odd, and k = n–1 for n even. This formula reduces the number of expensive cross-products by a factor of two (replacing them with vector subtractions). In total there are 3n+3 multiplications and 5n+1 additions making this formula roughly twice as fast as the classical one.

[Van Gelder, 1995] also states that this method can be applied to 2D polygons, but he does not write down the details. Working this out produces a formula that uses n multiplications and 3n–1 additions, which is not as fast as the prior 2D formula we have given that use only n multiplications and (2n–1) additions. We simply note this here, and do not pursue it further.

Projection to 2D

However, one can further significantly speed up the computation of 3D planar polygon area by projecting the polygon onto a 2D plane [Snyder & Barr, 1987]. Then, the area can be computed in 2D using our fastest formula, and the 3D area is recovered by using an area scaling factor. This method projects the 3D-embedded polygon onto an axis-aligned 2D subplane, by ignoring one of the three coordinates. Further, to get the correct area sign for the orientation of the projected polygons, the basis vectors of the 2D subplanes have to be ordered so that their orientation is compatible with the 3D basis vectors. (I’d like to thank Nathan Lay for pointing this out, which has resulted in a sign-correction in both the formula and the implemented code). In particular, if (x,y,z) are the 3D coordinates for the standard xyz basis, then the projections that ignore x, y, and z are onto the yz, zx, and xy subplanes respectively. That is, the projections are given by: Proj_x(x,y,z) = (y,z), Proj_y(x,y,z) = (z,x), and Proj_z(x,y,z) = (x,y). Then, these projections are “orientation-preserving”, in that the sign of the projected polygon’s area matches it’s orientation in the projection subplane. And because of this, the scaling factor uses the sign of the component that is ignored.

Next, to select the projection that best avoids degeneracy and optimizes robustness, we look at the polygon's normal vector n, and choose the component with the greatest absolute value as the one to ignore. Let Proj_c() be the projection that ignores the selected coordinate c = x, y, or z. Then, the ratio of areas for the projected polygon Proj_c() and original polygon with normal is given by:

This corrects an error in the previously published formula that was noted by [Nathan Lay, 2013].

An interesting consequence of this formula, is that the “area vector” a of is now given by:

This area vector a is normal to the plane of the 3D polygon, and it’s length is the area of the 3D polygon. So conversely, if the 3D polygon area A() is already known, and n is the unit normal vector for ’s plane, then A()n is the area vector, whose components are thus the areas of ’s projections. In effect, by computing the area of one polygon projection, one can quickly get the areas of the other projections.

As a result, the 3D planar area can be recovered by a single extra multiplication in addition to the 2D area computation, and in total this algorithm uses n+5 multiplications, 2n+1 additions, 1 square root (when n is not a unit normal), plus a small overhead choosing the coordinate to ignore. This is a very significant improvement over the classical 3D formula, achieving a 6 times speed-up!! We give an efficient implementation below in the routine area3D_Polygon().

 
Implementations

Here are some sample "C++" implementations of these formulas as algorithms. We just give the 2D case with integer coordinates, and use the simplest structures for a point, a triangle, and a polygon which may differ in your application. We represent a polygon as an array of points, but it is often more convenient to have it as a linked list of vertices (to allow insertion or deletion during drawing operations), and the polygon routines can be easily modified to scan through the linked list (see [O'Rourke, 1998] for an example of this approach).

// Copyright 2000 softSurfer, 2012 Dan Sunday
// This code may be freely used and modified for any purpose
// providing that this copyright notice is included with it.
// iSurfer.org makes no warranty for this code, and cannot be held
// liable for any real or imagined damage resulting from its use.
// Users of this code must verify correctness for their application.
 

// a Point (or vector) is defined by its coordinates
typedef struct {int x, y, z;} Point;  // set z=0 for a 2D Point

// a Triangle is given by three points: Point V0, V1, V2

// a Polygon is given by:
//       int n = number of vertex points
//       Point* V[] = an array of n+1 vertex points with V[n]=V[0]
 

// Note: for efficiency low-level short functions are declared to be inline.

// isLeft(): test if a point is Left|On|Right of an infinite 2D line.
//    Input:  three points P0, P1, and P2
//    Return: >0 for P2 left of the line through P0 to P1
//          =0 for P2 on the line
//          <0 for P2 right of the line
inline int
isLeft( Point P0, Point P1, Point P2 )
{
    return ( (P1.x - P0.x) * (P2.y - P0.y)
           - (P2.x - P0.x) * (P1.y - P0.y) );
}
//===================================================================

// orientation2D_Triangle(): test the orientation of a 2D triangle
//  Input:  three vertex points V0, V1, V2
//  Return: >0 for counterclockwise
//          =0 for none (degenerate)
//          <0 for clockwise
inline int
orientation2D_Triangle( Point V0, Point V1, Point V2 )
{
    return isLeft(V0, V1, V2);
}
//===================================================================

// area2D_Triangle(): compute the area of a 2D triangle
//  Input:  three vertex points V0, V1, V2
//  Return: the (float) area of triangle T
inline float
area2D_Triangle( Point V0, Point V1, Point V2 )
{
    return (float)isLeft(V0, V1, V2) / 2.0;
}
//===================================================================

// orientation2D_Polygon(): test the orientation of a simple 2D polygon
//  Input:  int n = the number of vertices in the polygon
//          Point* V = an array of n+1 vertex points with V[n]=V[0]
//  Return: >0 for counterclockwise
//          =0 for none (degenerate)
//          <0 for clockwise
//  Note: this algorithm is faster than computing the signed area.
int
orientation2D_Polygon( int n, Point* V )
{
    // first find rightmost lowest vertex of the polygon
    int rmin = 0;
    int xmin = V[0].x;
    int ymin = V[0].y;

    for (int i=1; i<n; i++) {
        if (V[i].y > ymin)
            continue;
        if (V[i].y == ymin) {   // just as low
            if (V[i].x < xmin)  // and to left
                continue;
        }
        rmin = i;      // a new rightmost lowest vertex
        xmin = V[i].x;
        ymin = V[i].y;
    }

    // test orientation at the rmin vertex
    // ccw <=> the edge leaving V[rmin] is left of the entering edge
    if (rmin == 0)
        return isLeft( V[n-1], V[0], V[1] );
    else
        return isLeft( V[rmin-1], V[rmin], V[rmin+1] );
}
 //===================================================================

// area2D_Polygon(): compute the area of a 2D polygon
//  Input:  int n = the number of vertices in the polygon
//          Point* V = an array of n+1 vertex points with V[n]=V[0]
//  Return: the (float) area of the polygon
float
area2D_Polygon( int n, Point* V )
{
    float area = 0;
    int  i, j, k;   // indices

    if (n < 3) return 0;  // a degenerate polygon

    for (i=1, j=2, k=0; i<n; i++, j++, k++) {
        area += V[i].x * (V[j].y - V[k].y);
    }
    area += V[n].x * (V[1].y - V[n-1].y);  // wrap-around term
    return area / 2.0;
}
//===================================================================

// area3D_Polygon(): compute the area of a 3D planar polygon
//  Input:  int n = the number of vertices in the polygon
//          Point* V = an array of n+1 points in a 2D plane with V[n]=V[0]
//          Point N = a normal vector of the polygon's plane
//  Return: the (float) area of the polygon
float
area3D_Polygon( int n, Point* V, Point N )
{
    float area = 0;
    float an, ax, ay, az; // abs value of normal and its coords
    int  coord;           // coord to ignore: 1=x, 2=y, 3=z
    int  i, j, k;         // loop indices

    if (n < 3) return 0;  // a degenerate polygon

    // select largest abs coordinate to ignore for projection
    ax = (N.x>0 ? N.x : -N.x);    // abs x-coord
    ay = (N.y>0 ? N.y : -N.y);    // abs y-coord
    az = (N.z>0 ? N.z : -N.z);    // abs z-coord

    coord = 3;                    // ignore z-coord
    if (ax > ay) {
        if (ax > az) coord = 1;   // ignore x-coord
    }
    else if (ay > az) coord = 2;  // ignore y-coord

    // compute area of the 2D projection
    switch (coord) {
      case 1:
        for (i=1, j=2, k=0; i<n; i++, j++, k++)
            area += (V[i].y * (V[j].z - V[k].z));
        break;
      case 2:
        for (i=1, j=2, k=0; i<n; i++, j++, k++)
            area += (V[i].z * (V[j].x - V[k].x));
        break;
      case 3:
        for (i=1, j=2, k=0; i<n; i++, j++, k++)
            area += (V[i].x * (V[j].y - V[k].y));
        break;
    }
    switch (coord) {    // wrap-around term
      case 1:
        area += (V[n].y * (V[1].z - V[n-1].z));
        break;
      case 2:
        area += (V[n].z * (V[1].x - V[n-1].x));
        break;
      case 3:
        area += (V[n].x * (V[1].y - V[n-1].y));
        break;
    }

    // scale to get area before projection
    an = sqrt( ax*ax + ay*ay + az*az); // length of normal vector
    switch (coord) {
      case 1:
        area *= (an / (2 * N.x));
        break;
      case 2:
        area *= (an / (2 * N.y));
        break;
      case 3:
        area *= (an / (2 * N.z));
    }
    return area;
}
//===================================================================

 
References

Donald Coxeter & Samuel Greitzer, Geometry Revisited (1967)

Ronald Goldman, "Area of Planar Polygons and Volume of Polyhedra" in Graphics Gems II (1994)

Nathan Lay, personal communication (2013).

Joseph O'Rourke, Computational Geometry in C (2nd Edition), Section 1.3 "Area of a Polygon" (1998)

J.P. Snyder and A.H. Barr, "Ray Tracing Complex Models Containing Surface Tessellations", ACM Comp Graphics 21, (1987)

Dan Sunday, "Fast Polygon Area and Newell Normal Computation", journal of graphics tools (jgt) Vol 7(2), (2002)

Allen Van Gelder, "Efficient Computation of Polygon Area and Polyhedron Volume" in Graphics Gems V (1995)

Eric Weisstein, MathWorld Site: Triangle or in The CRC Concise Encyclopedia of Mathematics (1998)
 

© Copyright 2012 Dan Sunday, 2001 softSurfer