三维空间刚体运动4-5:四元数多点离散数值解插值方法:Sping
三维空间刚体运动4-5:四元数多点离散数值解插值方法:Sping
- 1. 正切曲率κ(γ,t)\kappa(\gamma, t)κ(γ,t)在H1H_{1}H1上的离散数值解——Sping
- 1.1 离散版解决方案
- 1.2 梯度下降法求解
- 1.3 Sping算法步骤
- 2. 优化策略——计算Sping插值范数的替代方法及算法扩展
- 2.1 计算Sping插值范数的三种替代方法
- 2.2 角速度惩罚项
- 3. 实践Sping算法时的改进
- 3.1 多步放置多次迭代
- 3.2 简化参数宽度——l(qi)l(q_{i})l(qi)
- 3.3 对关键帧曲率加权
- 3.4 改进后的图像
- 3.5 其它细节说明
- 4. Squad与Sping对比
- 4.1 例1:准圆
- 4.2 例2:柔和曲线
- 4.3 例3:带尖点的插值曲线
- 4.4 例4:钟摆
- 4.5 例5:扰动摆
- 4.6 例6:全局属性
- 4.7 结论
- 5. 插值曲线总结
序:本篇系列文章参照高翔老师《视觉SLAM十四讲从理论到实践》,讲解三维空间刚体运动,为读者打下坚实的数学基础。博文将原第三讲分为五部分来讲解,其中四元数部分较多较复杂,又分为四部分。如果读者急于实践,可直接阅读第五部分的机器人运动轨迹,此部分详细讲解了安装准备工作。此系列总体目录如下:
- 旋转矩阵和变换矩阵;
- 旋转向量表示旋转;
- 欧拉角表示旋转;
- 四元数包括以下部分:
4-1. 四元数表示变换;
4-2. 四元数线性插值方法:LinEuler/LinMat/Lerp/Nlerp/Slerp;
4-3. 四元数多点插值方法:Squad;
4-4. 四元数多点连续解析解插值方法:Spicv;
4-5. 四元数多点离散数值解插值方法:Sping。 - 实践:SLAM中显示机器人运动轨迹及相机位姿。
摘要:四元数插值四个方法Slerp、Squad、Spicv和Sping既复杂又很重要,为了详细阐述,故每个方法独立成一篇博文讲解。没有插值的四元数是没有灵魂的,插值的重要性不言而喻。Slerp是经典的两点间一阶连续可导插值方法,Squad方法在Slerp的基础上实现多点间的一阶连续可导,Spicv是多点间连续解析解插值方法,而Sping则是多点间离散数值解插值方法,更适合复杂曲线,此博文也是国内首篇介绍Spicv和Sping的中文资料。
1. 正切曲率κ(γ,t)\kappa(\gamma, t)κ(γ,t)在H1H_{1}H1上的离散数值解——Sping
由《四元数多点连续解析解插值方法:Spicv》可知:虽然连续解析解过程推导严密,但遗憾的是,我们不能给出最优插值曲线的解析表达式。因此,我们将尝试提出一个用数值方法来解决问题的离散数值解版本。
首先我们把这个问题写成等价的离散形式,然后尝试梯度下降解决这个新版本的问题。由于插值算法使用数值解并用到梯度,所以命名为Sping(Spherical Interpolation using Nmerical Gradient descent)。
注释:原论文将此算法命名为Spring,但这个名字太普遍,容易混淆,所以笔者将其精简为Sping,更容易记忆,也减少重名。
1.1 离散版解决方案
现在回顾一下曲率平方的积分K(γ)K(\gamma)K(γ)的概念:给定控制点Q1,...,Qk∈H1Q_{1},...,Q_{k}\in H_{1}Q1,...,Qk∈H1,求γ(t)∈C2(I,H1)\gamma(t)\in C^{2}(I,H_{1})γ(t)∈C2(I,H1),使存在t1,...,tk∈(I)t_{1},...,t_{k}\in (I)t1,...,tk∈(I),满足γ(ti)=Qi\gamma(t_{i})=Q_{i}γ(ti)=Qi,并使下列表达式最小化:K(γ)=∫t1tN∥κ(γ,t)∥2dt(1.1)K(\gamma)=\int_{t_{1}}^{t_{N}}\left \| \kappa (\gamma,t) \right \|^{2}dt\tag{1.1}K(γ)=∫t1tN∥κ(γ,t)∥2dt(1.1)据此在离散版本中进行改写,我们将尝试解决以下问题:
给定控制点Q1,...,Qk∈H1Q_{1},...,Q_{k}\in H_{1}Q1,...,Qk∈H1,寻找q1,...,qN∈H1(N≥k)q_{1},...,q_{N}\in H_{1}(N\geq k)q1,...,qN∈H1(N≥k),对于t=1,...,kt=1,...,kt=1,...,k有qit=Qtq_{i_{t}}=Q_{t}qit=Qt,并使下列表达式最小化:E=∑i=1Nl(qi)∥κ~(qi)∥2(1.2)E=\sum_{i=1}^{N}l(q_{i})\left \| \tilde{\kappa}(q_{i}) \right \|^{2}\tag{1.2}E=i=1∑Nl(qi)∥κ~(qi)∥2(1.2)可见积分被求和代替了,我们积分区间的参数宽度称为l(qi)l(q_{i})l(qi),它可以表示为第iii个四元数前后间隔内参数宽度的中值:l(qi)=∥qi−qi−1∥+∥qi−qi+1∥2(1.3)l(q_{i})=\frac{\left \| q_{i}-q_{i-1} \right \|+\left \| q_{i}-q_{i+1} \right \|}{2}\tag{1.3}l(qi)=2∥qi−qi−1∥+∥qi−qi+1∥(1.3)在l(qi)l(q_{i})l(qi)的近似中,qiq_{i}qi和qi−1q_{i-1}qi−1之间的参数化距离的另一种度量方法是角度θi\theta _{i}θi,它须满足cosθi=qi⋅qi+1\cos \theta _{i}=q_{i}\cdot q_{i+1}cosθi=qi⋅qi+1。
方程(1.2)(1.2)(1.2)中另一个参数κ~\tilde{\kappa}κ~,它表示局部曲率的离散版:k~(qi)=qi′′−qi′′⋅qiqi⋅qiqi(1.4)\tilde{k}(q_{i})=q_{i}^{''}-\frac{q_{i}^{''}\cdot q_{i}}{q_{i}\cdot q_{i}}q_{i}\tag{1.4}k~(qi)=qi′′−qi⋅qiqi′′⋅qiqi(1.4)请注意,在定义局部曲率时,分母qi⋅qiq_{i}\cdot q_{i}qi⋅qi没有被忽略(如上篇中方程(3.1)(3.1)(3.1))。这是因为插值的四元数不一定为单位四元数,此部分请参考第2.1节。因此分母的值一般不等于1,所以它不能被省略。
在方程(1.4)(1.4)(1.4)中,使用的是插值曲线的离散近似的二阶导数。二阶导数的一个很好的中心近似是1:qi′′=qi−1−2qi+qi+1l(qi)2(1.5)q_{i}^{''}=\frac{q_{i-1}-2q_{i}+q_{i+1}}{l(q_{i})^{2}}\tag{1.5}qi′′=l(qi)2qi−1−2qi+qi+1(1.5)
1.2 梯度下降法求解
本节尝试求解方程(1.2)(1.2)(1.2),它可以用梯度下降法最小化。一般来说,梯度下降法可以描述如下:将被最小化的函数图像(通常称为能量函数)视为有坡度的丘陵,其中的函数值是每组坐标上的丘陵高度,函数在某点的梯度表示它最陡的上坡方向。梯度下降是基于对解的初始估计,从最初的估计开始计算梯度,并在梯度的相反方向上迈出一小步(即下山),从而产生一个新的点。这个过程在新的点上不断重复,直到函数值不再通过走一步变得更小,在这个过程中,梯度下降法产生一个近似的局部最小值。梯度下降理论的思想大致如此,不再做详细解释,但是我们将对推导方法进行足够详细的描述,以便那些在该领域没有前置知识的读者仍然能理解其推导过程。
当使用如梯度下降法的数值近似方法时,通常需要对解空间的进行限制;相反,如果解位于期望解空间之外,则会添加一个使解付出更大代价的附加项。据此我们寻找一个函数,它可以确定离散版本的插值曲线是否在H1H_{1}H1中,因此此函数是由单位四元数组成的,可通过下式实现:g(q)=q⋅q−1(1.6)g(q)=q\cdot q-1\tag{1.6}g(q)=q⋅q−1(1.6)其中:H1={q∈H∣g(q)=0}H_{1}= \{ q\in H|g(q)=0 \}H1={q∈H∣g(q)=0},即当g(q)=0g(q)=0g(q)=0时,qqq是单位四元数。确定四元数是否为单位四元数的方法,可以与能量函数EEE结合成新的能量函数FFF:F=∑i=1Nl(qi)∥κ(qi)~∥2+cg(qi)2(1.7)F=\sum_{i=1}^{N}l(q_{i})\left \| \tilde{\kappa(q_{i})} \right \|^{2}+cg(q_{i})^{2}\tag{1.7}F=i=1∑Nl(qi)∥∥∥κ(qi)~∥∥∥2+cg(qi)2(1.7)假设c∈Rc\in \mathbb{R}c∈R为相匹配的实数,当EEE在某处取最小值时,能量函数FFF也近似有一个最小值,在这里所有的qiq_{i}qi为近似单位四元数。
我们想用梯度下降法找到方程(1.7)(1.7)(1.7)中FFF的最小值,因此我们必须求出梯度。下面,我们将qi,xq_{i,x}qi,x记为离散插值曲线第iii个四元数中的第xxx坐标,其中x∈1,2,3,4x\in{1,2,3,4}x∈1,2,3,4,将四元数视为四维向量,因此梯度是4N4N4N维的。FFF关于每个坐标的偏导可以写成:∂F∂qi,x=∂∂qi,x(∑j=1Nl(qj)∥κ~(qj)∥2+cg(qj)2)\frac{\partial F}{\partial q_{i,x}}=\frac{\partial }{\partial q_{i,x}}\left ( \sum_{j=1}^{N} l(q_{j}) \left \| \tilde{\kappa}(q_{j}) \right \|^{2} + cg(q_{j})^{2} \right )∂qi,x∂F=∂qi,x∂(j=1∑Nl(qj)∥κ~(qj)∥2+cg(qj)2)在方程(1.3),(1.4),(1.5)(1.3),(1.4),(1.5)(1.3),(1.4),(1.5)中,qiq_{i}qi存在于算法中的l(qi−1),l(qi),l(qi+1),κ~(qi−1),κ~(qi),κ~(qi+1)l(q_{i-1}),l(q_{i}),l(q_{i+1}),\tilde{\kappa}(q_{i-1}),\tilde{\kappa}(q_{i}),\tilde{\kappa}(q_{i+1})l(qi−1),l(qi),l(qi+1),κ~(qi−1),κ~(qi),κ~(qi+1),因此相应的项必然出现在偏导中,因此有:∂F∂qi,x=∂∂qi,x(l(qi−1)∥κ~(qi−1)∥2+l(qi)∥κ~(qi)∥2+l(qi+1)∥κ~(qi+1)∥2+cg(qi)2)=∂∂qi,x(l(qi−1)κ~(qi−1)⋅κ~(qi−1)+l(qi)κ~(qi)⋅κ~(qi)+l(qi+1)κ~(qi+1)⋅κ~(qi+1)+cg(qi)2)=∂l(qi−1)∂qi,xκ~(qi−1)⋅κ~(qi−1)+∂l(qi)∂qi,xκ~(qi)⋅κ~(qi)+∂l(qi+1)∂qi,xκ~(qi+1)⋅κ~(qi+1)+2l(qi−1)κ~(qi−1)∂κ~(qi−1)∂qi,x+2l(qi)κ~(qi)∂κ~(qi)∂qi,x+2l(qi+1)κ~(qi+1)∂κ~(qi+1)∂qi,x+2cg(qi)∂g(qi)∂qi,x(1.8)\begin{aligned} \frac{\partial F}{\partial q_{i,x}} &= \frac{\partial }{\partial q_{i,x}}\left ( l(q_{i-1})\left \| \tilde{\kappa}(q_{i-1}) \right \|^{2} + l(q_{i})\left \| \tilde{\kappa}(q_{i}) \right \|^{2} + l(q_{i+1})\left \| \tilde{\kappa}(q_{i+1}) \right \|^{2} + cg(q_{i})^{2} \right ) \\ &= \frac{\partial }{\partial q_{i,x}}\left ( l(q_{i-1})\tilde{\kappa}(q_{i-1})\cdot \tilde{\kappa}(q_{i-1}) + l(q_{i})\tilde{\kappa}(q_{i})\cdot \tilde{\kappa}(q_{i}) + l(q_{i+1})\tilde{\kappa}(q_{i+1})\cdot \tilde{\kappa}(q_{i+1}) + cg(q_{i})^{2} \right ) \\ &= \frac{\partial l(q_{i-1})}{\partial q_{i,x}}\tilde{\kappa}(q_{i-1})\cdot \tilde{\kappa}(q_{i-1}) + \frac{\partial l(q_{i})}{\partial q_{i,x}}\tilde{\kappa}(q_{i})\cdot \tilde{\kappa}(q_{i}) + \frac{\partial l(q_{i+1})}{\partial q_{i,x}}\tilde{\kappa}(q_{i+1})\cdot \tilde{\kappa}(q_{i+1}) + \\& \quad \quad 2l(q_{i-1})\tilde{\kappa}(q_{i-1})\frac{\partial \tilde{\kappa}(q_{i-1})}{\partial q_{i,x}} + 2l(q_{i})\tilde{\kappa}(q_{i})\frac{\partial \tilde{\kappa}(q_{i})}{\partial q_{i,x}} + 2l(q_{i+1})\tilde{\kappa}(q_{i+1})\frac{\partial \tilde{\kappa}(q_{i+1})}{\partial q_{i,x}} + \\ & \quad \quad 2c g(q_{i})\frac{\partial g(q_{i})}{\partial q_{i,x}} \end{aligned}\tag{1.8}∂qi,x∂F=∂qi,x∂(l(qi−1)∥κ~(qi−1)∥2+l(qi)∥κ~(qi)∥2+l(qi+1)∥κ~(qi+1)∥2+cg(qi)2)=∂qi,x∂(l(qi−1)κ~(qi−1)⋅κ~(qi−1)+l(qi)κ~(qi)⋅κ~(qi)+l(qi+1)κ~(qi+1)⋅κ~(qi+1)+cg(qi)2)=∂qi,x∂l(qi−1)κ~(qi−1)⋅κ~(qi−1)+∂qi,x∂l(qi)κ~(qi)⋅κ~(qi)+∂qi,x∂l(qi+1)κ~(qi+1)⋅κ~(qi+1)+2l(qi−1)κ~(qi−1)∂qi,x∂κ~(qi−1)+2l(qi)κ~(qi)∂qi,x∂κ~(qi)+2l(qi+1)κ~(qi+1)∂qi,x∂κ~(qi+1)+2cg(qi)∂qi,x∂g(qi)(1.8)下面我们将推导方程(1.8)(1.8)(1.8)中子表达式的偏导数。
我们引入符号:1x1_{x}1x,它表示某个向量在第xxx坐标是1,在其他坐标是0。现在我们可以求出g(qi),l(qi−1),l(qi),l(qi+1),qi−1′′,qi′′,qi+1′′g(q_{i}),l(q_{i-1}),l(q_{i}),l(q_{i+1}),q_{i-1}^{''},q_{i}^{''},q_{i+1}^{''}g(qi),l(qi−1),l(qi),l(qi+1),qi−1′′,qi′′,qi+1′′以及κ~(qi−1),κ~(qi),κ~(qi+1)\tilde{\kappa}(q_{i-1}),\tilde{\kappa}(q_{i}),\tilde{\kappa}(q_{i+1})κ~(qi−1),κ~(qi),κ~(qi+1)的偏导了:∂g(qi)∂qi,x=∂∂qi,x(qi⋅qi−1)=2qi⋅∂qi∂qi,x=2qi,x(1.9)\begin{aligned} \frac{\partial g(q_{i})}{\partial q_{i,x}} &= \frac{\partial }{\partial q_{i,x}}(q_{i}\cdot q_{i}-1) \\&= 2q_{i}\cdot \frac{\partial q_{i}}{\partial q_{i,x}} = 2q_{i,x} \end{aligned}\tag{1.9}∂qi,x∂g(qi)=∂qi,x∂(qi⋅qi−1)=2qi⋅∂qi,x∂qi=2qi,x(1.9)∂l(qi−1)∂qi,x=∂∂qi,x∥qi−1−qi−2∥+∥qi−1−qi∥2=12∂∂qi,x((qi−1−qi−2)⋅(qi−1−qi−2)+(qi−1−qi)⋅(qi−1−qi))=12(−1x⋅(qi−1−qi)(qi−1−qi)⋅(qi−1−qi))=1x⋅(qi−qi−1)2∥qi−1−qi∥(1.10)\begin{aligned} \frac{\partial l(q_{i-1})}{\partial q_{i,x}} &= \frac{\partial }{\partial q_{i,x}}\frac{\left \| q_{i-1}-q_{i-2} \right \|+\left \| q_{i-1}-q_{i} \right \|}{2} \\&= \frac{1}{2}\frac{\partial }{\partial q_{i,x}}\left ( \sqrt{(q_{i-1}-q_{i-2})\cdot (q_{i-1}-q_{i-2})} + \sqrt{(q_{i-1}-q_{i})\cdot(q_{i-1}-q_{i})} \right ) \\&= \frac{1}{2}\left ( -\frac{1_{x}\cdot(q_{i-1}-q_{i})}{\sqrt{(q_{i-1}-q_{i})\cdot(q_{i-1}-q_{i})}} \right ) \\&= \frac{1_{x}\cdot(q_{i}-q_{i-1})}{2\left \| q_{i-1}-q_{i} \right \|} \end{aligned}\tag{1.10}∂qi,x∂l(qi−1)=∂qi,x∂2∥qi−1−qi−2∥+∥qi−1−qi∥=21∂qi,x∂((qi−1−qi−2)⋅(qi−1−qi−2)+(qi−1−qi)⋅(qi−1−qi))=21(−(qi−1−qi)⋅(qi−1−qi)1x⋅(qi−1−qi))=2∥qi−1−qi∥1x⋅(qi−qi−1)(1.10)∂l(qi)∂qi,x=∂∂qi,x∥qi−qi−1∥+∥qi−qi+1∥2=12∂∂qi,x((qi−qi−1)⋅(qi−qi−1)+(qi−qi+1)⋅(qi−qi+1))=12(1x⋅(qi−qi−1)(qi−qi−1)⋅(qi−qi−1)+1x⋅(qi−qi+1)(qi−qi+1)⋅(qi−qi+1))=1x⋅(qi−qi−1)2∥qi−qi−1∥+1x⋅(qi−qi+1)2∥qi+1−qi∥(1.11)\begin{aligned} \frac{\partial l(q_{i})}{\partial q_{i,x}} &= \frac{\partial }{\partial q_{i,x}}\frac{\left \| q_{i}-q_{i-1} \right \|+\left \| q_{i}-q_{i+1} \right \|}{2} \\&= \frac{1}{2}\frac{\partial }{\partial q_{i,x}}\left ( \sqrt{(q_{i}-q_{i-1})\cdot (q_{i}-q_{i-1})} + \sqrt{(q_{i}-q_{i+1})\cdot(q_{i}-q_{i+1})} \right ) \\&= \frac{1}{2}\left ( \frac{1_{x}\cdot(q_{i}-q_{i-1})}{\sqrt{(q_{i}-q_{i-1})\cdot(q_{i}-q_{i-1})}} + \frac{1_{x}\cdot(q_{i}-q_{i+1})}{\sqrt{(q_{i}-q_{i+1})\cdot(q_{i}-q_{i+1})}} \right ) \\&= \frac{1_{x}\cdot(q_{i}-q_{i-1})}{2\left \| q_{i}-q_{i-1} \right \|} + \frac{1_{x}\cdot(q_{i}-q_{i+1})}{2\left \| q_{i+1}-q_{i} \right \|} \end{aligned}\tag{1.11}∂qi,x∂l(qi)=∂qi,x∂2∥qi−qi−1∥+∥qi−qi+1∥=21∂qi,x∂((qi−qi−1)⋅(qi−qi−1)+(qi−qi+1)⋅(qi−qi+1))=21((qi−qi−1)⋅(qi−qi−1)1x⋅(qi−qi−1)+(qi−qi+1)⋅(qi−qi+1)1x⋅(qi−qi+1))=2∥qi−qi−1∥1x⋅(qi−qi−1)+2∥qi+1−qi∥1x⋅(qi−qi+1)(1.11)∂l(qi+1)∂qi,x=∂∂qi,x∥qi+1−qi∥+∥qi+1−qi+2∥2=12∂∂qi,x((qi+1−qi)⋅(qi+1−qi)+(qi+1−qi+2)⋅(qi+1−qi+2))=12(−1x⋅(qi+1−qi)(qi+1−qi)⋅(qi+1−qi))=1x⋅(qi−qi+1)2∥qi+1−qi∥(1.12)\begin{aligned} \frac{\partial l(q_{i+1})}{\partial q_{i,x}} &= \frac{\partial }{\partial q_{i,x}}\frac{\left \| q_{i+1}-q_{i} \right \|+\left \| q_{i+1}-q_{i+2} \right \|}{2} \\&= \frac{1}{2}\frac{\partial }{\partial q_{i,x}}\left ( \sqrt{(q_{i+1}-q_{i})\cdot (q_{i+1}-q_{i})} + \sqrt{(q_{i+1}-q_{i+2})\cdot(q_{i+1}-q_{i+2})} \right ) \\&= \frac{1}{2}\left ( -\frac{1_{x}\cdot(q_{i+1}-q_{i})}{\sqrt{(q_{i+1}-q_{i})\cdot(q_{i+1}-q_{i})}} \right ) \\&= \frac{1_{x}\cdot(q_{i}-q_{i+1})}{2\left \| q_{i+1}-q_{i} \right \|} \end{aligned}\tag{1.12}∂qi,x∂l(qi+1)=∂qi,x∂2∥qi+1−qi∥+∥qi+1−qi+2∥=21∂qi,x∂((qi+1−qi)⋅(qi+1−qi)+(qi+1−qi+2)⋅(qi+1−qi+2))=21(−(qi+1−qi)⋅(qi+1−qi)1x⋅(qi+1−qi))=2∥qi+1−qi∥1x⋅(qi−qi+1)(1.12)∂qi−1′′∂qi,x=∂∂qi,xqi−2−2qi−1+qil(qi−1)2=l(qi−1)2∂∂qi,x(qi−2−2qi−1+qi)−(qi−2−2qi−1+qi)∂∂qi,xl(qi−1)2l(qi−1)4=l(qi−1)2∂∂qi,xqi−(qi−2−2qi−1+qi)2l(qi−1)∂∂qi,xl(qi−1)l(qi−1)4=l(qi−1)21x−2(qi−2−2qi−1+qi)l(qi−1)∂∂qi,xl(qi−1)l(qi−1)4(1.13)\begin{aligned} \frac{\partial q_{i-1}^{''}}{\partial q_{i,x}} &= \frac{\partial }{\partial q_{i,x}}\frac{q_{i-2}-2q_{i-1}+q_{i}}{l(q_{i-1})^{2}} \\&= \frac{l(q_{i-1})^{2}\frac{\partial }{\partial q_{i,x}}(q_{i-2}-2q_{i-1}+q_{i})-(q_{i-2}-2q_{i-1}+q_{i})\frac{\partial }{\partial q_{i,x}}l(q_{i-1})^{2}}{l(q_{i-1})^{4}} \\&= \frac{l(q_{i-1})^{2}\frac{\partial }{\partial q_{i,x}}q_{i}-(q_{i-2}-2q_{i-1}+q_{i})2l(q_{i-1})\frac{\partial }{\partial q_{i,x}}l(q_{i-1})}{l(q_{i-1})^{4}} \\&= \frac{l(q_{i-1})^{2}1_{x}-2(q_{i-2}-2q_{i-1}+q_{i})l(q_{i-1})\frac{\partial }{\partial q_{i,x}}l(q_{i-1})}{l(q_{i-1})^{4}} \end{aligned}\tag{1.13}∂qi,x∂qi−1′′=∂qi,x∂l(qi−1)2qi−2−2qi−1+qi=l(qi−1)4l(qi−1)2∂qi,x∂(qi−2−2qi−1+qi)−(qi−2−2qi−1+qi)∂qi,x∂l(qi−1)2=l(qi−1)4l(qi−1)2∂qi,x∂qi−(qi−2−2qi−1+qi)2l(qi−1)∂qi,x∂l(qi−1)=l(qi−1)4l(qi−1)21x−2(qi−2−2qi−1+qi)l(qi−1)∂qi,x∂l(qi−1)(1.13)∂qi′′∂qi,x=∂∂qi,xqi−1−2qi+qi+1l(qi)2=l(qi)2∂∂qi,x(qi−1−2qi+qi+1)−(qi−1−2qi+qi+1)∂∂qi,xl(qi)2l(qi)4=l(qi)2∂∂qi,x(−2qi)−(qi−1−2qi+qi+1)2l(qi)∂∂qi,xl(qi)l(qi)4=−2l(qi)21x+2(qi−1−2qi+qi+1)l(qi)∂∂qi,xl(qi)l(qi)4(1.14)\begin{aligned} \frac{\partial q_{i}^{''}}{\partial q_{i,x}} &= \frac{\partial }{\partial q_{i,x}}\frac{q_{i-1}-2q_{i}+q_{i+1}}{l(q_{i})^{2}} \\&= \frac{l(q_{i})^{2}\frac{\partial }{\partial q_{i,x}}(q_{i-1}-2q_{i}+q_{i+1})-(q_{i-1}-2q_{i}+q_{i+1})\frac{\partial }{\partial q_{i,x}}l(q_{i})^{2}}{l(q_{i})^{4}} \\&= \frac{l(q_{i})^{2}\frac{\partial }{\partial q_{i,x}}(-2q_{i})-(q_{i-1}-2q_{i}+q_{i+1})2l(q_{i})\frac{\partial }{\partial q_{i,x}}l(q_{i})}{l(q_{i})^{4}} \\&= -\frac{2l(q_{i})^{2}1_{x}+2(q_{i-1}-2q_{i}+q_{i+1})l(q_{i})\frac{\partial }{\partial q_{i,x}}l(q_{i})}{l(q_{i})^{4}} \end{aligned}\tag{1.14}∂qi,x∂qi′′=∂qi,x∂l(qi)2qi−1−2qi+qi+1=l(qi)4l(qi)2∂qi,x∂(qi−1−2qi+qi+1)−(qi−1−2qi+qi+1)∂qi,x∂l(qi)2=l(qi)4l(qi)2∂qi,x∂(−2qi)−(qi−1−2qi+qi+1)2l(qi)∂qi,x∂l(qi)=−l(qi)42l(qi)21x+2(qi−1−2qi+qi+1)l(qi)∂qi,x∂l(qi)(1.14)∂qi+1′′∂qi,x=∂∂qi,xqi−2qi+1+qi+2l(qi+1)2=l(qi+1)2∂∂qi,x(qi−2qi+1+qi+2)−(qi−2qi+1+qi+2)∂∂qi,xl(qi+1)2l(qi+1)4=l(qi+1)2∂∂qi,xqi−(qi−2qi+1+qi+2)2l(qi+1)∂∂qi,xl(qi+1)l(qi+1)4=2l(qi+1)21x−2(qi−2qi+1+qi+2)l(qi+1)∂∂qi,xl(qi+1)l(qi+1)4(1.15)\begin{aligned} \frac{\partial q_{i+1}^{''}}{\partial q_{i,x}} &= \frac{\partial }{\partial q_{i,x}}\frac{q_{i}-2q_{i+1}+q_{i+2}}{l(q_{i+1})^{2}} \\&= \frac{l(q_{i+1})^{2}\frac{\partial }{\partial q_{i,x}}(q_{i}-2q_{i+1}+q_{i+2})-(q_{i}-2q_{i+1}+q_{i+2})\frac{\partial }{\partial q_{i,x}}l(q_{i+1})^{2}}{l(q_{i+1})^{4}} \\&= \frac{l(q_{i+1})^{2}\frac{\partial }{\partial q_{i,x}}q_{i}-(q_{i}-2q_{i+1}+q_{i+2})2l(q_{i+1})\frac{\partial }{\partial q_{i,x}}l(q_{i+1})}{l(q_{i+1})^{4}} \\&= \frac{2l(q_{i+1})^{2}1_{x} - 2(q_{i}-2q_{i+1}+q_{i+2})l(q_{i+1})\frac{\partial }{\partial q_{i,x}}l(q_{i+1})}{l(q_{i+1})^{4}} \end{aligned}\tag{1.15}∂qi,x∂qi+1′′=∂qi,x∂l(qi+1)2qi−2qi+1+qi+2=l(qi+1)4l(qi+1)2∂qi,x∂(qi−2qi+1+qi+2)−(qi−2qi+1+qi+2)∂qi,x∂l(qi+1)2=l(qi+1)4l(qi+1)2∂qi,x∂qi−(qi−2qi+1+qi+2)2l(qi+1)∂qi,x∂l(qi+1)=l(qi+1)42l(qi+1)21x−2(qi−2qi+1+qi+2)l(qi+1)∂qi,x∂l(qi+1)(1.15)∂κ~(qi−1)∂qi,x=∂∂qi,x(qi−1′′−qi−1′′⋅qi−1qi−1⋅qi−1qi−1)=∂qi−1′′∂qi,x−∂qi−1′′∂qi,xqi−1qi−1⋅qi−1qi−1(1.16)\begin{aligned} \frac{\partial \tilde{\kappa}(q_{i-1})}{\partial q_{i,x}} &= \frac{\partial }{\partial q_{i,x}} \left ( q_{i-1}^{''}-\frac{q_{i-1}^{''}\cdot q_{i-1}}{q_{i-1}\cdot q_{i-1}}q_{i-1} \right ) \\&= \frac{\partial q_{i-1}^{''}}{\partial q_{i,x}}-\frac{\frac{\partial q_{i-1}^{''}}{\partial q_{i,x}}q_{i-1}}{q_{i-1}\cdot q_{i-1}}q_{i-1} \end{aligned}\tag{1.16}∂qi,x∂κ~(qi−1)=∂qi,x∂(qi−1′′−qi−1⋅qi−1qi−1′′⋅qi−1qi−1)=∂qi,x∂qi−1′′−qi−1⋅qi−1∂qi,x∂qi−1′′qi−1qi−1(1.16)∂κ~(qi)∂qi,x=∂∂qi,x(qi′′−qi′′⋅qiqi⋅qiqi)=∂qi′′∂qi,x−qi′′⋅qiqi⋅qi∂qi∂qi,x−∂∂qi,x(qi′′⋅qiqi⋅qi)qi=∂qi′′∂qi,x−qi′′⋅qiqi⋅qi1x−∂∂qi,x(qi′′⋅qiqi⋅qi)qi=∂qi′′∂qi,x−qi′′⋅qiqi⋅qi1x−(∂qi′′⋅qi∂qi,xqi⋅qi−∂qi⋅qi∂qi,xqi′′⋅qi(qi⋅qi)2)qi=∂qi′′∂qi,x−qi′′⋅qiqi⋅qi1x−((∂qi′′∂qi,x⋅qi+qi′′⋅∂qi∂qi,x)qi⋅qi−2(qi⋅∂qi∂qi,x)qi′′⋅qi(qi⋅qi)2)qi=∂qi′′∂qi,x−qi′′⋅qiqi⋅qi1x−((∂qi′′∂qi,x⋅qi+qi′′⋅1x)qi⋅qi−2(qi⋅1x)qi′′⋅qi(qi⋅qi)2)qi(1.17)\begin{aligned} \frac{\partial \tilde{\kappa}(q_{i})}{\partial q_{i,x}} &= \frac{\partial }{\partial q_{i,x}} \left ( q_{i}^{''}-\frac{q_{i}^{''}\cdot q_{i}}{q_{i}\cdot q_{i}}q_{i} \right ) \\&= \frac{\partial q_{i}^{''}}{\partial q_{i,x}}-\frac{q_{i}^{''}\cdot q_{i}}{q_{i}\cdot q_{i}}\frac{\partial q_{i}}{\partial q_{i,x}}-\frac{\partial }{\partial q_{i,x}}\left ( \frac{q_{i}^{''}\cdot q_{i}}{q_{i}\cdot q_{i}} \right )q_{i} \\&= \frac{\partial q_{i}^{''}}{\partial q_{i,x}}-\frac{q_{i}^{''}\cdot q_{i}}{q_{i}\cdot q_{i}}1_{x}-\frac{\partial }{\partial q_{i,x}}\left ( \frac{q_{i}^{''}\cdot q_{i}}{q_{i}\cdot q_{i}} \right )q_{i} \\&= \frac{\partial q_{i}^{''}}{\partial q_{i,x}}-\frac{q_{i}^{''}\cdot q_{i}}{q_{i}\cdot q_{i}}1_{x}-\left ( \frac{\frac{\partial q_{i}^{''}\cdot q_{i}}{\partial q_{i,x}}q_{i}\cdot q_{i}-\frac{\partial q_{i}\cdot q_{i}}{\partial q_{i,x}}q_{i}^{''}\cdot q_{i}}{(q_{i}\cdot q_{i})^{2}} \right )q_{i} \\&= \frac{\partial q_{i}^{''}}{\partial q_{i,x}}-\frac{q_{i}^{''}\cdot q_{i}}{q_{i}\cdot q_{i}}1_{x}-\left ( \frac{\left ( \frac{ \partial q_{i}^{''}}{\partial q_{i,x}}\cdot q_{i}+q_{i}^{''}\cdot \frac{\partial q_{i}}{\partial q_{i,x}} \right )q_{i}\cdot q_{i} - 2\left ( q_{i}\cdot\frac{\partial q_{i}}{\partial q_{i,x}} \right )q_{i}^{''}\cdot q_{i}}{(q_{i}\cdot q_{i})^{2}} \right )q_{i} \\&= \frac{\partial q_{i}^{''}}{\partial q_{i,x}}-\frac{q_{i}^{''}\cdot q_{i}}{q_{i}\cdot q_{i}}1_{x}-\left ( \frac{\left ( \frac{ \partial q_{i}^{''}}{\partial q_{i,x}}\cdot q_{i}+q_{i}^{''}\cdot 1_{x} \right )q_{i}\cdot q_{i} - 2\left ( q_{i}\cdot 1_{x}\right )q_{i}^{''}\cdot q_{i}}{(q_{i}\cdot q_{i})^{2}} \right )q_{i} \end{aligned}\tag{1.17}∂qi,x∂κ~(qi)=∂qi,x∂(qi′′−qi⋅qiqi′′⋅qiqi)=∂qi,x∂qi′′−qi⋅qiqi′′⋅qi∂qi,x∂qi−∂qi,x∂(qi⋅qiqi′′⋅qi)qi=∂qi,x∂qi′′−qi⋅qiqi′′⋅qi1x−∂qi,x∂(qi⋅qiqi′′⋅qi)qi=∂qi,x∂qi′′−qi⋅qiqi′′⋅qi1x−⎝⎛(qi⋅qi)2∂qi,x∂qi′′⋅qiqi⋅qi−∂qi,x∂qi⋅qiqi′′⋅qi⎠⎞qi=∂qi,x∂qi′′−qi⋅qiqi′′⋅qi1x−⎝⎜⎜⎛(qi⋅qi)2(∂qi,x∂qi′′⋅qi+qi′′⋅∂qi,x∂qi)qi⋅qi−2(qi⋅∂qi,x∂qi)qi′′⋅qi⎠⎟⎟⎞qi=∂qi,x∂qi′′−qi⋅qiqi′′⋅qi1x−⎝⎜⎜⎛(qi⋅qi)2(∂qi,x∂qi′′⋅qi+qi′′⋅1x)qi⋅qi−2(qi⋅1x)qi′′⋅qi⎠⎟⎟⎞qi(1.17)∂κ~(qi+1)∂qi,x=∂∂qi,x(qi+1′′−qi+1′′⋅qi+1qi+1⋅qi+1qi+1)=∂qi+1′′∂qi,x−∂qi+1′′∂qi,xqi+1qi+1⋅qi+1qi+1(1.18)\begin{aligned} \frac{\partial \tilde{\kappa}(q_{i+1})}{\partial q_{i,x}} &= \frac{\partial }{\partial q_{i,x}} \left ( q_{i+1}^{''}-\frac{q_{i+1}^{''}\cdot q_{i+1}}{q_{i+1}\cdot q_{i+1}}q_{i+1} \right ) \\&= \frac{\partial q_{i+1}^{''}}{\partial q_{i,x}}-\frac{\frac{\partial q_{i+1}^{''}}{\partial q_{i,x}}q_{i+1}}{q_{i+1}\cdot q_{i+1}}q_{i+1} \end{aligned}\tag{1.18}∂qi,x∂κ~(qi+1)=∂qi,x∂(qi+1′′−qi+1⋅qi+1qi+1′′⋅qi+1qi+1)=∂qi,x∂qi+1′′−qi+1⋅qi+1∂qi,x∂qi+1′′qi+1qi+1(1.18)
现在我们对梯度中的所有项都有了简单且长的表达式。我们回到梯度方程(1.8)(1.8)(1.8),可看到其所有组成项均已导出,将方程(1.4)(1.4)(1.4)中κ~(qi−1),κ~(qi),κ~(qi+1)\tilde{\kappa}(q_{i-1}),\tilde{\kappa}(q_{i}),\tilde{\kappa}(q_{i+1})κ~(qi−1),κ~(qi),κ~(qi+1)的导出项和偏导数方程(1.16),(1.17),(1.18)(1.16),(1.17),(1.18)(1.16),(1.17),(1.18)中∂κ~(qi−1)∂qi,x,∂κ~(qi)∂qi,x,∂κ~(qi+1)∂qi,x\frac{\partial \tilde{\kappa}(q_{i-1})}{\partial q_{i,x}},\frac{\partial \tilde{\kappa}(q_{i})}{\partial q_{i,x}},\frac{\partial \tilde{\kappa}(q_{i+1})}{\partial q_{i,x}}∂qi,x∂κ~(qi−1),∂qi,x∂κ~(qi),∂qi,x∂κ~(qi+1)的方程代入梯度方程(1.8)(1.8)(1.8),即可显式计算出梯度。
下面我们看算法的具体步骤。
1.3 Sping算法步骤
上面我们推导出了Sping算法中梯度的显式离散表达式,现在我们可以根据梯度下降法给出Sping算法的详细步骤:
- 初始化
给出一个合理的初始估计值q0=(q1,...,qN)\mathbf{q^{0}}=(q_{1},...,q_{N})q0=(q1,...,qN)。容易想到的是使用已介绍过的方法SlerpSlerpSlerp或SquadSquadSquad。初始估计值越好,方法收敛越快。 - 迭代
根据公式(1.8)(1.8)(1.8)的推导,计算出梯度:▽F=((∂F∂q1,1,∂F∂q1,2,∂F∂q1,3,∂F∂q1,4),...,(∂F∂qN,1,∂F∂qN,2,∂F∂qN,3,∂F∂qN,4))\triangledown F=((\frac{\partial F}{\partial q_{1,1}},\frac{\partial F}{\partial q_{1,2}},\frac{\partial F}{\partial q_{1,3}},\frac{\partial F}{\partial q_{1,4}}),...,(\frac{\partial F}{\partial q_{N,1}},\frac{\partial F}{\partial q_{N,2}},\frac{\partial F}{\partial q_{N,3}},\frac{\partial F}{\partial q_{N,4}}))▽F=((∂q1,1∂F,∂q1,2∂F,∂q1,3∂F,∂q1,4∂F),...,(∂qN,1∂F,∂qN,2∂F,∂qN,3∂F,∂qN,4∂F))注:所有关键帧梯度初始估计值均设为零。
对初始估计值:qi+1=qi−e▽F\mathbf{q^{i+1}}=\mathbf{q^{i}}-e\triangledown Fqi+1=qi−e▽F进行迭代,其中eee定义为步长。
或者:在初始估计值上执行迭代:qi+1=qi−e▽F∥▽F∥\mathbf{q^{i+1}}=\mathbf{q^{i}}-e\frac{\triangledown F}{\left \| \triangledown F \right \|}qi+1=qi−e∥▽F∥▽F。该方法可以提供更大的数值稳定性。 - 终止条件
重复第二步,直到达到预定的终止条件。终止条件可以使用迭代次数,总能量FFF,与控制点数量有关的能量FFF,梯度的大小,或者更高级的策略。
这就是算法的基本步骤,下面我们看一下算法中对范数的改进及算法的扩展。
2. 优化策略——计算Sping插值范数的替代方法及算法扩展
2.1 计算Sping插值范数的三种替代方法
如上所述,我们希望插值曲线位于单位四元数的H1H_{1}H1空间,这是通过使用方程(1.7)(1.7)(1.7)中的惩罚函数g(qi)g(q_{i})g(qi)做到的。但是很难确定方程(1.7)(1.7)(1.7)中惩罚函数的权重因子c∈Rc\in \mathbb{R}c∈R 多大才能与EEE“相匹配”。理论上,ccc必须是无穷大,才能保证qi∈H1q_{i}\in H_{1}qi∈H1,有几种方法可以从不同角度解决这个问题,下面分别介绍:
- 投影
在解空间内和解空间外的点之间有一个简单且定义明确的连接,即投影。如前所述,通过原点的直线上的所有四元数都执行相同的旋转。因此,通过对生成的四元数进行归一化,可以简单地将解估计值投影到解空间中。这可以在每次迭代中执行一次投影,也可以在最后一次迭代之后再投影。 - 拉格朗日算子2
惩罚函数g(qi)g(q_{i})g(qi)也可以用拉格朗日算子λi\lambda_{i}λi来引入,方法如下:F=E+∑i=1Nλig(qi)F=E+\sum_{i=1}^{N}\lambda_{i}g(q_{i})F=E+i=1∑Nλig(qi)方程(1.1)(1.1)(1.1)的解是FFF的一个奇异点。然而奇异点并不是最小值,而是一个鞍点。因此梯度下降法对于此组合四元数FFF仍然适用,需要注意的是梯度上升必须在辅助变量λi\lambda_{i}λi上执行。
通过在方程(1.7)(1.7)(1.7)中加入惩罚函数和拉格朗日算法,该方法变得更具有数值鲁棒性,收敛速度更快。该方法的细节可以在附录中的[Platt & Barr, 1988]和[Barr et al.,1992]中找到。 - 极坐标
在上述的求解方法中,我们重新表述了问题,使解空间上的限制被整合到最小化的能量函数中(加入g(qi)g(q_{i})g(qi))。我们还可以重新定义对解空间的限制,降低最小化表达式的复杂度(投影)。这也可以通过用极坐标来表示四元数实现,此时四元数被写为:q=[r,(ρ,θ,ϕ)]q=[r,(\rho,\theta,\phi)]q=[r,(ρ,θ,ϕ)],其中rrr是半径,(ρ,θ,ϕ)(\rho,\theta,\phi)(ρ,θ,ϕ)是三个必要的旋转角度。
当使用这种表示方式时,可以很容易地保持在解空间的限制内H1H_{1}H1上求解。这可以通过不对半径坐标rrr执行梯度下降来实现,使该坐标保持恒定值1,这样就保持了解空间的限制。
从理论上看,我们认为对生成的四元数进行归一化是“作弊”行为,但也不是不可用。在上述表示中,极坐标将确保插值曲线停留在单位四元数的H1H_{1}H1空间中;使用拉格朗日算子将确保一个更鲁棒和更快的收敛算法。为了使算法尽可能简单,我们不会再进一步探索上述可能性,即我们将不再进一步讨论这三种替代方法中的任何一种。
2.2 角速度惩罚项
一般认为,梯度下降法是一种很容易推广的方法。但是对于插值曲线的预期性质,必须可以简单地描述为某些函数的零交叉或恒定状态。例如,我们可能想要确保整个插值曲线的角速度是恒定的,这可以使用另一个惩罚函数得到:w(qi)=∥qi−1−qi∥−∥qi−qi+1∥(2.1)w(q_{i}) =\left \| q_{i-1}-q_{i} \right \| - \left \| q_{i}-q_{i+1} \right \|\tag{2.1}w(qi)=∥qi−1−qi∥−∥qi−qi+1∥(2.1)可以看到,惩罚函数w(qi)w(q_{i})w(qi)随上步长qi−1q_{i-1}qi−1与下步长qi+1q_{i+1}qi+1之差而增加。该惩罚项可以简单地引入到算法中,比如将其引入到方程(1.7)(1.7)(1.7)的原始能量函数中:F=∑i=1Nl(qi)∥κ(qi)~∥2+cgg(qi)2+cww(qi)2(2.2)F=\sum_{i=1}^{N}l(q_{i})\left \| \tilde{\kappa(q_{i})} \right \|^{2}+c_{g}g(q_{i})^{2}+c_{w}w(q_{i})^{2}\tag{2.2}F=i=1∑Nl(qi)∥∥∥κ(qi)~∥∥∥2+cgg(qi)2+cww(qi)2(2.2)因此必须在梯度方程(1.8)(1.8)(1.8)上加上一项。这很容易推导,因为w(qi)=2l(qi)w(q_{i}) = 2l(q_{i})w(qi)=2l(qi),因此可以用方程(1.11)(1.11)(1.11)得到其导数:∂w(qi)∂qi,x=1x⋅(qi−qi−1)∥qi−qi−1∥+1x⋅(qi−qi+1)∥qi+1−qi∥(2.3)\frac{\partial w(q_{i})}{\partial q_{i,x}} = \frac{1_{x}\cdot(q_{i}-q_{i-1})}{\left \| q_{i}-q_{i-1} \right \|} + \frac{1_{x}\cdot(q_{i}-q_{i+1})}{\left \| q_{i+1}-q_{i} \right \|}\tag{2.3}∂qi,x∂w(qi)=∥qi−qi−1∥1x⋅(qi−qi−1)+∥qi+1−qi∥1x⋅(qi−qi+1)(2.3)添加上述惩罚项并不能保证插补曲线速度不变,只能使角速度平缓些。这是由于梯度包含了其他影响各个关键帧之间距离的项。如Spicv中第3.2节所述,局部曲率的定义包括几何曲率(一阶导)和插值曲线的角加速度(二阶导)的大小。在实践中,这意味着在关键帧周围,插值曲线的角速度将会更小(速度也小),而曲线曲率最大(关键帧通常为拐弯处)。
当梯度包含反作用的项时,算法的鲁棒性就会降低。因此,设置cwc_{w}cw“相匹配”并不一定保证插值帧之间的距离大致相同。相反,这可能导致该方法在数值上变得不稳定,产生不想要的结果。
3. 实践Sping算法时的改进
以上对算法的描述都是纯理论的,我们不希望仅仅从理论上分析算法的收敛性和稳定性,这样不能保证该方法在实践中有效。因此,我们将提出一些在实现算法时具体需要考虑的因素和改进措施。
3.1 多步放置多次迭代
在实践中,当初始估计给定很多帧时,该算法表现得很糟糕,因为每一帧只受它相邻帧影响,并且最初只有关键帧被正确定位。因此,插值曲线对关键帧的所有适配性必须通过关键帧和给定帧之间的帧进行传播,介于两者之间的帧越多,系统就需要更多的迭代才能达到能量最小值,这实际上给出了每个关键帧之间可接受的帧数的上限。
注意:每添加一次关键帧/给定帧称为放置一次,对估计值qi+1=qi−e▽F\mathbf{q^{i+1}}=\mathbf{q^{i}}-e\triangledown Fqi+1=qi−e▽F进行一次更新记为一次迭代。
解决这个问题的一个有效方法是将问题简化为多步放置多次迭代,类似m*n次循环。在第一步放置中使用相对较少的帧,这样系统在几次迭代后达到能量最小值。在第一步中计算出的帧在第二步中作为关键帧使用。在第二步放置中添加更多的中间帧,并且在几次迭代之后再次达到能量最小值。这个过程可以重复任意次数,直到达到所需的帧数。
在图3-1中可以看到算法的结果,在第一步中放置了所有帧。即使经过多次迭代,也会给出一个糟糕的近似曲线。如果用Slerp生成传递给优化算法的初始估计帧,那么该结果也与Slerp非常相似。
图3-1 一步放置多次迭代
说明:图3-1是一步放置多次迭代后插值曲线的图像。第一步插值350帧,迭代500次后的显示。
如果使用多步放置多次迭代算法,经过较少的迭代就可以得到更好的插值曲线。在图3.2中可以看到第一步插值帧数很少,这作为第二步的初始配置。最终的插值曲线如图3.3所示。当使用多步算法时,角速度图表现也更好。如果使用一步放置法(见图3.1),速度曲线类似于心电图,而多步算法(见图3.3)的速度曲线稍好一些,但仍然有些不均匀。这是我们算法的一个缺点,所以我们期待未来能出现一个具有更好收敛性更健壮的算法,能产生更好的速度图。
图3-2 多步放置多次迭代 第一步图像
说明:图3-2是多步放置多次迭代的第一步图像。第一步插值22帧,迭代200次后的显示。
总的来说,在多步放置方法中使用了550次迭代,当使用多步方法时,结果明显更好。
图3-3 多步放置多次迭代 最终图像
说明:图3-3是多步放置迭代的最终图像。第二步插值110帧,迭代150次;在第三步和最后一步,插值550帧并迭代100次。
3.2 简化参数宽度——l(qi)l(q_{i})l(qi)
公式(1.3)(1.3)(1.3)给出了插补曲线参数的“参数宽度”的数值近似。实践证明,该表达式对插值曲线的形状没有影响。然而,该表达式使算法在数值上变得不那么稳定。因此,我们使用更简单的表达式:l(qi)=1(3.1)l(q_{i})=1\tag{3.1}l(qi)=1(3.1)
对其求偏导:∂l(qi)∂qi,x=∂l(qi−1)∂qi,x=∂l(qi+1)∂qi,x=0(3.2)\frac{\partial l(q_{i})}{\partial q_{i,x}} = \frac{\partial l(q_{i-1})}{\partial q_{i,x}} = \frac{\partial l(q_{i+1})}{\partial q_{i,x}} = 0\tag{3.2}∂qi,x∂l(qi)=∂qi,x∂l(qi−1)=∂qi,x∂l(qi+1)=0(3.2)梯度的表达式方程(1.8)(1.8)(1.8)也相应简化为:∂F∂qi,x=2κ~(qi−1)∂κ~(qi−1)∂qi,x+2κ~(qi)∂κ~(qi)∂qi,x+2κ~(qi+1)∂κ~(qi+1)∂qi,x+2cg(qi)∂g(qi)∂qi,x(3.3)\frac{\partial F}{\partial q_{i,x}} = 2\tilde{\kappa}(q_{i-1})\frac{\partial \tilde{\kappa}(q_{i-1})}{\partial q_{i,x}} + 2\tilde{\kappa}(q_{i})\frac{\partial \tilde{\kappa}(q_{i})}{\partial q_{i,x}} + 2\tilde{\kappa}(q_{i+1})\frac{\partial \tilde{\kappa}(q_{i+1})}{\partial q_{i,x}} + 2c g(q_{i})\frac{\partial g(q_{i})}{\partial q_{i,x}}\tag{3.3}∂qi,x∂F=2κ~(qi−1)∂qi,x∂κ~(qi−1)+2κ~(qi)∂qi,x∂κ~(qi)+2κ~(qi+1)∂qi,x∂κ~(qi+1)+2cg(qi)∂qi,x∂g(qi)(3.3)
3.3 对关键帧曲率加权
曲率最小化的连续解析解Spicv算法可以确保曲率是最小的,而离散数值方法Sping算法只给出了最小曲率的近似解。而近似解的有效性取决于对其数学组成表达式的近似程度。例如,插值曲线的二阶导数近似表示为:qi′′=qi−1−2qi+qi+1l(qi)2q_{i}^{''}=\frac{q_{i-1}-2q_{i}+q_{i+1}}{l(q_{i})^{2}}qi′′=l(qi)2qi−1−2qi+qi+1上式即为方程(1.5)(1.5)(1.5)。这种近似是插值曲线局部曲率κ(γ,t)\kappa (\gamma,t)κ(γ,t)近似的核心。使用关键帧在特定曲线之间插值时,预测这种近似是否会给数值解带来“缺陷”是很有必要的。
比如尖锐曲线问题。在实际应用中,数值方法对尖锐曲线有些敏感。例如,图3-4中的曲线在关键帧的位置就太“尖锐”了:
图3-4 多步放置多次迭代 尖锐曲线图像
说明:图3-4是对尖锐曲线的关键帧使用多步放置算法的图像,插值200帧。
对曲率的近似不能充分地在关键帧上传播。让我们重新检查梯度(原版公式1.8的简化版,公式3.3):
∂F∂qi,x=2κ~(qi−1)∂κ~(qi−1)∂qi,x+2κ~(qi)∂κ~(qi)∂qi,x+2κ~(qi+1)∂κ~(qi+1)∂qi,x+2cg(qi)∂g(qi)∂qi,x\frac{\partial F}{\partial q_{i,x}} = 2\tilde{\kappa}(q_{i-1})\frac{\partial \tilde{\kappa}(q_{i-1})}{\partial q_{i,x}} + 2\tilde{\kappa}(q_{i})\frac{\partial \tilde{\kappa}(q_{i})}{\partial q_{i,x}} + 2\tilde{\kappa}(q_{i+1})\frac{\partial \tilde{\kappa}(q_{i+1})}{\partial q_{i,x}} + 2c g(q_{i})\frac{\partial g(q_{i})}{\partial q_{i,x}}∂qi,x∂F=2κ~(qi−1)∂qi,x∂κ~(qi−1)+2κ~(qi)∂qi,x∂κ~(qi)+2κ~(qi+1)∂qi,x∂κ~(qi+1)+2cg(qi)∂qi,x∂g(qi)最后一项确保四元数仍然是单位四元数,与曲率无关,它可以忽略掉,这样在确定梯度以及迭代过程中单个帧的移动时,就只留下了三个因素需要考虑:2κ~(qi−1)∂κ~(qi−1)∂qi,x+2κ~(qi)∂κ~(qi)∂qi,x+2κ~(qi+1)∂κ~(qi+1)∂qi,x2\tilde{\kappa}(q_{i-1})\frac{\partial \tilde{\kappa}(q_{i-1})}{\partial q_{i,x}} + 2\tilde{\kappa}(q_{i})\frac{\partial \tilde{\kappa}(q_{i})}{\partial q_{i,x}} + 2\tilde{\kappa}(q_{i+1})\frac{\partial \tilde{\kappa}(q_{i+1})}{\partial q_{i,x}}2κ~(qi−1)∂qi,x∂κ~(qi−1)+2κ~(qi)∂qi,x∂κ~(qi)+2κ~(qi+1)∂qi,x∂κ~(qi+1),这三个元素代表四元数本身以及相邻四元数的曲率变化。
我们通常认为考虑四元数本身的曲率就足够了。然而,如果不考虑相邻帧,算法结果将退化为SlerpSlerpSlerp,SlerpSlerpSlerp中除了关键帧之外的每一帧都是零曲率的,它不会对关键帧中的大曲率进行补偿。
因此,每个关键帧的相邻帧决定了关键帧的曲率是否最小,所以必须“注意”它们相邻帧的曲率,而这是由方程中的项2κ~(qi−1)∂κ~(qi−1)∂qi,x2\tilde{\kappa}(q_{i-1})\frac{\partial \tilde{\kappa}(q_{i-1})}{\partial q_{i,x}}2κ~(qi−1)∂qi,x∂κ~(qi−1)和2κ~(qi+1)∂κ~(qi+1)∂qi,x2\tilde{\kappa}(q_{i+1})\frac{\partial \tilde{\kappa}(q_{i+1})}{\partial q_{i,x}}2κ~(qi+1)∂qi,x∂κ~(qi+1)来保证的。如图3-4所示,但这还不够,关键帧的图像还是太“尖锐”。
因此,插值曲线二阶导数的近似值在这种情况下不是一个好的选择,这个近似太局部化了,应该有更宽的范围。但这不是一篇关于数值计算方法的报告,因此不在深究最佳的近似表达式。
不过,我们所描述的问题也可以通过简单的方法解决:每个关键帧可以简单地“要求”它的相邻帧更加“贴合”,换句话说,可以在梯度的每个表达式上添加一个权重函数,权重值取决于一个给定的帧是否是某个关键帧的邻居。例如,如果qi−1q_{i-1}qi−1是关键帧,则将2κ~(qi−1)∂κ~(qi−1)∂qi,x2\tilde{\kappa}(q_{i-1})\frac{\partial \tilde{\kappa}(q_{i-1})}{\partial q_{i,x}}2κ~(qi−1)∂qi,x∂κ~(qi−1)替换为2ckκ~(qi−1)∂κ~(qi−1)∂qi,x2\mathbf{c_{k}}\tilde{\kappa}(q_{i-1})\frac{\partial \tilde{\kappa}(q_{i-1})}{\partial q_{i,x}}2ckκ~(qi−1)∂qi,x∂κ~(qi−1)。在实际应用中,ckc_{k}ck取值约为1.2∼1.51.2\sim 1.51.2∼1.5,加权函数的效果如图3-5所示:
图3-5 多步放置多次迭代 加权尖锐曲线图像
说明:图3-5是对尖锐曲线的关键帧使用多步放置算法并加权的图像,对其插值200帧,关键帧周围的曲率加权系数为1.2。
然而“通过眼睛”(即2κ~(qi−1)∂κ~(qi−1)∂qi,x2\tilde{\kappa}(q_{i-1})\frac{\partial \tilde{\kappa}(q_{i-1})}{\partial q_{i,x}}2κ~(qi−1)∂qi,x∂κ~(qi−1))在算法中添加一个常数ckc_{k}ck,与我们声明的从客观标准推导插值曲线的目的是不一致的,但并不影响效果。另一种方法是分析不同数值近似法的性质,但这和收敛性和稳定性一样,不在这个项目的讨论范围之内。
3.4 改进后的图像
在图3-6中可以看到使用Sping并改进后的总体效果。结果只比在关键帧上不包含曲率特殊处理的版本的图3-3略好。在下一章中,我们将看到Sping更清晰地展示其价值的例子。
图3-6 Sping效果图
说明:图3-6是Sping效果图,对其插值550帧,关键帧周围的曲率加权系数为1.3。
3.5 其它细节说明
目前为止,我们还没有描述上述算法的终止要求。在实践中,我们选择了最简单的方法:迭代次数。相应地,我们也没有定义算法使用的产生初始估计值的方法。我们选择在每一对关键帧之间使用Slerp,这个起点显然不是最优的,使用Squad可以实现更快的收敛。但此处作为对比,我们的目的是获得比Squad更优秀的插值曲线,所以使用Slerp作为初始估计似乎更公平。当单独使用Sping时,采用Squad作为初始估计显然更好。
算法的其余常数(迭代中的步长eee和惩罚因子ccc)都可以根据鲁棒性和收敛性准则计算出来。正如前面提到的,我们将不再进一步描述这些属性,并将常量视为实现细节,详情请参见博主上传资源中的程序介绍——四元数插值绘图论文原始代码。
4. Squad与Sping对比
在四元数部分,我们讨论了一系列的插值算法。在已知的算法中,最有说服力的是Squad。在本章中,我们将比较Squad和我们自己的算法Sping,比较将以一些具有代表性的例子为基础进行说明。
4.1 例1:准圆
首先,我们检查一下Squad和Sping是否可以产生漂亮的圆形曲线。我们将关键帧放置在一个球面上正方形的直角上,围绕中心关键帧的插值曲线应该近似为一个准圆。图4-1和图4-2显示,这两条曲线都很好地满足了这一要求,相比较而言,Sping的曲线更接近正圆。
图4-1 带有柔和圆角的Squad简单曲线
说明:图4-1是带有柔和圆角的简单曲线,柔和圆角上使用Squad进行插值。插值中共使用550帧。
图4-2 带有柔和圆角的Sping简单曲线
说明:图4-2是带有柔和圆角的简单曲线,柔和圆角中使用Sping进行插值。生成步骤中总共使用了430帧和700次迭代,关键帧的曲率没有使用额外的权重。
4.2 例2:柔和曲线
这个例子对这两种算法都不是真正的挑战。预期的插值曲线是没有尖角的圆形曲线,曲线在交点处应该不成问题。
图4-3和图4-4的插值曲线都很好。然而,即使关键帧的曲率没有添加额外的权重,Sping的曲线依然比Squad的角更圆润,这意味着Sping的曲率更小。此外,Sping的速度图比Squad的速度图更稳定,这意味着Sping的行为更符合预期。
图4-3 带有交点的Squad柔和曲线
说明:图4-3是对一组有交点的柔和曲线进行Squad插值。在插值中总共使用了550帧。
图4-4 带有交点的Sping柔和曲线
说明:图4-4是对一组有交点的柔和曲线进行Sping插值,在插值步骤上共有550帧和700次迭代,关键帧周围的曲率没有增加额外权重。
4.3 例3:带尖点的插值曲线
Squad可以在关键帧处产生非常尖锐的插值曲线。插值曲线Squad(图4-5)呈现出一条尖尖的曲线。然而,Sping能够产生一个很光滑的圆角插值曲线(图4-6)。
另外,Squad的插值曲线有一个尖角并不与Squad可微分的事实相矛盾。在有尖角的关键帧,曲线的速度为零。因此,尽管曲线的几何外观不是直观上可微分的,但函数Squad仍然是可微分的。
图4-5 带有尖点的Squad插值曲线
说明:图4-5是用Squad插值了一个带有尖点的曲线,使用了200帧。
图4-6 带有尖点的Sping插值曲线
说明:图4-6是用Sping插值了一个带有尖点的曲线,总共有200帧和700次迭代分布在三个步骤中,关键帧周围曲率的相对权重为1.3。
4.4 例4:钟摆
我们继续研究曲率较大的曲线。在这个例子中,我们用的是曲率无穷大的曲线——钟摆运动,它通过三个关键帧实现的,其中第一帧和最后一帧是相等的。
期望的插值曲线的行为在视觉上并不明显,因为所有的关键帧都在一条弧线上,所以曲线应该保持在这个弧线上。
图4-7和图4-8显示了这两条曲线的相同行为——钟摆运动。注意:Lerp和Slerp将产生相同的曲线。
图4-7 Squad插值显示钟摆运动
说明:图4-7是Squad插值显示钟摆运动,其插值超过250帧
图4-8 Sping插值显示钟摆运动
说明:图4-8是Sping插值显示钟摆运动,它在三个步骤内插值超过200帧和700次迭代,关键帧周围曲率的相对权重为1.5。
4.5 例5:扰动摆
在钟摆运动中,尽管当关键帧都在同一弧线上时,插值曲线保持在同一弧线上是很合理的,但并不总是这样。在实际的钟摆曲线上,曲线在尖端关键帧处的曲率是无限大的。本期望Sping应该曲率更小,但实际上并不是,这多少有点令人失望。
而对于扰动摆,要理解它有必要记住这个算法:使用梯度下降,每一帧将在每一步轻微移动,以减少曲率。但是靠近中间关键帧的帧应该向哪个方向移动呢?因为曲线是对称的,所以每一帧的梯度都是零,但这并不表示曲率是一个局部极小值,恰恰相反,它是一个局部极大值。
因此,钟摆就是一个使Sping失去作用的例子。我们进一步研究这一点,通过略微扰乱第一帧和最后关键帧,使曲线不再是一个纯粹的钟摆,如图4-9所示,注意只是略微不同。
图4-10显示梯度下降算法Sping不再排斥纯圆弧的固定点,产生了较圆滑的曲线。因此这里我们以最小的扰动为代价,让插值曲线在中间关键帧处做到漂亮圆滑。
图4-9 Squad插值的扰动摆
说明:图4-9中Squad插值的扰动摆,使用250帧进行插值。
图4-10 Sping插值的扰动摆
说明:图4-10是用Sping插值的扰动摆,它在三个步骤中插入200帧和700次迭代,关键帧周围曲率的相对权重为1.5。
4.6 例6:全局属性
最后一个例子展示了Squad和Sping的根本区别。在每个间隔中,Squad的插值曲线都是由之前的两个关键帧和之后的两个关键帧定义的,比如局部曲率定义。相比之下,Sping的插值曲线是全局定义的。
图4-11显示了在一组共5个关键帧上的插值曲线。前三个关键帧近似地位于一条较平的弧线上,因此插值曲线在第一个间隔内是一条弧线。同样的,插值曲线在最后一个间隔内也形成一条弧线。
图4-11 Squad插值尖曲线
说明:图4-11是Squad使用550帧产生尖曲线。
相比之下,Sping的插值曲线(图4-12)是漂亮且平滑的。该算法的全局结构允许曲线在所有间隔上均匀分布曲率,因此在中间关键帧不存在过度曲率,而是将部分曲率传播到外部间隔。
图4-12 Sping插值尖曲线
说明:图4-12为Sping插值尖曲线,它使用430帧和900次迭代避免了产生尖端。关键帧周围曲率的相对权重为4。
需要注意的是,在这个例子中,有必要在关键帧周围的曲率上增加一个相对较大的权重。然而,毫无疑问,对组成表达式(特别是方程1.5中的qi′′q_{i}^{''}qi′′)有更好的数值近似的算法将产生相同的结果——而不必将程序的参数整合到示例中。
4.7 结论
下面显示了这两种方法的基本区别:
属性 | Squad特性 | Sping特性 |
---|---|---|
算法 | 简单 可解析 连续 区域化 | 复杂 数值化 离散化 全局化 |
插值曲线 | 简单曲线表现良好。关键帧的尖角处具有较大的曲率 | 在所有曲线中表现良好。程序的参数必须符合某些例子。 |
在Squad和Sping之间的做选择的标准并不明显。如果在大多数情况下需要一个简单的算法来产生良好的结果,那么简单的Squad就足够了。如果要求更好的插值算法,在所有情况下,更复杂的Sping将更合适。
5. 插值曲线总结
本篇的目标是客观地探索最优插值方法,从一般求解规则推导出插值曲线。我们首先研究了在文献中发现的或多或少的萌芽式插值曲线,其中包括初级的LinMat, LinEuler, Lerp, 还有较简单的Slerp,最后的结果是令人较满意的Squad。
以这些简单的插值曲线作为基础,我们试图定义插值曲线应该属于哪一类函数(上一篇Spicv的第2.2节),然而我们无法确定哪一种平滑性(CnC^{n}Cn)适合于定义期望的这类函数。
然后我们试图定义一条插值曲线,使插值曲线的局部曲率积分K(γ)K(\gamma)K(γ)最小化。这些推导要求插值曲线具有四阶可导的连续导函数(注意是导函数,不是导数):γ(t)∈C4(I,H1)\gamma(t)\in C^{4}(I,H^{1})γ(t)∈C4(I,H1)。不幸的是,这种推导产生了一个我们无法求解的四阶微分方程。因此,对于Spicv的第2.2节中的开放问题:“曲线应可导多少次,是否允许奇点?”,不再有讨论的必要。
因此,我们推导出了一个离散的数值解,并提出了一种基于梯度下降的方法,我们检查并改进了这种方法。最终的结果是一些非常令人满意的插值曲线。
作为我们最终的插值算法,我们将选择本篇的梯度下降算法,同时在子间隔内使用Squad产生相对应分布的插值帧,并在关键帧周围使用加权曲率,我们将这条插值曲线命名为Sping。
本文基于《视觉SLAM十四讲:从理论到实践》和《Quaternions, Interpolation and Animation》编写,但相对于原文会适当精简,同时为便于全面理解,会收集其他网络好文,根据作者理解,加入一些注解和扩展知识点,如果您觉得还不错,请一键四连(点赞关注收藏评论),让更多的人看到。
参考文献:
- 《视觉SLAM十四讲:从理论到实践》,高翔、张涛等著,中国工信出版社
- Quaternions, Interpolation and Animation
- [Kincaid & Cheney, 1991] David Kincaid & Ward Cheney. Numerical Analysis. Brooks/Cole Publishing Company, Pacific Grove, California, 1991.
- [Barr et al., 1992] Alan H. Barr, Bena Currin, Steven Gabriel, & John F. Hughes. Smooth interpolation of orientations with angular velocity constraints using quaternions. Computer Graphics, 26(2):313-320, July 1992.
- [Platt & Barr, 1988] John C. Platt & Alan H. Barr. Constraint methods for exible models.Computer Graphics, 22(4):279{288, August 1988.
推导请参考[Kincaid & Cheney, 1991], [Barr et al., 1992] ↩︎
关于拉格朗日算子请参考《如何理解拉格朗日乘子法?》。 ↩︎
三维空间刚体运动4-5:四元数多点离散数值解插值方法:Sping相关推荐
- 三维空间刚体运动4-4:四元数多点连续解析解插值方法:Spicv
三维空间刚体运动4-4:四元数多点连续解析解插值方法:Spicv 1. 总述:多点旋转插值的数学方法 2. 插值曲线及其连续性 2.1 插值曲线定义 2.2 插值曲线连续性的讨论 3. 最优插值曲线 ...
- 三维空间刚体运动4-1:四元数表示旋转(各形式相互转换加代码)
三维空间刚体运动4-1:四元数表示变换(各形式相互转换加代码) 1. 四元数的定义 1.1 为什么使用四元数 1.2 复数与四元数 1.3 四元数的形式 2. 四元数的运算 2.1 基础运算 2.2 ...
- 三维空间刚体运动4-3:四元数线性插值方法:Squad
三维空间刚体运动4-3:四元数线性插值方法:Squad Squad的引出 B e ˊ z i e r c u r v e B\acute{e}zier \space curveB e ˊ zier c ...
- 三维空间刚体运动5:详解SLAM中显示机器人运动轨迹及相机位姿(原理流程)
三维空间刚体运动5:详解SLAM中显示机器人运动轨迹及相机位姿(原理流程) 一.显示运动轨迹原理讲解 二.前期准备 三.git管理子模块及克隆源代码 1.学习使用Git Submodule 2.克隆源 ...
- 三维空间刚体运动3:欧拉角表示旋转(全面理解万向锁、RPY角和欧拉角)
三维空间刚体运动3:欧拉角表示旋转(全面理解万向锁.RPY角和欧拉角) 1. 欧拉角 1.1 定义 2.2 RPY角与Z-Y-X欧拉角 2. 欧拉角到旋转矩阵 3. 旋转矩阵到欧拉角 4. 万向锁 4 ...
- 三维空间刚体运动1:旋转矩阵与变换矩阵(详解加代码示例)
三维空间刚体运动1:旋转矩阵与变换矩阵(详解加代码示例) 1. 点.向量和坐标系 2.坐标系间的欧式变换 2.1 旋转 2.2 平移 3.齐次坐标和变换矩阵 4. 相似.仿射和射影变换 4.1 相似变 ...
- 三维空间刚体运动2:旋转向量与罗德里格斯公式(最详细推导)
三维空间刚体运动2:旋转向量与罗德里格斯公式(最详推导) 1.旋转向量定义 2.罗德里格斯公式-向量转换为矩阵 2.1 定义 2.2 推导 2.2.1 推导一 2.2.2 推导二 2.2.3 推导向量 ...
- 视觉SLAM十四讲:第3讲 三维空间刚体运动
第3讲:三维空间刚体运动 三维空间中刚体运动的描述方式:旋转矩阵.变换矩阵.四元数和欧拉角 3.1 旋转矩阵 3.1.1 点和向量,坐标系 三维空间中,给定线性空间基(e1,e2,e3)(\mathb ...
- 高博SLAM十四讲书本程序学习——第3讲 三维空间刚体运动
小白高博SLAM十四讲书本程序学习_1 第3讲 三维空间刚体运动 在高博原始注释上,针对我自己不明白的部分,做额外注释 如果有错误的地方,请大家指点指点 博文目录 一.P.48 eigenMatrix ...
最新文章
- Django介绍工程搭建
- 美国部分Android手机竟将用户隐私数据回传至上海服务器!
- xml 文本转json java_java将XML文档转换成json格式数据
- python文件读read()、readline()、readlines()对比
- 根据数据库表字段删除所有相关信息(删库)
- MySQL高级 - 锁 - MyISAM表锁 - 读锁
- 判断扫码是支付宝还是微信
- Takeown--夺取文件or文件夹所有权
- 数据结构——最小生成树之克鲁斯卡尔算法(Kruskal)
- vue项目中z-index不起作用(将vue实例挂在到window上面)
- 1.php查询数据,数据查询 · thinkphp5 · 看云
- GEE开发之Landsat8计算NDWI和数据分析
- 机器视觉应该先看什么书?
- TCP粘包和拆包问题
- 【渝粤教育】电大中专中药学基础 (2)作业 题库
- 我的2014--菜鸟慢慢在长大
- HTB打靶(Active Directory 101 Forest)
- springboot之shiro
- 2018医学考博英语阅读理解解题技巧
- 计算机二级操作题相关笔记
热门文章
- sob攻略超详细攻略_2020云南旅游超详细必看攻略(附带云南美食景点攻略)
- day27 java的集合(5) HashMap集合和与Hashtable的区别
- python为源文件指定系统默认_Python 设置系统默认编码
- config.class.php,The EventConfig class - PHP 7 中文文档
- ROS入门笔记(十):编写与测试简单的消息发布器和订阅器(C++)
- linux gret 文件内容,DataX插件开发指南.docx
- python 列表 元祖 字典,Python 列表、元组、字典
- uitextfield 键盘类型_以编程方式更改UITextField键盘类型
- 工作分流是什么意思_【嘉陵特装要闻】重庆嘉陵召开持续推进职工分流安置工作布置会...
- Struts 体系结构与工作原理 图