4.3 四元数
尽管四元数早在1843年就由William Rowan Hamilton爵士发明,作为复数的扩展,但直到1985年Shoemake[1633]才将它们引入计算机图形领域1。四元数用于表示旋转和方位。它们在几个方面都优于欧拉角和矩阵。任何三维方向都可以表示为围绕特定轴的单次旋转。给定轴和角度表示,与四元数转换相互转换很简单,然后任一方向的欧拉角转换则具有挑战性。四元数可用于稳定和恒定的方向插值,这是欧拉角无法很好完成的。
复数有实部和虚部。每个复数由两个实数表示,第二个实数乘以√−1。同样,四元数有四个部分。前三个值与旋转轴密切相关,旋转角度影响所有四个部分(更多可参考第4.3.2节)。每个四元数由四个实数表示,每个实数与不同的部分相关联。由于四元数有四个分量,我们选择将它们表示为向量,但为了区分它们,我们给它们戴上帽子:^qq。我们从四元数的一些数学背景开始,然后用它来构造各种有用的转换。
4.3.1 数学背景
我们从四元数的定义开始。
定义:四元数^qq可以通过以下方式定义,都是等价的。
^qq=(qqv,qw)=iqx+jqy+kqz+qw=qqv+qw,qqv=iqx+jqy+kqz=(qx,qy,qz),i2=j2=k2=−1,jk=−kj=i,ki=−ik=j,ij=−ji=k.(4.31)
变量qw称为四元数^qq的实部。虚部是qqv,而i、j和k称为虚数单位。
对于虚部qqv,我们可以使用所有法向量运算,例如加法、缩放、点积、叉积等。使用四元数的定义,推导出两个四元数^qq和^rr之间的乘法运算,如下所示。请注意,虚数单位的乘法是不可交换的。
乘法:乘法:^qq^rr=(iqx+jqy+kqz+qw)(irx+jry+krz+rw)=i(qyrz−qzry+rwqx+qwrx)+j(qzrx−qxrz+rwqy+qwry)+k(qxry−qyrx+rwqz+qwrz)+qwrw−qxrx−qyry−qzrz=(qqv×rrv+rwqqv+qwrrv,qwrw−qqv⋅rrv)(4.32)
从这个方程可以看出,我们同时使用叉积和点积来计算两个四元数的乘法。
除了四元数的定义,还需要加法、共轭、范数和单位元素的定义:
加法:加法:^qq+^rr=(qqv,qw)+(rrv,rw)=(qqv+rrv,qw+rw)共轭:共轭:^qq∗=(qqv,qw)∗=(−qqv,qw)范数:范数:n(^qq)=√^qq^qq∗=√^qq∗^qq=√qqv⋅qqv+q2w=√q2x+q2y+q2z+q2w单位元素:单位元素:^ii=(00,1)(4.33)
当n(^qq)=√^qq^qq∗被简化时(结果如上所示),虚部抵消,只剩下实部。范数有时表示为||^qq|=n(^qq)[1105]。上面的结果是可以导出一个乘法逆,用^qq−1表示。方程^qq−1^qq=^qq−1^qq=1对逆必须成立(这对于乘法逆是常见的)。我们从范数的定义推导出一个公式:
n(^qq)2=^qq^qq∗⟺^qq^qq∗n(^qq)2=1(4.34)
这给出了乘法逆,如下所示:
逆:逆:^qq−1=1n(^qq)2^qq∗
逆的公式使用了标量乘法,可以从方程4.3.1中看到的乘法推导出这个运算:s^qq=(00,s)(qqv,qw)=(sqqv,sqw), 并且^qqs=(qqv,qw)(00,s)=(sqqv,sqw)。这意味着标量乘法是可交换的:s^qq=^qqs=(sqqv,sqw)。
以下规则集合很容易从定义中推导出来:
共轭规则:共轭规则:(^qq∗)∗=^qq(^qq+^rr)∗=^qq∗+^rr∗(^qq^rr)∗=^rr∗^qq∗(4.36)
范数规则:范数规则:n(^qq∗)=n(^qq)n(^qq^rr)=n(^qq)n(^rr)(4.37)
乘法法则:乘法法则:线性度:线性度:^pp(s^qq+t^rr)=s^pp^qq+t^pp^rr(s^pp+t^qq)^rr=s^pp^rr+t^qq^rr结合性:结合性:^pp(^qq^rr)=(^pp^qq)^rr(4.38)
单位四元数^qq=(qqv,qw)使得n(^qq)=1。由此可知^qq可写作:
^qq=(sinϕuuq,cosϕ)=sinϕuuq+cosϕ(4.39)
对于某个三维向量uuq,使得||uuq||=1,因为
n(^qq)=n(sinϕuuq,cosϕ)=√sin2ϕ(uuq⋅uuq)+cos2ϕ=√sin2ϕ+cos2ϕ=1(4.40)
当且仅当uuq⋅uuq=1=||uuq||2. 正如将在下一节中看到的,单位四元数非常适合以最有效的方式创建旋转和方向。但在此之前,会为单位四元数引入一些额外的操作。
对于复数,二维单位向量可以写成cosϕ+isinϕ=eiϕ。 四元数的等价物是
^qq=sinϕuuq+cosϕ=eϕuuq(4.41)
根据公式4.41,单位四元数的对数函数和幂函数为:
对数:对数:log(^qq)=log(eϕuuq)=ϕuuq指数:指数:^qqt=(sinϕuuq+cosϕ)t=eϕtuuq=sin(ϕt)uuq+cos(ϕt)(4.42)
4.3.2 四元数变换
我们现在将研究四元数集的一个子类,即单位长度的子类,称为单位四元数。关于单位四元数的最重要的事实是它们可以表示任何三维旋转,而且这种表示非常紧凑和简单。
现在我们将描述是什么使单位四元数对旋转和方向如此有用。首先,将点或向量的四个坐标pp=(px py pz pw)T代入四元数^pp的分量,并假设我们有一个单位四元数^qq=(sinϕuuq,cosϕ)。可以证明
^qq^pp^qq−1(4.43)
将^pp(以及点pp)绕轴uuq旋转角度2ϕ。注意,由于^qq是一个单位四元数,^qq−1=^qq∗。参见图4.9。
图4.9. 由单位四元数表示的旋转变换的图示,^qq=(sinϕuuq,cosϕ)。变换围绕轴uuq旋转 2ϕ弧度。
^qq的任何非零实数倍数也表示相同的变换,这意味着^qq和−^qq表示相同的旋转。也就是说,取反轴uuq和实部qw会创建一个与原始四元数完全一样旋转的四元数。这也意味着从矩阵中提取四元数可以返回^qq或−^qq。
给定两个单位四元数^qq和^rr,首先应用^qq再应用^rr到四元数^pp(可以解释为点pp)的级联由公式 4.44 给出:
^rr(^qq^pp^qq∗)^rr∗=(^rr^qq)^pp(^rr^qq)∗=^cc^pp^cc∗(4.44)
这里,^cc=^rr^qq是单位四元数,表示单位四元数^qq和^rr的级联。
矩阵转换
由于经常需要将几种不同的变换组合起来,而且大部分都是矩阵形式,因此需要一种方法将式4.43转化为矩阵。四元数 ^qq可以转换为矩阵MMq,如公式4.45[1633,1634]所示:
MMq=⎛⎜
⎜
⎜
⎜⎝1−s(q2y+q2z)s(qxqy−qwqz)s(qxqz+qwqy)0s(qxqy+qwqz)1−s(q2x+q2z)s(qyqz−qwqx)0s(qxqz−qwqy)s(qyqz+qwqx)1−s(q2x+q2y)00001⎞⎟
⎟
⎟
⎟⎠(4.45)
在这里,标量是s=2/(n(^qq))2。 对于单位四元数,上式可简化为:
MMq=⎛⎜
⎜
⎜
⎜⎝1−2(q2y+q2z)2(qxqy−qwqz)2(qxqz+qwqy)02(qxqy+qwqz)1−2(q2x+q2z)2(qyqz−qwqx)02(qxqz−qwqy)2(qyqz+qwqx)1−2(q2x+q2y)00001⎞⎟
⎟
⎟
⎟⎠(4.46)
一旦构造了四元数,就不需要计算三角函数,因此转换过程在实践中很有效率。
从正交矩阵MMq到单位四元数^qq的逆转换要复杂一些。此过程的关键是与公式4.46中的矩阵的以下差异:
mq21−mq12=4qwqxmq02−mq20=4qwqymq10−mq01=4qwqz(4.47)
这些方程的含义是,如果qw已知,则可以计算向量vvq的值,从而推导出^qq。MMq的迹由下式计算:
tr(MMq)=4−2s(q2x+q2y+q2z)=4(1−q2x+q2y+q2zq2x+q2y+q2z+q2w)=4q2wq2x+q2y+q2z+q2w=4q2w(n(^qq))2(4.48)
对于单位四元数来说,可得出以下转换:
qw=12√tr(MMq)qx=(mq21−mq12)4qwqy=(mq02−mq20)4qwqz=(mq10−mq01)4qw(4.49)
为了有一个数值稳定的程序[1634],小数除法应该被避免。因此,首先设置t=q2w−q2x−q2y−q2z,由此得出:
m00=t+2q2xm11=t+2q2ym22=t+2q2zu=m00+m11+m22=t+2q2w(4.50)
这反过来意味着m00、m11、m22和u中的最大值,确定了qx、qy、qz和qw中最大的。如果qw最大,则使用公式4.49导出四元数。否则,我们注意到以下内容成立:
4q2x=+m00−m11−m22+m334q2y=−m00+m11−m22+m334q2z=−m00−m11+m22+m334q2w=tr(MMq)(4.51)
然后使用上述适当方程计算qx、qy和qz中的最大值,然后使用方程4.47计算^qq的其余分量。Schüler[1588]提出了一种无分支但使用四个平方根的变体。
球面线性插值
球面线性插值是一种运算,在给定两个单位四元数^qq和^rr以及参数t∈[0,1]的情况下,计算插值四元数。例如,这对于动画对象很有用。它对于插值相机方向没有用,因为相机的“向上”矢量在插值过程中可能会倾斜,通常是一种干扰效果。
此运算的代数形式由复合四元数^ss表示,如下所示:
^ss(^qq,^rr,t)=(^rr^qq−1)t^qq(4.52)
但是,对于软件实现,以下形式更合适,其中slerp代表球面线性插值:
^ss(^qq,^rr,t)=slerp(^qq,^rr,t)=sin(ϕ(1−t))sinϕ^qq+sin(ϕt)sinϕ^rr(4.53)
为了计算这个方程所需的φ,可以使用以下式子:cosϕ=qxrx+qyry+qzrz+qwrw[325]。对于 t∈[0,1],slerp函数计算(唯一的2)插值四元数,它们共同构成四维单位球面上从^qq(t=0)到^rr(t=1)的最短弧。圆弧位于由^qq,^rr和原点给出的平面与四维单位球面的交点形成的圆上。如图 4.10 所示。计算出的旋转四元数以恒定速度绕固定轴旋转。像这样的曲线具有恒定的速度,因此加速度为零,称为测地曲线[229]。过原点的平面与球面的交点在球面上形成一个大圆,这个圆的一部分称为大圆弧。
图4.10. 单位四元数表示为单位球面上的点。函数slerp用于四元数之间的插值,插值的路径是球体上的一个大圆弧。请注意,从^qq1到^qq2和从^qq1到^qq3到^qq2的插值不是一回事,即使它们到达相同的方向。
slerp函数非常适合在两个方位之间进行插值,并且有良好的效果(固定轴,恒速)。使用多个欧拉角进行插值时,情况并非如此。实际上,直接计算slerp是一项涉及调用三角函数的昂贵操作。Malysau[1114]讨论了将四元数集成到渲染管线中。他指出,当不使用slerp而只是在像素着色器中对四元数进行归一化时,对于90度角,三角形方向的误差最大为 4度。光栅化三角形时,此误差率是可以接受的。Li[1039, 1040]提供了更快的增量方法来计算slerps,并且不会牺牲任何准确性。Eberly[406]提出了一种仅使用加法和乘法计算slerps的快速技术。
当有两个以上的方向时,比如^qq0,^qq1,...,^qqn,并且我们想要从^qq0到^qq1到^qq2,依此类推直到^qqn−1,可以直接使用slerp。现在,当我们接近^qqi时,我们将使用^qqi−1和^qqi作为 slerp 的参数。通过^qqi后,我们将使用^qqi和^qqi+1作为 slerp 的参数。这将导致方向插值中出现突然的抖动,如图4.10所示。这类似于当点被线性插值时发生的情况;参见第720页图17.3的右上部分。有些读者可能希望在阅读第17章中的样条曲线后重新阅读以下段落。
更好的插值方法是使用某种样条。我们在^qqi和^qqi+1之间引入了四元数^aai和^aai+1。球面三次插值可以在四元数^qqi,^aai,^aai+1和^qqi+1的集合内定义。令人惊讶的是,这些额外的四元数计算如下[404]3:
^aai=^qqiexp[−log(^qq−1i^qqi−1)+log(^qq−1i^qqi+1)4](4.54)
通过三次样条平滑算法,^qqi和^aai用来对四元数进行球面插值,如公式4.55所示:
squad(^qqi,^qqi+1,^aai,^aai+1,t)=slerp(slerp(^qqi,^qqi+1,t),slerp(^aai,^aai+1,t),2t(1−t))(4.55)
从上面可以看出,squad函数是使用slerp从同样的球面插值构建的(第17.1.1节有关点的重复线性插值的信息)。插值将通过初始方向^qqi,i∈[0,...,n−1],但不通过^aai——这些用于指示初始方向的切线方向。
从一个向量到另一个向量的旋转
一个常见的操作是通过尽可能最短的路径从一个方向ss转换到另一个方向tt。四元数的数学运算极大地简化了这个过程,并显示了四元数与这种表示的密切关系。首先,将ss和tt归一化。然后计算单位旋转轴,称为uu,计算为uu=(ss×tt)/||ss×tt||。 接下来,e=ss⋅tt=cos(2ϕ)和||ss×tt||=sin(2ϕ),其中2ϕ是ss和tt之间的角度。表示从ss到tt的旋转的四元数是^qq=(sinϕuu,cosϕ)。事实上,使用半角关系和三角恒等式简化^qq=(sinϕsin2ϕ(ss×tt),cosϕ),可得出[1197]:
^qq=(qqv,qw)=(1√2(1+e)(ss×tt),√2(1+e)2)(4.56)
以这种方式直接生成四元数(相对于标准化叉积ss×tt)避免了当ss和tt指向几乎相同的方向时的数值不稳定[1197]。当ss和tt指向相反的方向时,两种方法都会出现稳定性问题,因为发生了除以零的情况。当检测到这种特殊情况时,任何垂直于ss的旋转轴都可以用来旋转到tt。
有时我们需要从ss到tt旋转的矩阵表示。在对等式4.46进行一些代数和三角简化后,旋转矩阵变为[1233]:
RR(ss,tt)=⎛⎜
⎜
⎜
⎜⎝e+hv2xhvxvy−vzhvxvz+vy0hvxvy+vze+hv2yhvyvz−vx0hvxvz−vyhvyvz+vxe+hv2z00001⎞⎟
⎟
⎟
⎟⎠(4.57)
在这个等式中,我们使用了以下中间计算:
vv=ss×tte=cos(2ϕ)=ss⋅tth=1−cos(2ϕ)sin2(2ϕ)=1−evv⋅vv=11+e(4.58)
可以看出,由于简化,所有平方根和三角函数都消失了,因此这是创建矩阵的有效方法。请注意,公式4.57的结构类似于公式4.30的结构,并注意后一种形式如何不需要三角函数。
请注意,当ss和tt平行或接近平行时必须小心,因为这时||ss×tt||≈0。如果ϕ≈0,那么我们可以返回单位矩阵。然而,如果2ϕ≈π,那么我们可以绕任意轴旋转π弧度。该轴可以是ss和任何其他不平行于ss的向量的叉积(第4.2.4节)。Möller和Hughes使用Householder矩阵以不同的方式处理这种特殊情况[1233]。
1公平地说,Robinson[1502]在1958年使用四元数进行刚体模拟。
2当且仅当^qq和^rr不相反。
3Shoemake[1633]给出了另一个推导。
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 一个奇形怪状的面试题:Bean中的CHM要不要加volatile?
· [.NET]调用本地 Deepseek 模型
· 一个费力不讨好的项目,让我损失了近一半的绩效!
· .NET Core 托管堆内存泄露/CPU异常的常见思路
· PostgreSQL 和 SQL Server 在统计信息维护中的关键差异
· CSnakes vs Python.NET:高效嵌入与灵活互通的跨语言方案对比
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· Plotly.NET 一个为 .NET 打造的强大开源交互式图表库
· 上周热点回顾(2.17-2.23)