Inertia tensor in Matrix Form from InertiaTensor and InertiaTensorRotation

Hello,

Inertia tensor is normally defined as a 3x3 Matrix

Unity and Physx give something equivalent, an IntertiaTensor, and an InertiaTensorRotation. I think they do it like this because it is cheaper (3+4 = 7 values) against (3x3=9 values).

My question is, having an InertiaTensor and an InertiaTensorRotation, how can I calculate the 3x3 Inertia Tensor?

I reckon this should be quite straightforward but I can't get the maths right

I am tempted to do:

Code (CSharp):

1. var intertiaMatrix = Matrix4x4.TRS(Vector3.One, inertiaTensorRotation, inertiaTensor)

Because TRS transforms the quaternion in a canonical form and it is scaled by the inertia tensor

Thank you

 

Did you ever figure this out? I'm new to Unity and have a similar requirement.

 

No, not really. But this is something I would like to study as soon as I can

 

I also tried to get more information on the inertia tensor rotation, but even the PhysX devs seem reluctant to talk about this.

I've been able to convert an inertia tensor matrix to the equivalent inertia tensor vector and inertia tensor rotation by translating the PhysX code, but haven't found the formula that does the opposite.

 

Yes actually it is a very opaque topic from the PhysX's side

I could understand that PhysX decided to go to two structures (vector and quaternion) instead of a 3x3 matrix in order to save a little of memory, although this change could increase the number of calculations and operations.

Anyway, I thought the maths from vector + quaternion to 3x3 should be straightforward, but when I did it manually, the results where different, so I think I am missing something here.

 

It's more than that. Actually it seems part of a heavily simplified rigidbody dynamics model.

In reality (and realistic rb dynamics models) a rotating object has angular momentum. The angular velocity results of dividing the angular momentum by the inertia [matrix]. Therefore, modifying the inertia actually changes the angular velocity.

But in PhysX's simplified model this doesn't happen. Instead of integrating angular momentum, they integrate angular velocity. So when a torque is applied, the resulting change in the angular velocity is calculated by dividing the torque by the inertia. Further changes in the inertia don't modify the angular velocity. Only when a new torque is applied the angular velocity changes based on the current inertia value. Therefore, freely rotating objects (no forces/torques applied) exhibit different behaviors in PhysX compared with the real world (and with realistic models).

While working in Project 424 we found severe differences when applying torques in PhysX using a non-identity inertia tensor rotation, compared with the same exact situation in an high-end dynamics simulation model. We haven't investigated further yet, but the preliminary results point that the combination of inertia tensor vector/rotation provide very different effects than the expected from the original inertia matrix.

 

I suspect the values of the Vector3 are the principal moments of inertia and the quaternion describing the differential rotation between the local coordinate system and the axes of the principal moments of inertia. There I would also guess that creating the rotation matrix with Matrix4x4.TRS is the right way to go.

Did you dig deeper into this @Edy? I'm not well practised in doing matrix transforms, but you should do something like T.inverse * X * T (or the other way round?). It's also not clear to me in which coordinate system the matrix you suggested in your initial post lives.

As for the differences between PhysX and more sophisticated models: I assume from a small test a while ago that Euler's equations are not respected. https://en.wikipedia.org/wiki/Euler's_equations_(rigid_body_dynamics) meaning that the behaviour in the following video will not show in PhysX.

(skip to 0:20)

 

That where I started to think about this.
Simulating the bicycle equilibrium, I got very different outputs. Simulating in Matlab with a 3x3 Inertia Matrix, I got the bicycle stabilized, but I couldn't reproduce that behavior with the vector/quaternion model in Unity

So the classic example about a figure skater increasing his angular velocity decreasing his moment of inertia is not possible in Physx, isn't? Very important point, I didn't know that

 

Exactly, that doesn't work. Modifying the inertia alone doesn't cause any effect in the rigidbody.

Maybe, just maybe, that example could be simulated explicitly by computing the angular momentum, modifying the inertia, then compute and apply a new angular speed out of the angular momentum and the new inertia.

 

We already thought about that, but no, the results were different. We don't know how to reconstruct the original inertia matrix out of the vector + rotation.

