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).

$\iiint_V\left(\mathbf{\nabla}\cdot\mathbf{F}\right)\,dV=$ $\oiint$ $\scriptstyle 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.x², 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).x²+(B+C).x²+(C+A).x²) for the centroid in x. (Typically in math the notation to generalize across the coordinate axis uses unit coordinate vectors e₁, e₂, and e₃, in general e_d, so we have F(P)=1/2 (P.e_d)².)

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".