基于四元数的点旋转在GLM中的应用

Question

基于四元数的点旋转在GLM中的应用

3

我正在尝试使用GLM中实现的四元数来旋转一个点。最终目标是使用这段代码创建轨道相机，但这只是一个辅助说明代码背后的动机。

为了更好地理解基于四元数的旋转，我编写了一些包含两个循环的代码。第一个循环将通过绕X轴旋转步骤逐渐改变四元数的方向，直到90度；第二个循环将继续在Z轴周围以步骤方式应用旋转，直到90度。每个循环执行4个步骤。因此，每个循环分别围绕它们各自的轴增量旋转90/4 = 22.5度。使用四元数乘法应用方向变化，并使用欧拉角跟踪。循环应以一个四元数结束，该四元数将把点(0, 0, 3)旋转到(3, 0, 0)。请注意，我不仅仅想确定将执行此旋转的四元数。目标是执行一系列增量旋转。

如果我们看下面的图片，从C到I的过渡发生在第一个循环中，然后从I到R的过渡发生在第二个循环中（请原谅稀疏的点命名）。

一个点的旋转被定义为（参见这里和这里）:

v' = q * v * q^-1

在这里，v应被视为一个纯四元数（其标量项w为零），而q需要是长度为1的单位四元数。据我所知，需要用四元数的逆进行右手乘法，以保持结果v'在三维空间中，并且不会得到一个4D向量。因此，v'也需要是一个纯四元数。

然后是旋转的加倍效应，左手乘法与q贡献了所需旋转的一半，右手乘法与逆元素添加了所需旋转的另一半。

这里有一个由Ben Eater和Grant Sanderson提供的出色交互式可视化和四元数解释，我将其用作交叉参考。它可以在这里找到。

因此，我们首先需要使用绕X轴旋转11.25度的四元数，GLM为欧拉角返回了此四元数（使用四元数表示[w，[x，y，z]]）：

Rotation of [ 11.25, 0.00,  0.00] deg => Q: [ 0.9952, [ 0.0980,  0.0000,  0.0000]]

根据 this，由于我们只是围绕X轴纯粹旋转，因此可以通过对四元数的w分量执行acos来验证GLM计算的四元数的旋转量：

float angle = acosf(q.w)

然后：

acos(0.9952) = 0.0980 rad / 5.6 degrees

这是所需角度的一半...在与交互式动画进行交叉验证时也得到了确认（对舍入表示歉意）。

GLM返回的四元数实际上只旋转了所需角度的一半... 如果我们查看GLM代码，从欧拉角计算w分量的过程会更加复杂，因为旋转可以发生在任意轴周围... 但是欧拉角确实减半了：

template <typename T, precision P>
GLM_FUNC_QUALIFIER tquat<T, P>::tquat(tvec3<T, P> const & eulerAngle)
{
    tvec3<T, P> c = glm::cos(eulerAngle * T(0.5));
    tvec3<T, P> s = glm::sin(eulerAngle * T(0.5));
    
    this->w = c.x * c.y * c.z + s.x * s.y * s.z;
    this->x = s.x * c.y * c.z - c.x * s.y * s.z;
    this->y = c.x * s.y * c.z + s.x * c.y * s.z;
    this->z = c.x * c.y * s.z - s.x * s.y * c.z;
}

我的第一个问题是为什么GLM要将角度减半？

尽管所需旋转角度不同，我还是进行了两个循环的旋转结果检查。而结果是...出乎意料。

如果我使用“不正确的”旋转形式（一些OpenGL在线教程建议的形式）并且仅通过左手乘法旋转点（但为22.5度的完整步骤）：

v' = q * v

我得到了我期望的结果。重点是正确地遵循所有中间步骤，并从（0，0，3）到达（3，0，0）。此外，在所有中间步骤中，w分量均为0。

但是，如果我使用“正确的”旋转形式，并通过左手乘以q和右手乘以q的逆来旋转该点（为了解决旋转加倍的问题而进行半步旋转，角度为11.25度）：

