Creating Quaternions 创建四元数

整理从角轴、from-to shortest arc、Look At 以及旋转矩阵创建四元数的方法。

May 11, 2025约 12 分钟Quaternionsysysimon

Creating Quaternions 创建四元数

我们一般直接使用 角度-转轴 信息来创建四元数，我们想要表示的是绕 给定轴向 $\hat{n}$ 进行 逆时针 旋转 $\theta^\circ$ 的行为

根据我们前面章节的推导，我们知道用于表示旋转的 四元数 是一个特殊的 单位四元数，它的标量分量 $w$ 并不直接等于旋转的角度 $\theta^\circ$ ，但却与之相关，准确来说是

w = \cos{\frac{\theta}{2}}

旋转发生在与转轴正交的平面内，也就是说平面的 法向量 就是转轴方向，为了保持 四元数 的长度为 $1$ ，需要对转轴 (单位向量) 进行大小为 $\sin{\frac{\theta}{2}}$ 的缩放才能保持 四元数 整体为 单位四元数，这也在我们之前的推导中有所体现

\sqrt{\sin^2{\theta} + \cos^2{\theta}} = 1

正的角度表示 逆时针 旋转

Quaternion AngleAxis(float degrees, Vec3 axis) {
    float radians = degrees * 0.0174533f;
 
    if (MagnitudeSq(axis) != 1) { // Do epsilon check here!
        axis = Normalize(axis);
    }
 
    Quaternion result;
    result.x = axis.x * sinf(radians * 0.5f);
    result.y = axis.y * sinf(radians * 0.5f);
    result.z = axis.z * sinf(radians * 0.5f);
    result.w = cosf(radians * 0.5f);
 
    return result;
}

从一个方向转向另一个方向的最短弧 shortest arc of from-to

如果我们已知两个 归一化 的方向向量 $\vec{p}_0$ 和 $\vec{p}_1$ ，如何如何求出 “最短”、“唯一” 的那一次旋转，把 $\vec{p}_0$ 转到 $\vec{p}_1$ ？实际上存在无数种旋转，但是我一般来说想要的是代表 最短弧旋转 的那个 四元数

从几何上看，可以将 单位向量 看成位于 单位球面 上的 “点”，连接 $\vec{p}_0$ , $\vec{p}_1$ 以及球心形成的平面，球面与该平面的交线是一条大圆弧；在这条弧上的弧长就是 “最短弧”

Shortest arc quaternion construction between two unit direction vectors on a sphere. — 把两个单位方向看成单位球面上的点，连接它们的大圆短弧给出 From-To 旋转的最短路径。

把两个单位方向看成单位球面上的点，连接它们的大圆短弧给出 From-To 旋转的最短路径。

该弧所在的平面垂直于一条轴，这条轴的方向就是 转轴方向，它可以借助叉乘得出，我们知道叉乘是

\vec{p}_0 \times \vec{p}_1 = \| \vec{p}_0 \| \| \vec{p}_1 \| \sin{\theta} \ \cdot \hat{n}

向量叉乘给出的结果天然正交于两个输入向量

因为我们假设 $\vec{p}_0$ 和 $\vec{p}_1$ 都是 单位向量，所以 $\| \vec{p}_0 \| = \| \vec{p}_1 \| = 1$ ，那么

\vec{p}_0 \times \vec{p}_1 = \sin{\theta} \ \cdot \hat{n}

我们可以求得转轴

\mathrm{axis} = \hat{n} = \frac{ \vec{p}_0 \times \vec{p}_1 }{ \sin{\theta} }

而两向量的夹角，也就是应当旋转的角度，则可以由点乘得出

\vec{p}_0 \cdot \vec{p}_1 = \cos{\theta} \implies \theta = \cos^{-1}{(\vec{p}_0 \cdot \vec{p}_1)}

图中的 $\frac{\theta}{2}$ 是指如果使用四元数，那么应该使用半角

如此一来我们就能得到 角度-转轴 信息来构造表达这个旋转的 四元数，但是由于 $\cos^{-1}$ 或者说 acos 在运行时计算的很慢，我们可以用一些技巧来避免计算它，我们知道

\begin{align*} \vec{p}_0 \times \vec{p}_1 &= \sin{\theta} \ \cdot \mathrm{axis} \\ \vec{p}_0 \cdot \vec{p}_1 &= \cos{\theta} \end{align*}

而要构造从 $\vec{p}_0$ 旋转到 $\vec{p}_1$ 的 四元数，根据我们前面章节的推导可得

