在两个四元数之间插值很长的路

Question

我们可以像这样在两个四元数之间进行 slerp 插值：

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;
}

这将以最短的方式在两个四元数之间进行插值。然而，在四元数之间进行插值还有很长的路要走。如下图所示（来源Maya）。

我们如何进行长距离内插？

Answer 1

单位四元数的性质以及它们映射到 3D 旋转的方式意味着它们可以用两种方式描述每个 3D 旋转值 - 如q(r, v') 和 q(-r, -v')（将它们想象为轴角旋转 - 反转两个轴并且角度会导致相同的 3D 旋转）。

四元数实际上是 4D 单位球面上的点，这两个值表示该球面上的对跖点。

对于两个四元数的 slerp（或 nlerp）要遵循最短路径，相应的 4D 点必须位于 4D 球体的同一半球（这也是为什么加权平均超过 2 个四元数没有唯一解）。这映射到非负点积，通常在插值代码中进行测试。

简单地否定一个源四元数会给你一个“在 4D 球体的另一侧”的点，并导致插值“很长的路”（并解释为什么否定插值参数会导致相同的结果） .

Answer 2

这可以通过改变球面插值的角度来完成。

在通常的 SLERP(q1, q2, t) 中，当 t=0 时得到 q1，当 t=1 时得到 q2。行进的测地线距离实际上是 q1 和 q2 之间的角度，我们将其命名为 theta。

我们在这里要做的是移动补码距离，即 2PI - theta，但旋转方向相反。我们将其称为补集 theta。

我们想找到一个四元值函数 Q(t) 使得：

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

当 t = 0 时

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

当 t=1

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

所以我们知道Q(0) = 1（恒等四元数）和Q(1) = (q1^-1 q2)。

事实证明，我们可以根据四元数的指数映射和主对数来定义这样的函数 Q(t)：

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

您可以通过为 t 赋予值（例如 t=0 和 t=1）来检查它是否有效。

到目前为止还不错，但是 Q(t) 将导致我们拥有常规的 SLERP，而不是我们想要的。让我们仔细看看对数：

Log(q1^-1 q2) = theta V

其中V是一个单位向量（实际上是纯单位四元数），它是四元数q1^-1 q2的旋转轴。而 theta 基本上是 q1 和 q2 之间的角度。

我们实际上需要改变这个对数，这样 Q(t) 才能走得更远，这就是补码 theta 距离：

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

其中 CompTheta = theta - 2PI。

回想一下指数函数是：

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

现在，我们如何找到对数，即 theta V？

当你将 q1^-1 q2 相乘时，你会得到一个新的四元数，我们称之为 q12。

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

q12 = w + V'

地点：

w = cos(theta/2)

V' = sin(theta/2) V

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

V = V'/|V'|

所以最后你的 SLERP2(q1,q2, t) 等于：

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

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

Discraimler：我还没有测试过这个。如果您可以测试它，请在此处发表评论。

Answer 3

因此，为了避免其他可能遇到此问题的人感到困惑，我拼凑了一些 SDL/OpenGL 代码，以使 Mauricio 或 Martin 的答案能够正常工作。我找到了马丁的答案，因为它在实施时有点模糊，即使它陈述了事实。不幸的是，即使在毛里西奥的帮助下，我也没有设法得到他的答案。

我也犯了一些错误，我从不同的地方尝试了许多不同的 slerp 函数来进行健全性检查，结果导致我有些困惑，所以我最终从头开始实现自己的 slerp（代码中的 SlerpIam()）不检查最近的路径。

在代码 Slerp1() 和 Slerp2() 中，我认为在没有选择最短路径时会被破坏，这是我困惑的一部分 - 从我发现的无数 slerp 中，我认为它们被错误地修改为尝试和支持最长的路径，但他们没有。所以我最初是想像 Martin 提到的那样破解它们，但结果大错特错。

我的测试用例显示一个点围绕 Z 轴旋转/滚动到 270 度。

我在 Windows 上使用 SDL2 编译了代码，您需要包含 SDL 标头和链接等：

#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;
}

所以从头开始使用我自己的 SlerpIam()，我认为我的理智已经恢复，马丁的回答本质上是正确的。我得到以下我认为正确的函数（注意它们目前不处理小角度 lerp 回退）：

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;
}