Sign in to follow this  

Rigid Body Physics - Inertia Tensor Space

This topic is 850 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

I implemented simple rigid body system (simulating a simple box in 3D) and am confused about the space the inertia tensor should be, or to be more precise, how to compute the proper space.

I'll just dive into pretty simple code (Python).

Force addition to rigid body:

def AddForce(self, force, point):

    self.totalForce = self.totalForce + force

    oriMat = quaternion.ToMatrix(self.ori)
    r = point - self.pos
    torque = vector.Cross(r, force)
    self.totalTorque = self.totalTorque + torque
Note here that ori is orientation quaternion. Point, self.pos and forces are all provided in world space.

Now integration:

def Integrate(self, deltaTime):

    linAcc = 1.0/self.mass * self.totalForce

    oriMat = quaternion.ToMatrix(self.ori)
    inertia_world = self.inertia_local * oriMat    # HERE
    matrix.Invert(inertia_world)
    angAcc = vector.Transform(self.totalTorque, inertia_world)

    self.pos = self.pos + self.linVel*deltaTime
    self.ori = quaternion.AddVector(self.ori, self.angVel*deltaTime)
    self.ori = quaternion.Normalized(self.ori)        

    self.linVel = self.linVel + linAcc*deltaTime
    self.angVel = self.angVel + angAcc*deltaTime

    self.totalForce = vector.Vector3()
    self.totalTorque = vector.Vector3()
As you can see, in line "HERE" I calculate the inertia tensor in world space, by simply multiplying the local inertia matrix by the rigid body's orientation matrix (which is effectively the local-to-world matrix without the translation part).

Note that I use row-major matrix order.

The way I transform the local inertia tensor seemed natural to me this way. I have local space, I have local-to-world transform, so I just mul by that matrix to get inertia in world space. And it works just fine. I even made a simple simulation in Unity using their rigid body to make sure that my behaves more or less the same. And it does.

However, I googled a bit about the inertia tensor transform and found that people do it differently. For instance here:
http://gamedev.stackexchange.com/questions/70355/inertia-tensor-and-world-coordinate-conversion
or here:

http://www.pixar.com/companyinfo/research/pbm2001/pdf/notesg.pdf (code on page G22).

They both suggest "embedding" the local tensor matrix in between local-to-world and transpose(local-to-world). But when I do this in my simulation it simply works incorrectly.

Any idea where this difference comes from?

Share this post


Link to post
Share on other sites

Well, the embedding is correct. You can visualize this using parentheses when computing the world space angular momentum:

 

L = I * omega = R * I' * RT * omega = R * ( I' * ( RT * omega ) )

 

- First you basically transform the world space angular velocity to body space: RT * omega

- Then you apply the body space inertia (hence the prime to indicate body space): I' * ( RT * omega )

