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?

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

You can check the derivation Dirk described here or at David Eberly Game Physics book.

The integrals for simple shapes (such blocks and spheres) are given in tables. Based on the volume, density, mass of the object you can get the correct local inertia.

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.

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.

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 .

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

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 ;).

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.

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 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.