在两个四元数之间进行插值

问题描述 投票：0回答：3

我们可以像这样在两个四元数之间进行 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）。

我们如何插值长距离？

c++ c interpolation linear-algebra quaternions
3个回答
3
投票

单位四元数的性质以及它们映射到 3D 旋转的方式意味着它们可以用两种方式描述每个 3D 旋转值 - 作为 

q(r, v')
和 
q(-r, -v')
（将它们想象为轴角旋转 - 反转轴和角度导致相同的 3D 旋转）。

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

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

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

1
投票

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

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

我们想要做的是移动补距离，即 2PI - 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) = θV

其中V是单位向量（实际上是纯单位四元数），它是四元数q1^-1 q2的旋转轴。 θ 基本上是 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(θ/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：我还没有测试过这个。如果你可以测试一下，请在这里评论。

1
投票

因此，为了避免其他可能遇到此问题的人感到困惑，我拼凑了一些 SDL/OpenGL 代码，努力让 Mauricio 或 Martin 的答案发挥作用。我发现马丁的答案是因为它在实施时有点模糊，尽管它陈述了事实。不幸的是，即使在毛里西奥的帮助下，我也没能得到他的答案。

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

在代码 Slerp1() 和 Slerp2() 中，我认为当未选择最短路径时就被破坏了，这是我困惑的一部分 - 从无数的 slerps 中我发现我认为它们被错误地修改为尝试支持最长路径但他们没有。所以我最初试图像马丁提到的那样破解它们，但结果出了严重的错误。

我的测试用例显示一个点绕 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;
}
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.