三维空间刚体运动4-1:四元数表示变换(各形式相互转换加代码——下篇)

三维空间刚体运动4-1:四元数表示变换(各形式相互转换加代码——下篇)

  • 4. 四元数到其它旋转表示的相互转换
    • 4.1 旋转向量
    • 4.2 旋转矩阵
    • 4.3 欧拉角
      • 4.3.1 转换关系
      • 4.3.2 转换中的万象锁问题
  • 5. 四元数的其他性质
    • 5.1 旋转的复合
    • 5.2 双倍覆盖
    • 5.3 指数形式
    • 5.4 幂函数性质
  • 6. 实践
    • 6.1 四元数旋转运算
    • 6.2 坐标变换

序:本篇系列文章参照高翔老师《视觉SLAM十四讲从理论到实践》,讲解三维空间刚体运动,为读者打下坚实的数学基础。博文将原第三讲分为五部分来讲解,其中四元数部分较多较复杂,又分为四部分。如果读者急于实践,可直接阅读第五部分的机器人运动轨迹,此部分详细讲解了安装准备工作。此系列总体目录如下:

  1. 旋转矩阵和变换矩阵;
  2. 旋转向量表示旋转;
  3. 欧拉角表示旋转;
  4. 四元数包括以下部分:
    4-1.1 四元数表示变换——上篇;
    4-1.2 四元数表示变换——下篇;
    4-2. 四元数线性插值方法:LinEuler/LinMat/Lerp/Nlerp/Slerp;
    4-3. 四元数多点插值方法:Squad;
    4-4. 四元数多点连续解析解插值方法:Spicv;
    4-5. 四元数多点离散数值解插值方法:Sping。
  5. 实践:SLAM中显示机器人运动轨迹及相机位姿。

4. 四元数到其它旋转表示的相互转换

任意单位四元数描述了一个旋转,该旋转也可用旋转向量、旋转矩阵或欧拉角描述。现在来考察四元数与旋转向量、旋转矩阵以及欧拉角的相互转换关系。

4.1 旋转向量

