Interpolate between two quaternions the long way

Question

We can do a slerp interpolation between two quaternions like this:

quat slerp(quat q1, quat q2, float t) {
    float angle = acos(dotProduct(q1, q2));
    float denom = sin(angle);
    //check if denom is zero
    return (q1*sin((1-t)*angle)+q2*sin(t*angle))/denom;
}

This will interpolate between the two quaternions the shortest way. However there's also a long way to interpolate between quaternions. As seen in the image below (source Maya).

How do we interpolate the long way?

Answer 1

If acos(dot(a/|a|, b/|b|)) = 1 is the angle between vectors a and b, then the rest of the circle is 2pi - acos(dot(a/|a|, b/|b|)), but the measure is 1/acos(dot(a/|a|, b/|b|)) from the first eq; so we get (2pi - acos(dot(a/|a|, b/|b|)))/acos(dot(a/|a|, b/|b|)), which simplifies to 2pi/acos(dot(a/|a|, b/|b|)) - 1.

So slerp with t = [0; -2pi/acos(dot(a/|a|, b/|b|)) + 1] would go the long way.

Answer 2

Idk if this is unnecessarily complicating things but:

Take the midpoint of the short path. The opposite of that will be the midpoint of the long path.

So we short slerp from q1 to the opposite in the first half of the time, then short slerp from the opposite to q2 in the second half of the time.

quat longSlerp(quat q1, quat q2, float t) {
    quat midpoint = slerp(q1, q2, 0.5);
    quat opposite = -midpoint;
    if(t <= 0.5)
        return slerp(q1, opposite, 2t);
    else
        return slerp(opposite, q2, 2(t-0.5));
}

Answer 3

So in an effort to save some confusion from others who may go through this I knocked together some SDL/OpenGL code in an effort to get either Mauricio's or Martin's answer to work. I found Martin's answer as it stands a little bit nebulous when implementing even though it states the truth. I haven't managed to get Mauricio's answer to work even with his help unfortunately.

I made a few mistakes too, and the number of different slerp functions I tried from various places to sanity check things ended up causing me some confusion so I ended up implementing my own slerp (SlerpIam() in the code) from scratch without checking for the nearest path.

In the code Slerp1() & Slerp2() I think are broken for when the shortest path is not selected which was part of my confusion - as from the myriad of slerps I found I think they were erroneously modified to try and support longest paths but they don't. So I was originally trying to hack them as Martin mentioned but it went horribly wrong.

My test case shows a point being rotated/rolled around the Z axis to 270 degrees.

I compiled the code with SDL2 on Windows and you will need to include SDL headers and link etc:

#include <cmath>


constexpr float PI = 3.14159265358979323846;


struct Quat         { float x, y, z, w; };
struct Vec3         { float x, y, z; };
struct AxisAngle    { Vec3 axis; float angle; };


float ToRadian(float degree)    { return degree * PI / 180.0f; }

Quat operator+(Quat a, Quat b)  { return { a.x + b.x, a.y + b.y, a.z + b.z, a.w + b.w }; }
Quat operator*(Quat q, float s) { return { q.x * s, q.y * s, q.z * s, q.w * s }; }
Quat operator*(float s, Quat q) { return { q.x * s, q.y * s, q.z * s, q.w * s }; }
Quat operator*(Quat second, Quat first)
{
    return Quat
    {
        second.w*first.x + second.x*first.w + second.y*first.z - second.z*first.y,
        second.w*first.y - second.x*first.z + second.y*first.w + second.z*first.x,
        second.w*first.z + second.x*first.y - second.y*first.x + second.z*first.w,
        second.w*first.w - second.x*first.x - second.y*first.y - second.z*first.z
    };
}

