How can I calculate the inertia tensor of a hollow object defined by a triangle mesh?

Question

I want to calculate the mass, center of mass, and inertia tensor of player created objects. Some of the objects will be hollow instead of solid. I am creating a closed triangle mesh for their object. Note that this is for a physics building game so I want the values to be as accurate as possible (within reason).

Based on the references listed below, I can get the properties for a solid object by creating a tetrahedron from each triangle of the mesh and calculating the signed volume which gives the mass (volume * density) and the center of mass. I even got the inertia tensor calculation working.

How do I do the same for a hollow object?

1. For the mass and center of mass, I am iterating through the triangles of the mesh and totaling their area and calculating the area-weighted average of their positions, then multiplying the surface area by a "density" value to get the mass.

    public static (float, Vector3) CalculateSurfaceArea(Mesh mesh)
    {
        Vector3[] vertices = mesh.vertices;
        int[] triangles = mesh.triangles;
   
        float totalArea = 0f;
        Vector3 temp = Vector3.zero;
   
        for (int i = 0; i <= triangles.Length - 3; i += 3)
        {
            Vector3 v1 = vertices[triangles[i]];
            Vector3 v2 = vertices[triangles[i + 1]];
            Vector3 v3 = vertices[triangles[i + 2]];
            Vector3 edge21 = v2 - v1;
            Vector3 edge31 = v3 - v1;
   
            float area = Vector3.Cross(edge21, edge31).magnitude / 2f; // area of the triangle
            totalArea += area;
   
            Vector3 triCenter = (v1 + v2 + v3) / 3f; // center of the triangle
            temp += triCenter * area;
        }
   
        Vector3 areaCenter = temp / totalArea;
   
        return (totalArea, areaCenter);
    }
   

2. For the inertia tensor, I am trying a similar approach where I iterate through all the triangles and total up their moments of inertia and using the parallel-axis theorem to account for their positions but (a) I am not sure this is correct and (b) how do I calculate the products of inertia (Ixy, Ixz, Iyz)?

    public static (Vector3, Quaternion) CalculateHollowInertiaTensor(Mesh mesh)
    {
        Vector3[] vertices = mesh.vertices;
        int[] triangles = mesh.triangles;
   
        double Ixx = 0f;
        double Iyy = 0f;
        double Izz = 0f;
        double Ixy = 0f;
        double Ixz = 0f;
        double Iyz = 0f;
   
        for (int i = 0; i <= triangles.Length - 3; i += 3)
        {
            Vector3 v1 = vertices[triangles[i]];
            Vector3 v2 = vertices[triangles[i + 1]];
            Vector3 v3 = vertices[triangles[i + 2]];
            Vector3 edge21 = v2 - v1;
            Vector3 edge31 = v3 - v1;
            Vector3 center = (v1 + v2 + v3) / 3f; // center of the triangle
            Vector3 offset = center - v1;
   
            float area = Vector3.Cross(edge21, edge31).magnitude / 2f; // area of the triangle
   
            // Moment of inertia of triangle rotating around the first vertex
            // https://en.wikipedia.org/wiki/List_of_moments_of_inertia
            // I = (1/6)m(P.P + P.Q + Q.Q)
            float triIxx = (edge21.y * edge21.y + edge21.z * edge21.z + edge21.y * edge31.y + edge21.z * edge31.z + edge31.y * edge31.y + edge31.z * edge31.z) / 6f;
            float triIyy = (edge21.x * edge21.x + edge21.z * edge21.z + edge21.x * edge31.x + edge21.z * edge31.z + edge31.x * edge31.x + edge31.z * edge31.z) / 6f;
            float triIzz = (edge21.x * edge21.x + edge21.y * edge21.y + edge21.x * edge31.x + edge21.y * edge31.y + edge31.x * edge31.x + edge31.y * edge31.y) / 6f;
   
            // Shift to the center of the triangle
            triIxx -= offset.y * offset.y + offset.z * offset.z;
            triIyy -= offset.x * offset.x + offset.z * offset.z;
            triIzz -= offset.x * offset.x + offset.y * offset.y;
   
            // Shift to the origin (using parallel-axis theorem)
            triIxx += center.y * center.y + center.z * center.z;
            triIyy += center.x * center.x + center.z * center.z;
            triIzz += center.x * center.x + center.y * center.y;
   
            Ixx += triIxx * area;
            Iyy += triIyy * area;
            Izz += triIzz * area;
            //Ixy += area * center.x * center.y;
            //Ixz += area * center.x * center.z;
            //Iyz += area * center.y * center.z;
        }
   
        Matrix<double> inertiaTensor = Matrix<double>.Build.Dense(3, 3);
        inertiaTensor[0, 0] = Ixx;
        inertiaTensor[1, 1] = Iyy;
        inertiaTensor[2, 2] = Izz;
        inertiaTensor[0, 1] = inertiaTensor[1, 0] = -Ixy;
        inertiaTensor[0, 2] = inertiaTensor[2, 0] = -Ixz;
        inertiaTensor[1, 2] = inertiaTensor[2, 1] = -Iyz;
   
        Debug.Log(inertiaTensor);
   
        //Matrix<double> inertiaTensorInverse = inertiaTensor.Inverse();
   
        // Find the principal axes and simplified inertia tensor
        Evd<double> evd = inertiaTensor.Evd();
        Vector3 inertiaTensorDiagonal = MakeVector3(evd.EigenValues.Real());
        Quaternion inertiaTensorRotation = Quaternion.LookRotation(
            MakeVector3(evd.EigenVectors.Column(2)),
            MakeVector3(evd.EigenVectors.Column(1))
            );
   
        return (inertiaTensorDiagonal, inertiaTensorRotation);
    }
   