v' = q * v * q^-1

当第二个循环开始绕Z轴旋转点时，我开始得到错误的结果。一个小但明显的Z分量开始悄悄地出现，旋转仅短了22.5度的全步骤。这在下面图像中的绿色点中可见。

对于旋转的两种方法，旋转点的w分量保持为0...

有人能解释一下为什么使用GLM旋转可以通过单个左乘获得正确结果吗？

这是一种将操作次数减少到最少的优化吗？

我能否在GLM中使用v' = q * v旋转以获得所有旋转的一致和正确结果？

代码：

const int rotSteps = 4;
// Rotate around X axis in steps to 90deg
vec3 eulerState = vec3(0.0f);
// point we want to rotate (use vec4 to track the w component during rotations)
vec4 v = vec4(0.0f, 0.0f, 3.0f, 0.0f);

// Full Euler steps for q * v rotation
quat orientF   = quat(1.0f, 0.0f, 0.0f, 0.0f);
vec3 euler     = vec3(RAD(90.0f), RAD(0.0f), RAD(0.0f));
vec3 eulerStep = euler / (float)rotSteps;
quat qEulerF   = quat(eulerStep); // GetRotQuat(eulerStep);

vec4 qa          = ToAngularForm(qEulerF);
vec3 orientEuler = eulerAngles(qEulerF);
CLogD(TAG, "Rot Full Step    Q [W, X, Y, Z]: " FMT_Q(4)  " / " FMT_V3(2) "deg / " FMT_QA(2), PAR_Q(qEulerF), PAR_V3(degrees(orientEuler)), PAR_QA(qa));

// Half Euler steps for q * v * q^-1 rotation
quat orientH    = quat(1.0f, 0.0f, 0.0f, 0.0f);
vec3 eulerStepH = eulerStep / 2.0f;
quat qEulerH    = quat(eulerStepH); // GetRotQuat(eulerStepH);

qa          = ToAngularForm(qEulerH);
orientEuler = eulerAngles(qEulerH);
CLogD(TAG, "Rot Half Step    Q [W, X, Y, Z]: " FMT_Q(4) " / " FMT_V3(2) "deg / " FMT_QA(2), PAR_Q(qEulerH), PAR_V3(degrees(orientEuler)), PAR_QA(qa));

quat qEulerHI = inverse(qEulerH);
vec4 qai      = ToAngularForm(qEulerHI);
orientEuler   = eulerAngles(qEulerHI);
CLogD(TAG, "Rot Half Step Q^-1 [W, X, Y, Z]: " FMT_Q(4) " / " FMT_V3(2) "deg / " FMT_QA(2), PAR_Q(qEulerHI), PAR_V3(degrees(orientEuler)), PAR_QA(qai));


for (int rotStep = 1; rotStep <= rotSteps; ++rotStep)
{
    // Track the absolute Euler rotation
    eulerState += eulerStep;
    // Rotate by incremental rotation as defined by Euler angles
    orientH = qEulerH * orientH;
    orientEuler = eulerAngles(orientH);
    CLogI(TAG, "Rot Step %d. Curr Abs Q: " FMT_Q(4) "/" FMT_V3(2) "deg, Abs Euler: " FMT_V3(2) "deg",
          rotStep, PAR_Q(orientH), PAR_V3(degrees(orientEuler)), PAR_V3(degrees(eulerState)));

    // Transform the point using the correct q * v * q^-1 rotation and multiply from Left and Right
    quat orientHI = inverse(orientH);
    qa  = ToAngularForm(orientH);
    qai = ToAngularForm(orientHI);

    vec4 rotV = orientH * v * orientHI;
    CLogD(TAG, "Rot      QL: " FMT_Q(4) " / " FMT_QA(1), PAR_Q(orientH), PAR_QA(qa));
    CLogD(TAG, "Rot      QR: " FMT_Q(4) " / " FMT_QA(1), PAR_Q(orientHI), PAR_QA(qai));
    CLogD(TAG, "Rot LR   -> " FMT_V4(1), PAR_V4(rotV));

    // Transform the point using the incorrect q * v rotation and multiply from Left only
    orientF = qEulerF * orientF;
    qa      = ToAngularForm(orientF);

    rotV = orientF * v;
    CLogD(TAG, "Rot      QR: " FMT_Q(4) " / " FMT_QA(1), PAR_Q(orientF), PAR_QA(qa));
    CLogD(TAG, "Rot L    -> " FMT_V4(1), PAR_V4(rotV));
}