The gyroscopic effect is a different thing, but no, it's not simulated in PhysX either.

 

The page 67 of this presentation from Stan Melax at GDC2014 describes the actual equivalency:

​

In our case, S would be the inertia tensor matrix, D a diagonal matrix with the inertia tensor vector, and R the matrix-equivalent of the inertia tensor rotation.

 

@Edy does this imply that there is no reversal process? Or is there an algorithm for working it back? I'm currently trying to calculate the dynamic effective mass of a rigidbody when applying forces at a specific point like detailed in this answer, but I think I need the matrix inertia tensor. Or would there be a way of using the diagonal tensor + rotation to obtain the result to this formula?

 

This topic has gone back and forward for a while but I think with the slide shared by @Edy we can put all the pieces together.

The Inertia Matrix is a Symmetric 3x3 matrix. In Unity, instead of a matrix, we have a quaternion and a vector.

Physx diagonalizes the Inertia Matrix.This means that a symmetric matrix is decomposed in M = A*D*inv(A)

This corresponds to the eigenvectors and d with the eigenvalues.

When the symmetric matrix corresponds to a inertia matrix, I = R*D*inv(R). In the rigidbody, InertiaTensor vector is the values of the diagonal of D. inertiaTensorRotation is the equivalent quaternion of the rotation matrix R

Getting inertiaTensor and inertiaTensorRotation from the Inertia Matrix:
We need to decomposed the matrix to the form I = R*D*inv(R). R corresponds to the eigenvectors and D to the eigenvalues so any algorithm could work. Jacobi algorithm, for example. When we have R, transform R to a quaternion to get the inertiaTensorRotation, and get the diagonal values of D to get inertiaTensor

Getting the inertia tensor matrix from inertiaTensor and inertiaTensorRotation:
Transform inertiaTensorRotation from a quaternion to a rotation matrix. Create a diagonal matrix with inertia tensor in the diagonal. Apply the formula I = R*D*inv(R).

I created this helper class. I hope it helps:

Code (CSharp):

