First the code then the discussion (if you care).  From the source file volint.cpp.   Note that for the code on this page here, the input parameters here are list of vertices and indexed list of triangles. So every int3 represents a triangle on the mesh. 

Here is how to compute Volume:

 
 float Volume(const float3 *verts, const int3 *tris, const int count) 
 {
 // count is the number of triangles (tris) 
 float vol=0;
 for(int i=0; i < count; i++) // for each triangle
 {
 vol+=Determinant(float3x3(verts[tris[i][0]],verts[tris[i][1]],verts[tris[i][2]])); 
 }
 return vol/6.0f; // since the determinant give 6 times tetra volume
 }
 
 
 

Here is how to compute Center Of Mass:

 
 float3 CenterOfMass(const float3 *vertices, const int3 *tris, const int count) 
 {
 // count is the number of triangles (tris) 
 float3 com(0,0,0);
 float volume=0; // actually accumulates the volume*6
 for(int i=0; i < count; i++) // for each triangle
 {
 float3x3 A(vertices[tris[i][0]],vertices[tris[i][1]],vertices[tris[i][2]]); 
 float vol=Determinant(A); // dont bother to divide by 6 
 com += vol * (A.x+A.y+A.z); // divide by 4 at end
 volume+=vol;
 }
 com /= volume*4.0f; 
 return com;
 }

 
 
 

Here is how to compute Inertial Tensor:

 
 float3x3 Inertia(const float3 *vertices, const int3 *tris, const int count, 

 const float3& com=float3(0,0,0)) 
 {
 // count is the number of triangles (tris) 
 // The moments are calculated based on the center of rotation com 

 // com if not provided is assumed to be [0,0,0] by default
 // assume mass==1.0 you can multiply by mass later.
 // for improved accuracy the next 3 variables, the determinant d, 

 // and its calculation should be changed to double
 float volume=0; // technically this variable accumulates the volume times 6
 float3 diag(0,0,0); // accumulate matrix main diagonal integrals [x*x, y*y, z*z]
 float3 offd(0,0,0); // accumulate matrix off-diagonal integrals [y*z, x*z, x*y]
 for(int i=0; i < count; i++) // for each triangle
 {
 // matrix trick for volume calc by taking determinant
 float3x3 A(vertices[tris[i][0]]-com,

 vertices[tris[i][1]]-com,

 vertices[tris[i][2]]-com); 
 // vol of tiny parallelapiped= d * dr * ds * dt 

 // where dr,ds,dt are the 3 partials of my tetral triple integral equasion
 float d = Determinant(A); 
 volume +=d; // add vol of current tetra (note it could be negative - that's ok)
 for(int j=0;j < 3;j++)
 {
 int j1=(j+1)%3; 
 int j2=(j+2)%3; 

 // defer the division of diag by 60 and offd by 120 until after loop
 diag[j] += (A[0][j]*A[1][j] + A[1][j]*A[2][j] + A[2][j]*A[0][j] + 
 A[0][j]*A[0][j] + A[1][j]*A[1][j] + A[2][j]*A[2][j] ) *d; 
 offd[j] += (A[0][j1]*A[1][j2] + A[1][j1]*A[2][j2] + A[2][j1]*A[0][j2] +
 A[0][j1]*A[2][j2] + A[1][j1]*A[0][j2] + A[2][j1]*A[1][j2] +
 A[0][j1]*A[0][j2]*2+ A[1][j1]*A[1][j2]*2+ A[2][j1]*A[2][j2]*2 ) *d; 

 }
 }
 // divide by total volume (vol/6) since density=1/volume
 diag /= volume*(60.0f /6.0f); 

 offd /= volume*(120.0f/6.0f);
 return float3x3(diag.y+diag.z , -offd.z , -offd.y,
 -offd.z , diag.x+diag.z, -offd.x,
 -offd.y , -offd.x , diag.x+diag.y );
 }
 
 

The above code is useful when preprocessing geometry for physical simulation and rigid body dynamics. It shouldn’t be too hard to cut and paste the above. A couple of things to keep in mind:

• The inertial properties are calculated based on an object of uniform mass of 1. Note this is not the same as assuming a density of 1. All you have to do is multiply the resulting matrix by your objects mass to scale it appropriately. Alternatively, to automatically make big things heavier and small things lighter, your application might specify things by density instead of mass. In that case, calculate the volume and multiply by that if you want.
• Volume calculations can be done from any reference point. However, the contribution to inertia by a small volume is proportional to the distance squared between that volume and the axis of rotation. Hence the importance of computing the summation about the right point. The inertial properties are based on a center of rotation which by default the routine assumes is 0,0,0 (the local origin). There is a way to translate and rotate inertia tensors but its not the same as transforming a regular matrix. Doing this is useful when combining a number of shapes into an aggregate rigidbody. In the case of a single shape, its easiest to just compute the center of mass first and use that when computing the inertia.

The implementations here are simple, fast, and have been written in such a way that hopefully will be easy for you to quickly grab and adapt to your own math library and get used properly by your game code.

The standard set of routines for computing polyhedral mass properties has been based on Brian Mirtich’s JGT paper and the associated source code provided. His approach, based on Green’s lemma, solves the volume integration by solving the boundary integral. The implementation presented here takes a more direct approach by just doing the volume integration over the mesh. By solving the equations of interest for the general tetrahedron, we can then sum results from the implicit tetrahedron for each triangle to integrate the entire mesh volume, provided it’s a proper manifold.

The formulas in the code are all based on tetrahedral volume summation. The vertices of triangle along with a reference point (origin) make up a tetrahedron. When the winding is opposite, then the tetrahedron is inside out and thus makes a negative contribution to the sum. Otherwise, this code couldn't work on anything other than convex polyhedra with the origin inside the solid.

Its not likely a problem, but since this is likely to be preprocessing, it is best to use double precision instead of single. Numerical precision is an issue when you are summing lots of numbers whose individual magnitudes can be much higher than the final result. For example calculating 1000000 – 999999 requires that you have at least 7 decimals of precision so it would be bad to use short integers here. Similarly, if you have many tiny triangles far away from the origin, then it might be wise to use double precision. To give a ballpark where things can start to break, a unit cube (1x1x1) placed 10000 units away from the origin, will not calculate the center of mass correctly (last decimal place off by one) when 32 bit precision is used (regular floats).