使用 GLM 的基于四元数的点旋转答案

【问题标题】：Quaternion based point rotations using GLM使用 GLM 的基于四元数的点旋转
【发布时间】：2020-11-06 07:19:02
【问题描述】：

我正在尝试使用 GLM 中实现的四元数来旋转一个点。最终目标是使用此代码创建一个轨道相机，但这是帮助理解代码背后的动机的旁注。

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

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

点的旋转定义为（参见here 和here）：

v' = q * v * q^-1

其中 v 应被视为纯四元数（标量项 w 为零），q 需要是单位四元数（长度为 1）。根据我的理解，需要与四元数的逆相乘来将结果 v' 保持在 3D 空间中，而不是最终得到 4D 向量。所以 v' 也需要是一个纯四元数。

然后是旋转的加倍效果，其中左手乘以 q 贡献一半的期望旋转，右手乘以逆增加另一半期望旋转。

Ben Eater 和 Grant Sanderson 对四元数进行了出色的交互式可视化和解释，我将其用作交叉参考。可以找到here。

所以我们首先需要使用一个围绕 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 返回的 11.25 度的四元数实际上旋转了所需角度的一半...如果我们查看 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]

【问题讨论】：

标签： c++ math rotation quaternions glm-math

【解决方案1】：

我有我的问题的答案和一个正常工作的轨道相机，但没有时间仔细检查示例代码现在是否可以正常工作 - 应该。

第一个问题是为什么 GLM 在四元数转换期间将角度减半，而根据 extended Euler's formula... 看起来它必须这样做。这部分可能需要更多调查，但由于时间不够，我不得不接受它。

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

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 代码而不是假设它做了什么......

【讨论】：