Quaternions and hardware skinning

This is a question about a GLSL implementation, but the main problem is purely mathematical.

I’m implementing hardware skinning via a vertex shader in my engine. In order to support a reasonable amount of bones and not to be limited by the number of uniforms, I’m using 2 arrays: 1 for bone positions, the other for rotation quaternions. Shader source is attached at the end of the post.

Here’s the problem. I have a fish model that should look like screenshot 1 (Blender screenshot). If I convert the quaternions to matrices, I get correct results (see screenshot 2; well, the TBN matrix seems to be inverted - thus the lighting against the sun, but it’s a minor problem). However, this is not very efficient and I would prefer to use the quaternions directly to rotate the vectors. Unfortunately, my implementation is obviously wrong and produces incorrect results (see screenshot 3).

I admit quaternions are a little bit obscure to me. I mean, I understand the general concept more or less, but despite extensive web search, I haven’t been able to find anything more than a very general description of rotating a vector by a quaternion (v’ = q * v * q^-1), and I’m not even sure how should the multiplication of quaternion components look like. :o I’m basing my code on the Bullet physics engine’s math library.

So, if anyone could point out what I’m doing wrong, I’d be most grateful.

Also, do the TBN vectors need to be normalized after the rotation even if the untransformed ones are guaranteed to be of unit length? I’ve seen several skinning implementations (both hardware and software) and they seem to be inconsistent in this matter.

The screenshots:

The shader (xyz are the quaternion vector components, w is the scalar; DST_MAX_JOINTS is a definition prepended to all shader sources by the engine at shader compilation time):

``````// TBN matrix components
attribute vec3 dstTangent;
attribute vec3 dstBinormal;
//attribute vec3 dstNormal;	// normal is uploaded as gl_Normal

// skinning
attribute vec4 dstJointIndices;
attribute vec3 dstJointOffset1;
attribute vec3 dstJointOffset2;
attribute vec3 dstJointOffset3;
attribute vec3 dstJointOffset4;
//attribute vec4 dstJointWeights;	// weights are uploaded as gl_Vertex

// skeleton pose information
uniform vec4 dstSkeletonRot[DST_MAX_JOINTS];
uniform vec4 dstSkeletonPos[DST_MAX_JOINTS];

varying mat3 invTBN;

// comment this line to fall back to converting quaternions to matrices
#define ROTATE_USING_QUATS

#ifdef ROTATE_USING_QUATS
// v' = q * v * q^-1
vec4 multQuat(vec4 q1, vec4 q2) {
return vec4(
q1.w * q2.x + q1.x * q2.w + q1.y * q2.z - q1.z * q2.y,
q1.w * q2.y + q1.y * q2.w + q1.z * q2.x - q1.x * q2.z,
q1.w * q2.z + q1.z * q2.w + q1.x * q2.y - q1.y * q2.x,
q1.w * q2.w - q1.x * q2.x - q1.y * q2.y - q1.z * q2.z
);
}

vec3 rotateVectorByQuat(vec4 q, vec3 v) {
#if 0
vec4 qv = vec4(
q.w * v.x + q.y * v.z - q.z * v.y,
q.w * v.y + q.z * v.x - q.x * v.z,
q.w * v.z + q.x * v.y - q.y * v.x,
-q.x * v.x - q.y * v.y - q.z * v.z
);
/*return vec3(
qv.x * q.w + qv.y * -q.z - qv.z * -q.y,
qv.y * q.w + qv.z * -q.x - qv.x * -q.z,
qv.z * q.w + qv.x * -q.y - qv.y * -q.x
);*/
return vec3(
qv.w * -q.x + qv.x * q.w + qv.y * -q.z - qv.z * -q.y,
qv.w * -q.y + qv.y * q.w + qv.z * -q.x - qv.x * -q.z,
qv.w * -q.z + qv.z * q.w + qv.x * -q.y - qv.y * -q.x
);
#else
vec4 qv = multQuat(q, vec4(v, 0.0));
return multQuat(qv, vec4(-q.x, -q.y, -q.z, q.w)).xyz;
#endif
}
#else
mat3 quat2mat(vec4 q) {
float x2 = q.x * q.x;
float y2 = q.y * q.y;
float z2 = q.z * q.z;
float xx = q.x * x2;
float xy = q.x * y2;
float xz = q.x * z2;
float yy = q.y * y2;
float yz = q.y * z2;
float zz = q.z * z2;
float wx = q.w * x2;
float wy = q.w * y2;
float wz = q.w * z2;
return mat3(
1.0 - (yy + zz),	xy - wz,		xz + wy,
xy + wz,		1.0 - (xx + zz),	yz - wx,
xz - wy,		yz + wx,		1.0 - (xx + yy)
);
}
#endif