// Rotate for 90 degrees around the Z axis
// Full Euler steps for q * v rotation
euler = vec3(RAD(0.0f), RAD(0.0f), RAD(90.0f));
eulerStep = euler / (float)rotSteps;
qEulerF = quat(eulerStep); // GetRotQuat(eulerStep);

qa = ToAngularForm(qEulerF);
orientEuler = eulerAngles(qEulerF);
CLogD(TAG, "Rot Full Step    Q [W, X, Y, Z]: " FMT_Q(4)  " / " FMT_V3(2) "deg / " FMT_QA(2), PAR_Q(qEulerF), PAR_V3(degrees(orientEuler)), PAR_QA(qa));

// Half Euler steps for q * v * q^-1 rotation
eulerStepH = eulerStep / 2.0f;
qEulerH = quat(eulerStepH); // GetRotQuat(eulerStepH);

qa = ToAngularForm(qEulerH);
orientEuler = eulerAngles(qEulerH);
CLogD(TAG, "Rot Half Step    Q [W, X, Y, Z]: " FMT_Q(4) " / " FMT_V3(2) "deg / " FMT_QA(2), PAR_Q(qEulerH), PAR_V3(degrees(orientEuler)), PAR_QA(qa));

qEulerHI = inverse(qEulerH);
qai = ToAngularForm(qEulerHI);
orientEuler = eulerAngles(qEulerHI);
CLogD(TAG, "Rot Half Step Q^-1 [W, X, Y, Z]: " FMT_Q(4) " / " FMT_V3(2) "deg / " FMT_QA(2), PAR_Q(qEulerHI), PAR_V3(degrees(orientEuler)), PAR_QA(qai));


for (int rotStep = 1; rotStep <= rotSteps; ++rotStep)
{
    // Track the absolute Euler rotation
    eulerState += eulerStep;
    // Rotate by incremental rotation as defined by Euler angles
    orientH = qEulerH * orientH;
    orientEuler = eulerAngles(orientH);
    CLogI(TAG, "Rot Step %d. Curr Abs Q: " FMT_Q(4) "/" FMT_V3(2) "deg, Abs Euler: " FMT_V3(2) "deg",
        rotStep, PAR_Q(orientH), PAR_V3(degrees(orientEuler)), PAR_V3(degrees(eulerState)));

    // Transform the point using the correct q * v * q^-1 rotation and multiply from Left and Right
    quat orientHI = inverse(orientH);
    qa = ToAngularForm(orientH);
    qai = ToAngularForm(orientHI);

    vec4 rotV = orientH * v * orientHI;
    CLogD(TAG, "Rot      QL: " FMT_Q(4) " / " FMT_QA(1), PAR_Q(orientH), PAR_QA(qa));
    CLogD(TAG, "Rot      QR: " FMT_Q(4) " / " FMT_QA(1), PAR_Q(orientHI), PAR_QA(qai));
    CLogD(TAG, "Rot LR   -> " FMT_V4(1), PAR_V4(rotV));

    // Transform the point using the incorrect q * v rotation and multiply from Left only
    orientF = qEulerF * orientF;
    qa = ToAngularForm(orientF);

    rotV = orientF * v;
    CLogD(TAG, "Rot      QR: " FMT_Q(4) " / " FMT_QA(1), PAR_Q(orientF), PAR_QA(qa));
    CLogD(TAG, "Rot L    -> " FMT_V4(1), PAR_V4(rotV));
}

输出：