q = (x,y,z,w) \implies \begin{cases} x = \sin{\frac{\theta}{2}} \ \mathrm{axis}_x \\ y = \sin{\frac{\theta}{2}} \ \mathrm{axis}_y \\ z = \sin{\frac{\theta}{2}} \ \mathrm{axis}_z \\ w = \cos{\frac{\theta}{2}} \end{cases}

此时向量部分的 $x,y,z$ 近乎是叉乘 $\vec{p}_0 \times \vec{p}_1$ ，而标量部分的 $w$ 则近乎是点乘 $\vec{p}_0 \cdot \vec{p}_1$ ，之所以说是 “近乎” 是因为构造这个 四元数 需要使用半角 $\frac{\theta}{2}$ ，而 $\vec{p}_0 \times \vec{p}_1$ 和 $\vec{p}_0 \cdot \vec{p}_1$ 中使用的则是原始的全角 $\theta$ ，我们有两种办法来等到等价的半角结果

半程向量
将 归一化 的 $\vec{p}_0$ 和 $\vec{p}_1$ 相加之后再 归一化，得到 半程向量 $\vec{h}$
$\vec{h} = \frac{\vec{p}_0 + \vec{p}_1}{\|\vec{p}_0 + \vec{p}_1\|}$
然后只求 $\vec{p}_0$ 旋转到 $\vec{h}$ 的 四元数 $q$ ，这本质上就是 “半角” $\frac{\theta}{2}$ 了
$q= \begin{cases} x = \sin{\frac{\theta}{2}} \ \mathrm{axis}_x \\ y = \sin{\frac{\theta}{2}} \ \mathrm{axis}_y \\ z = \sin{\frac{\theta}{2}} \ \mathrm{axis}_z \\ w = \cos{\frac{\theta}{2}} \end{cases} = \begin{cases} x = (\vec{p_0} \times \vec{h})_x \\ y = (\vec{p_0} \times \vec{h})_y \\ z = (\vec{p_0} \times \vec{h})_z \\ w = \vec{p_0} \cdot \vec{h} \end{cases}$
```
Quaternion FromToRotation(Vector3 from, Vector3 to) {
Vector3 p0 = Normalize(from);
Vector3 p1 = Normalize(to);
 
if (p0 == -p1) {
    Vector3 mostOrthogonal = Vector3(1, 0, 0);
 
    // 分量最小的那个轴就是最 正交 于 p0 的轴
    if (abs(p0.y) < abs(p0.x)) {
        mostOrthogonal = Vector3(0, 1, 0);
    }
 
    if (abs(p0.z) < abs(p0.y) && abs(p0.z) < abs(p0.x)) {
        mostOrthogonal = Vector3(0, 0, 1)
    }
 
    Vector3 axis = Normalize(Cross(p0, mostOrthogonal));
    return Quaternion(axis.x, axis.y, axis.z, 0);
}
 
Vector3 half = Normalize(p0 + p1);
Vector3 axis = Cross(p0, half);
 
Quaternion result;
result.x = axis.x;
result.y = axis.y;
result.z = axis.z;
result.w = Dot(p0, half);
 
return result;
}
```
一个需要处理的问题是，如果两个向量方向完全相反 (即 $\vec{p}_0 = -\vec{p}_1$ ) 这意味着旋转角度为 $180^\circ$ ，在这种情况下，直接通过这两个向量的叉积无法计算旋转轴，因此必须通过一个与 $\vec{p}_0$ 正交的向量来构建旋转轴，做法是挑一个与 $\vec{p}_0$ 最不平行的坐标轴 (分量最小的那个轴)，再取 axis = normalize(p0 × mostOrthogonal) 生成 $180^\circ$ 四元数 $(\mathrm{axis},0)$
平方根整旋转
算出全角的 完整旋转四元数 $q$ ，也就是说直接使用 $\vec{p}_0 \times \vec{p}_1$ 和 $\vec{p}_0 \cdot \vec{p}_1$ 来计算
$q = \begin{cases} x = (\vec{p}_0 \times \vec{p}_1)_x \\ y = (\vec{p}_0 \times \vec{p}_1)_y \\ z = (\vec{p}_0 \times \vec{p}_1)_z \\ w = \vec{p}_0 \cdot \vec{p}_1 \end{cases}$
然后，想要 “旋转减半” 就是要一个新 四元数 $\sqrt{q}$ ，满足 $q = \sqrt{q}\sqrt{q}$ ，对 单位四元数 来说，这等同于把 向量部分 乘以单位化的系数，使角度减半，再正规化

