Severe artifact when interpolating between dual quaternions

Question

I'm having trouble with my implementation of dual quaternion skinning. I'm still learning about the subject, so for the moment I'm converting from the bone matrix to a dual quaternion CPU side, and back to a matrix in the shader.

The conversion does apparently work correctly for single bones, but if I try to linearly blend between dual quaternions, I get this artifact: http://imagizer.imageshack.us/a/img838/8671/nun.gif I don't know what's causing this. Maybe it's related to how I normalize the dual quaternion, maybe it's in how I convert from dual quat to matrix. I've tried searching for actual dual quaternion code, but all I find is a bunch of hard to read mathematical definitions.

I am including pieces of the shader code, as I'm pretty sure that's where the problem is. Hopefully somebody proficient in quaternion math can look through it!

Blending the dual quaternions. Boneweight2 = (1.0 - boneweight1), so they'll always sum up to 1.

vec4 blendReal = boneReal[bone1] * boneWeight1 + boneReal[bone2] * boneWeight2;
vec4 blendDual = boneDual[bone1] * boneWeight1 + boneDual[bone2] * boneWeight2;
float blend_norm_real = length(blendReal);
blendReal /= blend_norm_real;
blendDual /= blend_norm_real;

Create matrix from dual quaternion:

mat4 MatFromDualQuat(vec4 rq, vec4 dq) 
{
    //Source: Section 3.4 http://www.seas.upenn.edu/~ladislav/papers/sdq-i3d07/sdq-i3d07.pdf
    //rq = real quaternion
    //dq = dual quaternion
    mat4 M;
    M[0][0] = 1.0 - 2.0 * (rq.y * rq.y + rq.z * rq.z); //
    M[1][0] = 2.0 * (rq.x * rq.y + rq.w * rq.z);//
    M[2][0] = 2.0 * (rq.w * rq.y - rq.x * rq.z);//
    M[3][0] = 0.0;
    M[0][1] = 2.0 * (rq.x * rq.y - rq.w * rq.z);
    M[1][1] = 1.0 - 2.0 * (rq.x * rq.x + rq.z * rq.z);
    M[2][1] = 2.0 * (rq.y * rq.z + rq.w * rq.x);
    M[3][1] = 0.0;
    M[0][2] = - 2.0 * (rq.x * rq.z + rq.w * rq.y);
    M[1][2] = 2.0 * (rq.y * rq.z - rq.w * rq.x);
    M[2][2] = 1.0 - 2.0 * (rq.x * rq.x + rq.y * rq.y);
    M[3][2] = 0.0;
    M[0][3] = 2.0 * (-dq.w * rq.x + dq.x * rq.w + dq.z * rq.y - dq.y * rq.z);
    M[1][3] = 2.0 * (-dq.w * rq.y + dq.y * rq.w + dq.x * rq.z - dq.z * rq.x);
    M[2][3] = 2.0 * (-dq.w * rq.z + dq.z * rq.w + dq.y * rq.x - dq.x * rq.y);
    M[3][3] = 1.0;
    return M;
}

And then I multiply that with the bind pose vertex position.