Rot Full Step    Q [W, X, Y, Z]: [ 0.9808, [ 0.1951,  0.0000,  0.0000]] / [ 22.50, -0.00,  0.00]deg / cos( 11.25) + sin( 11.25)( 1.00i +  0.00j +  0.00k)
Rot Half Step    Q [W, X, Y, Z]: [ 0.9952, [ 0.0980,  0.0000,  0.0000]] / [ 11.25, -0.00,  0.00]deg / cos( 5.63) + sin( 5.63)( 1.00i +  0.00j +  0.00k)
Rot Half Step Q^-1 [W, X, Y, Z]: [ 0.9952, [-0.0980, -0.0000, -0.0000]] / [-11.25, -0.00,  0.00]deg / cos( 5.63) + sin( 5.63)(-1.00i + -0.00j + -0.00k)
Rot Step 1. Curr Abs Q: [ 0.9952, [ 0.0980,  0.0000,  0.0000]]/[ 11.25, -0.00,  0.00]deg, Abs Euler: [ 22.50,  0.00,  0.00]deg
Rot      QL: [ 0.9952, [ 0.0980,  0.0000,  0.0000]] / cos( 5.6) + sin( 5.6)( 1.0i +  0.0j +  0.0k)
Rot      QR: [ 0.9952, [-0.0980, -0.0000, -0.0000]] / cos( 5.6) + sin( 5.6)(-1.0i + -0.0j + -0.0k)
Rot LR   -> [ 0.0, -1.1,  2.8,  0.0]
Rot      QR: [ 0.9808, [ 0.1951,  0.0000,  0.0000]] / cos( 11.3) + sin( 11.3)( 1.0i +  0.0j +  0.0k)
Rot L    -> [ 0.0, -1.1,  2.8,  0.0]
Rot Step 2. Curr Abs Q: [ 0.9808, [ 0.1951,  0.0000,  0.0000]]/[ 22.50, -0.00,  0.00]deg, Abs Euler: [ 45.00,  0.00,  0.00]deg
Rot      QL: [ 0.9808, [ 0.1951,  0.0000,  0.0000]] / cos( 11.3) + sin( 11.3)( 1.0i +  0.0j +  0.0k)
Rot      QR: [ 0.9808, [-0.1951, -0.0000, -0.0000]] / cos( 11.2) + sin( 11.2)(-1.0i + -0.0j + -0.0k)
Rot LR   -> [ 0.0, -2.1,  2.1,  0.0]
Rot      QR: [ 0.9239, [ 0.3827,  0.0000,  0.0000]] / cos( 22.5) + sin( 22.5)( 1.0i +  0.0j +  0.0k)
Rot L    -> [ 0.0, -2.1,  2.1,  0.0]
Rot Step 3. Curr Abs Q: [ 0.9569, [ 0.2903,  0.0000,  0.0000]]/[ 33.75, -0.00,  0.00]deg, Abs Euler: [ 67.50,  0.00,  0.00]deg
Rot      QL: [ 0.9569, [ 0.2903,  0.0000,  0.0000]] / cos( 16.9) + sin( 16.9)( 1.0i +  0.0j +  0.0k)
Rot      QR: [ 0.9569, [-0.2903, -0.0000, -0.0000]] / cos( 16.9) + sin( 16.9)(-1.0i + -0.0j + -0.0k)
Rot LR   -> [ 0.0, -2.8,  1.1,  0.0]
Rot      QR: [ 0.8315, [ 0.5556,  0.0000,  0.0000]] / cos( 33.8) + sin( 33.8)( 1.0i +  0.0j +  0.0k)
Rot L    -> [ 0.0, -2.8,  1.1,  0.0]
Rot Step 4. Curr Abs Q: [ 0.9239, [ 0.3827,  0.0000,  0.0000]]/[ 45.00, -0.00,  0.00]deg, Abs Euler: [ 90.00,  0.00,  0.00]deg
Rot      QL: [ 0.9239, [ 0.3827,  0.0000,  0.0000]] / cos( 22.5) + sin( 22.5)( 1.0i +  0.0j +  0.0k)
Rot      QR: [ 0.9239, [-0.3827, -0.0000, -0.0000]] / cos( 22.5) + sin( 22.5)(-1.0i + -0.0j + -0.0k)
Rot LR   -> [ 0.0, -3.0,  0.0,  0.0]
Rot      QR: [ 0.7071, [ 0.7071,  0.0000,  0.0000]] / cos( 45.0) + sin( 45.0)( 1.0i +  0.0j +  0.0k)
Rot L    -> [ 0.0, -3.0,  0.0,  0.0]