1. ///   Author: Manuel Espino.
2. ///   github: https://github.com/Maesla
3. ///   unity forum: https://forum.unity.com/members/maeslezo.863356/
4. ///   Utilities to get the inertia tensor matrix from a rigidbody
5. ///   and to get the inertia tensor and inertia tensor rotation from a inertia tensor matrix
6. ///-----------------------------------------------------------------
7.  
8. using UnityEngine;
9.  
10. public static class InertiaTensorUtils
11. {
12.     public static Matrix4x4 CalculateInertiaTensorMatrix(Rigidbody rb)
13.     {
14.         return CalculateInertiaTensorMatrix(rb.inertiaTensor, rb.inertiaTensorRotation);
15.     }
16.  
17.     // Inertia Tensor Matrix can be decomposed in M = transpose(R)*D*R
18.     // M is the original matrix
19.     // R is a rotation matrix, stored in the rigidbody as a quaternion in inertiaTensorRotation
20.     // D is a diagonal matrix, stored in the rigidbody as a vector3 in inertiaTensor
21.     // D are the eigenvalues and R are the eigenvectors
22.     // Inertia Tensor Matrix is a 3x3 Matrix, so it will appear in the first 3x3 positions of the 4x4 Unity Matrix used here
23.     public static Matrix4x4 CalculateInertiaTensorMatrix(Vector3 inertiaTensor, Quaternion inertiaTensorRotation)
24.     {
25.         Matrix4x4 R = Matrix4x4.Rotate(inertiaTensorRotation); //rotation matrix created
26.         Matrix4x4 S = Matrix4x4.Scale(inertiaTensor); // diagonal matrix created
27.         return R * S * R.transpose; // R is orthogonal, so R.transpose == R.inverse
28.     }
29.  
30.     private const float epsilon = 1e-10f;
31.     private const int maxSweeps = 32;
32.  
33.     /// <summary>
34.     /// Diagonalization of M
35.     /// </summary>
36.     /// <remarks>
37.     /// M will be decomposed by M = transpose(R)*D*R.
38.     /// InertiaTensorQuaternion is the quaternion equivalent to R
39.     /// InertiaTensor is the diagonal values of D
40.     /// InertiaTensor stores the eigenvalues and R stores the eigenvectors. Since the eigenvectors are the rotation axis, the quaternion representing R is the rotation axis
41.     /// </remarks>
42.     /// <param name="m"></param>
43.     /// <param name="inertiaTensor"></param>
44.     /// <param name="inertiaTensorRotation"></param>
45.     public static void DiagonalizeInertiaTensor(Matrix4x4 m, out Vector3 inertiaTensor, out Quaternion inertiaTensorRotation)
46.     {
47.         float m11 = m[0, 0];
48.         float m12 = m[0, 1];
49.         float m13 = m[0, 2];
50.         float m22 = m[1, 1];
51.         float m23 = m[1, 2];
52.         float m33 = m[2, 2];
53.  
54.         Matrix4x4 r = Matrix4x4.identity;
55.         for (int a = 0; a < maxSweeps; a++)
56.         {
57.             // Exit if off.diagonal entries small enough
58.             if ((fabs(m12) < epsilon) && (fabs(m13) < epsilon) && (fabs(m23) < epsilon))
59.                 break;
60.  
61.             // Annihilate (1,2) entry.
62.             if (m12 != 0.0F)
63.             {
64.                 float u = (m22 - m11) * 0.5F / m12;
65.                 float u2 = u * u;
66.                 float u2p1 = u2 + 1.0F;
67.                 float t = (u2p1 != u2) ?
68.                 ((u < 0.0F) ? -1.0F : 1.0F) * (sqrt(u2p1) - fabs(u))
69.                 : 0.5F / u;
70.                 float c = 1.0F / sqrt(t * t + 1.0F);
71.                 float s = c * t;
72.                 m11 -= t * m12;
73.                 m22 += t * m12;
74.                 m12 = 0.0F;
75.                 float temp = c * m13 - s * m23;
76.                 m23 = s * m13 + c * m23;
77.                 m13 = temp;
78.                 for (int i = 0; i < 3; i++)
79.                 {
80.                     float tempInner = c * r[i, 0] - s * r[i, 1];
81.                     r[i, 1] = s * r[i, 0] + c * r[i, 1];
82.                     r[i, 0] = tempInner;
83.                 }
84.             }
85.  
86.             // Annihilate (1,3) entry.
87.             if (m13 != 0.0F)
88.             {
89.                 float u = (m33 - m11) * 0.5F / m13;
90.                 float u2 = u * u;
91.                 float u2p1 = u2 + 1.0F;
92.                 float t = (u2p1 != u2) ?
93.                 ((u < 0.0F) ? -1.0F : 1.0F) * (sqrt(u2p1) - fabs(u))
94.                 : 0.5F / u;
95.                 float c = 1.0F / sqrt(t * t + 1.0F);
96.                 float s = c * t;
97.                 m11 -= t * m13;
98.                 m33 += t * m13;
99.                 m13 = 0.0F;
100.                 float temp = c * m12 - s * m23;
101.                 m23 = s * m12 + c * m23;
102.                 m12 = temp;
103.                 for (int i = 0; i < 3; i++)
104.                 {
105.                     float tempInner = c * r[i, 0] - s * r[i, 2];
106.                     r[i, 2] = s * r[i, 0] + c * r[i, 2];
107.                     r[i, 0] = tempInner;
108.                 }
109.             }
110.  
111.             // Annihilate (2,3) entry.
112.             if (m23 != 0.0F)
113.             {
114.                 float u = (m33 - m22) * 0.5F / m23;
115.                 float u2 = u * u;
116.                 float u2p1 = u2 + 1.0F;
117.                 float t = (u2p1 != u2) ?
118.                 ((u < 0.0F) ? -1.0F : 1.0F) * (sqrt(u2p1) - fabs(u))
119.                 : 0.5F / u;
120.                 float c = 1.0F / sqrt(t * t + 1.0F);
121.                 float s = c * t;
122.                 m22 -= t * m23;
123.                 m33 += t * m23;
124.                 m23 = 0.0F;
125.                 float temp = c * m12 - s * m13;
126.                 m13 = s * m12 + c * m13;
127.                 m12 = temp;
128.                 for (int i = 0; i < 3; i++)
129.                 {
130.                     float tempInner = c * r[i, 1] - s * r[i, 2];
131.                     r[i, 2] = s * r[i, 1] + c * r[i, 2];
132.                     r[i, 1] = tempInner;
133.                 }
134.             }
135.         }
136.  
137.         inertiaTensor.x = m11;
138.         inertiaTensor.y = m22;
139.         inertiaTensor.z = m33;
140.  
141.         inertiaTensorRotation = r.rotation;
142.     }
143.  
144.     private static float fabs(float f)
145.     {
146.         return Mathf.Abs(f);
147.     }
148.  
149.     private static float sqrt(float f)
150.     {
151.         return Mathf.Sqrt(f);
152.     }
153. }
154.  

