Interpolation over a sequence of quaternions

Question

I am trying to interpolate over a sequence of rotations represented by quaternions. I am currently employing Squad (Spherical and Quadrangle Interpolation). I successfully applied the function to 4 rotations and the results corresponded to the expected one.

I would like to generalize the technique to an arbitrary (t>4) number of rotations. In Quaternions, Interpolation and Animation, on page 54, it is hinted that Squad can be applied to any sequence of rotations.

Squad is defined as: $$Squad(q_i, q_{i+1}, s_i, s_{i+1}, t) := Slerp(Slerp(q_i, q_{i+1}, t), Slerp(s_i, s_{i+1}, t), 2t(1-t))$$ Support rotations are defined as:$$s_i := q_i\cdot exp{\Big(-{\frac{\log(q_i^{-1}\cdot\log(q_{i-1})) + \log(q_i^{-1}\cdot q_{i+1})}{4}}\Big)}$$

The report continues stating that supports should be enhanced as follows: $$s_0 = q_0; s_n = q_n$$

If I understood correctly, computing support rotations and invoking Squad with the corresponding four quaternions is sufficient.

For each quaternion i in the range [1..n-2] I compute the corresponding s_i. At interpolation time I map the t parameter to the sequence of rotation and identify an i index, I remodulate the t parameter and eventually invoke $$Squad(q_i, q_{i+1}, s_i, s_{i+1}, t)$$ No need to say I implemented this technique but the results are wrong.

I double-checked my implementation of Slerp, I tested log(), exp(), and other quaternion methods, but they all seem to be correct. I concluded I didn't understand how to employ the technique, so I am asking for advice here.

Do you see any flaw in the reasoning I have just shown?

Edit: I am sharing part of the code, perhaps the problem is strictly practical.

Slerp:

static Quaternion Slerp(Quaternion &lhs, Quaternion &rhs, float t, TM_BOOL invert = true) {
    t = Unit(t);    // Clamps in range 0..1.
    if (t < 1e-20)
        return lhs;
    if (t > 1.0 - 1e-20)
        return rhs;
    if (lhs == rhs)
        return lhs;
    if (lhs == -rhs)
        return lhs;
    float dot = lhs.dot(rhs);
    if (abs(dot) >= 0.9999999) {
        if (dot >= 0.0)
            return Lerp(lhs, rhs, t).normalize();
        else
            return Lerp(lhs, -rhs, t).normalize();
    } else {
        float tetha = acos(dot);
        if (!invert)
            return (lhs * sin(tetha * (1.0 - t)) + rhs * sin(tetha * t)) / sin(tetha);
        else {
            if (dot >= 0.0)
                return (lhs * sin(tetha * (1.0 - t)) + rhs * sin(tetha * t)) / sin(tetha);
            else
                return (lhs * sin(tetha * (1.0 - t)) + (-rhs) * sin(tetha * t)) / sin(tetha);
        }
    }
}

Squad:

static Quaternion Squad(
    Quaternion &start,
    Quaternion &end,
    Quaternion &firstSupport,
    Quaternion &secondSupport,
    float t
) {
    t = Unit(t);
    return Slerp(
        Slerp(start, end, t, false),
        Slerp(firstSupport, secondSupport, t, false),
        2.0 * t * (1.0 - t),
        false
    );
}

Squad with support:

static Quaternion Squad(
    Quaternion &q0,
    Quaternion &q1,
    Quaternion &q2,
    Quaternion &q3,
    float t
) {
    q1 = (q0 + q1).squaredMagnitude() < (q0 - q1).squaredMagnitude() ? -q1 : q1;
    q2 = (q1 + q2).squaredMagnitude() < (q1 - q2).squaredMagnitude() ? -q2 : q2;
    q3 = (q2 + q3).squaredMagnitude() < (q2 - q3).squaredMagnitude() ? -q3 : q3;

    Quaternion
        firstSupport = Support(q0, q1, q2),
        secondSupport = Support(q1, q2, q3);

    return Squad(q0, q3, firstSupport, secondSupport, t);
}

Support:

inline static Quaternion Support(
    Quaternion &previous,
    Quaternion &current,
    Quaternion &next
) {
    Quaternion
        inverseOfQuat = current.inverse(),
        leftProduct = inverseOfQuat * previous,
        rightProduct = inverseOfQuat * next,
        leftLog = leftProduct.log(),
        rightLog = rightProduct.log(),
        sum = leftLog + rightLog,
        param = -sum / 4.0,
        result = current * param.exp();
    return result;
}

Rotation Sequence C'tor:

RotationSequence(
    std::vector<Quaternion> rotations
)  {
    if (rotations.size() < 4)
        throw std::runtime_error("Rotation sequence requires at least four rotations.");
    for (Quaternion rotation : rotations)
        this->rotations.push_back(rotation);
    /* Copy the last position to force the interpolation to behave correctly. */
    this->rotations.push_back(rotations[rotations.size() - 1]);
    /* Preprocesses the sequence, selectively negating all quaternions. */
    for (size_t ui = 1; ui < this->rotations.size(); ui++)
        if (this->rotations[ui].dot(this->rotations[ui - 1]) < 0.0)
            this->rotations[ui].negate();
    /* Compute supports for the whole sequence, 0-th and n-th correspond to quaternions 0-th and n-th. */
    supports.push_back(this->rotations[0]);
    for (size_t ui = 1; ui <= this->rotations.size()-2; ui++)
        supports.push_back(
            Support(
                this->rotations[ui - 1],
                this->rotations[ui + 0],
                this->rotations[ui + 1]
            )
        );
    supports.push_back(rotations[rotations.size() - 1]);

Interpolation:

Quaternion Interpolate(
    float t
) {
    t = Unit(t);
    /* Here, I identify the corresponding rotation by the $t$ factor. */
    float progress = float(rotations.size() - 2) * t;
    size_t ui = size_t(floor(progress));
    t = progress - float(ui);
    assert(ui + 0 < rotations.size());
    assert(ui + 1 < rotations.size());
    return Squad(
        rotations[ui],
        rotations[ui + 1],
        supports[ui],
        supports[ui + 1],
        t
    );
}

Answer 1

The mod expression looked a bit funny to me — do you observe any differences if you use this instead?

Quaternion Interpolate(
    float t
) {
    t = Unit(t);

    float progress = float(rotations.size() - 1) * t;
    size_t ui = size_t(progress);
    assert(ui + 1 < rotations.size() - 1);

    t = progress - float(ui);

    return Squad(
        rotations[ui],
        rotations[ui + 1],
        supports[ui],
        supports[ui + 1],
        t
    );
}

I think you also need some handling for the case when the input t == 1, since that should be a valid place to query the interpolation (it should just return the final quaternion), but with this code it will trip your assert.