float Dot(Quat a, Quat b)   { return a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w; }
float Length(Quat q)        { return sqrtf(Dot(q, q)); }
Quat Normalise(Quat q)      { return q * (1.0f / sqrtf(Dot(q, q))); }
Quat Conjugate(Quat q)      { return{ -q.x, -q.y, -q.z, q.w }; }
Quat Reciprocal(Quat q)     { return Conjugate(q) * (1.0f / Dot(q, q)); }
Vec3 Rotate(Quat q, Vec3 v) { Quat vr = q * Quat{ v.x, v.y, v.z, 0.0f } *Conjugate(q); return { vr.x, vr.y, vr.z }; }

Quat ToQuat(AxisAngle r)
{
    float halfAngle = 0.5f * r.angle;
    float sinHalfAngle = sinf(halfAngle);


    return{ r.axis.x * sinHalfAngle, r.axis.y * sinHalfAngle, r.axis.z * sinHalfAngle, cosf(halfAngle) };
}


AxisAngle ToAxisAngle(Quat q)
{
    float s = 1.0f / sqrtf(1.0f - q.w * q.w);


    return { { q.x * s, q.y * s, q.z * s }, acosf(q.w) * 2.0f };
}


Quat Exp(Quat q)
{
    double b = sqrt(q.x * q.x + q.y * q.y + q.z * q.z);


    if (fabs(b) <= 1.0e-14 * fabs(q.w))
        return { 0.0f, 0.0f, 0.0f, expf(q.w) };
    else
    {
        float e = expf(q.w);
        float f = sinf(b) / b;


        return { e * f * q.x, e * f * q.y,  e * f * q.z, e * cosf(b) };
    }
}



Quat SlerpIam(Quat a, Quat b, float t)
{
    float dotAB = Dot(a, b);
    float theta = acosf(dotAB);
    float sinTheta = sinf(theta);
    float af = sinf((1.0f - t) * theta) / sinTheta;
    float bf = sinf(t * theta) / sinTheta;


    return a * af + b * bf;
}


Quat Slerp1(Quat q0, Quat q1, float t, bool shortPath = true)
{
    float d = Dot(q0, q1);
    float s0, s1;
    float sd = shortPath ? (d > 0) - (d < 0) : 1.0f;


    d = fabs(d);

    if (d < 0.9995f)
    {
        float s = sqrtf(1 - d * d);    //   Sine of relative angle
        float a = atan2f(s, d);
        float c = cosf(t*a);


        s1 = sqrtf(1 - c * c) / s;
        s0 = c - d * s1;
    }
    else
    {
        s0 = 1.0f - t;
        s1 = t;
    }


    return q0 * s0 + q1 * sd * s1;
}


Quat Slerp2(Quat q0, Quat q1, float t, bool shortPath = true)
{
    float a = 1.0f - t;
    float b = t;
    float d = Dot(q0, q1);
    float c = fabsf(d);


    if (c < 0.9995f)
    {
        c = acosf(c);
        b = 1.0f / sinf(c);
        a = sinf(a * c) * b;
        b *= sinf(t * c);

        if (shortPath && d < 0)
            b = -b;
    }


    return q0 * a + q1 * b;
}


Quat FarSlerpMauricio(Quat q0, Quat q1, float t)
{
    Quat q01 = Reciprocal(q0) * q1;
    Quat Vdash{ q01.x, q01.y, q01.z, 0.0f };
    Quat V = Vdash * (1.0f / Length(Vdash));
    float theta = 2.0f * atan2f(Length(Vdash), q01.w);
    float CompTheta = theta - 2.0f * M_PI;


    return q1 * Exp(t * CompTheta * V * 0.5f);
}


