前言
之前对写了一篇关于刚体运动学相关知识博客：刚体运动学——欧拉角、四元数、旋转矩阵，本篇博客就举例来说明，如何在运动捕捉数据中进行四元数插值。
国际惯例，参考博客：
探讨：向量（方向）之间的插值-四元数法VS.旋转矩阵法的性能比较
书籍《3D数学基础：图形与游戏开发》
插值理论
问题：3D空间中，在等长度的两个交角为 θ \theta θ的向量 V 1 ( x 1 , y 1 , z 1 ) V_1(x_1,y_1,z_1) V1(x1,y1,z1)和 V 2 ( x 2 , y 2 , z 2 ) V_2(x_2,y_2,z_2) V2(x2,y2,z2)。
实例：行星绕太阳转动，找到旋转过程的两个位置 p 1 , p 2 p_1,p_2 p1,p2，现在模拟从 p 1 p_1 p1到 p 2 p_2 p2的过程。
思路：
1 .一般线性插值
线性插值方法：
这里可以看出，插值的部分就是向量 V 3 V_3 V3.下面来证明 V 3 V_3 V3与 t t t的关系 V 2 − V 1 V_2-V_1 V2−V1得到 V 1 V_1 V1指向 V 2 V_2 V2的向量，再乘以 t t t就是 V 1 V_1 V1指向 V 3 V_3 V3的向量了，最后加上向量 V 1 V_1 V1就是向量 V 3 V_3 V3了，公式为：
v ( t ) = v 1 + t ∗ ( v 2 − v 1 ) （ 0 ≤ t ≤ 1 ) v(t) = v_1 + t*(v_2-v_1)（0\leq t\leq1) v(t)=v1+t∗(v2−v1)（0≤t≤1)
【注：可以看出一般线性插值长度变化了，不满足要求，用球面线性插值就不会变化】
2.一般球面线性插值
将插值结果放大一个放大系数 k ( t ) k(t) k(t)，使其长度放大到 ∣ v 1 ∣ |v_1| ∣v1∣或者 ∣ v 2 ∣ |v_2| ∣v2∣(简单的说就是保持长度不变)。
v ( t ) = k ( t ) ( v 1 + t ( v 2 − v 1 ) ) v(t) = k(t)(v_1 + t(v_2-v_1)) v(t)=k(t)(v1+t(v2−v1))
其中 k ( t ) = ∣ v 1 ∣ ∣ v ( t ) ∣ = ∣ v 1 ∣ ∣ v 1 + t ∗ ( v 2 − v 1 ) ∣ k(t) =\frac{|v_1|}{|v(t)|}=\frac{|v_1|}{|v_1+t*(v_2-v_1)|} k(t)=∣v(t)∣∣v1∣=∣v1+t∗(v2−v1)∣∣v1∣.
这样，插值向量 v ( t ) v(t) v(t)的端点就会沿着 v 1 v_1 v1， v 2 v_2 v2端点构成的圆弧进行， v 1 v_1 v1和 v 2 v_2 v2是等长的，圆弧实际位于 v 1 和 v_1和 v1和 v 2 v_2 v2构成的曲面上的一段，所以又叫球面线性插值。
这个插值解决了3D空间中旋转的插值，在关键帧动画中可以用来计算两个关键帧之间的动画。但是，由于它的插值不是等角速度的，而是变速的，所以如果用来实现案例中的效果的话还需进一步处理。
【注】一般球面线性插值 v ( t ) v(t) v(t)与 v 1 v_1 v1的夹角 θ ( t ) \theta(t) θ(t)不是t的线性函数。
证明过程如下(我滴妈呀，我的字好丑o(╯□╰)o)：
3.改进的球面线性插值，有两种方法：
1> 四元数工具
变换方法：
构造四元数 q ( cos θ ， sin θ ∗ v 1 ’ ) ， r ( cos θ ， sin θ ∗ v 2 ’ ) q(\cos \theta，\sin\theta *v_1’)，r(\cos \theta，\sin \theta *v_2’) q(cosθ，sinθ∗v1’)，r(cosθ，sinθ∗v2’)( v 1 ’ ， v 2 ’ v_1’，v_2’ v1’，v2’为单位 v 1 , v 2 v_1,v_2 v1,v2向量)，以及参数 t ( 0 ≤ t ≤ 1 ) t(0\leq t\leq1) t(0≤t≤1),则构造四元数变换：
四元数 s ( w , v ’ ) = r ∗ ( q − 1 ) t ∗ q s(w,v’)=r*(q-1)t*q s(w,v’)=r∗(q−1)t∗q即为球面线性插值变换。其中，s的虚部 v 1 ’ v_1’ v1’和 v 2 ’ v_2’ v2’间的插值向量，乘以长度 x 2 + y 2 + z 2 \sqrt{x^2+y^2+z^2} x2+y2+z2 ，即得到 v 1 , v 2 v_1,v_2 v1,v2间插值向量 v v v
另一种变换形式是对四元数进行插值变换
s ( w , v ′ ) = a ∗ q + b ∗ r s(w,v')=a*q+b*r s(w,v′)=a∗q+b∗r
其中 a = sin ( α ( 1 − t ) ) sin α , b = s i n ( α t ) sin α , cos α = x 1 ∗ x 2 + y 1 ∗ y 2 + w 1 ∗ w 2 a=\frac{\sin(\alpha(1-t))}{\sin\alpha},b=\frac{\\sin(\alpha t)}{\sin\alpha},\cos\alpha=x_1*x_2+y_1*y_2+w_1*w_2 a=sinαsin(α(1−t)),b=sinαsin(αt),cosα=x1∗x2+y1∗y2+w1∗w2
S的虚部 v ′ v' v′即为 v 1 ′ v_1' v1′和 v 2 ′ v_2' v2′间的插值向量，乘以长度 x 2 + y 2 + z 2 \sqrt{x^2+y^2+z^2} x2+y2+z2 ，即得 v 1 , v 2 v_1,v_2 v1,v2间插值向量 v v v
2> 利用旋转矩阵
变换方法： v = v 1 ∗ T r o t v=v1*Trot v=v1∗Trot
其中， T r o t Trot Trot即绕任意轴旋转的矩阵变换矩阵，因为 v 1 v_1 v1到 v 2 v_2 v2间的插值可以看成是 v 1 v_1 v1绕垂直于 v 1 , v 2 v_1,v_2 v1,v2组成的平面的向量的旋转，所以实际上是绕轴旋转的问题，不过相应参数变成 θ = t ∗ θ \theta=t*\theta θ=t∗θ，轴 q ( q 1 , q 2 , q 3 ) q(q1,q2,q3) q(q1,q2,q3)变成向量 v 1 × v 2 ∣ v 1 × v 2 ∣ = y 1 z 2 − z 1 y 2 , z 1 x 2 − x 1 z 2 , x 1 y 2 − y 1 x 2 sin θ \frac{v_1\times v_2}{|v_1\times v_2|}=\frac{y_1z_2-z_1y_2,z_1x_2-x_1z_2,x_1y_2-y_1x_2}{\sin\theta} ∣v1×v2∣v1×v2=sinθy1z2−z1y2,z1x2−x1z2,x1y2−y1x2
四元数插值
##第一种插值方法
四元数比较重要的一个用途就是球面线性插值(Spherical Linear Interpolation)，可以在两个四元数之间平滑插值。
插值步骤：
① 计算两个值的差： q 0 q_0 q0到 q 1 q_1 q1的角位移由 Δ q = q 0 − 1 q 1 \Delta q=q_0^{-1}q_1 Δq=q0−1q1给出
② 计算差的一部分，四元数求幂可以做到，差的一部分由 Δ q t \Delta q_t Δqt给出
③ 在开始值上加上差的一部分，用四元数乘法组合角位移 q 0 Δ q t q_0\Delta q_t q0Δqt
这样就可以得到slerp公式:
s l e r p ( q 0 , q 1 , t ) = q 0 ( q 0 − 1 q 1 ) t slerp(q_0,q_1,t)=q_0(q_0^{-1}q_1)^t slerp(q0,q1,t)=q0(q0−1q1)t
看看matlab中的函数实现：
function y = slerp(q1, q2, t)% The third parameter, t, gives the 'distance' along the 'arc' between the% quaternions, 0 representing q1 and 1 representing q2. If q1 and q2 are% unit pure quaternions, the interpolation is along a great circle of the% sphere between the points represented by q1 and q2. If q1 and q2 are unit% full quaternions, the interpolation is along the 'arc' on the 4-sphere:% this means the result is a quaternion which represents a rotation% intermediate between the two rotations represented by q1 and q2. If the% first two parameters are not unit quaternions, then there is also% interpolation in modulus.error(nargchk(3, 3, nargin)), error(nargoutchk(0, 1, nargout))if ~isnumeric(t) || ~isreal(t)error('Third parameter must be real and numeric.');endif any(any(t < 0.0)) || any(any(t > 1.0))error('Third parameter must have values between 0 and 1 inclusive.');endif ~(all(size(q1) == size(q2)) || isscalar(q1) || isscalar(q2))error(['First two parameters cannot be of different sizes unless' ...' one is a scalar.']);endif ~isscalar(t) if ~(all(size(q1) == size(t)) || all(size(q2) == size(t)) || ...(isscalar(q1) && isscalar(q2)) ...)error(['Third parameter cannot be an array unless' ...' the first two are scalars, or it has the'...' same size as one of the first two parameters.']);endendy = q1 .* (q1.^-1 .* q2).^t;
然后使用此函数尝试在运动捕捉数据中进行插值：
%方法一：matlab自带函数slerpclearclcclose alladdpath(genpath('.'))%读取两个运动数据skel,A,Bload sample.mat% skelPlayDataA(skel,[A;B])%将欧拉角转换为四元数quatA=joint_euler2quat(skel,A);quatB=joint_euler2quat(skel,B);%执行四元数插值，插20帧internum=20;temp_quat=zeros(31,4);%31个关节，每个关节一个四元数newMotion=zeros(internum,62);%20帧，每帧62维for i=1:internumt=i/internum;%对于角度采用四元数插值for j=1:size(quatA,1)temp_quat(j,:)=slerp(quatA(j,:),quatB(j,:),t); endtemp_quat(find(isnan(temp_quat)))=0;temp_quat=real(temp_quat);newMotion(i,:)=joint_quat2euler(temp_quat);%对于位置采用线性插值posA=A(1,1:3);posB=B(1,1:3);newMotion(i,1:3)=(1-t)*posA+t*posB;endnewMotion(find(isnan(newMotion)))=0;skelPlayDataA(skel,[A;newMotion;B])
结果
第二种插值方法
Slerp的思想就是沿着 4 D 4D 4D球面上连接两个四元数的弧插值。
先看平面上的两个 2 D 2D 2D向量 v 0 v_0 v0和 v 1 v_1 v1都是单位向量，我们需要计算 v t v_t vt它是沿着 v 0 v_0 v0到 v 1 v_1 v1弧的平滑插值。设 w w w是 v 0 v_0 v0到 v t v_t vt弧所截的角，那么 v t v_t vt就是 v 1 v_1 v1沿弧旋转 t w tw tw的结果。
需要考虑两点问题：一是四元数 q q q和 − q -q −q代表同一方位，但是作为slerp的参数时，可能有不一样的结果，是因为 4 D 4D 4D球面不是欧式空间的直接扩展，而这种现象在 2 D 3 D 2D 3D 2D3D空间是不会发生的。解决方法是选择 q 0 q_0 q0和 q 1 q_1 q1的符号使得点乘 q 0 ⋅ q 1 q_0\cdot q_1 q0⋅q1的结果是非负。第二就是如果 q 0 q_0 q0和 q 1 q_1 q1非常接近， sin θ \sin\theta sinθ会非常小，这时除法会出现问题，解决方法是此时采用线性插值。
在论文《从运动捕获数据中提取关键帧》也有介绍到这种四元数插值方法，这里直接贴过来，有兴趣去看看论文：
若 q 1 = [ w 1 , x 1 , y 1 , z 1 ] q_1=[w_1,x_1,y_1,z_1] q1=[w1,x1,y1,z1]和 q 2 = [ w 2 , x 2 , y 2 , z 2 ] q_2=[w_2,x_2,y_2,z_2] q2=[w2,x2,y2,z2]为两个单位四元数，它们之间的球面线性插值为
s l e r p ( q 1 , q 2 ; t ) = sin ( 1 − t ) θ sin θ q 1 + sin t θ sin θ q 2 slerp(q_1,q_2;t)=\frac{\sin(1-t)\theta}{\sin\theta}q_1+\frac{\sin t\theta}{\sin\theta}q_2 slerp(q1,q2;t)=sinθsin(1−t)θq1+sinθsintθq2
其中 θ = arccos ( w 1 w 2 + x 1 x 2 + y 1 y 2 + z 1 z 2 ) \theta=\arccos(w_1w_2+x_1x_2+y_1y_2+z_1z_2) θ=arccos(w1w2+x1x2+y1y2+z1z2)
直接撸代码：
function [ q3 ] = jointslerp( q1, q2, t )%SLERP quaternion slerp% computes the slerp of value t between quaternions q1 and q2%/simonlynen/5349167q1 = q1 ./ norm(q1);q2 = q2 ./ norm(q2);one = 1.0 - eps;d = q1'*q2;absD = abs(d);if(absD >= one)scale0 = 1 - t;scale1 = t;else% theta is the angle between the 2 quaternionstheta = acos(absD);sinTheta = sin(theta);scale0 = sin( ( 1.0 - t ) * theta) / sinTheta;scale1 = sin( ( t * theta) ) / sinTheta;endif(d < 0)scale1 = -scale1;endq3 = scale0 * q1 + scale1 * q2;q3 = q3 ./ norm(q3);end
同样使用此算法对运动捕捉数据进行插值：
%第二个插值方法clearclcclose alladdpath(genpath('.'))%读取两个运动数据skel,A,Bload sample.mat% skelPlayDataA(skel,[A;B])%将欧拉角转换为四元数quatA=joint_euler2quat(skel,A);quatB=joint_euler2quat(skel,B);%执行四元数插值，插20帧internum=20;temp_quat=zeros(31,4);%31个关节，每个关节一个四元数newMotion=zeros(internum,62);%20帧，每帧62维for i=1:internumt=i/internum;%对于角度采用四元数插值for j=1:size(quatA,1)temp_quat(j,:)=jointslerp(quatA(j,:)',quatB(j,:)',t); endnewMotion(i,:)=joint_quat2euler(temp_quat);%对于位置采用线性插值posA=A(1,1:3);posB=B(1,1:3);newMotion(i,1:3)=(1-t)*posA+t*posB;endnewMotion(find(isnan(newMotion)))=0;skelPlayDataA(skel,[A;newMotion;B])
结果：
后记
其实之前写过类似博客，但是不是用markdown写的，排版真的好丑，我就把它们删掉，写到此博客了。代码连接：链接：/s/1uLadyPL8yPlQWdPpLSWVrw 密码：asph
代码也可以到我个人的CSDN上传空间去找，或者微信公众号个人简介中的GitHub。此博客已同步更新至微信公众号
