How to calculate inertia tensor and tensor rotation manually

Question

I'm trying to set the moment of inertia of a compound object (one object with a Rigidbody and multiple child objects with simple colliders). I need to set the Rigidbody.inertiaTensor and Rigidbody.inertiaTensorRotation manually because the child objects have different material densities.

I have calculated the inertia tensors for each child object around its center of mass as well as the center of mass of the parent object. I believe I can calculate the moment of inertia around the parent center of mass by adding all the individual moments of inertia of the child objects and then adding the moments of inertia of the children as point masses around the parent center of mass (might be wrong on this?).

What I don't understand is how to translate this into this Vector3 inertia tensor and Quaternion tensor rotation that Unity uses; the tensor rotation in particular I don't really understand what it represents. If I had the final moment of inertia of the compound object as a Matrix4x4, could I turn it into this Vector3 + Quaternion system? Is it possible to do these calculations while staying in the Vector3 + Quaternion space?

If nothing else, I want to know how the Rigidbody.inertiaTensorRotation is calculated and what it actually represents. The documentation says that "the inertia tensor is rotated by the inertiaTensorRotation", but I can't seem to find a relationship between the child transform rotation and it's effect on the inertiaTensorRotation. I've tried visualizing it by making a parent object with a Rigidbody with one child 1x1x1m cube and drawing the inertiaTensor vector rotated by the inertiaTensorRotation, but it just seems to jump around wildly as I slowly rotate the different axes of the box. The only pattern I've noticed is that if I only rotate 1 axis and all the others are 0, then the inertiaTensorRotation == Quaternion.identity alt text

Appreciate any help with this.

Answer 1

Unity relies on PhysX to calculate these values. You can inspect ExtRigidBodyExt.cpp and ExtInertiaTensor.h to find out exactly how it's done. The particular points of interest are:

computeMassAndInertia in ExtRigidBodyExt.cpp which orchestrates this process by iterating through all child shapes and combining theirs inertias.
Ext::InertiaTensorComputer::transform and Ext::InertiaTensorComputer::add methods in ExtInertiaTensor.h which are necessary in order to combine child inertias. This file also contains methods for calculating inertias of various primitive shapes.
Finally, inertiaTensor and inertiaTensorRotation are extracted from resulting matrix. This is done via physx::PxDiagonalize in PsMathUtils.cpp

Here is the relevant source for step 3. The return value of physx::PxDiagonalize is what's called inertiaTensor in Unity and inertiaTensorRotationis the massFrame which is the second argument.

 // indexed rotation around axis, with sine and cosine of half-angle
 PxQuat indexedRotation(PxU32 axis, PxReal s, PxReal c)
 {
     PxReal v[3] = { 0, 0, 0 };
     v[axis] = s;
     return PxQuat(v[0], v[1], v[2], c);
 }
 
 PX_FOUNDATION_API PxVec3 physx::PxDiagonalize(const PxMat33& m, PxQuat& massFrame)
 {
     // jacobi rotation using quaternions (from an idea of Stan Melax, with fix for precision issues)
 
     const PxU32 MAX_ITERS = 24;
 
     PxQuat q = PxQuat(PxIdentity);
 
     PxMat33 d;
     for(PxU32 i = 0; i < MAX_ITERS; i++)
     {
         PxMat33 axes(q);
         d = axes.getTranspose() * m * axes;
 
         PxReal d0 = PxAbs(d[1][2]), d1 = PxAbs(d[0][2]), d2 = PxAbs(d[0][1]);
         PxU32 a = PxU32(d0 > d1 && d0 > d2 ? 0 : d1 > d2 ? 1 : 2); // rotation axis index, from largest off-diagonal
         // element
 
         PxU32 a1 = shdfnd::getNextIndex3(a); a2 = shdfnd::getNextIndex3(a1);
         if(d[a1][a2] == 0.0f || PxAbs(d[a1][a1] - d[a2][a2]) > 2e6f * PxAbs(2.0f * d[a1][a2]))
             break;
 
         PxReal w = (d[a1][a1] - d[a2][a2]) / (2.0f * d[a1][a2]); // cot(2 * phi), where phi is the rotation angle
         PxReal absw = PxAbs(w);
 
         PxQuat r;
         if(absw > 1000)
             r = indexedRotation(a, 1 / (4 * w), 1.f); // h will be very close to 1, so use small angle approx instead
         else
         {
             PxReal t = 1 / (absw + PxSqrt(w * w + 1)); // absolute value of tan phi
             PxReal h = 1 / PxSqrt(t * t + 1);          // absolute value of cos phi
 
             PX_ASSERT(h != 1); // |w|<1000 guarantees this with typical IEEE754 machine eps (approx 6e-8)
             r = indexedRotation(a, PxSqrt((1 - h) / 2) * PxSign(w), PxSqrt((1 + h) / 2));
         }
 
         q = (q * r).getNormalized();
     }
 
     massFrame = q;
     return PxVec3(d.column0.x, d.column1.y, d.column2.z);
 }

Note that shdfnd::getNextIndex3 in the above code is simply (a+1)%3. If you perform all the above steps, you will get exact values of inertiaTensor and inertiaTensorRotation.

Answer 2

This is degree level math.

The simple answer is to learn about the backgrouns math of matrices, inertia tensors, and their properties. The thing that confuses this is unity are using a non standard (from math perspective) way of representing the data.

In a nutshell, an inertia tensor must be diagonalized before it can be used for calculation easily with physics, and this involves a complex algorithm for finding eigenvectors.

An Eigen vector could be described as a way to change the matrix so it is in a new basis of representation, in the case of inertia tensors, this is the principle axis of rotation. Once calculated, the eigen vectors provide the axis upon which the inertia tensor would "balance" well, and then enters the physics equations easier.

If you decided to take an object which balances well in the long axis, but represent its position diagonally across the centre of mass, such that one tip occupied the position 1,1 and the other tip -1, -1 then the diagonalization and eigenvector representation results in a new orientation to consider the motions from. It then enables the math to be performed on the object with its long axis nicely aligned to the axes in use. This is dependent on the shape and orientation of the object.

So, point is, an object may not wish to rotate about the axis we define for them, and, the result is to turn the axes on the object to be able to align it such that the axes are along the lines it does wish to rotate on. The physically dictated axes of natural balance are referred to as the principle axes of rotation.

It would appear that rather than exposing the complex and somewhat irrelevant math(from the perspective of users of unity) that is used to get to this point, they have nicely tidied up the result into a rotation quaternion, such that if you change the local axes of the object by the quaternion, you have the principle axes of rotation.

I assume your research into inertia tensors found enough to make an inertia tensor, but the rotation part will be very difficult. You may however, be able to define the inertia tensor and have unity create the quaternion for you.

The quaternion is best described as a reformatting of the standard output from this sequence of math, usually you get a new set of basis vectors in the form of a 3x3 matric, calculated from the inertia tensor, but, of course, a 3x3 basis matrix can be made to represent a rotation, and then translated into that quaternion. The quaternion is a seperate derivation of the same data.

How to calculate inertia tensor and tensor rotation manually

2 Replies

Your answer

Follow this Question

Related Questions