void main() {
gl_TexCoord[0] = gl_MultiTexCoord0;

// put this in an array for convenience
vec3 offset[4] = vec3[](
dstJointOffset1,
dstJointOffset2,
dstJointOffset3,
dstJointOffset4
);

// compute vertex position in object space
vec3 pos = vec3(0.0);
vec3 tangent = vec3(0.0);
vec3 binormal = vec3(0.0);
vec3 normal = vec3(0.0);
// gl_Vertex actually holds weights
for (int i = 0; i < 4; i++) {
int index = int(dstJointIndices[i]);
#ifdef ROTATE_USING_QUATS
pos += gl_Vertex[i] * (rotateVectorByQuat(dstSkeletonRot[index], offset[i]) + dstSkeletonPos[index].xyz);
tangent += gl_Vertex[i] * (rotateVectorByQuat(dstSkeletonRot[index], dstTangent));
binormal += gl_Vertex[i] * (rotateVectorByQuat(dstSkeletonRot[index], dstBinormal));
normal += gl_Vertex[i] * (rotateVectorByQuat(dstSkeletonRot[index], gl_Normal.xyz));
#else
mat3 jointMat = quat2mat(dstSkeletonRot[index]);
pos += gl_Vertex[i] * ((jointMat * offset[i]) + dstSkeletonPos[index].xyz);
tangent += gl_Vertex[i] * (jointMat * dstTangent);
binormal += gl_Vertex[i] * (jointMat * dstBinormal);
normal += gl_Vertex[i] * (jointMat * gl_Normal.xyz);
#endif
}
// transform to eye space
gl_Position = gl_ModelViewProjectionMatrix * vec4(pos, 1.0);

// find the inverese TBN matrix
// apply world transform
mat3 TBN = mat3(
gl_ModelViewMatrix[0].xyz,
gl_ModelViewMatrix[1].xyz,
gl_ModelViewMatrix[2].xyz
) * mat3(
tangent,
binormal,
normal
);
// invert it (orthogonal -> transpose)
invTBN = transpose(TBN);
}
``````

>> Also, do the TBN vectors need to be normalized after the rotation
No, rotation never changes the length of a vector. But the interpolation (in varyings) changes it.

I know that, I only wasn’t sure if the bone weighting doesn’t change it.

Can somebody take a look at the rotateVectorByQuat function in the shader?

OK, by googling around and reading some more articles on the subject I’ve found that my quaternion multiplication function is correct. The inversion (conjugation) seems to be correct as well. Any idea why this should produce erratic results?

I’ve wrote a Doom3 model loader and hence, I used Quaternion a little bit.:

( http://www.youtube.com/watch?v=utdoRSeKg64 if you are interested)

You may have an error in your multQuat, this is what I got in my Quaternion class:

``````
void operator*= (const Quaternion& otherQ)
{

w = (w * otherQ.w) - (x * otherQ.x) - (y * otherQ.y) - (z * otherQ.z);
x = (x * otherQ.w) + (w * otherQ.x) + (y * otherQ.z) - (z * otherQ.y);
y = (y * otherQ.w) + (w * otherQ.y) + (z * otherQ.x) - (x * otherQ.z);
z = (z * otherQ.w) + (w * otherQ.z) + (x * otherQ.y) - (y * otherQ.x);
}

``````

I compared and I am not sure we have a match, I can’t confirm as I am not sure what convention you used for your vec4 to match the x,y,z,w components.

In my implementation x, y, z are the imaginary part and w is the real component.

``````vec4 multQuat(vec4 q1, vec4 q2) {
return vec4(
q1.w * q2.x + q1.x * q2.w + q1.z * q2.y - q1.y * q2.z,
q1.w * q2.y + q1.y * q2.w + q1.x * q2.z - q1.z * q2.x,
q1.w * q2.z + q1.z * q2.w + q1.y * q2.x - q1.x * q2.y,
q1.w * q2.w - q1.x * q2.x - q1.y * q2.y - q1.z * q2.z
);
}
``````

It should work unless I misunderstood your convention (q1q2 != q2q1)

Nope, it doesn’t work, the fish is twisted in exactly the same way.

There’s a bug in your quaternion to matrix conversion. The x2, y2, and z2 variables should be initialized to the sum of corresponding quaternion coordinates, not their product. This bug could be masking the root cause of your problem since you say the matrix version “seems” to work.

Your quaternion inverse in rotateVectorByQuat only works on unit quaternions. Verify that this is true in your situation.

The TBN vectors don’t need to be normalized after a rotation if they’re normalized beforehand. This assumes that you are doing a pure rotation (using unit quaternions). They DO need to be renormalized after weighting.

Don’t be so sure that converting to matrix format for transformation is slower. The extra cost for conversion is offset by a faster transform operation.

Are you sure of this? To be honest, I just copied that code over from Bullet. Besides, I’ve been trying this code over and over again with different models and animations and they all seem to be all right.

It is, I designed the file format to imply that.

Which means I should normalize them because I’m weighting, right?

Ah. I guess I’ll have to benchmark the two approaches, then.