void Draw()
{
    float t = float(SDL_GetTicks() % 6000) / 6000.0f;
    Quat id{ 0.0f, 0.0f, 0.0f, 1.0f};
    Quat target = ToQuat({{0.0f, 0.0f, 1.0f}, ToRadian(270.0f)});

    //Quat r = FarSlerpMauricio(id, target, t);
    Quat r = SlerpIam(id, target, t);
    //Quat r = Slerp1(id, target, t);
    //Quat r = Slerp2(id, target, t);

    Vec3 p = Rotate(r, { 1.0f, 0.0f, 0.0f });
    

    glClearColor(0.2f, 0.2f, 0.2f, 1.0f);
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    glBegin(GL_LINES);

    //  Floor grid
    glColor3f(0.13f, 0.13f, 0.13f);
    for (int i = 0; i < 8; ++i)
    {
        float f = 2.0f * float(i) / 7.0f - 1.0f;
        glVertex3f(-1.0f, 0.0f, f);
        glVertex3f(+1.0f, 0.0f, f);
        glVertex3f(f, 0.0f, -1.0f);
        glVertex3f(f, 0.0f, +1.0f);
    }

    //  Axii
    glColor3f(0.8f, 0.0f, 0.0f);
    glVertex3f(0.0f, 0.0f, 0.0f);
    glVertex3f(1.0f, 0.0f, 0.0f);

    glColor3f(0.0f, 0.8f, 0.0f);
    glVertex3f(0.0f, 0.0f, 0.0f);
    glVertex3f(0.0f, 1.0f, 0.0f);

    glColor3f(0.0f, 0.0f, 0.8f);
    glVertex3f(0.0f, 0.0f, 0.0f);
    glVertex3f(0.0f, 0.0f, 1.0f);

    //  Ray to path
    glColor3f(1.0f, 1.0f, 1.0f);
    glVertex3f(0.0f, 0.0f, 0.0f);
    glVertex3fv(&p.x);

    glEnd();
}


int main()
{
    SDL_GLContext openGL;
    SDL_Window* window;
    bool run = true;


    if (SDL_Init(SDL_INIT_VIDEO) < 0)
        return -1;
    
    SDL_GL_SetAttribute(SDL_GL_CONTEXT_MAJOR_VERSION, 2);
    SDL_GL_SetAttribute(SDL_GL_CONTEXT_MINOR_VERSION, 1);
    SDL_GL_SetAttribute(SDL_GL_DOUBLEBUFFER, 1);
    SDL_GL_SetAttribute(SDL_GL_DEPTH_SIZE, 24);
    SDL_GL_SetAttribute(SDL_GL_MULTISAMPLEBUFFERS, 1);
    SDL_GL_SetAttribute(SDL_GL_MULTISAMPLESAMPLES, 8);

    if (!(window = SDL_CreateWindow("slerp", 100, 100, 800, 800, SDL_WINDOW_OPENGL | SDL_WINDOW_SHOWN)))
        return -1;

    openGL = SDL_GL_CreateContext(window);

    glViewport(0, 0, 800, 800);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    glOrtho(-2.0f, 2.0f, -2.0f, 2.0f, -2.0f, 2.0f);
    glRotatef(45.0f, 1.0f, 0.0f, 0.0f);
    glRotatef(45.0f, 0.0f, 1.0f, 0.0f);
    
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();

    glEnable(GL_DEPTH_TEST);
    glDepthFunc(GL_LEQUAL);
    glClearDepth(1.0f);
        
    glDisable(GL_CULL_FACE);
    glCullFace(GL_BACK);
    glFrontFace(GL_CCW);
    
    while (run)
    {
        SDL_Event event;


        while (SDL_PollEvent(&event) != 0)
        {
            if (event.type == SDL_MOUSEBUTTONDOWN || event.type == SDL_MOUSEMOTION)
                ;

            if (event.type == SDL_QUIT)
                run = false;
        }

        Draw();
        SDL_GL_SwapWindow(window);
    }

    SDL_GL_DeleteContext(openGL);
    SDL_DestroyWindow(window);


    return 0;
}

So working from my own SlerpIam() from scratch, I think my sanity is restored and Martin's answer is in essence correct. I get the following functions that I think are correct (Note they don't deal with small angle lerp fallback presently):

