前言
之前对写了一篇关于刚体运动学相关知识博客:刚体运动学——欧拉角、四元数、旋转矩阵,本篇博客就举例来说明,如何在运动捕捉数据中进行四元数插值。
国际惯例,参考博客:
探讨:向量(方向)之间的插值-四元数法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∗w2S的虚部 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.');
end
if any(any(t < 0.0)) || any(any(t > 1.0))
error('Third parameter must have values between 0 and 1 inclusive.');
end
if ~(all(size(q1) == size(q2)) || isscalar(q1) || isscalar(q2))
error(['First two parameters cannot be of different sizes unless' ...
' one is a scalar.']);
end
if ~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.']);
end
end
y = q1 .* (q1.^-1 .* q2).^t;
然后使用此函数尝试在运动捕捉数据中进行插值:
%方法一:matlab自带函数slerp
clear
clc
close all
addpath(genpath('.'))
%读取两个运动数据skel,A,B
load 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:internum
t=i/internum;
%对于角度采用四元数插值
for j=1:size(quatA,1)
temp_quat(j,:)=slerp(quatA(j,:),quatB(j,:),t);
end
temp_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;
end
newMotion(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
%https://gist.github.com/simonlynen/5349167
q1 = 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 quaternions
theta = acos(absD);
sinTheta = sin(theta);
scale0 = sin( ( 1.0 - t ) * theta) / sinTheta;
scale1 = sin( ( t * theta) ) / sinTheta;
end
if(d < 0)
scale1 = -scale1;
end
q3 = scale0 * q1 + scale1 * q2;
q3 = q3 ./ norm(q3);
end
同样使用此算法对运动捕捉数据进行插值:
%第二个插值方法
clear
clc
close all
addpath(genpath('.'))
%读取两个运动数据skel,A,B
load 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:internum
t=i/internum;
%对于角度采用四元数插值
for j=1:size(quatA,1)
temp_quat(j,:)=jointslerp(quatA(j,:)',quatB(j,:)',t);
end
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;
end
newMotion(find(isnan(newMotion)))=0;
skelPlayDataA(skel,[A;newMotion;B])
结果:
后记
其实之前写过类似博客,但是不是用markdown写的,排版真的好丑,我就把它们删掉,写到此博客了。代码连接:链接:https://pan.baidu.com/s/1uLadyPL8yPlQWdPpLSWVrw 密码:asph
代码也可以到我个人的CSDN上传空间去找,或者微信公众号个人简介中的GitHub。此博客已同步更新至微信公众号