Triangle Integrals

$ \def \spekhaken {\iint} \def \J {\Delta} \def \half {\frac{1}{2}} \def \hieruit {\quad \Longrightarrow \quad} \def \EN {\quad \mbox{and} \quad} $ Our goal is to calculate some (first, second, third order) moments of an arbitrary triangle. To be mathematically precise: $$ \mu_x = \spekhaken x \, dx dy \EN \mu_y = \spekhaken y \, dx dy $$ $$ \sigma_{xx} = \spekhaken x^2 \, dx dy \EN \sigma_{yy} = \spekhaken y^2 \, dx dy \EN \sigma_{xy} = \spekhaken x y \, dx dy $$ Or even: $$ M_{xxx} = \spekhaken x^3 \, dx dy \EN M_{xxy} = \spekhaken x^2 y \, dx dy $$ $$ M_{xyy} = \spekhaken x y^2 \, dx dy \EN M_{yyy} = \spekhaken y^3 \, dx dy $$ Following the theory in the previous paragraph, global coordinates $(x,y)$ can be expressed in their local counterparts $(\xi,\eta)$ : $$ \begin{array}{ll} x - x_0 = \xi .(x_1 - x_0) + \eta.(x_2 - x_0) \\ y - y_0 = \xi .(y_1 - y_0) + \eta.(y_2 - y_0) \end{array} $$ It makes no difference for the outcome of the integrals if a more handsome choice for the coordinate system is to be preferred. Therefore, let one of the vertices of the triangle, say $(x_0,y_0)$, be selected as the origin of our global coordinate system. Then: $$ x = \xi .x_1 + \eta.x_2 \EN y = \xi .y_1 + \eta.y_2 $$ And the Jacobian $\J$ of this transformation is involved with: $$ dx.dy = (x_1 y_2 - x_2 y_1) \, d\xi d\eta = \J \, d\xi d\eta $$ Limited use will be made of Newton's binomial formula: $$ (a + b)^n = \sum_{k=0}^n C(n,k)\, a^k b^{n-k} = \sum_{k=0}^n \frac{n!}{(n-k)!\,k!} \, a^k b^{n-k} $$ The formulas for any moment of a triangle take the following general form: $$ \spekhaken x^m y^n \, dx dy = \spekhaken ( \xi .x_1 + \eta.x_2 )^m ( \xi .y_1 + \eta.y_2 )^n \, \J \, d\xi d\eta = $$ $$ \spekhaken \sum_{i=0}^m C(m,i)\, \xi^i x_1^i \eta^{m-i} x_2^{m-i} \sum_{j=0}^n C(n,j)\, \xi^j y_1^j \eta^{n-j} y_2^{n-j} \, \J \, d\xi d\eta = $$ $$ \J \, \sum_{i=0}^m \sum_{j=0}^n C(m,i)\, x_1^i x_2^{m-i} C(n,j)\, y_1^j y_2^{n-j} \spekhaken \xi^{i+j} \eta^{m+n-i-j} \, d\xi d\eta $$ The above formula - being far too complicated - will not be used in the sequel. It turns out, however, that we have to calculate integrals like: $$ F(m,n) = \spekhaken \xi^m \eta^n \,d\xi d\eta $$ Herewith the integration is carried out over a rectangular equilateral triangle, with local coordinates $\xi$ and $\eta$, where $0 \le \xi \le 1$ and $0 \le \eta \le (1-\xi)$ . Working out a few steps: $$ F(m,n) = \spekhaken \xi^m \eta^n \,d\xi d\eta = \int_0^1 \xi^m \left[ \int_0^{1-\xi} \eta^n d\eta \right] d\xi $$ $$ = \int_0^1 \xi^m \left[ \frac{(1-\xi)^{n+1}}{n+1} \right] d\xi = \int_0^1 \frac{(1-\xi)^{n+1}}{n+1} d\left(\frac{\xi^{m+1}}{m+1}\right) $$ $$ = \left[ \frac{(1-\xi)^{n+1}}{n+1} \frac{\xi^{m+1}}{m+1} \right]_0^1 - \int_0^1 \frac{\xi^{m+1}}{m+1} d\left(\frac{(1-\xi)^{n+1}}{n+1}\right) = 0 \, + \int_0^1 \frac{\xi^{m+1}}{m+1} (1-\xi)^n \, d\xi $$ $$ = \frac{n}{m+1} \int_0^1 \xi^{m+1} \left[ \frac{(1-\xi)^n}{n} \right] d\xi = \frac{n}{m+1} \int_0^1 \xi^{m+1} \left[ \int_0^{1-\xi} \eta^{n-1} d\eta \right] d\xi $$ $$ = \frac{n}{m+1} \spekhaken \xi^{m+1} \eta^{n-1} d\xi d\eta = \frac{n}{m+1} F(m+1,n-1) $$ Now we can set up the following sequence of formulas: $$ F(m,n) = \frac{n}{m+1} F(m+1,n-1) = \frac{n}{m+1} \frac{n-1}{m+2} F(m+2,n-2) = ... $$ $$ = \frac{ n (n-1)\, ...\, 2.1}{ (m+1)(m+2)\, ...\, (m+n-1)(m+n) } F(m+n,0) $$ $$ = \frac{ n(n-1)\, ...\, 2.1\, .\, m (m-1)\, ..\, 2.1} { 1.2.\, ... \, (m-1) m (m+1) \, ...\, (m+n) } F(m+n,0) = \frac{ m!\,n! }{ (m+n)! } F(m+n,0) $$ So only integrals of the form $F(m+n,0)$ are left to be calculated: $$ \spekhaken \xi^{m+n} \, d\xi d\eta = \int_0^1 \xi^{m+n} \left[ \int_0^{1-\xi} d\eta \right] d\xi = \int_0^1 \xi^{m+n} (1-\xi) d\xi $$ $$ = \int_0^1 \xi^{m+n} d\xi - \int_0^1 \xi^{m+n+1} d\xi = \left[ \frac{\xi^{m+n+1}}{m+n+1} \right]_0^1 - \left[ \frac{\xi^{m+n+2}}{m+n+2} \right]_0^1 $$ $$ = \frac{1}{m+n+1} - \frac{1}{m+n+2} = \frac{(m+n+2)-(m+n+1)}{(m+n+1)(m+n+2)} $$ $$ = \frac{1}{(m+n+1)(m+n+2)} = F(m+n,0) $$ Hence: $$ F(m,n) = \frac{ m!\,n! }{ (m+n)! } \frac{1}{(m+n+1)(m+n+2)} = \frac{ m!\,n! }{ (m+n+2)! } $$ This is the final result: $$ \spekhaken \xi^m \eta^n \,d\xi d\eta = \frac{ m!\,n! }{ (m+n+2)! } $$ Now let's calculate a few of these triangle moments. $$ \mbox{Area} = \spekhaken \, dx dy = \spekhaken \, d\xi d\eta \: \J = \frac{ 0!\,0! }{ (0+0+2)! } \, \J = \J / 2 = \half (x_1 y_2 - x_2 y_1 ) $$ Since all (other) moments have to be divided by this area, the outcome of their integrals have to be multiplied with a factor $2/\J$ . A first order moment is: $$ \frac{\spekhaken x \, dx dy}{\mbox{Area}} = 2/\J \, \spekhaken (x_1 \xi + x_2 \eta ) \, d\xi d\eta \: \J = 2 x_1 \, \spekhaken \xi \, d\xi d\eta + 2 x_2 \, \spekhaken \eta \, d\xi d\eta $$ $$ = 2 x_1 \, \frac{ 1!\,0! }{ (1+0+2)! } + 2 x_2 \, \frac{ 0!\,1! }{ (0+1+2)! } = x_1 \, 2/6 + x_2 \, 2/6 = (x_1 + x_2) / 3 $$ In very much the same way (replace $x$ by $y$) we can prove that: $$ \frac{\spekhaken y \, dx dy}{\mbox{Area}} = (y_1 + y_2) / 3 $$ Different though it seems, this is the same as the familiar result that the coordinates of the midpoint of a triangle equal one-third of the coordinates of the vertices: $$ \begin{array}{l} \overline{x} = x_0 + 1/3 \left[ (x_1 - x_0) + (x_2 - x_0) \right] = (x_0 + x_1 + x_2) /3 \\ \overline{y} = y_0 + 1/3 \left[ (y_1 - y_0) + (y_2 - y_0) \right] = (y_0 + y_1 + y_2) / 3 \end{array} $$ Second order moments are: $$ 2 / \J \spekhaken x^2 \, dx dy \EN 2 / \J \spekhaken y^2 \, dx dy \EN 2 / \J \spekhaken x y \, dx dy $$ It is sufficient to calculate only the last integral. Proper substitutions in $\overline{xy}$ will take care of the other two later on. $$ 2 / \J \spekhaken x y \, dx dy = 2 \spekhaken (x_1 \xi + x_2 \eta )(y_1 \xi + y_2 \eta ) \, d\xi d\eta $$ $$ = 2 x_1 y_1 \spekhaken \xi^2 \, d\xi d\eta + 2 x_1 y_2 \spekhaken \xi \eta \, d\xi d\eta + 2 x_2 y_1 \spekhaken \xi \eta \, d\xi d\eta + 2 x_2 y_2 \spekhaken \eta^2 \, d\xi d\eta $$ $$ = 2 x_1 y_1 \frac{ 2!\,0! }{ (2+0+2)! } + 2 x_1 y_2 \frac{ 1!\,1! }{ (1+1+2)! } + 2 x_2 y_1 \frac{ 1!\,1! }{ (1+1+2)! } + 2 x_2 y_2 \frac{ 0!\,2! }{ (0+2+2)! } $$ $$ = 2 x_1 y_1 . 1/12 \, + \, 2 x_1 y_2 . 1/24 \, + \, 2 x_2 y_1 . 1/24 \, + \, 2 x_2 y_2 . 1/12 $$ The result is: $$ \overline{x y} = ( 2 x_1 y_1 + x_1 y_2 + x_2 y_1 + 2 x_2 y_2 ) / 12 $$ Substitute $x$ instead of $y$ herein: $$ \overline{x x} = ( 2 x_1 x_1 + x_1 x_2 + x_2 x_1 + 2 x_2 x_2 ) / 12 \hieruit \overline{x^2} = ( x_1 x_1 + x_1 x_2 + x_2 x_2 ) / 6 $$ Or the reverse: $y$ instead of $x$. Giving: $$ \overline{y^2} = ( y_1 y_1 + y_1 y_2 + y_2 y_2 ) / 6 $$ However, second order moments should be evaluated preferrably with respect to the midpoint: $$ \overline{x y} -\overline{x} \overline{y} = ( 2 x_1 y_1 + x_1 y_2 + x_2 y_1 + 2 x_2 y_2 ) / 12 - (x_1 + x_2) / 3 . (y_1 + y_2) / 3 $$ $$ = ( 6 x_1 y_1 + 3 x_1 y_2 + 3 x_2 y_1 + 6 x_2 y_2 - 4 x_1 x_1 - 4 x_1 y_2 - 4 x_2 y_1 - 4 x_2 y_2 ) / 36 $$ $$ = ( 2 x_1 x_1 - x_1 y_2 - x_2 y_1 + 2 x_2 y_2 ) / 36 $$ By proper substitution: $$ \overline{x^2} - \overline{x}^2 = ( x_1 x_1 - x_1 y_2 + x_2 y_2 ) / 18 \EN \overline{y^2} - \overline{y}^2 = ( y_1 y_1 - y_1 y_2 + y_2 y_2 ) / 18 $$ The question remains to be answered why these triangle integrals are supposed to be so important. Well, virtually every domain in the plane can be thought as being built up from (small) triangles. This means that every integral over such a domain can effectively be thought as a (weighted) sum of integrals over nothing else but triangles $(k)$: $$ \frac{\spekhaken x^m y^n \, dx dy}{\spekhaken dx dy} = \frac{\sum_k \left[ \spekhaken x^m y^n \, dx dy \right]_k} {\sum_k \left[ \spekhaken dx dy \right]_k } = \frac{\sum_k \left[ \spekhaken x^m y^n \, d\xi d\eta \: \J \right]_k} {\sum_k \left[ \spekhaken d\xi d\eta \: \J \right]_k } $$ $$ = \sum_k \left[ \spekhaken x^m y^n \, d\xi d\eta \right]_k . \: w_k \qquad \mbox{where} \quad w_k = \frac{ \J_k }{ \half \sum_i \J_i } $$ That is: the triangle moments are weighted with their individual areas, divided by the total area of the domain. This is the main reason why triangle moments are so useful: you can compose all other planar moments out of them (within certain accuracy bounds). Provided that you are not too punctilious, anything else will not be needed.

