Skip to content

4.4.3 四元数的应用 ​

4.4.3.1 计算机绘图学中的 3D 旋转 ​

为了刻画运动流,需要利用旋转的插值. 因为 3D 旋转可以用单位四元数表示, 所以在计算机绘图学中发展了对于旋转的插值算法. 最直接的想法是从与欧氏空间的线性插值类似的定义出发. 基本算法是 Lerp, Slerp 和 Squad.

1. Lerp(线性 (l) 插值 (erp)) ​

设 p,q∈H1 及 t∈[0,1] ,那么

(4.143)Lerp(p,q,t)=p(1−t)+qt.
  • 这是 R4 中的一条连接 p∈H1∼S3⊂R4 和 q∈H1∼S3⊂R4 的线段.

  • 这条线段是 R4 中单位球的内部,并且不表示单位球 H1∼S3 上的任何连通曲线.

  • 因此旋转由所求得的四元数的规范化确定.

这个简单的算法几乎是完美的. 仅有的问题是在通过给定点和规范化了的所求得的四元数之间的割线上求得插值点后, 所得到的单位四元数不是等距四元数. 这个问题可以用下列算法解决.

2. Slerp(球面 (s) 线性 (l) 插值 (erp)) ​

设 p,q∈H1,t∈[0,1],φ(0<φ<π) 是一个介于 p,q 之间的角. 那么

(4.144)Slerp(p,q,t)=p(p¯q)t=p1−tqt=p[sin⁡((1−t)φ)sin⁡φ]+q[sin⁡(tφ)sin⁡φ].
  • 沿着单位球 S3⊂R4 上的大圆插值,连接 p 和 q .

  • 选取最短连接, −Sc(p,q)=⟨p,q⟩>0 必须成立 (此处 ⟨⋅,⋅⟩ 表示 p 和 q∈R4 的点积).

在图 4.5 中比较了依据 Lerp(见图 4.5(a)) 和 Slerp(见图 4.5(b)) 的插值.

019363af-d8ae-7006-ac42-15a9aafbc2ce_43_603_1161_435_246_0.jpg

特殊情形 p=1 设 p=1=(1,0,0,0)T 及 q=cos⁡φ+nq―sin⁡φ ,那么

(4.145)Slerp(p,q,t)=cos⁡(tφ)+nq――sin⁡(tφ).

特殊情形等距网格 令 ψ=φn ,则

qk:=Slerp(p,q,kn)=1sin⁡φ(sin⁡(φ−kψ)p+sin⁡(kψ)q),k=0,1,⋯,n.

(4.146)

Slerp插值的解释 为证明 (4.144) 中两个表达式的等价性,首先算出 Q= p−1q=p¯|p|2q=p¯q . 因为 p,q∈H1 ,所以标量部分是

(4.147)Q0=ScQ=Sc(p¯q)=⟨p,q⟩=cos⁡φ.

因为 p=p⋅1 ,以及 q=pp−1q=pQ ,所以 1 和 Q 间的插值是乘以 p 以保持 p 和 q 间的插值.

Q(t)=sin⁡((1−t)φ)sin⁡φ+Qsin⁡(tφ)sin⁡φ=sin⁡((1−t)φ)sin⁡φ+cos⁡φsin⁡(tφ)sin⁡φ+n→Q→sin⁡(tφ)sin⁡φ=sin⁡φcos⁡(tφ)−sin⁡(tφ)cos⁡φ+sin⁡(tφ)cos⁡φsin⁡φ+n→Q→sin⁡(tφ)sin⁡φ(4.148)=cos⁡(tφ)+n→Q→sin⁡(tφ)=etn→Q→φ=etlog⁡Q=Qt.

由此得到

(4.149)q(t)=pQ(t)=psin⁡((1−t)φ)sin⁡φ+qsin⁡(tφ)sin⁡φ=pQt=p(p−1q)t=p1−tqt.

3. Squad(球面 (s) 和四角形 (quad) 插值) ​

对于 qi,qi+1∈H1 及 t∈[0,1] ,法则是

Squad(qi,qi+1,si,si+1,t)(4.150)=Slerp(Slerp(qi,qi+1,t),Slerp(si,si+1,t),2t(1−t)),

其中

