Slerp on 2 quaternions vs the multiplication of "s"lerp of per-axis quaternions

Question

I split a quaternion into 3 rotations, and interpolate those, and multiply to merge the result, I don't get the same as spherically interpolating the non-split quaternions.

var quat_A_full = quaternion.EulerXYZ(ax_rad, ay_rad, az_rad);//-0.3161251f, 0.7749172f, 0.0828281f, 0.5410246f
var quat_B_full = quaternion.EulerXYZ(bx_rad, by_rad, bz_rad);//-0.5004139f, 0.5069054f, 0.1151911f, 0.6923611f

var quat_A_x = quaternion.RotateX(ax_rad);//0.2424405f, 0f, 0f, 0.9701663f
var quat_A_y = quaternion.RotateY(ay_rad);//0f, 0.8527206f, 0f, 0.5223673f
var quat_A_z = quaternion.RotateZ(az_rad);//0f, 0f, 0.5134861f, 0.8580979f

var quat_B_x = quaternion.RotateX(bx_rad);//-0.7160885f, 0f, 0f, 0.6980096f
var quat_B_y = quaternion.RotateY(by_rad);//0f, 0.4602496f, 0f, 0.8877895f
var quat_B_z = quaternion.RotateZ(bz_rad);//0f, 0f, -0.3183014f, 0.9479896f

var quat_A_full_reconstructed = math.mul(quat_A_z, math.mul(quat_A_y, quat_A_x);
// quat_A_full == quat_A_full_reconstructed <-- this is TRUE

var quat_B_full_reconstructed = math.mul(quat_B_z, math.mul(quat_B_y, quat_B_x);
// quat_B_full == quat_B_full_reconstructed <-- this is also TRUE

// but if I start to interpolate, sometimes I get results drifting off quite lot. 

var result_full = math.slerp(quat_A_full, quat_B_full, animation);

var result_parts = math.mul(
    math.slerp(quat_A_z, quat_B_z, animation), 
    math.mul(math.slerp(quat_A_y, quat_B_y, animation), 
             math.slerp(quat_A_x, quat_B_x, animation))
);

result_full != result_parts // it drifts by a lot sometimes:
(-0.4088441f, 0.6605082f, 0.09956253f, 0.621822f)
[ x: -107.9114542, y: 47.7336353, z: 80.7913939 ] (euler degrees)
vs
(-0.2624312f, 0.6446596f, 0.2598185f, 0.6693567f)
[ x: -87.4064322, y: 46.6058942, z: 87.1785867 ] (euler degrees)

I would need to interpolate per-axis due to animation constraints. Perhaps it's not possible?

PS: This is (unity's) quaternion slerp function I use:

public static quaternion slerp(quaternion q1, quaternion q2, float t)
    {
        float dt = dot(q1, q2);
        if (dt < 0.0f)
        {
            dt = -dt;
            q2.value = -q2.value;
        }

        if (dt < 0.9995f)
        {
            float angle = acos(dt);
            float s = rsqrt(1.0f - dt * dt);    // 1.0f / sin(angle)
            float w1 = sin(angle * (1.0f - t)) * s;
            float w2 = sin(angle * t) * s;
            return quaternion(q1.value * w1 + q2.value * w2);
        }
        else
        {
            // if the angle is small, use linear interpolation
            return nlerp(q1, q2, t);
        }
    }

Answer 1

Sounds like STERP can be useful for you.

sterp quaternion interpolation