Moment of inertia for an arbitrary polygon Is there an easy way to calculate the moment of inertia of an arbitrary 2-dimensional polygon, with respect to an axis perpendicular to the plane of the polygon (passing through the center of mass)?  
What I intended to do:

*

*Derive a formula for the moment of inertia of a triangle, with respect to an axis through a vertex of the triangle.

*Apply this formula on every triangle center of mass - vertex - next vertex in the polygon and add the results. 
 
What I achieved:

*

*For a triangle with vertices $(0,0)$, $(a,0)$ and $(b,c)$, the moment of inertia with respect to the axis $x=0\wedge y=0$ is: 
$$\int_0^{b}\int_0^{\frac{c}{b}\cdot x}x^2+y^2\text{ }dy\text{ }\text{ }\text{ }dx - \int_a^b\int_0^{\frac{c}{b-a}\cdot(x-a)}x^2+y^2\text{ }dy\text{ }\text{ }\text{ }dx$$
$$=\int_0^{b}\dfrac{\left(c^3+3b^2c\right)x^3}{3b^3}dx - \int_a^b\dfrac{c\left(\left(c^2+3b^2-6ab+3a^2\right)x^3+\left(-3ac^2-3ab^2+6a^2b-3a^3\right)x^2+3a^2c^2x-a^3c^2\right)}{3\left(b^3-3ab^2+3a^2b-a^3\right)}dx$$
$$=\dfrac{b\left(c^3+3b^2c\right)}{12} - \dfrac{\left(b-a\right)c\left(c^2+3b^2+2ab+a^2\right)}{12}$$
$$=\frac{c^3a+ca^3+bca^2+b^2ca}{12}$$
$$=\frac{ac}{12}\cdot\left(a^2+b^2+c^2+ab\right)$$

*Now we can calculate the moment of inertia of a triangle with vertices $A(x_1,y_1)$, $B(x_2,y_2)$ and $C(x_3,y_3)$ with respect to the axis through A as follows:
The Euclidean transformation that takes $A$ to $(0,0)$ and $B$ to $(\sqrt{x_2^2+y_2^2},0)$, takes $C$ to $(x_3',y_3')$ (which can be calculated quite easily).
As a Euclidean transformation preserves the moment of inertia, we can apply the preceding formula.

*Now we can apply this to every triangle center of mass - vertex - next vertex in the polygon and add the results to get the total moment of inertia.
  
I haven't tested this yet (I'll do that later), but could someone please check my calculations and my method? 
Furthermore, if there is an easier way to calculate this, please let me know ;)
 A: Turning Eli's comment into an answer to my own post:  