Third order moments, at last. It is handsom to evaluate only the following general expression in $(a,b,c)$ and specialize afterwards for $(a,b,c) = (x,x,x)$ , $(a,b,c) = (x,x,y)$, $(a,b,c) = (x,y,y)$, $(a,b,c) = (y,y,y)$ . Let's go: $$ \overline{a b c} = \spekhaken a(x,y)\, b(x,y)\, c(x,y) \, dx dy / \mbox{Area} $$ $$ = 2 \spekhaken (a_1 \xi + a_2 \eta )(b_1 \xi + b_2 \eta )(c_1 \xi + c_2 \eta ) \, d\xi d\eta $$ $$ =\, a_1 b_1 c_1 . 2 \spekhaken \xi^3 \eta^0 \, d\xi d\eta + a_1 b_1 c_2 . 2 \spekhaken \xi^2 \eta^1 \, d\xi d\eta + a_1 b_2 c_1 . 2 \spekhaken \xi^2 \eta^1 \, d\xi d\eta $$ $$ +\, a_1 b_2 c_2 . 2 \spekhaken \xi^1 \eta^2 \, d\xi d\eta + a_2 b_1 c_1 . 2 \spekhaken \xi^2 \eta^1 \, d\xi d\eta + a_2 b_1 c_2 . 2 \spekhaken \xi^1 \eta^2 \, d\xi d\eta $$ $$ +\, a_2 b_2 c_1 . 2 \spekhaken \xi^1 \eta^2 \, d\xi d\eta + a_2 b_2 c_2 . 2 \spekhaken \xi^0 \eta^3 \, d\xi d\eta $$ Using the formula $F(m,n) = m!\,n! / (m+n+2)!$ gives: $$ =\, a_1 b_1 c_1 . 2 . 1/20 + a_1 b_1 c_2 . 2 . 1/60 + a_1 b_2 c_1 . 2 . 1/60 + a_1 b_2 c_2 . 2 . 1/60 $$ $$ +\, a_2 b_1 c_1 . 2 . 1/60 + a_2 b_1 c_2 . 2 . 1/60 + a_2 b_2 c_1 . 2 . 1/60 + a_2 b_2 c_2 . 2 . 1/20 $$ Conclusion: $$ \half \; \overline{a b c} = (a_1 b_1 c_1 + a_2 b_2 c_2) / 20 $$ $$ + (a_1 b_1 c_2 + a_1 b_2 c_1 + a_1 b_2 c_2 + a_2 b_1 c_1 + a_2 b_1 c_2 + a_2 b_2 c_1)/60 $$ This has been implemented in ogenblik.pas as follows:

  function formule(a1,a2,b1,b2,c1,c2 : double) : double;
  var
    f : double;
  begin
    f := (a1*b1*c1 + a2*b2*c2)/20
       + (a2*b1*c1 + a1*b2*c1 + a1*b1*c2
        + a1*b2*c2 + a2*b1*c2 + a2*b2*c1)/60;
    formule := f;
  end;

A general remark is in place here. In order to prevent (too) large numbers to be involved with the calculation of higher order moments, we think it's a good idea to introduce, in advance, weighting factors for the coordinates. This can be accomplished by translating and rotating the coordinate system into its preferrred position, followed by dividing the global coordinates by their spreads in the $x$- and $y$-direction respectively: $$ x := \frac{x}{\sqrt{\sigma'_{xx}}} \EN y := \frac{y}{\sqrt{\sigma'_{yy}}} $$ Here $\sigma'_{xx}$ and $\sigma'_{yy}$ are the main moments of inertia.