我们已经知道两条已 归一化 的向量
$\vec{p}_0, \ \vec{p}_1 \in \mathbb{R}^3$
那么可以定义
$\begin{align*} d &= \vec{p}_0 \cdot \vec{p}_1 = \cos{\theta} \\ c &= \vec{p}_0 \times \vec{p}_1 \\ \|c\| &= \| \vec{p}_0 \| \| \vec{p}_1 \| \sin{\theta} = \sin{\theta} \end{align*}$
- $\| \vec{p}_0 \| = \| \vec{p}_1 \| = 1$
其中 $\theta$ 是它们的夹角，目标是构造把 $\vec{p}_0$ 旋转到 $\vec{p}_1$ 的 最短弧四元数 而不显式求 $\theta$

我们知道 半角恒等式
$\begin{align*} \cos{ \frac{\theta}{2} } &= \sqrt{ \frac{1+\cos{\theta} }{2} } = \sqrt{\frac{1+d}{2} } \\ \sin{ \frac{\theta}{2} } &= \sqrt{ \frac{1 - \cos{\theta} }{2} } = \sqrt{\frac{1 - d}{2} } \end{align*}$
以及 标准轴–角 到 四元数 的公式是
$q = (\mathbf{a}\sin{\frac{\theta}{2} }, \ \cos{\frac{\theta}{2} })$
- $\mathbf{a}$ 是 单位旋转轴
标量部分
$\cos{ \frac{\theta}{2} } = \sqrt{\frac{1+d}{2} }$
向量部分
$\mathbf{a}\sin{\frac{\theta}{2} } = \frac{c}{\|c\|}\sqrt{\frac{1-d}{2} }$
所以
$q = ( \frac{c}{\|c\|}\sqrt{\frac{1-d}{2} }, \ \ \sqrt{\frac{1+d}{2} } )$
虽然这已经不需要显式计算 $\theta$ ，只需要计算点乘和叉乘，但是我们仍然有优化的空间

根据 正弦的倍角公式 $\sin{2\theta} = 2\sin{\theta}\cos{\theta}$ ，我们知道
$\|c\| = \sin{\theta} = 2\sin{\frac{\theta}{2} }\cos{\frac{\theta}{2} }$
再加上之前已经提过的 半角恒等式
$\sin{\frac{\theta}{2} } = \sqrt{\frac{1-d}{2} }, \quad \cos{\frac{\theta}{2} } = \sqrt{\frac{1+d}{2} }$
我们可以将 向量部分 改写为
$\begin{align*} \mathbf{v} = \frac{c}{\|c\|}\sqrt{\frac{1-d}{2} } = c \ \frac{ \sqrt{ \frac{1-d}{2} } }{ 2\sin{\frac{\theta}{2} }\cos{\frac{\theta}{2} } } = c \ \frac{ \sqrt{ \frac{1-d}{2} } }{ 2\sqrt{\frac{1-d}{2} } \sqrt{\frac{1+d}{2} } } &= \frac{ c }{ 2 \sqrt{\frac{1+d}{2} } } \\ &= \frac{ c }{ \sqrt{2(1+d)} } \end{align*}$
相应的，我们也将 标量部分 写为 $\sqrt{2(1+d)}$ 的形式
$w = \sqrt{\frac{1+d}{2} } = \frac{\sqrt{2(1+d)} }{2}$
所以我们从头到尾 只需要计算一次平方根 即可
```
Quaternion FromToRotationFast(vec3 from, vec3 to)
{
    vec3  p0 = normalize(from);
    vec3  p1 = normalize(to);
 
    const float EPS = 1e-6f;
    float d = dot(p0, p1);
 
    // (a) 同向
    if (d > 1.0f - EPS)         // θ ≈ 0
        return quat(0,0,0,1);
 
    // (b) 反向
    if (d < -1.0f + EPS)        // θ ≈ π，需要任选正交轴
    {
        vec3 ortho = (fabs(p0.x) < 0.6f) ? vec3(1,0,0)
                    : (fabs(p0.y) < 0.6f) ? vec3(0,1,0)
                                          : vec3(0,0,1);
        vec3 axis = normalize(cross(p0, ortho));
        return quat(axis.x, axis.y, axis.z, 0);    // 180°
    }
 
    // (c) 一般情况：只用一次 sqrt
    vec3  c   = cross(p0, p1);             // 3 mul, 3 sub
    float s   = sqrtf(2.f * (1.f + d));    // 唯一一次 sqrt
    float inv = 1.f / s;                   // 或用 rsqrt 快速倒数
 
    return quat(c.x * inv,                // = c / √(2(1+d))
                c.y * inv,
                c.z * inv,
                0.5f * s);                // = √((1+d)/2)
}
```
- ALU 统计（常规路径）
  - 1 × dot + 1 × cross
  - 1 × sqrt + 1 × rsqrt/recip ：若 rsqrt 有硬件指令，可把 sqrt+div 折成 rsqrt + 1 乘