- Finally you transform everything back to word space to find get the world space angular momentum:  R * ( I' * ( RT * omega ) )

 

If you use row vectors you have to revert the whole equation. E.g.

 

L = omega * I = omega * RT * I' * R => I = RT * I' * R

 

HTH,

-Dirk

Edited by Dirk Gregorius

Share this post


Link to post
Share on other sites

I think I understand Dirk's explanation but it is described in terms of angular momentum whereas what I need is angular acceleration. I also looked again at those Siggraph notes and it appears they compute angular velocity using the inertia (1) and angular momentum. In my code however I compute angular velocity by integrating angular acceleration, which is computed from torque and intertia (2). Any "my" inertia (2) is different than inertia (1), that is, I compute world inertia using just one local-to-world matrix whereas world inertia (1) is computed using local-to-world and world-to-local.

I can't see where the Siggraph notes make use of the torque.

Share this post


Link to post
Share on other sites

Well, torque is just the change of angular momentum. The tensor transform stays the same. Assuming column vector you transform the local inertia like this:

 

I_world = R * I_local * RT

 

If you use row vectors it works like this:

 

I_world = RT * I_local * RT

 

My guess is that you might *not* have transformed all your equations to proper row vector form. I know row vectors are very popular in graphics, but for physics it is an easy way to shoot yourself into the foot. 

Share this post


Link to post
Share on other sites

Just my 2 cents: in the case you're using row major storage you can write your matrix multiplication in the reverse order. Like this:

inline CMatrix4x4 operator*(const CMatrix4x4& _mOther) const {
		CMatrix4x4 mOut;

		mOut._11 = _11 * _mOther._11 + _21 * _mOther._12 + _31 * _mOther._13 + _41 * _mOther._14;
		mOut._12 = _12 * _mOther._11 + _22 * _mOther._12 + _32 * _mOther._13 + _42 * _mOther._14;
		mOut._13 = _13 * _mOther._11 + _23 * _mOther._12 + _33 * _mOther._13 + _43 * _mOther._14;
		mOut._14 = _14 * _mOther._11 + _24 * _mOther._12 + _34 * _mOther._13 + _44 * _mOther._14;

		mOut._21 = _11 * _mOther._21 + _21 * _mOther._22 + _31 * _mOther._23 + _41 * _mOther._24;
		mOut._22 = _12 * _mOther._21 + _22 * _mOther._22 + _32 * _mOther._23 + _42 * _mOther._24;
		mOut._23 = _13 * _mOther._21 + _23 * _mOther._22 + _33 * _mOther._23 + _43 * _mOther._24;
		mOut._24 = _14 * _mOther._21 + _24 * _mOther._22 + _34 * _mOther._23 + _44 * _mOther._24;

		mOut._31 = _11 * _mOther._31 + _21 * _mOther._32 + _31 * _mOther._33 + _41 * _mOther._34;
		mOut._32 = _12 * _mOther._31 + _22 * _mOther._32 + _32 * _mOther._33 + _42 * _mOther._34;
		mOut._33 = _13 * _mOther._31 + _23 * _mOther._32 + _33 * _mOther._33 + _43 * _mOther._34;
		mOut._34 = _14 * _mOther._31 + _24 * _mOther._32 + _34 * _mOther._33 + _44 * _mOther._34;

		mOut._41 = _11 * _mOther._41 + _21 * _mOther._42 + _31 * _mOther._43 + _41 * _mOther._44;
		mOut._42 = _12 * _mOther._41 + _22 * _mOther._42 + _32 * _mOther._43 + _42 * _mOther._44;
		mOut._43 = _13 * _mOther._41 + _23 * _mOther._42 + _33 * _mOther._43 + _43 * _mOther._44;
		mOut._44 = _14 * _mOther._41 + _24 * _mOther._42 + _34 * _mOther._43 + _44 * _mOther._44;

		return mOut;
	}

and use 

 


I_world = R * I_local * RT

 

to get the correct inertia. Most of the books uses column major multiplication, so it's better to keep like that smile.png .

Edited by Irlan

Share this post


Link to post
Share on other sites

I think we're misunderstanding each other a bit here smile.png.

First of all, the code I posted works perfectly fine. And I transform the local inertia tensor to world inertia tensor (in column-major order) like this:

I_world = R * I_local

which is these parts of my code (note row-major order):

    oriMat = quaternion.ToMatrix(self.ori)
    inertia_world = self.inertia_local * oriMat

And I think this the "real" world inertia tensor. In general case, when you have a vector V in some local space, have local-to-world matrix transform, and want to get vector V in world space you do this:

V_world = local-to-world * V_local

not this:

V_world = local-to-world * V_local * world-to-local

I think your solution:

I_world = R * I_local * transpose(R)

is only "valid" in this specific formula:

L = R * I_local * transpose(R) * omega_world

but it's not because of this:

L = ( R * I_local * transpose(R) ) * omega_world

but because of this:

L = R * I_local  * ( transpose(R) * omega_world )

The "interpretation" for me here is that we take omega_world, transform it to local space (hence mul by transpose of R), then we can mul by local inertia (because omega is now in local space) and finally we take all this to world space by mul with R.

 

Sorry guys but I'm a proponent of row-major order and I'm going to use it ;).

Edited by maxest

Share this post


Link to post
Share on other sites

Yes, that is the correct interpretation. Column or row major order is purely a notational entity. In memory pretty much all matrices should be stored in row major format since graphics rendering hardware expects this. As far as notation is concerned, column major or row major form can be used due to operator overloading. So it really doesn't matter as long as the user is comfortable with the style.

 

Since I do a lot of programming with physics and graphics I find it's really nice to be comfortable with both notations.

Share this post


Link to post
Share on other sites

I am glad it works, but personally I would be suspicious and assume I might potentially miss a transform somewhere.

 

@Randy:

I never talked about row or column major order, but row and column vectors. These are two different things. The first is about whether the row or columns are continuous in memory. The other one is about whether you transform your vector like vT' = vT * MT or like v' = M * v. You can still have row major layout and use column vectors, but this would be less efficient. Actually a naive matrix and vector implementation would use exactly this.

 

Also the default on graphics HW is column major (first sentence):

https://msdn.microsoft.com/en-us/library/windows/desktop/bb509634(v=vs.85).aspx#Matrix_Ordering

Edited by Dirk Gregorius

Share this post


Link to post
Share on other sites

I am glad it works, but personally I would be suspicious and assume I might potentially miss a transform somewhere.

 

@Randy:

I never talked about row or column major order, but row and column vectors. These are two different things. The first is about whether the row or columns are continuous in memory. The other one is about whether you transform your vector like vT' = vT * MT or like v' = M * v. You can still have row major layout and use column vectors, but this would be less efficient. Actually a naive matrix and vector implementation would use exactly this.

 

Also the default on graphics HW is column major (first sentence):

https://msdn.microsoft.com/en-us/library/windows/desktop/bb509634(v=vs.85).aspx#Matrix_Ordering

 

I might be mistaken but that looks it's just stating HLSL notation, not the layout of the matrix in memory. So in memory you have a pre-transposed matrix, such that the notational columns lay within the rows in memory.

Share this post


Link to post
Share on other sites


I am glad it works, but personally I would be suspicious and assume I might potentially miss a transform somewhere.

I'm almost  sure this code works fine. The rigid body behaves quite predictably and very similar to simulation I made in Unity.

I could share my code if anyone wanted to play with it.

 

I'd also be happy if anyone could look again at the original code, which works, and its other variant that I think should work as well but for some reason doesn't.

Original code:

        oriMat = quaternion.ToMatrix(self.ori)
        inertia_world = self.inertia_local * oriMat
        matrix.Invert(inertia_world)
        angAcc = vector.Transform(self.totalTorque, inertia_world)

To recap: oriMat is effectively local-to-world transform. self.totalTorque is torque in world space. So this code transforms local inertia to world inertia and uses its inverse along with world torque to compute world angular acceleration.

 

Modified version:

        oriMat = quaternion.ToMatrix(self.ori)
        self.totalTorque = vector.Transform(self.totalTorque, matrix.Inverted(oriMat))
        angAcc = vector.Transform(self.totalTorque, matrix.Inverted(self.inertia_local))
        angAcc = vector.Transform(angAcc, oriMat)

Here I wanted to do the torque-inertia transform in local space. To do so, I transform self.totalTorque to rigid body's local space using inverse of oriMat. Then I mul that by inverted local inertia to get angular acceleration in local space. Finally, I transform angular acceleration to world space by transforming it with oriMat.

For some reason this variant doesn't yield correct behaviour of my rigid body and I really can't tell why. Anyone?

Share this post


Link to post
Share on other sites

I was pointed out by a friend that in the modified version I should skip the last line. So the coude should look this:

        oriMat = quaternion.ToMatrix(self.ori)
        self.totalTorque = vector.Transform(self.totalTorque, matrix.Inverted(oriMat))
        angAcc = vector.Transform(self.totalTorque, matrix.Inverted(self.inertia_local))

And indeed that works. I wrote down on paper all these transforms and this code gives the exact same series of transforms as the original one.

 

Original:

angAcc = torque * inertia^-1 = torque_world * (inertia_local * localToWorld)^-1 = torque_world * localToWorld^-1 * inertia_local^-1

Clearly this gives angAcc in world-space what I need.

 

The modified new:

angAcc = torque * inertia^-1 = (torque_world * localToWorld^-1) * inertia_local^-1 = torque_world * localToWorld^-1 * inertia_local^-1

This I would expect to give me angAcc in local-space (which further needs to be mul'led by localToWorld to be in world space) but math clearly shows this is the same as the previous formula.

 

Anyone sees what is wrong with my understanding? :)

Share this post


Link to post
Share on other sites

The mistake is that inertiaWorld != localToWorld  * inertiaLocal 

 

The reason is that inertia is a tensor and not a matrix and that they transform differently. Hence the physics definition of a tensor: "A tensor is what transforms like a tensor". smile.png

 

You have to build your example the other way around:

localAngularAcc = localInertia^-1 * localTorque

We can, of course, compute the body space torque like:

localTorque = localToWorld^T * worldTorque

Now we substitute equation (2) into equation (1) which yields:

localAngularAcc = localInertia^-1 * ( localToWorld^T * worldTorque )