References I have found so far:

• How to calculate the volume of a 3D mesh object the surface of which is made up triangles How to calculate the volume of a 3D mesh object the surface of which is made up triangles

• List of moments of inertia https://en.wikipedia.org/wiki/List_of_moments_of_inertia

• Polyhedral Mass Properties, David Eberly https://www.geometrictools.com/Documentation/PolyhedralMassProperties.pdf

• Explicit Exact Formulas for the 3-D Tetrahedron Inertia Tensor in Terms of its Vertex Coordinates, F. Tonon http://docsdrive.com/pdfs/sciencepublications/jmssp/2005/8-11.pdf

• Efficient Feature Extraction for 2D/3D Objects in Mesh Representation, Cha Zhang and Tsuhan Chen http://chenlab.ece.cornell.edu/Publication/Cha/icip01_Cha.pdf

===== EDIT =====

Based on John Alexiou's answer, I reimplemented the algorithm for hollow objects:

    public static (float, Vector3, Vector3, Quaternion) CalculateHollowMassCenterInertiaTensor(Mesh mesh)
    {
        Vector3[] vertices = mesh.vertices;
        int[] triangles = mesh.triangles;

        double Ixx = 0f;
        double Iyy = 0f;
        double Izz = 0f;
        double Ixy = 0f;
        double Ixz = 0f;
        double Iyz = 0f;

        float totalMass = 0f;
        Vector3 temp = Vector3.zero;

        for (int i = 0; i <= triangles.Length - 3; i += 3)
        {
            Vector3 v1 = vertices[triangles[i]];
            Vector3 v2 = vertices[triangles[i + 1]];
            Vector3 v3 = vertices[triangles[i + 2]];

            Vector3 center = (v1 + v2 + v3) / 3f; // center of the triangle
            float area = (Vector3.Cross(v1, v2) + Vector3.Cross(v2, v3) + Vector3.Cross(v3, v1)).magnitude / 2f; // area of the triangle

            totalMass += area;
            temp += center * area;

            Ixx += area * (v1.y * v1.y + v1.z * v1.z + v2.y * v2.y + v2.z * v2.z + v3.y * v3.y + v3.z * v3.z) / 3f;
            Iyy += area * (v1.x * v1.x + v1.z * v1.z + v2.x * v2.x + v2.z * v2.z + v3.x * v3.x + v3.z * v3.z) / 3f;
            Izz += area * (v1.x * v1.x + v1.y * v1.y + v2.x * v2.x + v2.y * v2.y + v3.x * v3.x + v3.y * v3.y) / 3f;
            Ixy += area * (v1.x * v1.y + v2.x * v2.y + v3.x * v3.y) / 3f;
            Ixz += area * (v1.x * v1.z + v2.x * v2.z + v3.x * v3.z) / 3f;
            Iyz += area * (v1.y * v1.z + v2.y * v2.z + v3.y * v3.z) / 3f;
        }

        Vector3 centerOfMass = temp / totalMass;

        Matrix<double> inertiaTensor = Matrix<double>.Build.Dense(3, 3);
        inertiaTensor[0, 0] = Ixx;
        inertiaTensor[1, 1] = Iyy;
        inertiaTensor[2, 2] = Izz;
        inertiaTensor[0, 1] = inertiaTensor[1, 0] = -Ixy;
        inertiaTensor[0, 2] = inertiaTensor[2, 0] = -Ixz;
        inertiaTensor[1, 2] = inertiaTensor[2, 1] = -Iyz;

        Debug.Log(inertiaTensor);

        // Find the principal axes and simplified inertia tensor
        Evd<double> evd = inertiaTensor.Evd();
        Vector3 inertiaTensorDiagonal = MakeVector3(evd.EigenValues.Real());
        Quaternion inertiaTensorRotation = Quaternion.LookRotation(
                MakeVector3(evd.EigenVectors.Column(2)),
                MakeVector3(evd.EigenVectors.Column(1))
            );

        return (totalMass, centerOfMass, inertiaTensorDiagonal, inertiaTensorRotation);
    }

