Introduction
========================================================================
Sometimes you stumble upon a problem that you think might be incredibly
complicated but then turns out to have a surprisingly easy solution,
of course the opposite happens perhaps even more frequently but this post
is about the first kind where want to calculate the area of a polygon.
The area of any shape is together with the length of its circumference or
boundary the most basic metric we could want to investigate and gives
one answer to the question of _how large_ an object is.
![Figure [areas]: Some simple polygonal shapes, ordered and colored by their relative area.](img/areas.svg width=100%)
Suppose you wanted to measure the size of a lake but you don't have a map,
unlikely nowadays but try to imagine such a scenario. One idea might be to
try to divide the lake into square patches that you could measure
by putting out marks at certain points.
It could be done but you would have to make sure the squares are large or small
enough to cover the lake while not overlapping and you would have to keep
track of which marks belong to which patch.
Or you could maybe _sweep_ the lake using small slices, rowing back and
forth, and measuring how wide the lake is at that point.
Certainly feasible for some lakes but if we have something with a shape
like the one in Figure [introshape] this could mean using a lot of small square
patches or a lot of rowing across.
![Figure [introshape]: To find the area of this shape we only have
to walk along the black outline.](img/intro_shape.svg width=60%)
In this post I will show how we we only have to follow the waterside and
walk around the lake, while keeping track of our relative position, to
find out its surface area.
All thanks to Greens theorem!
!!!
**I wasn't magic, it was multivariate calculus!**
When I took my first course in C++ (a distressing number of years ago)
one of the assignments contained a quite common example of a heirarchy
of classes that represented various geometric shapes such as circle,
square, triangle and polygon (some people have strong oppinions about
object oriented programming, I think it's fine in moderation) and they
all came with a method for calulating their area. The formulas for all
the common primitive shapes are known by most people but for the polygon
we where provided a function that computed some sum over all edges and
I just accepted it as some mathematical magic that did its job and
left it at that.
I simply did not expect anything I had seen in multivariate calculus
such as Green's Theorem to show up in an introductory course in C++
since those two subjects still lived in different worlds for me!
Green's theorem and the area of a simple connected domain
========================================================================
The fundamental theorem of calculus states that when evaluating the integral of
a function defined as $\frac{d}{dx} f(x)$ on some interval $[a, b]$ in fact
we only have to know the value of the function $f(x)$ at the endpoints of the
interval.
\begin{equation}
\int_{a}^{b} \frac{d}{dx} f(x) dx = f(b) - f(a)
\end{equation}
A version of this in 2D is [Green's theorem](https://en.wikipedia.org/wiki/Green%27s_theorem) [#Green].
It states that the integral of the $\unitk$-component of the curl
of a vector field, $\textbf{F}(x, y)$, over some simple connected domain, $\Omega$,
is equal to the line integral of the tangential, $\textbf{T}$, components of $\textbf{F}$ along
the boundary $\bnd$ such that we have the following equation.
\begin{equation}
\iint_{\Omega} \nabla \times \textbf{F} \cdot \textbf{k} dA = \oint_{\bnd} \textbf{F} \cdot d \textbf{r}
\end{equation}
The domain and its boundary is illustrated in Figure [domain] and we should note
that the $\unitk$-component will point _out_ of the screen in such a way that
$\textbf{N} = \textbf{T} \times \unitk$ along the boundary.
![Figure [domain]: An oriented bounded domain with tangent $\textbf{T}$ and normal $\textbf{N}$ indicated at one point on the boundary.](img/oriented_domain.svg width=60%)
The use of the _nabla_ operator $\nabla \times \textbf{F}$ is called the _curl_
of $\textbf{F}$ and is defined as
\begin{equation}
\nabla \times \textbf{F} = \begin{vmatrix}
\uniti & \unitj & \unitk \\
\frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\
F_1 & F_2 & F_3
\end{vmatrix}
= (\frac{\partial F_3}{\partial y} -\frac{\partial F_2}{\partial z}) \uniti +
(\frac{\partial F_1}{\partial z} -\frac{\partial F_3}{\partial x}) \unitj +
(\frac{\partial F_2}{\partial x} -\frac{\partial F_1}{\partial y}) \unitk \\
\end{equation}
Where $F_i$ denotes the $i$ -th component of $\textbf{F}$.
We use this and expand the integral expression.
\begin{equation}
\iint_{\Omega} \left( \frac{\partial F_2}{\partial x} - \frac{\partial F_1}{\partial y} \right) dx dy =
\oint_{\bnd} F_1(x,y) dx + F_2(x,y) dy
\end{equation}
So just as we only have to evaluate a function at the endpoints of the interval
in 1D we only have to evaluate a function along the boundary of the domain in
2D!
In order to evaluate the area we want to have a vector field such that
$\frac{\partial F_2}{\partial x} - \frac{\partial F_1}{\partial y} = 1$
so that we get the integral
\begin{equation}
\iint_{\Omega} \frac{\partial F_2}{\partial x} - \frac{\partial F_1}{\partial y} dx dy = \iint_{\Omega} 1 dx dy = \text{Area of } \Omega
\end{equation}
Some such vector fields are $\textbf{F} = x \unitj$, $\textbf{F} = - y \uniti$
and $\textbf{F} = \frac{1}{2} (-y \uniti + x \unitj)$.
All of these give us integrals that evaluate to the area of $\Omega$.
\begin{equation}
\oint_{\bnd} x dy = -\oint_{\bnd} y dx = \oint_{\bnd} x dy - y dx = \iint_{\Omega} 1 dA = \text{Area of } \Omega
\end{equation}
A simple polygon in 2D
========================================================================
A polygon is an ordered set of points, also called vertices, $v_i \in V = \{v_0, ..., v_N\}$
that have an x and a y-coordinate, $v_i \in \mathbb{R}^2$.
Two consecutive vertices are connected by an edge, defined as the indices of
its two edges $\{i, i+1\}$, where the list loops around so that the first and the
last vertices are also connected by an edge.
A simple polygon is one where every vertex is connected to two edges only and
edges only intersect at vertices.
See my previous blog post on [simple polygons and a method for finding a triangulation](http://all-systems-phenomenal.com/articles/ear_clipping_triangulation/index.php).
We also want the polygon to be consistently oriented, so we can say that it has
an inside and an outside we pick the same counter clockwise orientation as
we did for the boundaries of general simple connected domains.
Any simple connected domain in 2D may be approximated by a simple polygon and
here we see a polygonal representation of the shape in Figure [domain].
![Figure [polygon]: A simple polygon with 15 vertices that make up the boundary of a simple connected domain.](img/polygon.svg width=60%)
Area of the polygon
========================================================================
Calculating the integral along one edge
-------------------------------------------------------------------------------
We begin with the integral expression directly from Green's theorem using
$\mathbf{F} = -y \uniti$ as our vector field. We then consider
the integral along the line segment between two consecutive vertices,
$v_0 = (x_0, y_0)$ and $v_1 = (x_1, y_1)$, on the polygon and directly
substitute $y$ for the expression of the line.
\begin{align}
-\int_{x_0}^{x_1} y dx &= -\int_{x_0}^{x_1} \frac{y_1 - y_0}{x_1 - x_0}(x_0 - x_1) + y_0 dx \\
&= - \int_{x_0}^{x_1} B x - B x_0 + y_0 dx
\end{align}
We evaluate the integral expression and simplify by gathering and
factorization of the terms with $B$.
\begin{equation}
-\left[ \frac{1}{2} B x^2 - B x_0 x + y_0 x \right]_{x_0}^{x_1}
\end{equation}
\begin{align}
=& - \left( \frac{1}{2} B x_{1}^{2} - \frac{1}{2} B x_{0}^{2} - A x_0 x_1 + B x_{0}^{2} + y_0 x_1 - y_0 x_0 \right) \\
=& - \left( \frac{1}{2} B x_{0}^{2} + \frac{1}{2} B x_{1}^{2} - A x_0 x_1 + y_0 x_1 - y_0 x_0 \right) \\
=& - \left( \frac{1}{2} B (x_{0} - x_{1})^{2} + y_0 x_1 - y_0 x_0 \right)
\end{align}
We then expand $B$ and simplify by eliminating terms.
\begin{align}
& - \left( \frac{1}{2} \frac{y_1 - y_0}{x_1 - x_0} (x_0 - x_1)^{2} + y_0 x_1 - y_0 x_0 \right) \\
=& - \left( \frac{1}{2} \frac{y_0 - y_1}{x_0 - x_1} (x_0 - x_1)^{2} + y_0 x_1 - y_0 x_0 \right) \\
=& - \left( \frac{1}{2} (y_0 - y_1) (x_0 - x_1) + y_0 x_1 - y_0 x_0 \right) \\
=& - \left( \frac{1}{2} y_0 x_0 - \frac{1}{2} y_1 x_0 - \frac{1}{2} y_0 x_1 + \frac{1}{2} y_1 x_1 + y_0 x_1 - y_0 x_0 \right) \\
=& - \left(-\frac{1}{2} y_0 x_0 - \frac{1}{2} y_1 x_0 + \frac{1}{2} y_0 x_1 + \frac{1}{2} y_1 x_1 \right)
\end{align}
And finally we end up with:
\begin{equation}
\frac{1}{2} \left( y_0 x_0 + y_1 x_0 - y_0 x_1 - y_1 x_1 \right)
\end{equation}
The integral along all edges of the polygon
-------------------------------------------------------------------------------
Now let's go back and plug this into the sum we got from integrating over the
polygon. So we get one contribution from each edge of the polygon.
Let's begin with the edge $\{0, 1\}$ and then continue all around to calculate
the area $A$.
\begin{align}
A =& \frac{1}{2} \left( {\color{magenta}\cancel{y_0 x_0}} + y_1 x_0 - y_0 x_1 - {\color{red}\cancel{y_1 x_1}} \right) \\
+& \frac{1}{2} \left( {\color{red}\cancel{y_1 x_1}} + y_2 x_1 - y_1 x_2 - {\color{green}\cancel{y_2 x_2}} \right) \\
+& \frac{1}{2} \left( {\color{green}\cancel{y_2 x_2}} + y_3 x_2 - y_2 x_3 - {\color{blue}\cancel{y_2 x_2}} \right) \\
+& \frac{1}{2} \left( {\color{blue}\cancel{y_3 x_3}} + y_4 x_3 - y_3 x_4 - {\color{purple}\cancel{y_4 x_4}} \right) \\
+& \cdots \\
+& \frac{1}{2} \left( {\color{orange}\cancel{y_N x_N}} + y_0 x_N - y_N x_0 - {\color{magenta}\cancel{y_0 x_0}} \right) \\
\end{align}
We see that we get a sum where the terms to the sides cancel out so that the
area can be calculated as a sum over all edges $E = \{\{0, 1\}, \{1, 2\}, \dotsc , \{N, 0\}\}$.
\begin{equation}
A = \frac{1}{2} \sum_{\{i, j\} \in E} y_j x_i - y_i x_j
\end{equation}
This criss cross cancelation has caused this expression to be called the
[_shoelace formula_](https://en.wikipedia.org/wiki/Shoelace_formula) [#Wiki].
Python implementation
========================================================================
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Python
def computeAreaOfPolygon(vertices):
"""
Compute the area of a polygon in 2D.
Assumption that the polygon is simple, i.e is closed and has no
crossings and also that its vertex order is counter clockwise.
"""
area = 0.0
n = len(vertices)
for i in range(n):
j = i + 1
j %= n
p0 = vertices[i]
p1 = vertices[j]
area += p1[1]*p0[0] - p0[1]*p1[0]
area *= 0.5
return area
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!!! Tip
What really made me raise my eyebrows wasn't actually seeing this formula
but it turns out that **it has a natural generalization in 3D** as well, using
**Stoke's theorem**.
When I first started working with triangulated surface meshes I quite early
found myself wanting to calculate or estimate their volumes and at first
I thought that maybe I would have to first compute some tetrahedalized
version of them wich is not an easy task. However when I looked online I
found, just as for polygons and their edges in 2D, that it could be
computed by some magical sum over all the triangles.
It was only much later that I understood from where this sum actually
originated and why it worked!
Polygons with holes
========================================================================
One very nice property of this way of calculating the area is that it
works just as well for domains, or polygons, with holes.
The holes are defined as other simple polygons that lie completely
inside the first one. The polygon inside should be oriented in clockwise
direction so that our notion of inside and outside is consistent.
We can actually see that we should be able to calculate the area of the
shaded domain by simply calculating it first for the one without the
hole, then calculating it for the hole itself and subtracting this from
the first, total area. And this is exactly what will happen when we just
apply the method itself. Since the orientation of the hole is the opposite
of that of the outer boundary we could put a negative sign in front of its
integral, this applies for any directed integrals, and will get a negative
contribution to the area.
So we don't care if the polygon has holes or not.
The method for calculating the area works just the same!
![Figure [polygon-with-hole]: A simple domain with a hole inside.
The hole is oriented in opposite direction to the outer boundary.](img/domain_with_hole.svg width=60%)
Conclusion
========================================================================
I have presented one way of deriving an expression for calculating the area
of a polygon in 2D.
Another very elegant way of doing this is shown in Chapter 5 of [#Crane21].
The area of a polygon is shown to be computed by the signed area of every
triangle formed by taking two vertices of every edge and some other point
placed at any arbitrary position in space, inside or outside the polygon.
The illustration is very intuitive and I highly recommend giving it a look.
The proof is based on the exterior calculus version of Stokes theorem.
There is nothing wrong with taking something for granted, after all we don't
have unlimited time on out hands to understand _everything_. But it's always
interesting, and often a good exercise, to try to understand where certain
formulas, methods or algorithms come from in order to, step by step, deepen
your own knowledge in whatever field you are interested in.
+++++
**Bibliography**:
[#Crane21]: K. Crane. (2021).
Discrete Differential Geometry: An Applied Introduction
http://www.cs.cmu.edu/~kmcrane/Projects/DGPDEC/paper.pdf
[#Wiki]: Wikipedia
Shoelace formula
https://en.wikipedia.org/wiki/Shoelace_formula
[#Green]: Wikipedia
Green's theorem
https://en.wikipedia.org/wiki/Green%27s_theorem
Code for all figures and examples:
https://bitbucket.org/nils_olovsson/area_of_polygon/src/master/