- 访存
  只读 2 × vec3，写 1 × quat；无额外临时向量存储
这里，不需再做第二次 sqrt，因为 s 已经是那条根号；也不需再做一次除法——乘常量 0.5 会跟最终归一化常数折叠到同一条 FMA（乘-加）指令里，CPU/GPU 都很便宜，这就完成了 最短弧旋转四元数 的构造，且再无额外根号或除法

Look At

要把物体 局部坐标 的 前向(forward) 转到 目标方向(direction)，同时保证物体的新 up 方向 与期望的 up 向量 一致，我们需要

步骤	说明	关键运算
① 找旋转	求一个四元数，把默认前向 (0, 0, 1) 转到目标方向 `direction`	from-to rotation
② 正交化参考坐标	先算右向 `right = cross(direction, desiredUp)`，再算新的 up：`desiredUp = cross(right, direction)`，保证 up ⟂ direction	叉积
③ 计算物体原 up	用步骤①得到的四元数把参考 up (0, 1, 0) 旋转过去，得到 “物体眼中的 up”	向量×四元数
④ 再建四元数	求一个四元数，把步骤③得到的 up 转到步骤②的 `desiredUp`	from-to rotation
⑤ 组合旋转	把①④两个四元数相乘（先用①，再用④，所以乘法顺序反过来）得到最终旋转	四元数乘

四元数复合后仍然只是 一次旋转，但它等价于“先绕球面走到目标前向，再绕同一球面走到正确的 up”。因此组合后的结果角轴既唯一又没有万向锁

Quaternion LookAt(Vector3 direction, Vector3 desiredUp) {
    // Normalize input data
    direction = Normalize(direction);
    desiredUp = Normalize(desiredUp)
 
    // 共线判定
    const float EPS = 1e-5f;
    if (length(direction) < EPS) return Quaternion::Identity();
    if (abs(dot(direction, desiredUp)) > 1.f - EPS) {
        // 选一个与 direction 最不平行的世界轴替代 up
        desiredUp = abs(direction.y) < 0.99f ? Vector3(0,1,0) : Vector3(1,0,0);
    }
 
 
    // Step 1, Find quaternion that rotates from forward to direction
    Quaternion fromForwardToDirection = FromToRotation(Vector3(0, 0, 1), direction);
 
    // Step 2, Make sure up is perpendicular to desired direction
    Vector3 right = Cross(direction, desiredUp);
    desiredUp = Cross(right, direction); 
 
    // Step 3, Find the up vector of the suaternion from Step 1
    // Quaternion-vector multiplication (will be covered later)
    Vector3 objectUp = Mul(Vector3(0, 1, 0), fromForwardToDirection);
    
    // Step 4, Create quaternion from object up to desired up
    Quaternion fromObjectUpToDesiredUp = FromToRotation(objectUp, desiredUp);
 
    // Step 5, Combine rotations (in reverse! forward applied first, then up)
    // Quaternion-quaternion multiplication (will be covered later)
    Quaternion result = Mul(fromObjectUpToDesiredUp, fromForwardToDirection);
 
    // Should not be needed, but normalize output data
    return Q_Normalize(result);
}

从四元数提取旋转轴与旋转角

轴(axis) = 取四元数的向量部分 (x,y,z) 再归一化

Vector3 GetAxis(Quaternion q) { return normalize(Vector3(q.x, q.y, q.z)); }

角(angle) = 利用标量部分 w = cos(θ/2)

float GetAngleRadians(Quaternion q) { return 2 * acos(q.w); }
 
float GetAngleDegrees(Quaternion quat) { return 2.0 * acos(quat.w) * 57.2958;}

Creating Quaternions 创建四元数

从一个方向转向另一个方向的最短弧 shortest arc of from-to

Look At

从四元数提取 旋转轴 与 旋转角

从四元数提取旋转轴与旋转角