Rot Full Step    Q [W, X, Y, Z]: [ 0.9808, [ 0.0000,  0.0000,  0.1951]] / [ 0.00, -0.00,  22.50]deg / cos( 11.25) + sin( 11.25)( 0.00i +  0.00j +  1.00k)
Rot Half Step    Q [W, X, Y, Z]: [ 0.9952, [ 0.0000,  0.0000,  0.0980]] / [ 0.00, -0.00,  11.25]deg / cos( 5.63) + sin( 5.63)( 0.00i +  0.00j +  1.00k)
Rot Half Step Q^-1 [W, X, Y, Z]: [ 0.9952, [-0.0000, -0.0000, -0.0980]] / [ 0.00, -0.00, -11.25]deg / cos( 5.63) + sin( 5.63)(-0.00i + -0.00j + -1.00k)
Rot Step 1. Curr Abs Q: [ 0.9194, [ 0.3808,  0.0375,  0.0906]]/[ 45.00,  0.00,  11.25]deg, Abs Euler: [ 90.00,  0.00,  22.50]deg
Rot      QL: [ 0.9194, [ 0.3808,  0.0375,  0.0906]] / cos( 23.2) + sin( 23.2)( 1.0i +  0.1j +  0.2k)
Rot      QR: [ 0.9194, [-0.3808, -0.0375, -0.0906]] / cos( 23.2) + sin( 23.2)(-1.0i + -0.1j + -0.2k)
Rot LR   -> [ 1.0, -2.8,  0.0,  0.0]
Rot      QR: [ 0.6935, [ 0.6935,  0.1379,  0.1379]] / cos( 46.1) + sin( 46.1)( 1.0i +  0.2j +  0.2k)
Rot L    -> [ 1.1, -2.8,  0.0,  0.0]
Rot Step 2. Curr Abs Q: [ 0.9061, [ 0.3753,  0.0747,  0.1802]]/[ 45.00, -0.00,  22.50]deg, Abs Euler: [ 90.00,  0.00,  45.00]deg
Rot      QL: [ 0.9061, [ 0.3753,  0.0747,  0.1802]] / cos( 25.0) + sin( 25.0)( 0.9i +  0.2j +  0.4k)
Rot      QR: [ 0.9061, [-0.3753, -0.0747, -0.1802]] / cos( 25.0) + sin( 25.0)(-0.9i + -0.2j + -0.4k)
Rot LR   -> [ 1.9, -2.4,  0.1,  0.0]
Rot      QR: [ 0.6533, [ 0.6533,  0.2706,  0.2706]] / cos( 49.2) + sin( 49.2)( 0.9i +  0.4j +  0.4k)
Rot L    -> [ 2.1, -2.1,  0.0,  0.0]
Rot Step 3. Curr Abs Q: [ 0.8841, [ 0.3662,  0.1111,  0.2682]]/[ 45.00,  0.00,  33.75]deg, Abs Euler: [ 90.00,  0.00,  67.50]deg
Rot      QL: [ 0.8841, [ 0.3662,  0.1111,  0.2682]] / cos( 27.9) + sin( 27.9)( 0.8i +  0.2j +  0.6k)
Rot      QR: [ 0.8841, [-0.3662, -0.1111, -0.2682]] / cos( 27.9) + sin( 27.9)(-0.8i + -0.2j + -0.6k)
Rot LR   -> [ 2.5, -1.6,  0.3,  0.0]
Rot      QR: [ 0.5879, [ 0.5879,  0.3928,  0.3928]] / cos( 54.0) + sin( 54.0)( 0.7i +  0.5j +  0.5k)
Rot L    -> [ 2.8, -1.1,  0.0,  0.0]
Rot Step 4. Curr Abs Q: [ 0.8536, [ 0.3536,  0.1464,  0.3536]]/[ 45.00,  0.00,  45.00]deg, Abs Euler: [ 90.00,  0.00,  90.00]deg
Rot      QL: [ 0.8536, [ 0.3536,  0.1464,  0.3536]] / cos( 31.4) + sin( 31.4)( 0.7i +  0.3j +  0.7k)
Rot      QR: [ 0.8536, [-0.3536, -0.1464, -0.3536]] / cos( 31.4) + sin( 31.4)(-0.7i + -0.3j + -0.7k)
Rot LR   -> [ 2.9, -0.7,  0.4,  0.0]
Rot      QR: [ 0.5000, [ 0.5000,  0.5000,  0.5000]] / cos( 60.0) + sin( 60.0)( 0.6i +  0.6j +  0.6k)
Rot L    -> [ 3.0,  0.0,  0.0,  0.0]