Quat SlerpNear(Quat a, Quat b, float t)
{
    float dotAB = Dot(a, b);

    if (dotAB < 0.0f)
    {
        dotAB = -dotAB;
        b = -b;
    }

    float theta = acosf(dotAB);
    float sinTheta = sinf(theta);
    float af = sinf((1.0f - t) * theta) / sinTheta;
    float bf = sinf(t * theta) / sinTheta;


    return a * af + b * bf;
}


Quat SlerpFar(Quat a, Quat b, float t)
{
    float dotAB = Dot(a, b);

    if (dotAB > 0.0f)
    {
        dotAB = -dotAB;
        b = -b;
    }

    float theta = acosf(dotAB);
    float sinTheta = sinf(theta);
    float af = sinf((1.0f - t) * theta) / sinTheta;
    float bf = sinf(t * theta) / sinTheta;


    return a * af + b * bf;
}

Answer 4

The nature of unit quaternions and the way they map to 3D rotations means they can describe each 3D rotation value in two ways - as q(r, v') and as q(-r, -v') (imagine them as axis-angle rotations - inverting both the axis and the angle leads to the same 3D rotation).

Quaternions are actually points on a 4D unit spherical surface, and these two values represent anti-podal points on that sphere.

For a slerp (or nlerp) of two quaternions to follow the shortest path, the corresponding 4D points have to lie on the same hemisphere of the 4D sphere (this is also the reason why a weighted average of more than 2 quaternions doesn't have a unique solution). This maps to a non-negative dot product, and is usually something tested for in the interpolation code.

Simply negating one of the source quaternions will give you a point "on the opposite side of the 4D sphere", and lead to interpolation "the long way around" (and explains why negating the interpolation parameter leads to the same result).

Answer 5

This can be done by changing the angle in the spherical interpolation.

In the usual SLERP(q1, q2, t) you get q1 when t=0 and q2 when t=1. The geodesic distance traveled is actually the angle between q1 and q2, which we name theta.

What we want to do here is to travel the complement distance which is 2PI - theta, but in the opposite sense of rotation. We will call this the complement theta.

We want to find a quaternion-valued function Q(t) such that:

SLERP2(q1, q2, t) = q1 Q(t)

when t = 0

SLERP2(q1, q2, 0) = q1 Q(0) = q1

and when t=1

SLERP2(q1, q2, 1) = q1 Q(1) = q2.

So we know that Q(0) = 1 (identity quaternion) and Q(1) = (q1^-1 q2).

It turns out that we can define such a function Q(t) in terms of the exponential map and principal logarithm of quaternions:

Q(t) = Exp(t Log(q1^-1 q2)/2)

You can check that it works by giving values to t like t=0 and t=1.

So far ao good, but that Q(t) will lead us to have the regular SLERP, not the one that we want. Lets take a close look at the logarithm:

Log(q1^-1 q2) = theta V

Where V is a unit vector (actually pure unit quaternion) which is the axis of rotation of quaternion q1^-1 q2. And theta is basically the angle between q1 and q2.

We actually need to change that logarithm so that Q(t) will go the long way, which is the complement theta distance:

Q(t) = Exp(t CompTheta V/2)

Where CompTheta = theta - 2PI.

Recall that exponential function is:

Exp(t CompTheta V/2) = cos(t CompTheta/2) - sin(t CompTheta/2) V

Now, how do we find the Logarithm i.e., theta V?

When you multiply q1^-1 q2 you get a new quaternion, lets call it q12.

q12 = cos(theta/2) - sin(theta/2) V

q12 = w + V'

Where:

w = cos(theta/2)

V' = sin(theta/2) V

theta = 2 atan2(|V'|, w)

V = V'/|V'|

So finally your SLERP2(q1,q2, t) is equal to:

SLERP2(q1,q2, t) = q1 Q(t)

SLERP2(q1,q2, t) = q1 Exp(t CompTheta V/2)

Discraimler: I haven't tested this. If you can test it please comment here.