关于四元数的数学概念:四元数原理学习笔记
四元数的旋转
四元数并不只能用来表示旋转,任何旋转量都是为了计算空间坐标服务的,无论是四元数,还是矩阵,欧拉角,你都可以找用他们来表示坐标的形式。四元数表示坐标的形式很自然,就是虚数部分对应三维坐标的三个维度,实部为零。即,表示坐标的四元数必须是一个纯四元数。
p
=
x
i
+
y
j
+
z
k
\textbf{p} = x\textbf{i} + y\textbf{j} +z\textbf{k}
p=xi+yj+zk
则绕Y轴旋转
θ
\theta
θ角的四元数记为(其他轴同理):
q
=
c
o
s
1
2
θ
+
s
i
n
1
2
θ
⋅
j
\textbf{q}=cos\frac{1}{2}\theta + sin\frac{1}{2}\theta \cdot \textbf{j}
q=cos21θ+sin21θ⋅j
记旋转之后的新顶点为
p
′
\textbf{p}'
p′,四元数的旋转运算规则为:
p
′
=
qp
q
−
1
\textbf{p}' = \textbf{q} \textbf{p} \textbf{q}^{-1}
p′=qpq−1
链式运算为:
p
′
=
q
n
q
n
−
1
.
.
.
q
1
p
q
1
−
1
.
.
.
q
n
−
1
−
1
q
n
−
1
\textbf{p}' = \textbf{q}_n\textbf{q}_{n-1}...\textbf{q}_1\ \ \textbf{p}\ \ \textbf{q}_1^{-1}...\textbf{q}_{n-1}^{-1}\textbf{q}_{n}^{-1}
p′=qnqn−1...q1 p q1−1...qn−1−1qn−1
举例:比如点
(
0
,
1
,
0
)
(0, 1, 0)
(0,1,0),将其绕Z轴旋转60°,我们可以很自然按照初中的几何知识计算得出最终坐标为
(
±
3
2
,
1
2
,
0
)
(\pm\frac{\sqrt{3}}{2},\frac{1}{2}, 0)
(±23,21,0)(正负号取决于你如何定义三个轴的方向和顺序),用四元数验证一下:此时坐标四元数为
p
=
j
\textbf{p} = j
p=j,旋转四元数为
q
=
c
o
s
30
°
+
s
i
n
30
°
⋅
k
\textbf{q} = cos30° + sin30°\cdot\textbf{k}
q=cos30°+sin30°⋅k:
p
′
=
qp
q
−
1
=
qp
q
∗
∣
q
∣
2
=
(
c
o
s
30
°
+
s
i
n
30
°
⋅
k
)
⋅
j
⋅
(
c
o
s
30
°
−
s
i
n
30
°
⋅
k
)
1
=
1
2
j
−
3
2
i
=
(
−
3
2
,
1
2
,
0
)
\begin{aligned} \textbf{p}' &= \textbf{q} \textbf{p} \textbf{q}^{-1}\\ &= \frac{\textbf{q} \textbf{p} \textbf{q}^{*}}{|\textbf{q}|^{2}}\\ &= \frac{(cos30° + sin30°\cdot\textbf{k})\cdot \textbf{j} \cdot (cos30° - sin30°\cdot\textbf{k})}{1}\\ &= \frac{1}{2}\textbf{j} - \frac{\sqrt{3}}{2}\textbf{i}\\ &= (-\frac{\sqrt{3}}{2},\frac{1}{2},0) \end{aligned}
p′=qpq−1=∣q∣2qpq∗=1(cos30°+sin30°⋅k)⋅j⋅(cos30°−sin30°⋅k)=21j−23i=(−23,21,0)
Unity中的四元数和旋转
从上面可以知道,对任何一个点,我们用 p \textbf{p} p记录点的初始坐标,用 q \textbf{q} q记录点的旋转状态,然后用公式 p ′ = qp q − 1 \textbf{p}' = \textbf{q} \textbf{p} \textbf{q}^{-1} p′=qpq−1 计算点经过旋转之后的坐标。注意,如果有新的旋转 q 1 \textbf{q}_1 q1,根据运算规则,它是加在原旋转状态的左边: q 1 q \textbf{q}_1\textbf{q} q1q。
对于Unity中的模型上的每一个点都共享同一个旋转状态,点的初始坐标为制作模型的时候我们定义的点的坐标(这取决于物体的pivot点),这个坐标都被记录在Mesh.vertices
里。对于模型的任意一个点
i
i
i,其在世界坐标系下的顶点坐标为
o
b
j
P
=
t
r
a
n
s
f
o
r
m
.
p
o
s
i
t
i
o
n
v
P
i
=
m
e
s
h
.
v
e
r
t
i
c
e
s
[
i
]
q
=
t
r
a
n
s
f
o
r
m
.
r
o
t
a
t
i
o
n
v
G
P
i
=
o
b
j
P
+
q
⋅
v
P
i
⋅
q
−
1
\begin{aligned} objP &= transform.position \\ vP_i &= mesh.vertices[i] \\ q &= transform.rotation \\ vGP_i &= objP + q \cdot vP_i \cdot q^{-1} \end{aligned}
objPvPiqvGPi=transform.position=mesh.vertices[i]=transform.rotation=objP+q⋅vPi⋅q−1
所以transform.rotation
记录的是从世界坐标系下的旋转状态。或者说,模型坐标系相对于世界坐标系的旋转状态。
- 对于任何一个新的旋转
q
,q * transform.rotation
表示在当前旋转状态之上叠加旋转,是相对于世界坐标系的。 - 如果是
transform.rotation * q
呢,它有点违反四元数旋转运算的规则,因为新的旋转四元数应该加在左边,但是它也是有意义的,表示此时是相对模型坐标系进行旋转。这其实也很好理解:我们知道运算transform.rotation * p
是先计算旋转p
,现在我们把物体A经过旋转p
之后的物体当成另一个物体B,之后他们都经历了世界坐标下的旋转transform.rotation
,请问,B和A的差别差在哪里?很明显就是B就是围绕自身模型坐标旋转之后的结果。打个比方,初始状态下,A是一个人,现在又有一个人,记为B,但是头朝下,显然B在A的基础上先经历了旋转p
。现在把A和B当成一个整体,经历旋转transform.rotation
,那B相对于A还是头朝下么?当然是了!
那transform.localrotation
记录的是什么呢?记录的是相对于父元素的旋转状态。与其这么解释,不如说rotation
和localrotation
都是记录相对于全局坐标系的旋转状态,只是一个认为世界空间是全局,一个认为父物体空间是全局。
欧拉角和四元数
老生常谈的话题。为什么用四元数不用欧拉角教科书已经写得很清楚了,不外乎下面三点:
- 欧拉角表示旋转不唯一
- 欧拉角旋转插值非线性
- 欧拉角的万向锁问题
旋转不唯一是说有对于任意一个旋转,有无数中欧拉角可以表示这个旋转。根源其实在于三角函数的周期性。插值非线性我没有仔细思考过,但是原因估计跟三角函数的非线性有关,即对于角度本身是线性插值,但是映射到三角函数上就不是了。
重点说一下万向锁的问题,在此之前,先说明一个重要的结论,即欧拉角的旋转是有顺序的。这是说,先绕着x轴旋转 θ \theta θ角,再绕着y轴旋转 ϕ \phi ϕ角和先转y轴再转x轴得到的结果不同。这很好验证,拿自己的手机做就好了。
之后在Unity中做一个实验:新建一个Cube,reset Transform组件。然后将Rotation的x置为90°。现在你任意改变Rotation的y和z,你会发现,效果竟然是一样的!!!他们都绕着同一个轴旋转,这就是万向锁现象,但这是为何呢?
具体指明一下结果,就是,修改Rotation.y的时候,物体绕着世界坐标系的y轴旋转;而修改Rotation.z的时候,物体绕着模型坐标系的z轴旋转。这是为什么呢?
基于以下几个原因:
- Unity的旋转是有顺序的(为什么有顺序已经说了,是因为同一个欧拉角,不同的旋转顺序可以有不同的结果),顺序为 z → x →y
- 旋转都是绕着世界坐标系旋转
- Unity的旋转是一个状态不是一个过程
先理解第三句话,很重要:Unity的旋转是一个状态不是一个过程。这句话是说,即使你先指定的x,再指定的z,但是对于Unity来说,每当它要渲染一帧的时候,它拿到的都是三个数字。它不关心你先指定的哪个分量。一旦它拿到三个旋转参数,它会按照 z → x →y 的顺序旋转并计算模型每个顶点的最终位置。下面具体说明一下问题:
- 第一种情况就是你先在Inspector上先改变z,再改变x,物体都是绕着世界坐标轴旋转,没有毛病
- 第二种情况,你先改变x为90°(不是90°也没有问题,只是为了方便讲解和你的观察),然后再改变z,你会发现你在改变z的时候,物体是绕着自身的z轴旋转的。这是因为,其实你每一次改变z轴的瞬间,新的旋转参数发送给Unity,Unity都只拿到新的x分量和z分量,然后先绕着世界坐标的z轴旋转z分量,再绕着世界坐标的x轴旋转x分量。只不过这个顺序是计算的顺序,不是渲染的顺序,最终渲染出来的结果以及呈现在屏幕上的动画让你误以为物体是在绕着自身z轴旋转90°
换句话说,Unity的计算每一帧的顺序是这样的:
p
′
=
q
y
q
x
q
z
p
q
z
−
1
q
x
−
1
q
y
−
1
\textbf{p}' = \textbf{q}_y\textbf{q}_x\textbf{q}_z\ \ \textbf{p}\ \ \textbf{q}_z^{-1}\textbf{q}_x^{-1}\textbf{q}_y^{-1}
p′=qyqxqz p qz−1qx−1qy−1
其中无论
p
\textbf{p}
p是点在世界坐标下的初始位置,即旋转状态为0的位置,其次内部的运算有可能是用矩阵而不是四元数,但是其实四元数和矩阵的本质是一样的,可以相互转换,无碍。因为我们的例子没有旋转y轴,我们先把y去掉:
p
′
=
q
x
q
z
p
q
z
−
1
q
x
−
1
\textbf{p}' = \textbf{q}_x\textbf{q}_z\ \ \textbf{p}\ \ \textbf{q}_z^{-1}\textbf{q}_x^{-1}
p′=qxqz p qz−1qx−1
因为你上一帧已经改变了x分量和z分量,这一帧你想单独改变z分量,你的设想是将新的旋转在原有的旋转上面继续加,like this:
p
′
=
q
z
2
q
x
q
z
p
q
z
−
1
q
x
−
1
q
z
2
−
1
\textbf{p}' = \textbf{q}_{z2}\textbf{q}_x\textbf{q}_z\ \ \textbf{p}\ \ \textbf{q}_z^{-1}\textbf{q}_x^{-1}\textbf{q}_{z2}^{-1}
p′=qz2qxqz p qz−1qx−1qz2−1
但其实不是,Unity每一帧都从头开始,给你计算,like this:
p
′
=
q
x
q
z
′
p
q
z
′
−
1
q
x
−
1
\textbf{p}' = \textbf{q}_x\textbf{q}_{z}'\ \ \textbf{p}\ \ \textbf{q}_{z}'^{-1}\textbf{q}_x^{-1}
p′=qxqz′ p qz′−1qx−1
现在我们来解释一下,为什么呈现在你眼中的动画是绕着物体自身的z轴旋转,它的原理其实和上面解释transform.rotation * q
和q * transform.rotation
的时候几乎一样。
现在你改变z分量的时候,不同的帧之间,x分量是相同的,只有z分量不同。但是Unity计算的时候,是先计算z旋转对么?现在我们把不同帧的物体看成是不同的物体,那么在经历x轴旋转之前,他们的差别在哪里?答案是经历了不同的z旋转,而且是绕着世界坐标系的z轴旋转的。但是此时世界坐标系和模型坐标系重合,所以其实也可以说是绕着模型坐标系旋转的。现在这些物体共同经历了同样的x旋转,那旋转之后他们的差别在哪儿呢?答案很明显了,就是绕着自身z轴旋转了不同的角度。
至此,欧拉角的万向锁问题就讲完了。它出现的根源在于欧拉角是一个描述旋转状态的变量,同时欧拉角的旋转又关乎旋转顺序。
很多人,包括之前的我,一直没有理解万向锁不是以为不知道旋转是有顺序的,而是一直没有明白第一句话。网上关于解释万向锁有一个很著名的youtube视频,很多博客都有引用:Gimbal lock。看完这个视频我还是不明白,他演示了飞机旋转到一定角度之后,无论动最内层还是最外层,都无法使飞机仰起来,可是在Unity中,无论这个飞机朝向何处,我想让他仰起来不是轻而易举么?
这里的思维误区就在于,视频所说并不是你不能将飞机仰起来,而是当旋转顺序中间的分量为90°的时候,任意改变其他两个分量,这些欧拉角所描述的旋转状态都是在围绕一个轴旋转,丢失了一个维度,你想想你明明有两个维度的分量可以改变,但是飞机的旋转只在一个旋转维度上发生变更。这说明欧拉角到旋转状态并不是覆盖映射的(See also : Covering Space,覆盖空间)。
是欧拉角在描述旋转状态时出了问题,因为欧拉角的特性,它导致了万向锁,在出现万向锁的时候,尽管你有两个维度的变量,但是只能描述一个维度的旋转状态。这使得欧拉角不适合描述旋转状态,航天工程尤为如此。但是因为易于人理解,所以我们可以用欧拉角去施加旋转动作。