Boundary Integration and the Rotational Inertia Matrix

CS 482 Lecture, Dr. Lawlor

Torque, a vec3, directly affects angular momentum, but the conversion between the vec3 angular momentum (L) and the vec3 angular velocity (W) happens via the 3x3 rotational inertia matrix:
    L = I W

This 3x3 inertia matrix consists of 9 scalar values:
  [  Ixx  Ixy  Ixz  ]
[ Iyx Iyy Iyz ]
[ Izx Izy Izz ]
If we set the origin of our coordinate system to the center of mass, Ixz is the mass-weighted integral of x*z inside the object.

Calculating volume integrals on polygon meshes

We often want to compute a property of a volume, like the moment of inertia above, or mass, when all we have is a boundary representation such as a polygon mesh.  In these cases we often apply the Divergence Theorem (this is a 3D specialization of Stokes' Theorem, also known in 2D as Green's Theorem).

          S(\mathbf{F}\cdot\mathbf{n})\,dS .

The domain of the left integral V is any volume in 3D with a piecewise smooth boundary.  F is any continuously differentiable function returning a vector for any point in space.  ∇ . F gives the vector field divergence, the scalar dF.x/dx+dF.y/dy+dF.z/dz.  S is the surface of the volume.  F . n gives the dot product of the vector field with the outward-pointing surface normal n with unit length. 

Mass Example

For example, say we'd like to calculate the mass of a mesh with constant density r, by integrating r across the volume of the mesh V.  We pick F(P) = vec3(r*P.x,0,0) so that ∇ . F = r everywhere (we can pick any axis for F, or any combination of axes for F as long as we're consistent). 

Now to calculate the volume of the mesh, the divergence theorem tells us we can integrate
F . n across the surface of the mesh.  For polygons n is constant, so for a face with a normal lying along the x axis, integrating F . n is equivalent to multiplying F . n by the area of the face. 

For instance, to use this to calculate the mass of a cube with one vertex at the origin, and sides w,l, and h, we can ignore four of the six faces since their normal is orthogonal to F so
F . n is zero.  Of the remaining two faces, one is at the origin so F is zero, and the other has area l*h and F.x=r*P.x=r*w throughout.  This means the integral of F . n across the face is l*h*w*r, which is indeed the mass of a cube with density r.

Now let's calculate the mass of the volume whose boundary is given by a triangle mesh.  Each triangle with vertices A, B, C has this constant unit surface normal:
     n = normalize(cross(B-A,C-A))
Because n is constant, we can pull it out of the integral of F . n, giving n . integral of F over the triangle.  F is not constant over the triangle surface, but because F is linear, it is equivalent to evaluate it just once at the triangle's barycenter, which has coordinates (A+B+C)/3.  The integral is thus F((A+B+C)/3) = r*(A+B+C).x/3, times the surface normal n = normalize(cross(B-A,C-A)), times the area of the triangle, length(cross(B-A,C-A))/2.  Because the same vector is normalized and then multiplied by its own length, this is equivalent to:
   integral of F . n = 1/6 dot(vec3(
r*(A+B+C).x,0,0), cross(B-A,C-A))

If we sum this across each triangle, we indeed get the mass of the volume, including cases where the boundary includes convex and concave features crossing the X axis, and even for cases where the boundary has interior holes.  However, this will not calculate a consistent value for a mesh that is not closed, or has self-intersections.

Centroid and Inertia

We can compute the centroid in a similar fashion, using F(P)=1/2 P.x2, and similar formulas for y and z.  However, we must still calculate the integral of F across the triangle.  In class, we got lost in the details of how to do this via change of variables to barycentric coordinates.  But there is a useful theorem that states for a quadratic function like F, it is equivalent to evaluate the integral at the midpoints of the triangle edges, which yields a simple evaluation of F=1/24*n.x*((A+B).x2+(B+C).x2+(C+A).x2) for the centroid in x.  (Typically in math the notation to generalize across the coordinate axis uses unit coordinate vectors e1, e2, and e3, in general ed, so we have F(P)=1/2 (P.ed)2.)

For the details of one approach to generalize this calculation to the full inertia matrix, see Mirtich's 1996 paper "Fast and Accurate Computation of Polyhedral Mass Properties".