2.1 空间变换及其参数化
空间变换及其对应的参数化方法是配准任务进行形式化表达和求解的基础,如形式化表示和计算刚性配准涉及的欧式变换。不同优化方法进行变换参数估计时,都需要对变换进行参数化表示,以适应优化求解方法条件等要求。
2.1.1 什么是欧式变换与变换矩阵
众所周知,常见的物体存在于三维空间中,空间中的任意位置都可以由三维坐标表示。对于三维空间中的一个刚性物体,除了三维空间位置表示以外还具有自身的姿态表示。如图2-1所示,假设在该线性空间中存在一组基底(e1,e2,e3),则任意向量a在该基底下的坐标为(a1,a2,a3)T,我们通常选择e1=(1,0,0)T,e2=(0,1,0)T,e3=(0,0,1)T。向量a可以用基向量与坐标的线性组合得到,即:
图2-1 向量在空间坐标系中的表示
两个坐标系之间的运动(表示两个坐标系之间的空间变换)是由一个旋转和一个平移所组成,通常称这种运动为刚体运动。由于刚体运动保持了向量的长度和夹角,相当于把一个刚体原封不动地进行移动或旋转,不改变它自身的形状,因此它是一种欧式变换。一般而言,欧式变换由旋转和平移构成,假设单位正交基(e1,e2,e3)经过了一次旋转变成(,,),那么对于向量a而言,它在两个坐标系下(同一个点参考系不同)的坐标分别为(a1,a2,a3)T和(,,)T。由于a本身没有改变,因此:
式中左右两边同时左乘,则:
如果将式中3×3的矩阵写作R,则有a=Ra′,它描述了同一个向量坐标的变换关系,也就是旋转变换,因此称R为旋转矩阵(Rotation Matrix)。
旋转矩阵是一个行列式为1的正交矩阵,则R-1=RT,证明如下:
对于式(2-2)同时左乘,则有:
根据式(2-3)和式(2-4),二者分别表示为a=R1a′和a′=R2a,则。如果R1正交,则,也就是证明。根据式(2-3)可得:
式中,,以此类推可得,因此旋转矩阵是一个行列式为1的正交矩阵。将n维旋转矩阵的集合定义为:
式中SO(n)表示特殊正交群,旋转矩阵即描述了两个坐标的变换关系:
如果考虑平移,则a=Ra′可以更新表示为a=Ra′+t。因此两个坐标系的刚体运动可以由R和t进行描述。但是由于该变换不是一个线性变换,因此引入齐次坐标,对变换进行改写:
将原始的三维向量末尾增加1个维度,取值为1,该四维向量被称为齐次坐标(Homogene ous Coordinates),而矩阵写作T,称为变换矩阵(Transform Matrix),这种类型的矩阵又被称为特殊欧式群:
式中,R为3×3的矩阵,t为3×1的列向量,0T为1×3的行向量。同时T的逆表示进行反向的变换:
2.1.2 什么是轴角
旋转矩阵有9个变量,但是一次旋转只有3个自由度并且旋转矩阵自身存在正交约束和行列式值约束从而导致旋转矩阵表示旋转存在冗余变量,Rodrigues(罗德里格斯)最早提出使用一个三维向量来表示三维旋转变换,其方向与旋转轴一致,向量的模等于旋转角度,通常称这种向量为旋转向量。如果是用四个元素进行描述,即三个元素描述旋转轴,另外一个元素描述旋转角,则称这种形式的描述为轴角(Axis Angle),轴角共有三个量且无约束,因此具有三个自由度。表示为:
式中,n与旋转轴方向一致的单位向量,θ为角度。从轴角变换到旋转矩阵可以使用罗德里格斯公式进行转换:
式中,n^表示向量n=[nx,ny,nz]T的反对称矩阵,即:
使用轴角表示方式无法表达两次连续的旋转,由于两次旋转的旋转轴会不一样,因此旋转角度不能直接相加。当旋转角度为0°或者180°时,罗德里格斯公式失效,此时的旋转矩阵为0。可以理解为这种情况下旋转角具有无数多种情况。由于旋转矩阵到轴角表示法的映射关系是一对多的关系,因此存在多种轴角组合对应一种旋转矩阵的情况。
2.1.3 什么是欧拉角
旋转矩阵和轴角并不是最直观的表示形式,而欧拉角(Euler Angle)为描述刚体在三维欧式空间的取向提供了一种非常直观的表示方式。欧拉角有两种分类方法,第一种是按照旋转的轴的顺序,一共12种。三个轴只用两个的:Proper Euler angles(Z-X-Z,X-Y-X,Y-Z-Y,Z-Y-Z,X-Z-X,Y-X-Y)。三个轴全都用的:Tait-Bryan angles(X-Y-Z,Y-Z-X,Z-X-Y,X-Z-Y,Z-Y-X,Y-X-Z)。第二种是按照绕着不动的轴(初始的世界坐标系),还是按照转动后的坐标轴(一直在转动的本体坐标系)来旋转。也就是固定轴旋转和运动轴旋转两大类。因此一共有(Intrinsic rotations+Extrinsic rotations)×(Proper Euler angles+Tait-Bryan angles)=24种。不同领域采用的方式有所差异,但概念类似,比如经典力学中使用ZXZ,量子力学使用的是ZYZ,航空航天使用ZYX/ZXY。所以在跨行业或者跨模块协作的时候,一定要问清楚对方是哪一种欧拉角,本文只用一种来说明问题。从物体的初始状态(一般会选择和世界坐标系重合作为最初状态)绕着自身坐标系的XYZ三个轴进行旋转三个角度,来表示物体的朝向。欧拉角一般可以定义为静态和动态两种形式。首先定义绕各个旋转轴旋转相应角度时所对应的旋转矩阵。
当绕X轴旋转α角度时对应的旋转矩阵为:
当绕Y轴旋转β角度时对应的旋转矩阵为:
当绕Z轴旋转γ角度时对应的旋转矩阵为:
所谓静态的欧拉角指绕世界坐标系三个固定轴的旋转,旋转过程中世界坐标系的三个轴不发生变化,也称为外旋旋转方式。如图2-2所示,按照外旋方式,X-Y-Z旋转顺序(指先绕固定轴X旋转α,再绕固定轴Y旋转β,最后绕固定轴Z旋转γ),这种形式的欧拉角也称为RPY角,α、β和γ对应航空领域的Roll、Pitch和Yaw。最终可得旋转矩阵R1:
图2-2 欧拉角按照固定坐标轴旋转
而动态欧拉角是指绕物体坐标系三个轴的旋转,由于物体旋转过程中坐标轴随着旋转变换运动,也称为内旋旋转方式。如图2-3所示,此时若按照Z-Y-X的顺序依次旋转γ、β和α的角度,这种形式的欧拉角也称为Z-Y-X欧拉角,则最终可得旋转矩阵R2:
图2-3 欧拉角按照自身坐标轴旋转
这里要注意:内旋时绕自身坐标系旋转,R2右乘坐标向量,坐标系(基底)在变换,列变换;外旋时绕固定坐标系旋转,R1左乘坐标向量,坐标(向量)在变换,行变换。这种情况下如果使用RPY角进行描述,即绕物体的Z轴旋转γ,得到偏航角Yaw,绕旋转之后的Y轴旋转β,得到俯仰角Pitch,绕旋转之后的X轴旋转α,得到滚转角Roll。因此Yaw、Pitch和Roll分别对应γ、β和α。
由式(2-17)和式(2-18)可知,ZYX顺序的内旋等价于XYZ顺序的外旋,即:
欧拉角的一个显著缺点就是会碰到著名的万向锁(或万向节死锁)(Gimbal Lock)问题,比如对于外旋RPY角,当俯仰角Pitch为±90°时,第一次旋转与第三次旋转使用同一个轴,使得系统丢失了一个自由度(由3次旋转变成了2次旋转),如图2-4所示。这被称为一种奇异性问题。
图2-4 万向锁示意(第3次旋转变成与第1次相同)
2.1.4 什么是四元数
四元数最早于1843年由哈密顿(William Rowan Hamilton)发明,作为复数(Complex Numbers)的扩展。直到1985年才由Shoemake把四元数引入到计算机图形学中。四元数在一些方面优于欧拉角、轴角和旋转矩阵,比如解决万向锁问题。由于仅需存储4个浮点数,相比矩阵更加轻量,使得四元数解决求逆、串联(多个变换的叠加变换)等操作,相比矩阵更加高效,所以综合考虑,现在主流游戏或动画引擎都会以缩放向量加上旋转四元数和平移向量的形式进行存储角色的运动数据。任意一个在三维空间中的朝向都可以表示为一个绕某个特定轴的旋转。给定旋转轴及旋转角度,很容易把其他形式的旋转表示转化为四元数或者从四元数转化为其他形式。四元数具有紧凑性和非奇异性,包括3个虚部,1个实部,通常表示形式为:
式中,q0、q1、q2、q3均为实数,i、j、k为四元数的三个虚单位,它们之间有如下关系:
如果使用一个标量和一个向量来表示四元数,则:
式中s为四元数的实部,v为四元数的虚部。如果实部为0,则称为纯四元数。而单位四元数满足四元数的模为1,即:
假设有两个四元数qa,qb,它们的向量表示形式为[sa,va]T,[sb,vb]T,即:
则四元数的一些运算及其性质如下。
加减法:
乘法:
式(2-26)的简洁表达形式为:
共轭:
模:
逆:
数乘:
点乘:
任意单位四元数描述了一个转轴和绕该转轴的旋转角度,也就是说,任意单位向量v沿着该单位向量定义的旋转轴u旋转θ度之后的v′可以用四元数乘法表示。令两个四元数v=[0,v],,那么:
式中v′表示绕旋转轴u旋转θ后的v所构成的四元数,v′中实部为0,虚部为罗德里格斯公式的结果,即:
2.1.5 其他空间变换
已知欧式变换的变换矩阵为。除了欧式变换以外,三维空间还存在其他变换方式,包括相似变换、仿射变换和射影变换。
1)相似变换是等距变换和均匀缩放的一个复合,其矩阵表达为:
式中,S为缩放因子,,R为一个正交矩阵。
2)仿射变换的矩阵表达为:
式中,仿射变换的A是一个可逆矩阵,仿射变换是旋转、缩放、平移和错切这四种变换的组合表示。
3)射影变换的矩阵形式为:
射影变换是一种一般变换,射影变换可以理解为把理想点(平行直线在无穷远处的交点)变换到图像上,式中的aT表示缩放。