It is giving me the wrong values for the inertia tensor though.

	mass	Expected Inertial Tensor	Actual Inertial Tensor
Cube	6.00	1.6667 1.6667 1.6667	3.0000 3.0000 3.0000
Sphere	3.11	0.5236 0.5236 0.5236	0.5151 0.5162 0.5163
Cylinder	7.80	4.5488 1.7671 4.5488	8.7134 1.8217 8.7134

Cube has 1m sides, Sphere has 0.5m radius, Cylinder has 0.5m radius and 2m height. All are centered around the origin.

Note that I am using the total surface area as the mass. Assume that the thickness * density = 1 so mass = area * thickness * density becomes mass = area * 1. I will be picking better values for those later because I assume that will not have that big of an effect on the algorithm.

Also note that some rounding is to be expected because the shapes are being approximated by a triangular mesh.

Answer 1

For a shell body to have mass and mass moment of inertia, the sides must have some thickness ε>0.

This defines the mass of a triangle defined by the vectors A, B and C and thickness ε and density ρ as

area = 1/2*Math.Abs( Vector3.Cross(A,B) + Vector3.Cross(B,C) + Vector3.Cross(C,A) )
mass = ρ*area*ε

or the above can be used to find the density from the total mass of an object, once the total surface area is calculated.

To find the mass moment of inertia you need to define a function returning the MMOI matrix of a unit mass particle from the location vector r=(x,y,z).

Matrix3 Mmoi(Vector3 r)
{
    //     | y^2+z^2    -x y      -x z   |
    // I = |  -x y     x^2+z^2    -y z   |
    //     |  -x z      -y z     x^2+y^2 |
    return new Matrix3(
        r.y*r.y + r.z*r.z, -r.x*r.y, -r.x*r.y,
        -r.x*r.y, r.x*r.x + r.z*r.y, -r.y*r.z,
        -r.x*r.z, -r.y*r.z, r.x*r.x + r.y*r.y);
}

and then calculate the MMOI of a triangle from the vertices and the mass

Matrix3 Mmoi(double m, Vector3 A, Vector3 B, Vector3 C)
{
    return (m/3)*(Mmoi(A)+Mmoi(B)+Mmoi(C));
}

_{The above is derived from the surface integral over the triangle, and since [SO] does not support math formatting I am omitting the details here.}

Yes the above is true, the MMOI of a surface triangle is that of the average MMOI of the three vertices.

Update

In correlating the above with CAD I realized you have to integrate over both front and back surfaces of the triangle. The result that matches CAD mass properties is

Matrix3 Mmoi(double m, Vector3 A, Vector3 B, Vector3 C)
{
    return 2*(m/3)*(Mmoi(A)+Mmoi(B)+Mmoi(C));
}