Logo_iSurf_org

Subscribe
to eAlerts

for Geom Site
Updates

CLICK HERE

 

Area of Triangles and Polygons

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: 

Area_eqn0a
Area_eqn1

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 theta 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 Area_eqn1B, and thus the area is:

Area_eqn2

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:

Area_eqn3a

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:

Area_eqn3b2

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 theta and phi are both adjacent to the known base b. Then the formula for area is:

Area_eqn3c

 

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 vxw-abs=vw_sin where theta is the angle between the two vectors v and w. Thus for a 3D triangle with vertices V0V1V2 putting v=V1-V0 and w=V2-V0, one gets:

Area_pic4
Area_eqn4

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 Vi=(xi,yi)=(xi,yi,0) for i = 0,1,2, we can compute that:

Area_eqn5a

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.

Area_eqn5b2
Area_pic5

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 V0V1V2 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 V2 is to the left (positive) or the right (negative) of the directed line segment V0V1 . 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 THETA-cap-bold, 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 THETA-cap-bold is given by:

Quad_eqn0

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 Quad_eqn0_P=V0V1V2V3, we have:

Quad_pic1
Quad_eqn1


In 2D, with vertices Vi=(xi,yi)=(xi,yi,0) for i=0,3, this becomes:

Quad_eqn2

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 Quad_eqn3_Q=V0V1V2V3, let the midpoint vertices be Quad_eqn3_M0M1M2M3 as shown in the diagram:

Quadrangle

From elementary geometry, we know that in triangle V0V1V2 the midpoint line M0M1 is parallel to the base V0V2. In triangle V0V2V3, the line M2M3 is parallel to that same base V0V2. Thus, M0M1 and M2M3 are parallel to each other. Similarly, M0M3 and M1M2 are parallel, which shows that Quad_eqn3_M0M1M2M3 is a parallelogram. The area relation is also easy to demonstrate. And we can then compute the area as:

Quad_eqn3

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 Vi=(xi,yi), this becomes a formula for any quadrilateral’s signed area:

Quad_eqn4

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 Quad_eqn3_Q=V0V1V2V3, we have:

Quad_pic2
Quad_eqn5

 

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 OMEGA be defined by its vertices Vi=(xi,yi) for i=0,n with Vn=V0. Also, let P be any point; and for each edge ei=ViVi+1 of OMEGA, form the triangle DELTA-i=DELTA-PViVi+1. Then, the area of OMEGA is equal to the sum of the signed areas of all the triangles DELTA-i for i=0,n-1; and we have:

Polygon_pic1
Polygon_eqn1

Notice that, for a counterclockwise oriented polygon, when the point P is on the "inside" left side of an edge ViVi+1, then the area of DELTA-i is positive; whereas, when P is on the "outside" right side of an edge ViVi+1, then DELTA-i 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 DELTA-2 and DELTA-n-1 have positive area, and contribute positively to the total area of polygon OMEGA. However, as one can see, only part of DELTA-2 and DELTA-n-1 are actually inside OMEGA and there is a part of each triangle that is also exterior. On the other hand, the triangles DELTA-0 and DELTA-1 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 OMEGA.

One can make the formula more explicit by picking a specific point P and expanding the terms. By selecting P=(0,0), the area formula of each triangle reduces to 2A(DELTA-i)=(xiyi+1-xi+1yi). This yields:

Polygon_eqn2
Polygon_pic2


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 Vn+1=V1, 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 DELTA-V0V1V2 is given by half the magnitude of the cross product of two edge vectors; namely, Area_eqn4b.

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.

Polygon_pic3A general 3D planar polygon OMEGA has vertices Vi=(xi,yi,zi) for i=0,n with Vn=V0, 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 ei=ViVi+1 of OMEGA, form the 3D triangle DELTA-i=DELTA-PViVi+1. We would like to relate the sum of the areas of all these triangles to the area of the polygon OMEGA in the plane P . But what we have is a pyramidal cone with P as an apex over the polygon OMEGA 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.

Polygon_pic4To achieve this, start by associating to each triangle DELTA-i=DELTA-PViVi+1 an area vector Area_eqn4c, which is perpendicular to DELTA-i, and whose magnitude we know is equal to that triangle's area. Next, drop a perpendicular from P to a point P0 on P , and consider the projected triangle T-i=DELTA-P0ViVi+1. Then drop a perpendicular P0Bi from P0 to Bi on the edge ei=ViVi+1. Since PP0 is also perpendicular to ei, the three points PP0Bi define a plane that is perpendicular to ei, and thus PBi is a perpendicular from P to ei. Thus |PBi| is the height of DELTA-i, and |P0Bi| is the height of Ti. Further, the angle between these two altitudes = theta = the angle between n and alpha_vector-i since a 90 rotation (in the PP0Bi plane) results in congruence. This gives:

Polygon_eqn3a

This signed area computation is positive if the vertices of Ti 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 Ti to get the area of the polygon OMEGA. Writing this down, we have:

Polygon_eqn3b

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

Polygon_pic5
Polygon_eqn4

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 Quad_eqn3_Q=V0V1V2V3 can be computed in terms of the cross-product of its diagonals; namely as: Polygon_eqn5, which reduces four expensive cross-product computations to just one!

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

Polygon_eqn6

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: Projx(x,y,z) = (y,z), Projy(x,y,z) = (z,x), and Projz(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 Projc() be the projection that ignores the selected coordinate c = x, y, or z. Then, the ratio of areas for the projected polygon Projc(OMEGA) and original polygon OMEGA with normal n=(nx,ny,nz) is given by:

Polygon_eqn7

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 OMEGA is now given by:

Polygon_eqn8

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(OMEGA) is already known, and n is the unit normal vector for OMEGA’s plane, then A(OMEGA)n is the area vector, whose components are thus the areas of OMEGA’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