si=qiexp⁡(−log⁡(qi−1qi+1)+log⁡(qi−1qi−1)4).
  • 所得曲线类似于贝济埃 (Bézier) 曲线, 代替线性插值, 它保持了球面插值.

  • 算法产生了对于四元数序列 q0,q1,⋯,qN 的插值曲线.

  • 在第一个和最后一个区间表达式未被定义,因为 q−1 对于计算 s0 ,以及 qN+1 对于计算 sN 是必须的. 一个可能的方法是选取 s0=q0 以及 sN=qN (或者定义 q−1 和 qN+1 ). 还有其他的基于四元数的算法: nlerp, log-lerp, islerp,四元数德卡斯特里奥 (de Casteeljau) 样条.

4.4.3.2 旋转矩阵的插值 ​

可以借助旋转矩阵完全类似地刻画 Slerp 算法. 3×3 旋转矩阵 R (即群 SO(3) 的元素) 的对数是需要的,并且用群论语言它被定义为满足 er=R 的斜对称矩阵 r (即李群 so(3) 的元素). 于是 Slerp 算法可以用于旋转矩阵 R0 和 R1 间的插值, 这可以刻画为

(4.151)R(t)=R0(R0−1R1)t=R0exp⁡(tlog⁡(R0−1R1)).

一般地说, 应用四元数基本算法并且依据单位四元数转换为旋转矩阵的计算, 由四元数 q(t) 确定 R(t) 要简单些.

4.4.3.3 球极平面投影 ​

如果将 1∈H1∼S3 取作三维球 S3 的北极,那么单位四元数或三维球的元素可以用球极平面投影

H1∋q↦(1+q)(1−q)−1∈H0∼R3

分别映为纯四元数或 R3 的元素. 对应的逆映射是

(4.152)R3∼H0∋p↦(p−1)(p+1)−1∈H1∼S3.

4.4.3.4 卫星导航 ​

环绕地球运行的人造卫星的航向是已经确定的. 恒星看作在无穷远处, 于是它们相对于地球和相对于卫星的方向是相同的 (图 4.6). 任何测量上的差别都可从不同的坐标系推导出来, 因而也可由不同坐标系的相对旋转推出.

019363af-d8ae-7006-ac42-15a9aafbc2ce_45_619_1016_412_176_0.jpg

设 a→i 是在地球上的固定坐标系中指向第 i 颗恒星方向的单位向量, b→i 是在卫星上的固定坐标系中指向第 i 颗恒星方向的单位向量. 两个坐标系的相对旋转可以用单位四元数 h∈H1 刻画:

(4.153)b→i=ha→ih¯

如果考虑另一个恒星, 并且数据被测量误差覆盖, 那么解可以用最小二乘法, 即作为 (4.154) 的极小值确定,其中 h 是单位四元数, a―i=a→i 和 b―i=b→i 是单位向量:

Q2=∑i=1n⁡|b→i−ha→ih¯|2=∑i=1n⁡(b→i−ha→ih¯)⋅(b→i−ha→ih¯)(4.154)=∑i=1n⁡(b―i−ha―ih¯)(b―i−ha―ih¯)―=∑i=1n⁡(2−b―iha―ih¯−ha―ih¯b―i).

因为单位四元数群 H1 形成一个李群,所以 Q2 的临界点可以借助导数

(4.155)∂vh=limϑ→0⁡eϑvh−hϑ=vh(v,h 是四元数,ϑ 是实数 )

从

(4.156)∂vQ2=−∑i=1n⁡(b―iv―ha――ih¯+b―iha――i(v―h)―+v―ha―ih¯b――i+ha―i(v―h)―b――i)=0

确定. 这里 v―,b―i 以及 ha―ih¯ 是纯四元数,所以 v¯―=−v― ,因而 (4.156) 可以化简为

(4.157)∂v―Q2=−4v―⋅(∑i=1n⁡ha→ih¯×b―i)=0.

因为这里 v― 是任意的,所以这个表达式当

(4.158)∑i=1n⁡ha→ih¯×b―i=0―

时为零. 设 R 是由单位四元数 h 表示的旋转矩阵,即 ha→ih¯=Ra→i . 应用由向量 a→=(ax,ay,az)∈R3 定义的 3×3 矩阵

(4.159)K(a→)=(0−azayaz0−ax−ayax0),

对于任何向量 b→∈R3 我们得到

(4.160)K(a→)b→=a→×b→,K(K(a→)b→)=b→a→T−a→b→T.

由这个关系式可确定极小值问题的临界点:

(4.161)∑i=1n⁡K(Ra→i×b→i)=O⇔∑i=1n⁡(b→ia→iTRT−Ra→ib→iT)=O⇔RP=PTRT,