本节介绍旋转向量与四元数之间的转换关系。
绕坐标轴的多次旋转可以等效为绕某一转轴旋转一定的角度。假设已知等效旋转轴方向单位旋转向量 u u u的坐标为 u = [ u x u y u z ] u=\begin{bmatrix} u_{x}\\ u_{y}\\ u_{z} \end{bmatrix} u= uxuyuz ,旋转角为 θ \theta θ。根据3.2.3推导的定理3D旋转公式(一般情况四元数型),设其等效的单位四元数为 q = [ q 0 , q 1 , q 2 , q 3 ] q=[q_{0},q_{1},q_{2},q_{3}] q=[q0,q1,q2,q3],则有: { q 0 = c o s ( 1 2 θ ) q 1 = u x s i n ( 1 2 θ ) q 2 = u y s i n ( 1 2 θ ) q 3 = u z s i n ( 1 2 θ ) (4.1) \left\{\begin{matrix} q_{0}=cos(\frac{1}{2}\theta )\\ q_{1}=u_{x}sin(\frac{1}{2}\theta )\\ q_{2}=u_{y}sin(\frac{1}{2}\theta )\\ q_{3}=u_{z}sin(\frac{1}{2}\theta ) \end{matrix}\right.\tag{4.1} q0=cos(21θ)q1=uxsin(21θ)q2=uysin(21θ)q3=uzsin(21θ)(4.1)
同理,当已知单位四元数 q = [ q 0 , q 1 , q 2 , q 3 ] q=[q_{0},q_{1},q_{2},q_{3}] q=[q0,q1,q2,q3],其对应的旋转角 θ \theta θ和旋转向量 u u u为: { θ = 2 arccos ⁡ q 0 [ u x , u y , u z ] T = [ q 1 , q 2 , q 3 ] T / sin ⁡ ( θ 2 ) (4.2) \left\{\begin{matrix} \theta=2\arccos q_{0}\\ [u_{x},u_{y},u_{z}]^{T}=[q_{1},q_{2},q_{3}]^{T}/\sin(\frac{\theta}{2}) \end{matrix}\right.\tag{4.2} {θ=2arccosq0[ux,uy,uz]T=[q1,q2,q3]T/sin(2θ)(4.2)

4.2 旋转矩阵

已知单位四元数 q q q,根据博文《三维空间刚体运动2:旋转向量与罗德里格斯公式》中的罗德里格斯公式展开式(2.3)及本文公式(4.2),将单位旋转向量 u u u及旋转角 θ \theta θ替换为单位四元数 q q q,省去计算过程,得到单位四元数到旋转矩阵的转换公式为: R = [ 1 − 2 q 2 2 − 2 q 3 2 2 ( q 1 q 2 − q 0 q 3 ) 2 ( q 1 q 3 + q 0 q 2 ) 2 ( q 1 q 2 + q 0 q 3 ) 1 − 2 q 1 2 − 2 q 3 2 2 ( q 2 q 3 − q 0 q 1 ) 2 ( q 1 q 3 − q 0 q 2 ) 2 ( q 2 q 3 + q 0 q 1 ) 1 − 2 q 1 2 − 2 q 2 2 ] (4.3) R=\begin{bmatrix} 1-2q_{2}^{2}-2q_{3}^{2} & 2(q_{1}q_{2}-q_{0}q_{3}) & 2(q_{1}q_{3}+q_{0}q_{2})\\ 2(q_{1}q_{2}+q_{0}q_{3}) & 1-2q_{1}^{2}-2q_{3}^{2} & 2(q_{2}q_{3}-q_{0}q_{1})\\ 2(q_{1}q_{3}-q_{0}q_{2}) & 2(q_{2}q_{3}+q_{0}q_{1}) & 1-2q_{1}^{2}-2q_{2}^{2} \end{bmatrix}\tag{4.3} R= 12q222q322(q1q2+q0q3)2(q1q3q0q2)2(q1q2q0q3)12q122q322(q2q3+q0q1)2(q1q3+q0q2)2(q2q3q0q1)12q122q22 (4.3)
假设已知变换矩阵: R = [ r 11 r 12 r 13 r 21 r 22 r 23 r 31 r 32 r 33 ] R=\begin{bmatrix} r_{11} & r_{12} & r_{13}\\ r_{21} & r_{22} & r_{23}\\ r_{31} & r_{32} & r_{33} \end{bmatrix} R= r11r21r31r12r22r32r13r23r33 根据公式(4.3),推导对应的单位四元数为: { q 0 = 1 2 1 + r 11 + r 22 + r 33 q 1 = r 32 − r 23 4 q 0 q 2 = r 13 − r 31 4 q 0 q 3 = r 21 − r 12 4 q 0 (4.4) \left\{\begin{matrix} q_{0}=\frac{1}{2}\sqrt{1+r_{11}+r_{22}+r_{33}}\\ q_{1}=\frac{r_{32}-r_{23}}{4q_{0}}\\ q_{2}=\frac{r_{13}-r_{31}}{4q_{0}}\\ q_{3}=\frac{r_{21}-r_{12}}{4q_{0}} \end{matrix}\right.\tag{4.4} q0=211+r11+r22+r33 q1=4q0r32r23q2=4q0r13r31q3=4q0r21r12(4.4)

4.3 欧拉角

4.3.1 转换关系

那么将Z-Y-X欧拉角(或RPY角:绕固定坐标系的X-Y-Z依次旋转α,β,γ角)转换为四元数: q = [ cos ⁡ γ 2 0 0 sin ⁡ γ 2 ] [ cos ⁡ β 2 0 sin ⁡ β 2 0 ] [ cos ⁡ α 2 sin ⁡ α 2 0 0 ] = [ cos ⁡ α 2 cos ⁡ β 2 cos ⁡ γ 2 + sin ⁡ α 2 sin ⁡ β 2 sin ⁡ γ 2 sin ⁡ α 2 cos ⁡ β 2 cos ⁡ γ 2 − cos ⁡ α 2 sin ⁡ β 2 sin ⁡ γ 2 cos ⁡ α 2 sin ⁡ β 2 cos ⁡ γ 2 + sin ⁡ α 2 cos ⁡ β 2 sin ⁡ γ 2 cos ⁡ α 2 cos ⁡ β 2 sin ⁡ γ 2 − sin ⁡ α 2 sin ⁡ β 2 cos ⁡ γ 2 ] q=\begin{bmatrix} \cos\frac{\gamma}{2}\\ 0\\ 0\\ \sin\frac{\gamma}{2} \end{bmatrix} \begin{bmatrix} \cos\frac{\beta}{2}\\ 0\\ \sin\frac{\beta}{2}\\ 0 \end{bmatrix} \begin{bmatrix} \cos\frac{\alpha}{2}\\ \sin\frac{\alpha}{2}\\ 0\\ 0 \end{bmatrix} =\begin{bmatrix} \cos\frac{\alpha}{2}\cos\frac{\beta}{2}\cos\frac{\gamma}{2}+\sin\frac{\alpha}{2}\sin\frac{\beta}{2}\sin\frac{\gamma}{2}\\ \sin\frac{\alpha}{2}\cos\frac{\beta}{2}\cos\frac{\gamma}{2}-\cos\frac{\alpha}{2}\sin\frac{\beta}{2}\sin\frac{\gamma}{2}\\ \cos\frac{\alpha}{2}\sin\frac{\beta}{2}\cos\frac{\gamma}{2}+\sin\frac{\alpha}{2}\cos\frac{\beta}{2}\sin\frac{\gamma}{2}\\ \cos\frac{\alpha}{2}\cos\frac{\beta}{2}\sin\frac{\gamma}{2}-\sin\frac{\alpha}{2}\sin\frac{\beta}{2}\cos\frac{\gamma}{2} \end{bmatrix} q= cos2γ00sin2γ cos2β0sin2β0 cos2αsin2α00 = cos2αcos2βcos2γ+sin2αsin2βsin2γsin2αcos2βcos2γcos2αsin2βsin2γcos2αsin2βcos2γ+sin2αcos2βsin2γcos2αcos2βsin2γsin2αsin2βcos2γ
根据上面的公式可以求出逆解,即由四元数 q = ( q 0 , q 1 , q 2 , q 3 ) q=(q_{0},q_{1},q_{2},q_{3}) q=(q0,q1,q2,q3)到欧拉角的转换为: [ α β γ ] = [ a t a n 2 ( 2 ( q 0 q 1 + q 2 q 3 ) , 1 − 2 ( q 1 2 + q 2 2 ) ) arcsin ⁡ ( 2 ( q 0 q 2 − q 1 q 3 ) a t a n 2 ( 2 ( q 0 q 3 + q 1 q 2 ) , 1 − 2 ( q 2 2 + q 3 2 ) ) ] \begin{bmatrix} \alpha\\ \beta\\ \gamma \end{bmatrix} =\begin{bmatrix} atan2(2(q_{0}q_{1}+q_{2}q_{3}),1-2(q^{2}_{1}+q^{2}_{2}))\\ \arcsin (2(q_{0}q_{2}-q_{1}q_{3})\\ atan2(2(q_{0}q_{3}+q_{1}q_{2}),1-2(q^{2}_{2}+q^{2}_{3})) \end{bmatrix} αβγ = atan2(2(q0q1+q2q3),12(q12+q22))arcsin(2(q0q2q1q3)atan2(2(q0q3+q1q2),12(q22+q32)) 这里使用了 a t a n 2 ( y , x ) atan2(y,x) atan2(y,x)而不是 a r c t a n ( y x ) arctan(\frac{y}{x}) arctan(xy),因为 a r c t a n ( y x ) arctan(\frac{y}{x}) arctan(xy)的取值范围为 [ − π 2 , π 2 ] [-\frac{\pi}{2},\frac{\pi}{2}] [2π,2π],只有 180 ° 180° 180°,而绕某个轴旋转时范围是 360 ° 360° 360° a t a n 2 ( y , x ) atan2(y,x) atan2(y,x)正好满足需求。 a t a n 2 ( y , x ) atan2(y,x) atan2(y,x)是一个函数,在C语言里返回的是指方位角,函数原型为 d o u b l e a t a n 2 ( d o u b l e y , d o u b l e x ) double \space atan2(double y, double x) double atan2(doubley,doublex) ,返回以弧度表示的 y / x y/x y/x 的反正切。y 和 x 的值的符号决定了正确的象限。也可以理解为计算复数 x + y i x+yi x+yi 的辐角,计算时atan2 比 atan 稳定。

4.3.2 转换中的万象锁问题

另外需要注意的就是奇异性问题,即万向锁,下面分析这种情况。当刚体绕Y轴旋转了 90 ° 90° 90°(俯仰角pitch=90)时,如何计算横滚角roll和偏航角yaw?这时会发生自由度丢失的情况,即Yaw和Roll会变为一个自由度。此时再使用上面的公式根据四元数计算欧拉角会出现问题: arcsin ⁡ ( 2 ( q 0 q 2 − q 1 q 3 ) \arcsin (2(q_{0}q_{2}-q_{1}q_{3}) arcsin(2(q0q2q1q3)的定义域为 [ − 1 , 1 ] [-1,1] [1,1],当 2 ( q 0 q 2 − q 1 q 3 ) = ± 0.5 2(q_{0}q_{2}-q_{1}q_{3})=\pm 0.5 2(q0q2q1q3)=±0.5时(在程序中浮点数不能直接进行等于判断,要使用合理的阈值),俯仰角 β \beta β ± 90 ° \pm 90° ±90°,将其带入正向公式计算出四元数 ( q 0 , q 1 , q 2 , q 3 ) (q_{0},q_{1},q_{2},q_{3}) (q0,q1,q2,q3),然后可以发现逆向公式中atan2函数中的参数全部为0,即出现了 0 0 \frac{0}{0} 00的情况!无法计算。
β = 90 ° \beta=90° β=90°时, sin ⁡ β 2 = cos ⁡ β 2 = 0.707 \sin \frac{\beta}{2}=\cos \frac{\beta}{2}=0.707 sin2β=cos2β=0.707,将其带入公式中有: q = [ q 0 q 1 q 2 q 3 ] = [ 0.707 ( cos ⁡ α 2 cos ⁡ γ 2 + sin ⁡ α 2 sin ⁡ γ 2 ) 0.707 ( sin ⁡ α 2 cos ⁡ γ 2 − cos ⁡ α 2 sin ⁡ γ 2 ) 0.707 ( cos ⁡ α 2 cos ⁡ γ 2 + sin ⁡ α 2 sin ⁡ γ 2 ) 0.707 ( cos ⁡ α 2 sin ⁡ γ 2 − sin ⁡ α 2 cos ⁡ γ 2 ) ] = [ 0.707 cos ⁡ α − γ 2 0.707 sin ⁡ α − γ 2 0.707 cos ⁡ α − γ 2 − 0.707 sin ⁡ α − γ 2 ] q=\begin{bmatrix} q_{0}\\ q_{1}\\ q_{2}\\ q_{3} \end{bmatrix} =\begin{bmatrix} 0.707(\cos\frac{\alpha}{2}\cos\frac{\gamma}{2}+\sin \frac{\alpha}{2}\sin \frac{\gamma}{2})\\ 0.707(\sin\frac{\alpha}{2}\cos\frac{\gamma}{2}-\cos\frac{\alpha}{2}\sin \frac{\gamma}{2})\\ 0.707(\cos\frac{\alpha}{2}\cos\frac{\gamma}{2}+\sin \frac{\alpha}{2}\sin \frac{\gamma}{2})\\ 0.707(\cos\frac{\alpha}{2}\sin\frac{\gamma}{2}-\sin \frac{\alpha}{2}\cos\frac{\gamma}{2}) \end{bmatrix} =\begin{bmatrix} 0.707\cos\frac{\alpha - \gamma}{2}\\ 0.707\sin \frac{\alpha - \gamma}{2}\\ 0.707\cos \frac{\alpha - \gamma}{2}\\ -0.707\sin \frac{\alpha - \gamma}{2} \end{bmatrix} q= q0q1q2q3 = 0.707(cos2αcos2γ+sin2αsin2γ)0.707(sin2αcos2γcos2αsin2γ)0.707(cos2αcos2γ+sin2αsin2γ)0.707(cos2αsin2γsin2αcos2γ) = 0.707cos2αγ0.707sin2αγ0.707cos2αγ0.707sin2αγ q 1 q 0 = − q 3 q 2 = tan ⁡ α − γ 2 \frac{q_{1}}{q_{0}}=-\frac{q_{3}}{q_{2}}=\tan\frac{\alpha - \gamma}{2} q0q1=q2q3=tan2αγ,于是有: α − γ = 2 a t a n 2 ( q 1 , q 0 ) \alpha - \gamma = 2atan2(q_{1},q_{0}) αγ=2atan2(q1,q0)通常令 α = 0 \alpha=0 α=0,这时 γ = − 2 a t a n 2 ( q 1 , q 0 ) \gamma = -2atan2(q_{1},q_{0}) γ=2atan2(q1,q0)。当俯仰角 β \beta β为-90°,即刚体竖直向下时的情况也与之类似,可以推导出奇异姿态时的计算公式。
将四元数转换为欧拉角的代码可以参考附件。需要注意欧拉角有12种旋转次序,而上面推导的公式是按照Z-Y-X顺序进行的,所以有时会在网上看到不同的转换公式(因为对应着不同的旋转次序),在使用时一定要注意旋转次序是什么。比如ADAMS软件里就默认Body 3-1-3次序,即Z-X-Z欧拉角,而VREP中则按照X-Y-Z欧拉角旋转。另外万向锁问题代码的Z-Y-X欧拉角代码见类CameraSpacePoint。

5. 四元数的其他性质

为更全面理解四元数和方便引入slerp插值,这一节补充四元数的其它性质:旋转复合、双倍覆盖和指数形式。

5.1 旋转的复合

旋转的复合其实在我们之前证明 q 2 = q q = [ c o s ( 2 θ ) , s i n ( 2 θ ) u ] q^{2}=qq=[cos(2\theta),sin(2\theta)u] q2=qq=[cos(2θ),sin(2θ)u]的时候就已经涉及到一点了,但是这里我们考虑的是更一般的情况。
假设有两个表示沿着不同轴、不同角度旋转的四元数 q 1 、 q 2 q_{1}、q_{2} q1q2。首先,我们实施 q 1 q_{1} q1的变换,变换之后的 v ′ v^{'} v v ′ = q 1 v q 1 ∗ v^{'}=q_{1}vq^{*}_{1} v=q1vq1接下来,对 v ′ v^{'} v进行 q 2 q_{2} q2的变换,得到 v ′ ′ v^{''} v′′ v ′ ′ = q 2 v ′ q 2 ∗ = q 2 q 1 v q 1 ∗ q 2 ∗ = q n e t v q n e t ∗ v^{''}=q_{2}v^{'}q^{*}_{2}=q_{2}q_{1}vq^{*}_{1}q^{*}_{2}=q_{net}vq^{*}_{net} v′′=q2vq2=q2q1vq1q2=qnetvqnet其中 q n e t = q 2 q 1 q_{net}=q_{2}q_{1} qnet=q2q1,另外写成上面的形式,还用到性质 q 1 ∗ q 2 ∗ = ( q 2 q 1 ) ∗ q^{*}_{1}q^{*}_{2}=(q_{2}q_{1})^{*} q1q2=(q2q1)
注意四元数乘法的顺序,这和矩阵、复数的复合非常相似,都是从右向左叠加。另外要注意的是, q 1 q_{1} q1 q 2 q_{2} q2的等价旋转 q n e t q_{net} qnet并不是沿着 q 1 q_{1} q1 q 2 q_{2} q2的两个旋转轴进行的两次旋转,它是沿着一个全新的旋转轴进行的一个等价旋转,仅仅只有旋转的结果相同。

5.2 双倍覆盖

四元数与 3D 旋转的关系并不是一对一的,同一个 3D 旋转可以使用两个不同的四元数来表示.对任意的单位四元数 q = [ c o s ( 1 2 θ ) , s i n ( 1 2 θ ) u ] q=[cos(\frac{1}{2}\theta),sin(\frac{1}{2}\theta)u] q=[cos(21θ),sin(21θ)u] q q q − q −q q代表的是同一个旋转。如果 q q q表示的是沿着旋转轴 u u u旋转 θ θ θ度,那么 − q −q q代表的是沿着相反的旋转轴 − u −u u旋转 ( 2 π − θ ) (2π − θ) (2πθ),负负得正: − q = [ − c o s ( 1 2 θ ) , − s i n ( 1 2 θ ) u ] = [ c o s ( π − 1 2 θ ) , s i n ( π − 1 2 θ ) ( − u ) ] \begin{aligned} -q &= [-cos(\frac{1}{2}\theta),-sin(\frac{1}{2}\theta)u] \\ &= [cos(\pi - \frac{1}{2}\theta),sin(\pi - \frac{1}{2}\theta)(-u)]\end{aligned} q=[cos(21θ),sin(21θ)u]=[cos(π21θ),sin(π21θ)(u)]所以,这个四元数旋转的角度为 2 ( π − 1 2 θ ) = 2 π − θ 2(\pi - \frac{1}{2}\theta)=2\pi-\theta 2(π21θ)=2πθ,从下面的图中我们可以看到,这两个旋转是完全等价的:在这里插入图片描述
其实从四元数的旋转公式中也能推导出相同的结果: ( − q ) v ( − q ) ∗ = ( − 1 ) 2 q v q ∗ = q v q ∗ (-q)v(-q)^{*}=(-1)^{2}qvq^{*}=qvq^{*} (q)v(q)=(1)2qvq=qvq所以,我们经常说单位四元数与3D旋转有一个「2对1满射同态」(2-1 Surjective Homomorphism) 关系,或者说单位四元数双倍覆盖(Double Cover)了3D旋转。它的严格证明会用到一些李群的知识,这里不再拓展。
有一点需要注意的是,虽然 q q q − q −q q是两个不同的四元数,但是由于旋转矩阵中的每一项都包含了四元数两个分量的乘积,它们的旋转矩阵是完全相同的,即旋转矩阵并不会出现双倍覆盖的问题。

5.3 指数形式

在讨论复数的时候,我们利用欧拉公式将2D的旋转写成了 v ′ = e i θ v v^{'}=e^{i\theta v} v=eiθv这样的指数形式。实际上,我们也可以利用四元数将 3D 旋转写成类似的形式。
类似于复数的欧拉公式,四元数也有一个类似的公式。如果 u u u是一个单位向量,那么对于单位四元数 u q = [ 0 , u ] u_{q}= [0, u] uq=[0,u],有: e u θ = c o s ( θ ) + u q s i n ( θ ) = c o s ( θ ) + u s i n ( θ ) e^{u\theta}=cos(\theta)+u_{q}sin(\theta)=cos(\theta)+usin(\theta) euθ=cos(θ)+uqsin(θ)=cos(θ)+usin(θ)这也就是说, q = [ cos ⁡ ( θ ) , sin ⁡ ( θ ) u ] q = [\cos(\theta), \sin(\theta)u] q=[cos(θ),sin(θ)u]可以使用指数表示为 e u θ e^{u\theta} euθ。这个公式的证明与欧拉公式的证明非常类似,直接使用级数展开就可以了,这里不再扩展。
注意,因为 u u u是一个单位向量, u 2 = [ − u ⋅ u , 0 ] = − ∥ u ∥ 2 = − 1 u^{2}= [−u · u, 0] = −∥u∥^{2}= −1 u2=[uu,0]=u2=1.这与欧拉公式中的 i i i是非常类似的。有了指数型的表示方式,我们就能够将之前四元数的旋转公式改写为指数形式了,由此可得到定义:
3D旋转公式(指数型):任意向量 v v v沿着以单位向量定义的旋转轴 u u u旋转 θ θ θ角度之后的 v ′ v^{'} v可以使用四元数的指数表示。令 w = [ 0 , v ] , u q = [ 0 , u ] w= [0, v], u_{q}= [0, u] w=[0,v],uq=[0,u],那么: w ′ = e u q θ 2 v e − u q θ 2 w^{'}=e^{u_{q}\frac{\theta}{2}}ve^{-u_{q}\frac{\theta}{2}} w=euq2θveuq2θ有了四元数的指数定义,我们就能够定义四元数的更多运算了。首先是自然对数 l o g log log,对任意单位四元数 q = [ c o s ( θ ) , s i n ( θ ) u ] = e u q θ q = [cos(θ), sin(θ)u]=e^{u_{q}\theta} q=[cos(θ),sin(θ)u]=euqθ l o g ( q ) = l o g ( e u q θ ) = [ 0 , u θ ] log(q)=log(e^{u_{q}\theta})=[0,{u\theta}] log(q)=log(euqθ)=[0,uθ]接下来是单位四元数的幂运算: q t = ( e u q θ ) t = e u q ( t θ ) = [ c o s ( t θ ) , s i n ( t θ ) u ] q^{t}=(e^{u_{q}\theta})^{t}=e^{u_{q}(t\theta)}=[cos(tθ), sin(tθ)u] qt=(euqθ)t=euq()=[cos(),sin()u]可以看到,一个单位四元数的 t t t次幂等同于将它的旋转角度缩放至 t t t倍,并且不会改变它的旋转轴( u u u必须是单位向量,所以一般不能与 t t t结合)。这些运算会在之后讨论四元数插值时非常有用。限于篇幅,四元数插值的讲解见下一篇博客《三维空间刚体运动4-2:四元数插值(lerp,Nlerp,Slerp,Squad等)》

5.4 幂函数性质

上面给出四元数的指数形式后,可推导出一个重要的指函数性质,这一性质下一篇将用到。首先看推论:
四元数推论1 :设 q ∈ H 1 , p = [ a , b v ] q\in \mathbb{H}_{1},p=[a,b\mathbf{v}] qH1,p=[a,bv],其中 a , b , r ∈ R , v ∈ R 3 a,b,r\in \mathbb{R},\mathbf{v}\in\mathbb{R}^{3} a,b,rR,vR3,如果 q [ r , v ] q ∗ = [ r , v ′ ] q[r,\mathbf{v}]q^{*}=[r,\mathbf{v^{'}}] q[r,v]q=[r,v],那么 q [ a , b v ] q ∗ = [ a , b v ′ ] q[a,b\mathbf{v}]q^{*}=[a,b\mathbf{v^{'}}] q[a,bv]q=[a,bv]
推论证明: q p q ∗ = q [ a , b v ] q ∗ = q b [ a b , v ] q ∗ = b [ a b , v ′ ] = [ a , b v ′ ] \begin{aligned}qpq^{*}&=q[a,b\mathbf{v}]q^{*}\\&=qb[\frac{a}{b},\mathbf{v}]q^{*}\\&=b[\frac{a}{b},\mathbf{v^{'}}]\\&=[a,b\mathbf{v^{'}}]\end{aligned} qpq=q[a,bv]q=qb[ba,v]q=b[ba,v]=[a,bv]根据推论,得出以下定理:
四元数定理2 :设 q , p ∈ H 1 , p = [ cos ⁡ θ , sin ⁡ θ v ] , t ∈ R q,p\in \mathbb{H}_{1},p=[\cos\theta,\sin\theta\mathbf{v}],t\in\mathbb{R} q,pH1,p=[cosθ,sinθv],tR,那么 q p t q ∗ = ( q p q ∗ ) t qp^{t}q^{*}=(qpq^{*})^{t} qptq=(qpq)t
定理证明:
根据推论,存在 v ′ ∈ R 3 \mathbf{v^{'}}\in\mathbb{R}^{3} vR3,使得 q [ cos ⁡ θ , sin ⁡ θ v ] q ∗ = [ cos ⁡ θ , sin ⁡ θ v ′ ] q[\cos\theta,\sin\theta\mathbf{v}]q^{*}=[\cos\theta,\sin\theta\mathbf{v^{'}}] q[cosθ,sinθv]q=[cosθ,sinθv],因此得到: q p t q ∗ = q ( exp ⁡ ( t log ⁡ p ) ) q ∗ = q ( exp ⁡ ( t [ 0 , θ v ] ) ) q ∗ = q ( exp ⁡ [ 0 , t θ v ] ) q ∗ = q ( [ cos ⁡ t θ , sin ⁡ t θ v ] ) q ∗ = [ cos ⁡ t θ , sin ⁡ t θ v ′ ] = exp ⁡ ( t [ 0 , θ v ′ ] ) = exp ⁡ ( t log ⁡ [ cos ⁡ θ , sin ⁡ θ v ′ ] ) = exp ⁡ ( t log ⁡ ( q p q ∗ ) ) = ( q p q ∗ ) t \begin{aligned} qp^{t}q^{*} &=q(\exp(t\log p))q^{*} \\&= q(\exp(t[0,\theta\mathbf{v}]))q^{*} \\&= q(\exp[0,t\theta\mathbf{v}])q^{*} \\&= q([\cos t\theta,\sin t\theta\mathbf{v}])q^{*} \\ &= [\cos t\theta,\sin t\theta\mathbf{v^{'}}] \\&= \exp(t[0,\theta\mathbf{v^{'}}]) \\&= \exp(t\log[\cos \theta,\sin \theta\mathbf{v^{'}}]) \\&= \exp(t\log(qpq^{*})) \\&= (qpq^{*})^{t}\end{aligned} qptq=q(exp(tlogp))q=q(exp(t[0,θv]))q=q(exp[0,v])q=q([cos,sinv])q=[cos,sinv]=exp(t[0,θv])=exp(tlog[cosθ,sinθv])=exp(tlog(qpq))=(qpq)t
至此,四元数表示旋转的理论部分讲解完了。

6. 实践

现在,我们通过两个小程序实际演练四元数的运算。其中四元数的基础运算和高阶运算程序实现放在下一讲Slerp中。

6.1 四元数旋转运算

第一个小程序是演示四元数的常规运算,包括与旋转矩阵、旋转向量的转换,以及用四元数旋转一个向量,如下所示:

#include<iostream>
#include<cmath>
using namespace std;#include<eigen3/Eigen/Core>
#include<eigen3/Eigen/Geometry>using namespace Eigen;int main(int argc, char **argv){//Eigen/Geometry模块提供了各种旋转和平移的表示,3D旋转矩阵直接使用Matrix3d或Matrix3fMatrix3d rotation_matrix = Matrix3d::Identity();//旋转向量使用AngleAxis,它底层不直接是Matrix,但运算可以当做矩阵,因为重载了运算符AngleAxisd rotation_vector(M_PI/4, Vector3d(0, 0, 1));//设置输出精度cout.precision(3);cout<<"rotation matrix =\n"<<rotation_matrix<<endl;cout<<"rotation vector =\n"<<rotation_vector.matrix()<<endl;//旋转向量转换的矩阵可以直接赋值给旋转矩阵rotation_matrix = rotation_vector.toRotationMatrix();cout<<"rotation vector to Matrix =\n"<<rotation_matrix<<endl;//旋转向量Vector3d v(1, 0, 0);cout<<"v = "<<v.transpose()<<endl;//四元数,可以直接把AngleAxis赋值给四元数,反之亦然Quaterniond q = Quaterniond(rotation_vector);//coeffs:多项式系数(coefficients),其顺序为(x,y,z,w),w为实部,xyz为虚部cout<<"quaternion from rotation vector = "<<q.coeffs().transpose()<<endl;//也可以直接把旋转矩阵赋给它q = Quaterniond(rotation_matrix);cout<<"quaternion from rotation matrix = "<<q.coeffs().transpose()<<endl;//使用四元数旋转一个向量,使用重载的乘法即可//注意q*v在数学上是qvq^{-1}v_rotated = q*v;cout<<"(1,0,0) after rotation = "<<v_rotated.transpose()<<endl;//用常规向量乘法表示qvq^{-1},则计算如下:Quaterniond q_rotate_v = q*Quaterniond(0,1,0,0)*q.inverse();cout<<"should be equal to "<<q_rotate_v.coeffs().transpose()<<endl;return 0;

输出结果如下:

rotation matrix =
1 0 0
0 1 0
0 0 1
rotation vector =0.707 -0.707      00.707  0.707      00      0      1
rotation vector to Matrix =0.707 -0.707      00.707  0.707      00      0      1
v = 1 0 0
v transofrmed = 1.71 3.71    4
quaternion from rotation vector =     0     0 0.383 0.924
quaternion from rotation matrix =     0     0 0.383 0.924
(1,0,0) after rotation = 0.707 0.707     0
should be equal to 0.707 0.707     0     0

请读者注意,程序代码通常和数学表示有一些细微的差别。例如,通过运算符重载,四元数和三维向量可以直接计算乘法,但在数学上则需要先把向量转成虚四元数,再利用四元数乘法进行计算,同样的情况也适用于变换矩阵乘三维向量的情况。总体而言,程序中的用法会比数学公式更灵活。

6.2 坐标变换

下面我们举一个小例子来演示坐标变换。
例子设有小萝卜一号和小萝卜二号位于世界坐标系中。记世界坐标系为 W W W,小萝卜们的坐标系分别为 R 1 R_{1} R1 R 2 R_{2} R2。小萝卜一号的位姿为 q 1 = [ 0.35 , 0.2 , 0.3 , 0.1 ] T q_{1}=[0.35,0.2,0.3,0.1]^{T} q1=[0.35,0.2,0.3,0.1]T t 1 = [ 0.3 , 0.1 , 0.1 ] T t_{1}=[0.3,0.1,0.1]^{T} t1=[0.3,0.1,0.1]T。小萝卜二号的位姿为 q 2 = [ − 0.5 , 0.4 , − 0.1 , 0.2 ] T q_{2}=[-0.5,0.4,-0.1,0.2]^{T} q2=[0.5,0.4,0.1,0.2]T t 2 = [ − 0.1 , 0.5 , 0.3 ] T t_{2}=[-0.1,0.5,0.3]^{T} t2=[0.1,0.5,0.3]T。这里的 q q q t t t表达的是 T R k , W , k = 1 , 2 T_{R_{k},W},k=1,2 TRk,W,k=1,2,也就是世界坐标系到相机坐标系的变换关系。现在,小萝卜一号看到某个点在自身的坐标系下坐标为 p R 1 = [ 0.5 , 0 , 0.2 ] T p_{R_{1}}=[0.5,0,0.2]^{T} pR1=[0.5,0,0.2]T,求该向量在小萝卜二号坐标系下的坐标。
这是一个非常简单,但又具有代表性的例子。在实际场景中你经常需要在同一个机器人的不同部分,或者不同机器人之间转换坐标。计算过程也很简单,只需计算: p R 2 = T R 2 , W T W , R 1 p R 1 p_{R_{2}}=T_{R_{2},W}T_{W,R_{1}}p_{R_{1}} pR2=TR2,WTW,R1pR1下面请看程序:

#include<iostream>
#include<vector>
#include<algorithm>
#include<eigen3/Eigen/Core>
#include<eigen3/Eigen/Geometry>using namespace std;
using namespace Eigen;int main(int argc, char** argv){//位姿四元数q1,q2Quaterniond q1(0.35, 0.2, 0.3, 0.1), q2(-0.5, 0.4, -0.1, 0.2);cout<<"q1.coeffs=\n"<<q1.coeffs()<<endl;cout<<"q2.coeffs=\n"<<q2.coeffs()<<endl;//四元数转换为旋转矩阵cout<<"q1.matrix=\n"<<q1.matrix()<<endl;cout<<"q2.matrix=\n"<<q2.matrix()<<endl;//归一化,转换为单位四元数q1.normalize();q2.normalize();cout<<"q1.normalize="<<q1.matrix()<<endl;cout<<"q2.normalize="<<q2.matrix()<<endl;//位移向量t1,t2,小萝卜一号下的坐标pr1Vector3d t1(0.3, 0.1, 0.1), t2(-0.1, 0.5, 0.3);Vector3d pr1(0.5, 0, 0.2);//Eigen::Isometry3d:欧式变换矩阵4*4Isometry3d Twr1(q1), Tr2w(q2);cout<<"Twr1.matrix="<<Twr1.matrix()<<endl;cout<<"Tr2w.matrix="<<Tr2w.matrix()<<endl;//坐标转换,计算T=Ra+tTwr1.pretranslate(t1);Tr2w.pretranslate(t2);//先将pr1转换为世界坐标系,然后转换为小萝卜二号下的坐标pr2Vector3d pr2 = Tr2w * Twr1.inverse()*pr1;cout<<"pr2.transpose="<<pr2.transpose()<<endl;return 0;
}

输出结果如下:

q1.coeffs=0.20.30.1
0.35
q2.coeffs=0.4
-0.10.2
-0.5
q1.matrix=0.8  0.05  0.250.19   0.9 -0.08
-0.17   0.2  0.74
q2.matrix=0.9  0.12  0.26
-0.28   0.6  0.360.06 -0.44  0.66
q1.normalize=0.238095   0.190476   0.9523810.72381   0.619048  -0.304762-0.647619   0.761905 0.00952381
q2.normalize=0.782609   0.26087  0.565217
-0.608696  0.130435  0.7826090.130435 -0.956522   0.26087
Twr1.matrix=0.238095   0.190476   0.952381          00.72381   0.619048  -0.304762          0-0.647619   0.761905 0.00952381          00          0          0          1
Tr2w.matrix=0.782609   0.26087  0.565217         0
-0.608696  0.130435  0.782609         00.130435 -0.956522   0.26087         00         0         0         1
pr2.transpose=
-0.0309731    0.73499   0.296108

四元数的第一部分讲解到此。

本文基于《视觉SLAM十四讲:从理论到实践》和《Quaternions, Interpolation and Animation》编写,但相对于原文会适当精简,同时为便于全面理解,会收集其他网络好文,根据作者理解,加入一些注解和扩展知识点,如果您觉得还不错,请一键四连(点赞关注收藏评论),让更多的人看到。

参考文献:

  1. 《视觉SLAM十四讲:从理论到实践》,高翔、张涛等著,中国工信出版社
  2. Quaternions, Interpolation and Animation
  3. 四元数与三维旋转
  4. 四元数与欧拉角(RPY角)的相互转换
  5. 如何形象地理解四元数?

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/pingmian/64297.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

使用layui的table提示Could not parse as expression(踩坑记录)

踩坑记录 报错图如下 原因&#xff1a; 原来代码是下图这样 上下俩中括号都是连在一起的&#xff0c;可能导致解析问题 改成如下图这样 重新启动项目&#xff0c;运行正常&#xff01;

大模型的构建与部署(2)——数据清洗

版权声明 本文原创作者:谷哥的小弟作者博客地址:http://blog.csdn.net/lfdfhl1. 数据清洗的必要性与影响 1.1 数据清洗对模型性能的影响 数据清洗是数据预处理的关键步骤,对于模型训练的性能和准确性有着直接的影响。原始数据中的缺失值、重复值、异常值以及数据格式不一致…

【MySQL】--- 数据库基础

Welcome to 9ilks Code World (๑•́ ₃ •̀๑) 个人主页: 9ilk (๑•́ ₃ •̀๑) 文章专栏&#xff1a; MySQL 本篇博客我们来建立一下数据库的相关概念&#xff0c;主要理解什么是数据库以及mysql和mysqld&#xff0c;MySQL架构等问题。 &#x1f3e0; 登录…

Vue中纯前端实现导出简单Excel表格的功能

Vue 前端Excel导出 Vue中纯前端导出简单Excel表格的方法(使用vue-json-excel插件) 前言 在许多的后台系统中少不了导出Excel表格的功能&#xff0c;在项目中纯前端使用vue-json-excel插件来实现简单Excel表格的导出功能。 使用方法 1、安装依赖 npm install vue-json-exc…

KeyFormer:使用注意力分数压缩KV缓存

Keyformer: KV Cache Reduction through Key Tokens Selection for Efficient Generative Inference 202403&#xff0c;发表在Mlsys Introduction 优化KV cache的策略&#xff0c;主要是集中在系统级别的优化上&#xff0c;比如FlashAttention、PagedAttention&#xff0c;它…

Linux 权限管理实践:精确控制用户对 systemctl 和 journalctl 命令的使用

前言 在 Linux 系统管理中&#xff0c;精确控制用户对特定命令的访问权限是一项关键的安全实践。使用 systemctl 和 journalctl 命令时&#xff0c;不当的权限设置可能会导致不必要的风险。本篇博客将详细讨论如何通过 sudoers 文件和 Polkit 策略为不同用户配置 systemctl 和…

SSH连接成功,但VSCode连接不成功

环境 在实验室PC上连接服务器234 解决方案&#xff1a;在VSCode中重新添加远程主机 删除旧的VSCode Server 在远程主机上&#xff0c;VSCode会安装一个‘vscode-server’服务来支持远程开发&#xff0c;有时旧的‘vscode-server’文件可能会导致问题&#xff0c;删除旧的&am…

【Qt】qt安装

在工作一年之后&#xff0c;还是想做一个Qt的教程&#xff0c;遥想研一刚刚接触Qt&#xff0c;从0到1学习&#xff0c;没有什么参考书籍&#xff0c;网上的资料也不多&#xff0c;幸好Qt官方文档写得好&#xff0c;加上自己肯研究&#xff0c;才堪堪入门。 现在我想自己写一个…

Web开发 -前端部分-CSS

CSS CSS&#xff08;Cascading Style Sheet&#xff09;:层叠样式表&#xff0c;用于控制页面的样式&#xff08;表现&#xff09;。 一 基础知识 1 标题格式 标题格式一&#xff1a; 行内样式 <!DOCTYPE html> <html lang"en"><head><meta…

YOLOv8目标检测(六)_封装API接口

YOLOv8目标检测(一)_检测流程梳理&#xff1a;YOLOv8目标检测(一)_检测流程梳理_yolo检测流程-CSDN博客 YOLOv8目标检测(二)_准备数据集&#xff1a;YOLOv8目标检测(二)_准备数据集_yolov8 数据集准备-CSDN博客 YOLOv8目标检测(三)_训练模型&#xff1a;YOLOv8目标检测(三)_训…

51c视觉~YOLO~合集6~

我自己的原文哦~ https://blog.51cto.com/whaosoft/12830685 一、其他yolo 1.1 Spiking-YOLO​ 使用常规深度神经网络到脉冲神经网络转换方法应用于脉冲神经网络域时&#xff0c;性能下降的很多&#xff0c;深入分析后提出了可能的解释&#xff1a;一是来自逐层归一化的效率…

如何在 Ubuntu 22.04 上安装 Strapi CMS

简介 Strapi 是一个使用 JavaScript 构建的开源、无头内容管理系统 (CMS)。与其他无头 CMS 一样&#xff0c;Strapi 开箱即用不带前端。它使用 API 作为其前端&#xff0c;允许你使用流行的框架&#xff08;如 React 和 Next.js&#xff09;构建网站。Strapi 基于插件系统&…

数字IC后端零基础入门基础理论(Day1)

数字IC后端设计导入需要用到的input数据如下图所示。 数字后端零基础入门系列 | Innovus零基础LAB学习Day9 Netlist: 设计的Gate level&#xff08;门级&#xff09;网表。下图所示为一个计数器设计综合后的门级netlist。 从这个netlist中我们看到这个设计顶层的名字叫counte…

序列模型的使用示例

序列模型的使用示例 1 RNN原理1.1 序列模型的输入输出1.2 循环神经网络&#xff08;RNN&#xff09;1.3 RNN的公式表示2 数据的尺寸 3 PyTorch中查看RNN的参数4 PyTorch中实现RNN&#xff08;1&#xff09;RNN实例化&#xff08;2&#xff09;forward函数&#xff08;3&#xf…

WSL2内部的Ubuntu怎么设置网络内桥接模式,弄了好久老是不成功,怎么办?

环境: Win10专业版 WSL2 Ubuntu22.04 问题描述: WSL2内部的Ubuntu怎么设置网络内桥接模式 解决方案: 方法一 1.控制面板开启,Hyper-V 管理器 2.重启电脑 3…创建外部虚拟交换机 打开 Hyper-V 管理器,在右侧操作面板中点击“虚拟交换机管理器”。 选择“创建虚…

redis集群 服务器更换ip,怎么办,怎么更换redis集群的ip

redis集群 服务器更换ip&#xff0c;怎么办&#xff0c;怎么更换redis集群的ip 1、安装redis三主三从集群2、正常状态的redis集群3、更改redis集群服务器的ip 重启服务器 集群会down4、更改redis集群服务器的ip 重启服务器 集群down的原因5、更改redis集群服务器的ip后&#xf…

记录学习《手动学习深度学习》这本书的笔记(五)

这一章是循环神经网络&#xff0c;太难了太难了&#xff0c;有很多卡壳的地方理解了好久&#xff0c;比如隐藏层和隐状态的区别、代码的含义&#xff08;为此专门另写了一篇【笔记】记录对自主实现一个神经网络的步骤的理解&#xff09;、梯度计算相关&#xff08;【笔记】记录…

人大金仓数据linux安装注意事项

人大金仓数据linux安装注意事项 本次是个人搭建虚拟机安装centos7的环境下进行安装。 1、安装流程参照https://help.kingbase.com.cn/v9/install-updata/install-linux/preface.html。 2、mount安装文件报错 操作手册提供mount的命令如下&#xff1a; mount KingbaseES_V009R0…

【GIS教程】使用GDAL-Python将tif转为COG并在ArcGIS Js前端加载-附完整代码

目录 一、数据格式 二、COG特点 三、使用GDAL生成COG格式的数据 四、使用ArcGIS Maps SDK for JavaScript加载COG格式数据 一、数据格式 COG&#xff08;Cloud optimized GeoTIFF&#xff09;是一种GeoTiff格式的数据。托管在 HTTP 文件服务器上&#xff0c;可以代替geose…

探索智能时代:如何利用AI一键生成PPT改变演示文稿的制作方式

在这个科技飞速发展的时代&#xff0c;信息的传递方式发生了翻天覆地的变化。曾几何时&#xff0c;我们还在为制作PPT而熬夜&#xff0c;手动选择模板、调整布局&#xff0c;甚至为每一张幻灯片的内容苦思冥想。然而&#xff0c;随着人工智能技术的不断进步&#xff0c;制作PPT…