Search Unity

  1. Welcome to the Unity Forums! Please take the time to read our Code of Conduct to familiarize yourself with the forum rules and how to post constructively.
  2. We are updating our Terms of Service for all Unity subscription plans, effective October 13, 2022, to create a more streamlined, user-friendly set of terms. Please review them here: unity.com/legal/terms-of-service.
    Dismiss Notice
  3. Have a look at our Games Focus blog post series which will show what Unity is doing for all game developers – now, next year, and in the future.
    Dismiss Notice
  4. Join us on Thursday, September 29, for a day with Unity's SRP teams here on the forum or on Reddit, and discuss topics around URP, HDRP, and the Scriptable Render Pipeline in general.
    Dismiss Notice

Inertia tensor in Matrix Form from InertiaTensor and InertiaTensorRotation

Discussion in 'Physics' started by Maeslezo, Jul 28, 2020.

  1. Maeslezo

    Maeslezo

    Joined:
    Jun 16, 2015
    Posts:
    239
    Hello,

    Inertia tensor is normally defined as a 3x3 Matrix
    upload_2020-7-28_14-2-20.png
    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
     
    Last edited: Jul 28, 2020
  2. myelin_price

    myelin_price

    Joined:
    Oct 24, 2019
    Posts:
    1
    Did you ever figure this out? I'm new to Unity and have a similar requirement.
     
  3. Maeslezo

    Maeslezo

    Joined:
    Jun 16, 2015
    Posts:
    239
    No, not really. But this is something I would like to study as soon as I can
     
  4. Edy

    Edy

    Joined:
    Jun 3, 2010
    Posts:
    2,189
    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.
     
  5. Maeslezo

    Maeslezo

    Joined:
    Jun 16, 2015
    Posts:
    239
    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.
     
    Last edited: Dec 16, 2020
  6. Edy

    Edy

    Joined:
    Jun 3, 2010
    Posts:
    2,189
    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.
     
  7. tjmaul

    tjmaul

    Joined:
    Aug 29, 2018
    Posts:
    397
    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)
     
  8. Maeslezo

    Maeslezo

    Joined:
    Jun 16, 2015
    Posts:
    239
    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
     
  9. Edy

    Edy

    Joined:
    Jun 3, 2010
    Posts:
    2,189
    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.
     
  10. Edy

    Edy

    Joined:
    Jun 3, 2010
    Posts:
    2,189
    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.
     
  11. Edy

    Edy

    Joined:
    Jun 3, 2010
    Posts:
    2,189
    The page 67 of this presentation from Stan Melax at GDC2014 describes the actual equivalency:

    upload_2021-9-23_11-15-53.png

    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.
     
    Last edited: Sep 23, 2021
    PedroCoriAG, tjmaul and Maeslezo like this.
  12. PedroCoriAG

    PedroCoriAG

    Joined:
    Nov 1, 2016
    Posts:
    7
    @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?
     
  13. Maeslezo

    Maeslezo

    Joined:
    Jun 16, 2015
    Posts:
    239
    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:
    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)
     
    Nition and PedroCoriAG like this.
  14. PedroCoriAG

    PedroCoriAG

    Joined:
    Nov 1, 2016
    Posts:
    7
    @Maeslezo Thank you so much for putting that detailed reply together! It's been extremely helpful!
     
  15. Cloudwalker_

    Cloudwalker_

    Joined:
    Jan 3, 2014
    Posts:
    84
    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.
     
  16. Maeslezo

    Maeslezo

    Joined:
    Jun 16, 2015
    Posts:
    239
    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;
     
  17. Maeslezo

    Maeslezo

    Joined:
    Jun 16, 2015
    Posts:
    239
unityunity