Then we transform the local angular acceleration to world space

worldAngularAcc = localToWorld * localAngularAcc = localToWorld * localInertia^-1 * ( localToWorld^T * worldTorque )

Finally we define the world space inertia as:

worldSpaceInertia^-1 = localToWorld * localInertia^-1 * localToWorld^T

This yields the pleasing result:

worldAngularAcc = worldSpaceInertia^-1 * worldTorque 

The key point is really that inertia is a tensor and not a matrix. Note how all my operations are conceptually applied to a vector and not to the tensor. I find that confusing as well, but I am not aware of a better way to explain it! Hopefully the example makes clear how this works for the inertia tensor and how your math above is wrong missing a transformation.

 

Maybe you should think of it more as a change of basis. You transform the world torque into the local body space. Then you apply the local inertia tensor. And finally you transform the the result back to world space. 

 

HTH,

-Dirk

Edited by Dirk Gregorius

Share this post


Link to post
Share on other sites

Sorry for such a massive delay but I completely forgot about that thread...

 

Okay, I've gone through your equations and the link you provided but it's still not right.

You're saying (and the said link) that to transform tensor you go like transform * tensor * transform^T. But my code doesn't do it. It simply does transform * tensor and it works (in world-space). Look again at the code that computes world-space angular acc:

        localToWorld = quaternion.ToMatrix(self.ori)
        inertia_world = self.inertia_local * localToWorld
        angAcc_world = vector.Transform(self.totalTorque, matrix.Inverted(inertia_world))

Expanded and written in math notation it would be like:

angAcc_world =
self.totalTorque * inertia_world^-1 =
self.totalTorque * (self.inertia_local * localToWorld)^-1 =
self.totalTorque * localToWorld^-1 * self.inertia_local^-1

Again, localToWorld is only used once, not twice. And this code really works, he simulation behaves in predictable manner, as opposed to any other approach I try.

I don't really know where this:


Finally we define the world space inertia as:

worldSpaceInertia^-1 = localToWorld * localInertia^-1 * localToWorld^T

is taking place in my code.

 

One thing interesting here is that the last form of this equations could be written like this:

angAcc_world =
self.totalTorque * localToWorld^-1 * self.inertia_local^-1 =
(self.totalTorque * localToWorld^-1) * self.inertia_local^-1 =
totalTorque_local * self.inertia_local^-1 =
angAcc_local

Now that is darn confusing. Does it really mean that angAcc in world-space is the same as angAcc in local-space?

Edited by maxest

Share this post


Link to post
Share on other sites

Every math and physics book you will find will tell you that this is how it works, E.g. here:

https://en.wikipedia.org/wiki/Moment_of_inertia#The_inertia_matrix_in_different_reference_frames

 

Just because something looks correct, doesn't actually proof it is correct. Secondly, you essentially want to convince me here that you are right and the whole physics world is wrong (even though I have shown you the mistake several times backing it up with external resources). My suggestion to approach your problem would be to apply the correct transformation and now that you can rule that out as an error source find your real bug. Or maybe there is no bug and you have a wrong expectation what the outcome of your experiment is.

 

And obviously:

 

angAcc_world == angAcc_local  would only be correct iff local_to_world is the identity matrix. You have essentially proven yourself that your are wrong. 

Edited by Dirk Gregorius

Share this post


Link to post
Share on other sites

If you do Transform * Tensor the results will be *hugely* incorrect. Imagine if the body is almost symmetric (Tensor is almost a unit matrix times mass) - you will effectively multiply your acceleration by the transform matrix. So if you rotate the body 180 degrees around the X axis, all new torques will be completely invertex around that axis.

Share this post


Link to post
Share on other sites

That is a nice example! Let me quickly extend on it. Say we actually have a sphere inertia which is invariant under rotation.

 

I_world = R * I_local * RT

 

The sphere inertia I_local = 2/5 * mass * radius^2 * identity_matrix. Let's write this as I_local = s * E (where s = 2/5 * mass * radius^2 and E = identity matrix). Then

 

I_world = R * ( s * E ) * RT = s (R * E * RT) = s (R * RT) = s * E

 

Which is the correct result because of the symmetry. You formula breaks apart here.

Share this post


Link to post
Share on other sites

It's not that I'm stubborn and decided to proof that the whole world is wrong and I'm right. But I re-checked my simulation nubmer of times and know it behaves right. And I discovered why...

 

