Area of Triangles and Polygons (2D & 3D)

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) has 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 q between them, then Euclid says this is enough to determine the triangle and its area. Using trigonometry, the height of the triangle over the base b is given by h = a sin q, 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 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 Kevin Brown's Heron's Formula page and 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 q and j 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 parallelogram and triangle can be expressed as the magnitude of the cross-product of two edge vectors, since |v×w| = |v||w||sin q| where q is the angle between the two vectors v and w. Thus for a 3D triangle with vertices V₀V₁V₂ putting v=V₁-V₀ and w=V₂-V₀, one gets:

In 2 dimensional space (2D), a vector can be viewed as embedded in 3D by adding a third component which is set = 0. This lets one take the cross-product of 2D vectors, and use it to compute area. Given a triangle with vertices V_i=(x_i,y_i)=(x_i,y_i,0) for i=0,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 or 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 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. The 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 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 a 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 Q, 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 Q 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 V₀V₁V₂V₃, we have:

In 2D, with vertices V_i=(x_i,y_i)=(x_i,y_i,0) for i=0,3, this becomes:

which is again a signed area just as we had for triangles.

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]. So, for a quadrilateral Q=V₀V₁V₂V₃, let this parallelogram have midpoint vertices M₀M₁M₂M₃ 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₂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 M₀M₁M₂M₃ is a parallelogram. The area relation is also easy to demonstrate, and we can compute the quadrilateral's area as:

which equals 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 V_i=(x_i,y_i), it becomes a formula for a 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 Q=V₀V₁V₂V₃, 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 W be defined by its vertices V_i=(x_i,y_i) for i=0,n with V_n=V₀. Also, let P be any point; and for each edge V_iV_i+1 of W, form the triangle D_i=DPV_iV_i+1. Then, the area of W is equal to the sum of the signed areas of all the triangles D_ifor i=0,n-1; and we have:

Notice that, for a counterclockwise oriented polygon, when the point P is on the "inside" left side of an edge V_iV_i+1, then the area of D_i is positive; whereas, when P is on the "outside" right side of an edge V_iV_i+1, then D_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 D₂=DPV₂V₃ and D_n-1=DPV_n-1V_n have positive area, and contribute positively to the total area of polygon W. However, as one easily observes, only part of D₂ and D_n-1 are actually inside W and there is a part of each triangle that is also exterior. On the other hand, the triangles D₀ and D₁ have negative area and this cancels out the exterior excesses of positive area triangles. In the final analysis, exterior areas all get cancelled, and one is left with exactly the area of the polygon W.

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

A little algebra shows that the second and third summations are equal to the first. 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, but to avoid any overhead from computing the index i(mod n), one must extend the polygon array up to V_n+1=V₁.

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.

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 D=V₀V₁V₂ is given by half the magnitude of the cross product of two of its edge vectors; namely, ½ |V₀V₁× V₀V₂|.

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 W has vertices V_i=(x_i,y_i,z_i) for i=0,n, with V_n=V₀. All the vertices lie on the same 3D plane p which has a unit normal vector n. Now, the same as in the 2D case, let P be any 3D point (not generally on the plane p); and for each edge e_i=V_iV_i+1 of W, form the 3D triangle D_i=DPV_iV_i+1. We would like to relate the sum of the areas of all these triangles to the area of the polygon W in the plane p. But what we have is a pyramidal cone with P as an apex over the polygon W as a base. We need to project the triangular sides of this cone onto the plane p of the base polygon, and compute signed areas of the projected triangles. If we can do this, 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 D_i an area vector a_i=½(PV_i × PV_i+1), perpendicular to D_i, 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 T_i=DP₀V_iV_i+1. Drop a perpendicular P₀B_i from P₀ to B_i on the edge e_i=V_iV_i+1. Since PP₀ is also perpendicular to e_i, the three points PP₀B_i define a plane normal to e_i, and PB_i is a perpendicular from P to the edge e_i. Thus, |PB_i| is the height of D_i, and |P₀B_i| is the height of T_i. Further, the angle between these two altitudes = q = the angle between n and a_i since a 90° rotation 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. The same 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 W. Writing this down, we have:

Finally, by selecting P=(0,0,0), we have PV_i=V_i and this yields the simple 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. He first noted that the area of a 3D planar quadrilateral Q=V₀V₁V₂V₃ can be computed in terms of the cross-product of its diagonals; that is:

which reduces four expensive cross-product computations to just one!

Then, any large polygon W 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 £(n-1)/2. For n 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 classic standard 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 formulas we have given. We simply note this here, and do not pursue it further.

Projection to 2D

One can 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 using an area scaling factor. This method is implemented by projecting onto an axis-aligned plane by ignoring one of the three coordinates. To avoid degeneracy and optimize robustness, we look at the plane's normal vector n (see the April 2001 Algorithm About Planes), and choose the component with the greatest absolute as the one to ignore. Let Proj_c() be the projection that ignores the coordinate c = x, y, or z. Then, the ratio of areas for the projected polygon Proj_c(W) and original planar polygon W with normal n = (n_x,n_y,n_z) is:

Thus, the 3D planar area can be recovered by a single extra multiplication, 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 of the classic standard formula, achieving almost 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 (www.softsurfer.com)
// This code may be freely used and modified for any purpose
// providing that this copyright notice is included with it.
// SoftSurfer 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;    // exclude z for 2D
// 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 points with V[n]=V[0], V[n+1]=V[1]

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

// isLeft(): tests if a point is Left|On|Right of an infinite line.
//    Input: three points P0, P1, and P2
//    Return: >0 for P2 left of the line through P0 and 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 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 triangle
//    Input: three vertex points V0, V1, V2
//    Return: the (float) area of T
inline float
area2D_Triangle( Point V0, Point V1, Point V2 )
{
    return (float)isLeft(V0, V1, V2) / 2.0;
}
//===================================================================

// orientation2D_Polygon(): tests the orientation of a simple polygon
//    Input: int n = the number of vertices in the polygon
//            Point* V = an array of n+1 vertices 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 this rmin vertex
    // ccw <=> the edge leaving 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(): computes the area of a 2D polygon
//    Input: int n = the number of vertices in the polygon
//            Point* V = an array of n+2 vertices
//                       with V[n]=V[0] and V[n+1]=V[1]
//    Return: the (float) area of the polygon
float
area2D_Polygon( int n, Point* V )
{
    float area = 0;
    int   i, j, k;     // indices

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

// area3D_Polygon(): computes the area of a 3D planar polygon
//    Input: int n = the number of vertices in the polygon
//            Point* V = an array of n+2 vertices in a plane
//                       with V[n]=V[0] and V[n+1]=V[1]
//            Point N = unit 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

    // 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
    for (i=1, j=2, k=0; i<=n; i++, j++, k++)
        switch (coord) {
        case 1:
            area += (V[i].y * (V[j].z - V[k].z));
            continue;
        case 2:
            area += (V[i].x * (V[j].z - V[k].z));
            continue;
        case 3:
            area += (V[i].x * (V[j].y - V[k].y));
            continue;
        }

    // 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*ax));
        break;
    case 2:
        area *= (an / (2*ay));
        break;
    case 3:
        area *= (an / (2*az));
    }
    return area;
}
//===================================================================

References

Kevin Brown, MathPages : Heron's Formula

Donald Coxeter & Samuel Greitzer, Geometry Revisited (1967)

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

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)

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)