Apparently, there exists something called the second moment of area (see also: https://en.wikipedia.org/wiki/Second_moment_of_area). 
This is the mathematical definition of the physical moment of inertia.  

The second moment of area for an arbitrary shape $R$ with respect to an arbitrary axis $a$ is defined as:
$$J_{a} = \iint\limits_{R}\rho^2 dA$$
When the axis is the $x$-axis, the second moment of area can be calculated as follows:
$$J_x=\iint\limits_R y^2 dx dy$$
And when it's the $y$-axis:
$$J_y=\iint\limits_R x^2 dy dx$$  
For a polygon $A_0A_1A_2\cdots A_{n-1}$ with $A_i=(x_i,y_i)$ and $A_n=A_0$, those are equal to:
$$J_x = \frac{1}{12}\sum_{i=0}^{n-1}\left(x_iy_{i+1}-x_{i+1}y_i\right)\left(y_i^2+y_iy_{i+1}+y_{i+1}^2\right)$$
$$J_y = \frac{1}{12}\sum_{i=0}^{n-1}\left(x_iy_{i+1}-x_{i+1}y_i\right)\left(x_i^2+x_ix_{i+1}+x_{i+1}^2\right)$$
Now, denote by $z$ the line perpendicular to the $xy$-plane that passes through the origin. Then:
$$J_z = \iint\limits_R\rho^2 dA = \iint\limits_Rx^2+y^2\text{ }dA = \iint\limits_Ry^2\text{ }dxdy + \iint\limits_Rx^2\text{ }dydx = J_x+J_y$$
(this is the so-called perpendicular axis theorem) 
Thanks to this, we can calculate $J_z$!   
I strongly believe that I was going to end up with this as well, due to the $12$ in the denominator and the fact that only subsequent $x$- and $y$-values are multiplied.
A: We can also derive the formula for the second moment of area w.r.t to the axis through the origin directly. 
Let $A_0(=A_n), A_1\ldots, A_{n-1}$ be the vertices of your polygon $P$ where $A_i = (x_i, y_i)$.
For $\vec{r} = (x,y)$ notice that $\operatorname{div}(r^2\vec{r}) = 4r^2$ so using the two-dimensional divergence theorem we get
$$J_z = \int_{P} r^2\,dA = \frac14\int_P \operatorname{div}(r^2\vec{r})\,dA = \frac14\int_{\partial P}r^2\,\vec{r}\cdot \hat{n}\,d\gamma$$
where we integrate over the boundary $\gamma = \partial P$ of the polygon with $\hat{n}$ being the outward pointing normal.
Notice that for $i=0,\ldots, n-1$ you can parameterize the line segment $[A_i, A_{i+1}]$ with $\gamma_i : [-1,1] \to \mathbb{R}^2$ given as
$$\gamma_i(t) = \frac{A_{k+1}+A_k + t(A_{k-1}-A_k)}2, \quad t\in[-1,1]$$
and the outward pointing normal is $\hat{n} = \frac{(A_{k+1}-A_k)\times \hat{z}}{\|A_{k+1}-A_k\|}$ where $\hat{z} = (0,0,1)$.
We have $\gamma_i'(t) = \frac12(A_{k-1}-A_k) $ so $\hat{n}\,d\gamma_i = \frac12 (A_{k+1}-A_k)\times \hat{z}\,dt$. Therefore
\begin{align}
J_z &= \frac14\int_{\partial P}r^2\,\vec{r}\cdot \hat{n}\,d\gamma \\
&= \frac14\sum_{i=0}^{n-1}\int_{\gamma_i}r^2\,\vec{r}\cdot \hat{n}\,d\gamma_i \\
&= \frac1{64}\sum_{i=0}^{n-1}\int_{-1}^1 \big\|A_{k+1}+A_k + t(A_{k-1}-A_k)\big\|^2\,\big(A_{k+1}+A_k + t(A_{k-1}-A_k)\big)\cdot  (A_{k+1}-A_k)\times \hat{z}\,dt
\end{align}
Using the triple product rule we get
\begin{align}\big(A_{k+1}+A_k + t(A_{k-1}-A_k)\big)\cdot  \left((A_{k+1}-A_k)\times \hat{z}\right) &= \left(\big(A_{k+1}+A_k + t(A_{k-1}-A_k)\big)\times  (A_{k+1}-A_k)\right)\cdot \hat{z}\\
&= 2(A_{k+1}\times A_k)\cdot \hat{z}\\
&= 2(x_{k}y_{k+1}-x_{k+1}y_k)
\end{align}
Hence by integrating it follows
\begin{align}
J_z &= \frac1{32}\sum_{k=0}^{n-1}(x_{k}y_{k+1}-x_{k+1}y_k)\int_{-1}^1 \big\|A_{k+1}+A_k + t(A_{k-1}-A_k)\big\|^2\,dt\\
&=  \frac1{12}\sum_{k=0}^{n-1}(x_{k}y_{k+1}-x_{k+1}y_k)\left(x_{k+1}^2+x_{k+1}x_k+x_k^2 + y_{k+1}^2+y_{k+1}y_k+y_k^2\right)
\end{align}
which recovers your formula.
Inspired by this more general calculation.