其中 P=∑i=1n⁡a→ib→iT . 如果 P 分解为乘积 P=RpTS ,其中矩阵 S 是对称的,而 P=Rp 是旋转矩阵,那么由 (4.161) 推出

(4.162)RRpTS=SRpRT,

并且

(4.163)R=Rp

显然是一个解; 因为在这种情形,由于 RpRpT=E ,所以 RpRpTS=S=SRpRpT 但无论如何, 还有另外 3 个解, 即

(4.164)R=RjRp(j=1,2,3),

其中 Rj 表示绕 S 的第 j 个特征向量旋转 π ,即有 RjSRj=S . 我们可以从

RjRpRpTS=SRpRpTRjT⇔RjS=SRjT⇔RjSRj=S

看出 R=RjRp 是 (4.162) 的解.

使得 Q2 取极小的解是

(4.165)R={Rp,detP>0,Rj0Rp,detP<0,

其中 Rj0 表示绕与 S 的绝对值最小的特征值相应的特征向量旋转 π .

4.4.3.5 向量分析 ​

如果将算子 ∇(参见第 933 页 (13.67),13.2.6.1) 和向量)v→(参见第 919 页 13.1.3) 等同于四元数计算中的 ∇Q 和 v― ,即

(4.166)∇Q=i∂∂x+j∂∂y+k∂∂z(4.167)v―(x,y,z)=v1(x,y,z)i+v2(x,y,z)j+v3(x,y,z)k,

其中 i,j,k 依据 (4.107)(参见第 387 页),那么四元数的乘法法则 (参见第 389 页 (4.109b)) 给出

(4.168)∇Qv―=−∇⋅v→+∇×v→=−divv→+rotv→.

(还可参见第 389 页 4.4.1.1,4. 中 ◼ ). 将

(4.169a)D=∂∂t+i∂∂x+j∂∂y+k∂∂z

及

w(t,x,y,z)=w0(t,x,y,z)+w1(t,x,y,z)i+w2(t,x,y,z)j+w3(t,x,y,z)k(4.169b)=w0(t,x,y,z)+w―(t,x,y,z)

代入, 则有

(4.169c)Dw=∂∂tw0−divw―+rotw―+gradw0.

特别,对于任何二次连续可微函数 t(t,x,y,z) ,

(4.170a)∇Q∇¯Qf=∇¯Q∇Qf=∇∇f=∂2f∂x2+∂2f∂y2+∂2f∂z2=Δ3f,

以及

(4.170b)∇Q∇Qf=−∇∇f=−∂2f∂x2−∂2f∂y2−∂2f∂z2=−Δ3f,

其中 Δ3 表示 R3 中的拉普拉斯算子 (参见第 934 页 13.2.6.5 中 (13.75)).

(4.170c)DD¯f=D¯Df=∂2f∂t2+∂2f∂x2+∂2f∂y2+∂2f∂z2=Δ4f,

其中 Δ4 表示 R4 中的拉普拉斯算子. ∇Q ,恰如 D 一样,常称为狄拉克或柯西-黎曼算子.

4.4.3.6 规范化四元数和刚体运动 ​

1. 八元数 ​

一个八元数 hˇ 有形式

(4.171)hˇ=h0+ϵh1,h0,h1∈H.

这里 ϵ 是对偶单位,它与每个四元数交换,并且 ϵ2=0 . 乘法运算是通常的四元数乘法 (参见第 391 页 (4.115)).

2. 刚体运动 ​

在 R3 中可以借助单位八元数,即具有下列性质的八元数

(4.172)hˇhˇ―=(h0+ϵh1)(h0―+ϵh¯1)=1⇔{h0h¯0=1,h0h¯1+h1h¯0=0

刻画刚体运动 (旋转和平移相互交替) (表 4.1).

表 4.1 用八元数表示刚体运动

元素

表示

空间中点 p―=(px,py,pz)

pˇ=1+p―ϵ ,其中 p―=pxi+pyj+pzk

旋转

r∈H1 ,单位四元数

平移 t―=(tx,ty,tz)

1+12t―ϵ ,其中 t―=txi+tyj+tzk

因为 h 和 −h 刻画相同的刚体运动,所以单位八元数

(4.173)hˇ=h0+h1ϵ=(1+12t―ϵ)r=r+12t―rϵ,t―∈H0,r∈H1,

给出 R3 中刚体运动群 SE(3) 的双重覆盖.

version 1.24.0