References:

• Mathematics for 3D Game Programming and Computer Graphics. Eric Lengyel
  • 3.6 Diagonalization
  • 16.3 Eigenvalues and Eigenvectors
• https://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm
• https://en.wikipedia.org/wiki/Diagonalizable_matrix
• https://en.wikipedia.org/wiki/Moment_of_inertia#Principal_axes

Final notes:

• DiagonalizeInertiaTensor maybe is not the fastest algorithm. Anyway, it should be use just one time to cache the matrix. You can find another version, probably faster, here: https://answers.unity.com/questions/1484654/how-to-calculate-inertia-tensor-and-tensor-rotatio.html
• I am using the builtin Matrix4x4 matrix, but the InertiaTensor is 3x3
  
• Not sure why in the slide is written D = R*S*inv(R) instead of S = R*D*inv(R)
• Sometimes you may find transpose(R) instead of inverse(R).. They are equivalent because R is orthogonal, so inverse(R) = transpose(R)

 

@Maeslezo Thank you so much for putting that detailed reply together! It's been extremely helpful!

 

Is the resulting matrix in the reference frame of the rigidbody. Could this be used to combine the inertia tensors of two or more objects? I'd like to combine two objects that are jointed together to determine scaling to apply my pd controller that will torque about one of the combined bodies.

 

Yes, the matrix is in local reference.
If you want to use it in world reference, you have to rotate the matrix

Code (CSharp):

1. Matrix4x4 rotationMatrix = Matrix4x4.Rotate(rb.transform.rotation);
2. Matrix4x4 inertiaTensorLocal = rotationMatrix * InertiaTensor * rotationMatrix.transpose;

 

gist created:
https://gist.github.com/Maesla/7b3cffbda7d0a5b02aa7166d3eed5def

 

Hello, and thank you @Maeslezo and others for your work on this topic.

I'm looking to manually calculate the inertia tensor and rotation given any rigidbody. I see all the formulas you've come up with for going the other way, and I've also taken the time to translate PhysX's Diagonalization code (taken from here: https://answers.unity.com/questions/1484654/how-to-calculate-inertia-tensor-and-tensor-rotatio.html) into c# so I can do that much myself in Unity.

However, I lack the 3x3 matrix that needs to go into that code. I'm assuming this is the body frame inertia matrix for the rigidbody. Any idea how Unity is getting this? I assume, based on their documentation, that is has to do with the compound collider of the rigidbody in some way.

Any insight here would be GREATLY appreciated. Thank you in advance!

 

Here's a list of inertia tensors for multiple shapes:
https://en.wikipedia.org/wiki/List_of_moments_of_inertia

If you have a rigidbody with multiple colliders -each with their own inertia tensor- you can combine them into a single inertia tensor. This thread will help: https://gamedev.net/forums/topic/361031-adding-inertia-tensors-together/3372851/

For arbitrary convex hulls, you have to calculate it by hand: decompose the hull into a list of tetrahedrons, and then combine their inertia tensors.

And if for some reason you want the inertia tensor of an arbitrary (possibly concave) closed volume, you have to perform a convex decomposition, find the inertia tensor of each convex piece, and then combine them. Tetrahedralizing the volume and combining the inertia tensors of all tetrahedrons is one way of doing this.

 

Thankfully I'm just doing a disc-shaped object, so only one or two formulas were needed, and I got it working!! Thank you for the insight.