Oh, I had that too, but it caused everything to spin like crazy. Matrix and quat libraries are a bit weird when it comes to column-major stuff, you need to triple-check the order.

Wish you had mentioned that a few posts before ;)

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.

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

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?

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.

Oh by the way, the additive quaternion update isn't quite accurate, isn't it? Is it better to separate the angle from the angular velocity and multiply by a correct quat instead? I.e. instead of doing q += t/2 * (0,w) * q, do q = (cos phi/2, sin phi/2 * w) * q ?

Yes, I found this better for "small" angular velocities. What you describe is sometimes referred to as "exponential map" or "Grassia's" method. You need to take some extra care if Theta -> 0 if I remember correctly.

Here is the reference if you are interested:

http://www.cs.cmu.edu/~spiff/moedit99/expmap.pdf

Section 3.0 is the part of practical relevance. A good read in general though.

I've just re-checked my statement from this post (http://www.gamedev.net/topic/667652-rigid-body-physics-inertia-tensor-space/#entry5246272 at the bottom) by rendering the object twice, using both approaches. And it appears that the outcome is a little bit different. It is very close but under some rare circumstances the two can behave differently but eventually they diverge to a very similar result. I also think I noticed that the bigger the integration timestep the bigger the difference might be.

I'll just assume that the "most correct" solution is the one with properly transforming the tensor matrix, although I find it very interesting that the other gives yields similar results.

Oh by the way, the additive quaternion update isn't quite accurate, isn't it? Is it better to separate the angle from the angular velocity and multiply by a correct quat instead? I.e. instead of doing q += t/2 * (0,w) * q, do q = (cos phi/2, sin phi/2 * w) * q ?

What's phi here? Angle of angular velocity? How do you extract angle from a 3D vector?

And overwriting operator+ for the integration is a poor choice and totally misleading in my opinion.

Why do you think so? It's addition of a vector (angular velocity) to a quaternion (orientation) so it seemed natural to overload + operator.

You asked about my quaternion mul formula. Here it is:

def Mul(q1, q2):

    temp = Quaternion()

    temp.x = q1.w*q2.x + q1.x*q2.w + q1.z*q2.y - q1.y*q2.z
    temp.y = q1.w*q2.y + q1.y*q2.w + q1.x*q2.z - q1.z*q2.x
    temp.z = q1.w*q2.z + q1.z*q2.w + q1.y*q2.x - q1.x*q2.y
    temp.w = q1.w*q2.w - q1.x*q2.x - q1.y*q2.y - q1.z*q2.z

    return temp

I can't imagine why not transforming the torque proprely, and then multiplying the quats in the wrong order would produce reasonable results, maybe its something specific to your case.

Phi is the angle, yes. Angular velocity is the direction vector times angular speed, so all you need to do is compute its length. However, as Dirk mentioned, since the formula essentially becomes v*sin(len/2)/len, you need to be extra careful when len is small. The function sin(len)/len does behave well near zero, but you need a special way to compute it, most likely using taylor series (the function is called sinc and its probably included in boost math libraries so you could take the code from there). But then again, when the angle is truly small, there is very little difference between the 'precise' approach and the one you've been using, as sin(x)/x is very close to 1. Just make sure to re-normalize your quat every time you update it, or it will drift no matter which method you use.

I think I have different signs in my quat update function:

  quat.adjustVec3 = function(out, q, v, t) {
    var qx = q[0], qy = q[1], qz = q[2], qw = q[3], vx = v[0] * t * 0.5, vy = v[1] * t * 0.5, vz = v[2] * t * 0.5;
    out[0] = qx + qw * vx + qz * vy - qy * vz;
    out[1] = qy + qx * vz + qw * vy - qz * vx;
    out[2] = qz + qy * vx - qx * vy + qw * vz;
    out[3] = qw - qx * vx - qy * vy - qz * vz;
    quat.normalize(out, out);
    return out;
  };

If you multiply omega by dt the length of the resulting vector is Phi. It is kind of a simplification since you assume a constant direction of omega over the integration interval.

f(x) = sin(x) / x in the limit x -> 0 is of the form 0/0. Interestingly the actual limit is f( x -> 0 ) = 1 which can be shown with the l'hopital rule.

https://en.wikipedia.org/wiki/L%27H%C3%B4pital%27s_rule

I don't think that overloading operator+ is a good idea because if someone else reads your code it might not be immediately apparent to him what your are trying to do. Writing a simple function IntegrateOrientation() would be much more meaningful in this case (e.g. look for self-documenting code). I also like the term of 'navigation depth' in that context. It means if I read a function in some code how deep do I need to navigate into sub functions to understand the code. Ideally I can understand the overall function by just looking at the top level code. It makes you a lot of friends writing readable code.

Regarding the original problem I think we have been through all details and you should be able to figure out the problem yourself. Just take your Unity example and then step through both programs side by side until you find the problem. This shouldn't really be too hard!