• 9
• 16
• 15
• 12
• 9

# Interpolating rotation matrices

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

## Recommended Posts

The recent "Why do game programmers use quaternions to rotate a vector?" made me remember a related problem that i had (and unelegantly solved with quaternions IMO) a while ago. It's quite simple, really. Let's say that you have a 3D vector representing a position P, and a 3D vector representing the velocity V at a given time t. Then, to interpolate this position, you simply do:
P' = P + t.V
My question is, what's the equivalent for orientations ? Each "current" orientation at a given time is represented by a rotation matrix M. From the physics engine (i use ODE) i get back an angular velocity represented as a 3D vector W. How do i calculate M', roughly equivalent (bad notation incoming) to this: ?
M' = M + t.W
At the moment, my solution involves converting everything to quaternions, using the slerp operator and converting back to a matrix, but i don't feel it's good... Y.

##### Share on other sites
I'm not quite sure the word "interpolation" is the correct term for this... It's used when you approximate the change of a known state into another -also known- state. What you're asking, is -basically- how you calculate the next orientation from the current one and the angular velocity.
Consider a vector v representing a vertex on a mesh, and the angular velocity for a given frame. Both the vertex's position and the angular velocity are expression of time -in general. To see how the vertex will move, one has to differentiate its position vector with respect to time.

It turns out that the the derivative of this vector must be a vector perpendicular to both the normal vector to the plane of rotation (that one is merely the angular velocity), as well as to the vector v itsself. (Because you want v to actually rotate) and its magnitude is given by the magnitude of the rate of change of the angle it rotates by, with respect to time.

This sums up to the rate of change of the vector v, being ω X v, where ω is the vector of angular velocity, and X denotes the cross product.

The above holds for all vectors in the model, thus it should also hold for its axes too. Since they are the ones that define the model's orientation (through its world/model matrix) one simply has to apply the above to each column vector (or row depending on the convention) of that matrix.

At that point, a specific form of anti-symmetrix matrix is quite often used. It's called the "anti-symmetrix skew" matrix ("cross-product" matrix is usual too), and has the property that given a cross product of vectors: y X v, it creates a matrix Y which yields the exact same result when it transforms v, thus: y X v == Y*v.
That matrix takes a vector y =(y1,y2,y3) and returns the matrix:
[ 0   -y3   y2 ][y3    0   -y1 ][-y2   y1    0 ]
.

It turns out that if you multiply this matrix by another matrix, it performs the effect described above, on each axis-vector in the second matrix. This way, you can differentiate the current orientation, and use the result to integrate it.
E.g. with simple euler integration:
O = O + O'*dt, where O is the orientation matrix, and O'= AntiSymmetricSkewMatrix(angularVelocity)*O.

This way you can keep track of an object's orientation, when you know its angular velocity.

##### Share on other sites
Using quaternions will help deal with drift, matrix-skew issues (you can use the "Star" matrix method as shown by someusername, but should periodically orthogonalize the orientation matrix).

From my vector library:

/*Two methods to update rotation with velocity:(currentVel is in world coordinates).1:Quat4 adr = .5f*mulQuat4BwIs0(currentRot,currentVel);              // Note w is 0.f.currentRot += adr*dt;   // Add derivative.currentRot.normalize();2:flt hdt = .5f*dt;Quat4 rdr(hdt*currentVel.x,hdt*currentVel.y,hdt*currentVel.z,1.f); // Note w is 1.f.currentRot = rdr ^ currentRot; // Rotate by delta transform.currentRot.normalize();*/

The first method is discussed in Baraff's papers. I derived the 2nd method while experimenting with basic quaternion/velocity concepts (which is also well known, though I have not seen the method in a publication). The methods are effectively equivalent (for small time steps, within float error, etc.).

Method 1 above and the "Star-matrix" method are both explained in Baraff's paper. Note that method 1 above is useful for higher order integration methods (RK2, RK4, etc.).

##### Share on other sites
Indeed, quaternions do much better with numerical drift. Especially in -self driven- physical simulations, it will be much longer till one starts noticing visible distortion in the models, when compared to matrices. (if they don't normalize every once in a while)

Also, the method I described is what is often seen in tutorials/papers on physical simulations with matrices; probably because it's the standard approach: "differentiate and then integrate".
However, one doesn't have to integrate the current orientation. In fact, this is where most of the error accumulates.

Just like the 2nd method suggested by John Schultz, the angular velocity can be used directly to create an axis-angle matrix that represents the instantaneous change in orientation. Once you have that matrix, you simply multiply it with the current orientation, to concatenate their effects, just like you do with most matrices.

This way you can limit the error introduced by integration to a minimum; just in the angle of rotation: dt*|ω|. However, an occasional normalization wouldn't hurt either...