By
Simple easy code for:
First we give you the code, adapted from file volint.cpp, followed by the theory (if you even care). 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 *vertices, const int3 *tris, const int count)
{
// count is the number of triangles (tris)
float volume=0;
for(int i=0; i < count; i++) // for each triangle
{
volume += Determinant(float3x3(vertices[tris[i][0]],vertices[tris[i][1]],vertices[tris[i][2]])); //divide by 6 later for efficiency
}
return volume/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 which is [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
{
float3x3 A(vertices[tris[i][0]]-com,vertices[tris[i][1]]-com,vertices[tris[i][2]]-com); // matrix trick for volume calc by taking determinant
float d = Determinant(A); // vol of tiny parallelapiped= d * dr * ds * dt (the 3 partials of my tetral triple integral equasion)
volume +=d; // add vol of current tetra (note it could be negative - that's ok we need that sometimes)
for(int j=0;j < 3;j++)
{
int j1=(j+1)%3;
int j2=(j+2)%3;
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; // divide by 60.0f later;
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 120.0f later
}
}
diag /= volume*(60.0f /6.0f); // divide by total volume (vol/6) since density=1/volume
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:
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).
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.
There is also an algorithm with source code to diagonalize a symmetric 3x3 which can be used on the inertia tensor matrix if your physics engine expect that or if you just want to find the principle axes of your object.
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 implementations are simple, fast, and have been written in such a way that should make it easy for you to grab and incorporate in your own code quickly.
At the time of developing this source code, long time game developer Jonathan Blow had also had developed a similar implementation. It should be posted on the web soon. It too uses a volume tetrahedron integration approach. Instead of solving the general tetrahedron, Jonathan’s approach was to solve a unit right-angle tetrahedron and then apply transforms to reshape it into a general tetrahedron. Even though, the formulas were derived and constructed differently, they all produce the same numerical result. Jonathan has an intuitive explanation for his covariant mass matrix that is worth reading for anyone interested in trying to get a deeper understanding.