I use the XNA animation component to perform blending between animations. It uses spherical linear interpolation between matrices.

In most cases this works without problems. However, I have one model where all of the animations are blended wrong.
When I blend two animations where the character is standing, during the blend period he will face to the right and appear to be in the bind pose.

When one of the animations has the model laying on the ground, there appears to be warping of the entire model during the blend period.

I am not an expert on blending or animation and I didn't write the code, but I do have all the source code. Could anyone explain what might be going wrong with this particular model/set of animations?

Here is the code for the method SlerpMatrix that does the interpolation between matrices:

[source lang="csharp"]
private static Quaternion qStart, qEnd, qResult;
private static Vector3 curTrans, nextTrans, lerpedTrans;
private static Vector3 curScale, nextScale, lerpedScale;
private static Matrix startRotation, endRotation;
private static Matrix returnMatrix;

public static void SlerpMatrix(
ref Matrix start,
ref Matrix end,
float slerpAmount,
out Matrix result)
{
if (start == end)
{
result = start;
return;
}
// Get rotation components and interpolate (not completely accurate but I don't want
// to get into polar decomposition and this seems smooth enough)
Quaternion.CreateFromRotationMatrix(ref start, out qStart);
Quaternion.CreateFromRotationMatrix(ref end, out qEnd);
Quaternion.Lerp(ref qStart, ref qEnd, slerpAmount, out qResult);

// Get final translation components
curTrans.X = start.M41;
curTrans.Y = start.M42;
curTrans.Z = start.M43;
nextTrans.X = end.M41;
nextTrans.Y = end.M42;
nextTrans.Z = end.M43;
Vector3.Lerp(ref curTrans, ref nextTrans, slerpAmount, out lerpedTrans);

// Get final scale component
Matrix.CreateFromQuaternion(ref qStart, out startRotation);
Matrix.CreateFromQuaternion(ref qEnd, out endRotation);
curScale.X = start.M11 - startRotation.M11;
curScale.Y = start.M22 - startRotation.M22;
curScale.Z = start.M33 - startRotation.M33;
nextScale.X = end.M11 - endRotation.M11;
nextScale.Y = end.M22 - endRotation.M22;
nextScale.Z = end.M33 - endRotation.M33;
Vector3.Lerp(ref curScale, ref nextScale, slerpAmount, out lerpedScale);

// Create the rotation matrix from the slerped quaternions
Matrix.CreateFromQuaternion(ref qResult, out result);

// Set the translation
result.M41 = lerpedTrans.X;
result.M42 = lerpedTrans.Y;
result.M43 = lerpedTrans.Z;

// Add the lerped scale component
result.M11 += lerpedScale.X;
result.M22 += lerpedScale.Y;
result.M33 += lerpedScale.Z;
}[/source]

This kind of shortcut seems wrong to me. Just by assuming that the matrix is orthogonal (I guess that what CreateFromRotationMatrix does) doesn't make it so, so who knows what kind of quaternion you obtain from this operation.
Also, this kind of interpolation of scaling doesn't make sense to me either, spatial transformations aren't generally formed like that. It should work fine when there's no scaling because you'll be interpolating between zero matrices, but with scaling I'm not even sure what kind of output it'll produce.

The right way to do it would be to decompose A=UP (polar decomposition), obtain quaternion from U and lerp/slerp it, while interpolating P component-wise.
I'm not sure why the original author refrained from implementing the polar decomposition as its implementation is straightforward and easy.

Great, that was the sort of insight I was looking for.

If it is not too much to ask for, could you give (pseudo-) code for the changes I need to do?

Basically you'll need to this (in this matlab'ish notation):

function A = MatrixSlerp(A0,A1,t)
[U0,P0]=polarDec(A0);
[U1,P1]=polarDec(A1);

q0=dcm2quat(U0);
q1=dcm2quat(U1);
q = slerp(q0, q1, t);

P = (1-t)*P0 + t*P1;

A = quat2dcm(q)*P;
end


Now, polarDec is not a matlab function, but it means polar decomposition. slerp is also not a matlab function, but it's clear what it means.

You can implement polarDec like this:

function [U,P] = polarDec(A)
U=A;
for i=1:1:10
uInv=inv(U);
g=sqrt(sum(uInv(.^2)/sum(U(.^2));
U=(U*g+uInv'/g)/2;
end

P=U'*A;
end


If you have translation as well (and you do), just interpolate it separately.

Thank you for the help, but I'm at a stage where I'm not sure how to proceed. I guess I just don't have enough knowledge to translate the example you provided into working code.

I will try and search for the answer somewhere...