- omahena

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- omahena · Accepted Answer

我已经得到了我问题的答案和一个工作正常的轨道相机，但还没有时间再次检查样例代码是否正确——它应该是正确的。

第一个问题是为什么GLM在四元数转换过程中会将角度减半，根据扩展欧拉公式来看，它必须这样做。虽然这部分可能需要进一步调查，但由于时间不足，我不得不接受它。

GLM中的向量旋转是使用乘法运算符实现的。这意味着将vec3乘以四元数时不会将向量转换为四元数然后执行乘法，而是执行向量旋转：

template <typename T, precision P>
GLM_FUNC_QUALIFIER tvec3<T, P> operator*(tquat<T, P> const & q, tvec3<T, P> const & v)
{
    tvec3<T, P> const QuatVector(q.x, q.y, q.z);
    tvec3<T, P> const uv(glm::cross(QuatVector, v));
    tvec3<T, P> const uuv(glm::cross(QuatVector, uv));

    return v + ((uv * q.w) + uuv) * static_cast<T>(2);
}

因此，使用四元数旋转向量的正确方法是通过在四元数和向量之间使用乘法运算符，像这样：

v' = q * v

或者在C++中：

vec3 posOrigin;
quat rotQ;
...
vec3 posRot = rotQ * posOrigin;

这段代码实际上并没有直接执行四元数乘法，而是进行了旋转。个人认为GLM提供像 rotate(quat, vec) 这样的函数调用更为理想... 但我相信重载运算符也有其存在的必要。

另外需要注意的是操作数的顺序是有影响的，因为向量和四元数之间的乘法定义如下：

template <typename T, precision P>
GLM_FUNC_QUALIFIER tvec3<T, P> operator*(tvec3<T, P> const & v, tquat<T, P> const & q)
{
    return glm::inverse(q) * v;
}

因此，它将以相反的方式旋转向量。

请注意，GLM还实现了四元数之间的乘法，但是需要使用两个四元数之间的乘法运算符：

template <typename T, precision P>
template <typename U>
GLM_FUNC_QUALIFIER tquat<T, P> & tquat<T, P>::operator*=(tquat<U, P> const & r)
{
    tquat<T, P> const p(*this);
    tquat<T, P> const q(r);

    this->w = p.w * q.w - p.x * q.x - p.y * q.y - p.z * q.z;
    this->x = p.w * q.x + p.x * q.w + p.y * q.z - p.z * q.y;
    this->y = p.w * q.y + p.y * q.w + p.z * q.x - p.x * q.z;
    this->z = p.w * q.z + p.z * q.w + p.x * q.y - p.y * q.x;
    return *this;
}

由于GLM文档短缺，我无法找到很多资料。这种运算符重载会导致错误的假设和时间的显著损失。因此，我认为我应该阅读GLM代码，而不是假设它的功能...