I thought that I would proof myself right by implementing the simulaton in Unity so that it's easier to compare with Unity's built-in rigidbody implementation (that one is for sure working right). I typed out the same simulation code I did for my Python app. And... in Unity it didn't work right. So I used the actual formula you're suggesting for transforming the inertia tensor. And it worked! Now that was a discovery... It pushed me to write it in yet another language, this time C++ (I thougt maybe I had screwed and/or misunderstood how objects in Python worked). And my simulation worked identically to the one in Python.

 

Since I'm stubborn to keep insisting on using row-major matrices order (and Unity uses column-major) I took a quick look at the function that adds vector to quaternion (to alter orientation with new angular velocity). Here it comes:

    quat operator + (const vec3& v) const
    {
        quat temp = quat(0.0f, v.x, v.y, v.z);
        temp = temp * *this;
        temp.w *= 0.5f;
        temp.x *= 0.5f;
        temp.y *= 0.5f;
        temp.z *= 0.5f;
        
        return quat(w + temp.w, x + temp.x, y + temp.y, z + temp.z);
    }

I took a blind shot and swapped the multiplication order in this line:

        temp = temp * *this;

to this:

        temp = *this * temp;

and surprisingly it worked!

To quickly recap:

When I use my original add-vector-to-quaternion function and inertia_local*local_to_world the simulation works.

When I use the new add-vector-to-quaternion (with quaternions mul swapped) and the "valid" inertia transformation the simulation works.

 

Why the wrong ordering of quaternion multiplication could compensate for the lack of the localToWorld^T matrix?

Edited by maxest

Share this post


Link to post
Share on other sites

First I am glad you found your problem. Congrats!

 

Quaternions and row vectors are another topic. For row vectors you have v' = v * R1 * R2 and for column vectors you have v' = R2 * R1 * v. The matrix closest to the vector is applied first. For quaternion multiplication we have the same as for column vectors. E.g. q = q2 * q1. The right most transform is applied first

 

What some people do. E.g. XNA/DirectXMath and also the Maya SDK is two swap the quaternion multiplication order to be consistent with row vector multiplication. I personally don't like that at all and it is one of the reasons I use column vectors. Anyway, it is something to be aware off and to keep in mind when using your math library I guess.

 

E.g. from the Maya SDK (http://download.autodesk.com/us/maya/2011help/API/class_m_quaternion.html#935885f085838e5c5cdc2b2a0fffad82):

Quaternions in Maya multiply on the right (post-multiply) the same as matrices. Many popular quaternion papers (Shoemake) use pre-multiplication where quaternions pre-multiply on the left so you must be aware of this when using quaternions.

 

The same hold true for XNA/DirectXMath (https://msdn.microsoft.com/en-us/library/windows/desktop/microsoft.directx_sdk.quaternion.xmquaternionmultiply(v=vs.85).aspx):

The result represents the rotation Q1 followed by the rotation Q2 to be consistent with XMMatrixMulplity concatenation since this function is typically used to concatenate quaternions that represent rotations (i.e. it returns Q2*Q1).

 

HTH,

-Dirk

Edited by Dirk Gregorius

Share this post


Link to post
Share on other sites

Yeah, I'm totally aware of how row- and column-major work :). That function I almost copied from a book (which uses column-major order) and simply missed that I needed to change the order.

 

I wonder now why the wrong ordering of quaternion multiplication could compensate for the lack of the localToWorld^T matrix in the inertia transform equation?

Share this post


Link to post
Share on other sites

The quaternion multiplication and using body and not the world space angular velocity in the quaternion integrator are two things that come to mind here which might be causing your problem.  The first is difficult to tell without seeing your quaternion multiplication. 

 

The quaternion derivative is defined as q' = 1/2 * w * q. Here I am assuming that we use a quaternion multiplication such that q = q2 * q1 (e,g. Shoemake) and w is the pure quaternion holding the world space angular velocity. If you want to use the body space angular velocity the formula conveniently changes to q' = 1/2 * q * w. Here w now holds the body space angular velocity. This is very easy to prove. Just plug w_global = q * w_local * conjugate( q ) into the formula. 

 

Also note that you need to normalize the quaternion after integration. And overwriting operator+ for the integration is a poor choice and totally misleading in my opinion.

Edited by Dirk Gregorius

Share this post


Link to post
Share on other sites

This topic is 850 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this