Applying weights to matrixes and vertices (bone rotation)

Question

I'm rotating the bones of a skeleton inside a mesh for a low poly 3D figure. On the vertex shader its applied like this.
glsl:

    vec4 vert1 = (bone_matrix[index1]*vertex_in)*weight;
    vec4 vert2 = (bone_matrix[index2]*vertex_in)*(1-weight);
    gl_Position =  vert1+vert2;

bone_matrix[index1] is the matrix of one bone and bone_matrix[index2] is the matrix of the other. weight designates vertex_in's membership to these bones. The problem is the closer the weight is to .5, the more the diameter of the elbow shrinks when a rotation is applied. I've tested it with around a 10,000 vertex cylinder shape (with a gradient of weights). The result looked like bending a garden hose.

I got my weighting method from these sources. Its actually the only way I could find:
http://www.opengl.org/wiki/Skeletal_Animation
http://ogldev.atspace.co.uk/www/tutorial38/tutorial38.html
http://blenderecia.orgfree.com/blender/skinning_proposal.pdf

The left is how the shape starts, the middle is how the above equation rotates it, and the right is my goal. The mid points are weighted 0.5. It only gets worse the more bent it is, at 180 degrees it has zero diameter.

• I've tried assembling the matrix on the shader, so that I can apply the weights to the rotation instead of the resulting vertices. It looks perfect like the one in picture on the right, but it requires assembling the matrix for every single vertex (expensive)
• I've looked into quaternions, but glsl doesn't natively support them (correct me if I'm wrong) and they're confusing. Is that what I need to do?
• I've considered having three bones per joint, and add a "kneecap" between every bone. This wouldn't eliminate the problem but would mitigate it.
• I'm considering projecting the vertex its original distance from the axis after they're rotated. This would fail at 180 degrees but would be (relatively) cheap.

So considering the options, or other options that I may not have considered, How have others avoid this pinching effect?

EDIT: I've gotten SLERP to work using quaternions but I opted not to use it as GLSL does not natively support it. I couldn't get the geometric SLERP to work as described by Tom. I got NLERP working for the first 90 degrees, so I added an extra "bone" between each joint. So to bend the forearm 40 degrees I bend the elbow and the forearm by 20 degrees each. This eliminates the pinching effect at the expense of doubling the quantity of bones which is not an ideal solution.

Answer 1

The problem

The cause of what your seeing is illustrated by the drawing in Levans answer. However, to understand what's going on consider what's happening when you execute the code:

If the first point vert1 has coordinates (p, 0) the coordinates of vert2 will be (p cos(α), p sin(α)) where α is the angle between the two bones (this is always possible given an appropriate coordinate transform). Adding these together using the appropriate weights w and 1-w we get the following coordinates:

x = w p + (1-w) p cos(α)
y = (1-w) p sin(α)

The length of this vector is:

length^2 = x^2 + y^2
         = (w p + (1-w) p cos(α))^2 + (1-w)^2 p^2 sin(α)^2
         = p^2 [w^2 + 2 w (1-w) cos(α) + (1-w)^2 cos(α)^2 + (1-w)^2 sin(α)^2]
         = p^2 [w^2 + (1-w)^2 + 2 w (1-w) cos(α)]

As an example, when w = 1/2 this simplifies to:

length^2 = p^2 (1/2 + 1/2 cos(α)) = p^2 cos(α/2)^2

And length = p |cos(α/2)| whereas the length of the original vectors is p (see graph). The length of the new vector becomes smaller, this is the shrinking effect that you perceived. The reason for this is that we are actually interpolating the two vertices along a straight line. If we want the keep the same length p we actually need to interpolate along a circle around the center of the rotation. One possible approach is to renormalize the resulting vector, preserving the width at the joint.

This means we need to divide the resulting vertex coordinates by |cos(α/2)| (or the more general result for arbitrary weights). This has as a side effect of course, a divide by zero whenever the angle is exactly 180° (for the same reason the width at the joint is zero with your technique).

I'm no skeletal animation expert, but it seems to me the original solution as you described it, is an approximation to work with small bone angles (where the shrinking effect is minimal).

Alternative approaches

A different approach is to interpolate your rotations instead of your vertices. See for example the slerp wiki page and this paper.

SLERP

The slerp technique is similar to the technique I described above in the sense that it also preserves the width at the joint, however it interpolates directly along a circular path around the joint. The general formula is:

gl_Position = [sin((1-w)α)*vert1 + sin(wα)*vert2]/sin(α)

Given the points from above vert1 = (p, 0) and vert2 = (p cos(α), p sin(α)) applying the SLERP formula yields result = (x, y) with:

x = p [sin((1-w)α) + sin(wα) cos(α)]/sin(α)
y = p sin(wα) sin(α)/sin(α) = p sin(wα)

Calculating the cosine cos θ of the angle between vert1 and result yields:

cos(θ) = vert1*result/(|vert1| |result|) = vert1*result/p^2
       = p^2 [sin(wα) + sin((1-w)α) cos(α)]/sin(α)/p^2
       = [sin(α) cos((1-w)α) - cos(α) sin((1-w)α) + sin((1-w)α) cos(α)]/sin(α)
       = cos((1-w)α)

The angle between vert2 and result is:

cos(φ) = vert2*result/p^2
       = [sin(wα) cos(α) + sin((1-w)α) cos(α)^2 + sin((1-w)α) sin(α)^2]/sin(α)
       = [sin(wα) cos(α) + sin((1-w)α) cos(α)]/sin(α)
       = [sin(wα) cos(α) + sin(α) cos(wα) - cos(α) sin(wα)]/sin(α)
       = cos(wα)

This means that θ/φ = (1-w)/w which expresses the fact that SLERP interpolates with constant radial velocity. When working with 3D rotation matrices we can express the rotation transforming vert1 into vert2 as M = inverse(A)*B = transpose(A)*B so that we can express the rotation angle α as:

cos(α) = (tr(M) - 1)/2 = (tr(transpose(A)*B) - 1)/2
       = (A[0][0]*B[0][0] + A[0][1]*B[1][0] + A[0][2]*B[2][0] + 
          A[1][0]*B[0][1] + A[1][1]*B[1][1] + A[1][2]*B[2][1] + 
          A[2][0]*B[0][2] + A[2][1]*B[1][2] + A[2][2]*B[2][2] - 1)/2

Quaternion LERP

When working with quaternions a good approximation to the SLERP is to linearly interpolate the quaternions directly after which you renormalize the result. This gives an interpolation curve identical to the one in SLERP, however interpolation does not occur at constant radial velocity.

If you really want to avoid these problems altogether you can always split your meshes at the joint and rotate these separately.