网站开发找公司好还是个人,阳谷聊城做网站,一个微信可以做两个网站支付,威海教育行业网站建设前言
之前对写了一篇关于刚体运动学相关知识博客#xff1a;刚体运动学——欧拉角、四元数、旋转矩阵#xff0c;本篇博客就举例来说明#xff0c;如何在运动捕捉数据中进行四元数插值。
国际惯例#xff0c;参考博客#xff1a;
探讨#xff1a;向量#xff08;方向…前言
之前对写了一篇关于刚体运动学相关知识博客刚体运动学——欧拉角、四元数、旋转矩阵本篇博客就举例来说明如何在运动捕捉数据中进行四元数插值。
国际惯例参考博客
探讨向量方向之间的插值-四元数法VS.旋转矩阵法的性能比较
书籍《3D数学基础图形与游戏开发》
插值理论
问题3D空间中在等长度的两个交角为θ\thetaθ的向量V1(x1,y1,z1)V_1(x_1,y_1,z_1)V1(x1,y1,z1)和V2(x2,y2,z2)V_2(x_2,y_2,z_2)V2(x2,y2,z2)。
实例行星绕太阳转动找到旋转过程的两个位置p1,p2p_1,p_2p1,p2现在模拟从p1p_1p1到p2p_2p2的过程。
思路
1 .一般线性插值
线性插值方法 这里可以看出插值的部分就是向量V3V_3V3.下面来证明V3V_3V3与ttt的关系V2−V1V_2-V_1V2−V1得到V1V_1V1指向V2V_2V2的向量再乘以ttt就是V1V_1V1指向V3V_3V3的向量了最后加上向量V1V_1V1就是向量V3V_3V3了公式为 v(t)v1t∗(v2−v1)0≤t≤1)v(t) v_1 t*(v_2-v_1)0\leq t\leq1) v(t)v1t∗(v2−v1)0≤t≤1) 【注可以看出一般线性插值长度变化了不满足要求用球面线性插值就不会变化】
2.一般球面线性插值 将插值结果放大一个放大系数k(t)k(t)k(t)使其长度放大到∣v1∣|v_1|∣v1∣或者∣v2∣|v_2|∣v2∣(简单的说就是保持长度不变)。 v(t)k(t)(v1t(v2−v1))v(t) k(t)(v_1 t(v_2-v_1)) v(t)k(t)(v1t(v2−v1)) 其中k(t)∣v1∣∣v(t)∣∣v1∣∣v1t∗(v2−v1)∣k(t) \frac{|v_1|}{|v(t)|}\frac{|v_1|}{|v_1t*(v_2-v_1)|}k(t)∣v(t)∣∣v1∣∣v1t∗(v2−v1)∣∣v1∣.
这样插值向量v(t)v(t)v(t)的端点就会沿着v1v_1v1v2v_2v2端点构成的圆弧进行v1v_1v1和v2v_2v2是等长的圆弧实际位于v1和v_1和v1和v2v_2v2构成的曲面上的一段所以又叫球面线性插值。
这个插值解决了3D空间中旋转的插值在关键帧动画中可以用来计算两个关键帧之间的动画。但是由于它的插值不是等角速度的而是变速的所以如果用来实现案例中的效果的话还需进一步处理。
【注】一般球面线性插值v(t)v(t)v(t)与v1v_1v1的夹角θ(t)\theta(t)θ(t)不是t的线性函数。
证明过程如下(我滴妈呀我的字好丑o(╯□╰)o) 3.改进的球面线性插值有两种方法
1 四元数工具 变换方法
构造四元数q(cosθsinθ∗v1’)r(cosθsinθ∗v2’)q(\cos \theta\sin\theta *v_1’)r(\cos \theta\sin \theta *v_2’)q(cosθsinθ∗v1’)r(cosθsinθ∗v2’)(v1’v2’v_1’v_2’v1’v2’为单位v1,v2v_1,v_2v1,v2向量)以及参数t(0≤t≤1)t(0\leq t\leq1)t(0≤t≤1),则构造四元数变换 四元数s(w,v’)r∗(q−1)t∗qs(w,v’)r*(q-1)t*qs(w,v’)r∗(q−1)t∗q即为球面线性插值变换。其中s的虚部v1’v_1’v1’和v2’v_2’v2’间的插值向量乘以长度x2y2z2\sqrt{x^2y^2z^2}x2y2z2即得到v1,v2v_1,v_2v1,v2间插值向量vvv 另一种变换形式是对四元数进行插值变换 s(w,v′)a∗qb∗rs(w,v)a*qb*r s(w,v′)a∗qb∗r 其中asin(α(1−t))sinα,bsin(αt)sinα,cosαx1∗x2y1∗y2w1∗w2a\frac{\sin(\alpha(1-t))}{\sin\alpha},b\frac{\\sin(\alpha t)}{\sin\alpha},\cos\alphax_1*x_2y_1*y_2w_1*w_2asinαsin(α(1−t)),bsinαsin(αt),cosαx1∗x2y1∗y2w1∗w2 S的虚部v′vv′即为v1′v_1v1′和v2′v_2v2′间的插值向量乘以长度x2y2z2\sqrt{x^2y^2z^2}x2y2z2即得v1,v2v_1,v_2v1,v2间插值向量vvv
2 利用旋转矩阵
变换方法vv1∗Trotvv1*Trotvv1∗Trot
其中TrotTrotTrot即绕任意轴旋转的矩阵变换矩阵因为v1v_1v1到v2v_2v2间的插值可以看成是v1v_1v1绕垂直于v1,v2v_1,v_2v1,v2组成的平面的向量的旋转所以实际上是绕轴旋转的问题不过相应参数变成 θt∗θ\thetat*\thetaθt∗θ轴q(q1,q2,q3)q(q1,q2,q3)q(q1,q2,q3)变成向量v1×v2∣v1×v2∣y1z2−z1y2,z1x2−x1z2,x1y2−y1x2sinθ\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×v2sinθy1z2−z1y2,z1x2−x1z2,x1y2−y1x2
四元数插值
##第一种插值方法
四元数比较重要的一个用途就是球面线性插值(Spherical Linear Interpolation)可以在两个四元数之间平滑插值。
插值步骤
① 计算两个值的差q0q_0q0到q1q_1q1的角位移由Δqq0−1q1\Delta qq_0^{-1}q_1Δqq0−1q1给出
② 计算差的一部分四元数求幂可以做到差的一部分由Δqt\Delta q_tΔqt给出
③ 在开始值上加上差的一部分用四元数乘法组合角位移q0Δqtq_0\Delta q_tq0Δqt
这样就可以得到slerp公式: slerp(q0,q1,t)q0(q0−1q1)tslerp(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.]);end
endy q1 .* (q1.^-1 .* q2).^t;然后使用此函数尝试在运动捕捉数据中进行插值
%方法一matlab自带函数slerp
clear
clc
close all
addpath(genpath(.))%读取两个运动数据skel,A,B
load sample.mat
% skelPlayDataA(skel,[A;B])
%将欧拉角转换为四元数
quatAjoint_euler2quat(skel,A);
quatBjoint_euler2quat(skel,B);
%执行四元数插值插20帧
internum20;
temp_quatzeros(31,4);%31个关节每个关节一个四元数
newMotionzeros(internum,62);%20帧每帧62维
for i1:internumti/internum;%对于角度采用四元数插值for j1:size(quatA,1)temp_quat(j,:)slerp(quatA(j,:),quatB(j,:),t); endtemp_quat(find(isnan(temp_quat)))0;temp_quatreal(temp_quat);newMotion(i,:)joint_quat2euler(temp_quat);%对于位置采用线性插值posAA(1,1:3);posBB(1,1:3);newMotion(i,1:3)(1-t)*posAt*posB;
end
newMotion(find(isnan(newMotion)))0;
skelPlayDataA(skel,[A;newMotion;B])
结果 第二种插值方法
Slerp的思想就是沿着4D4D4D球面上连接两个四元数的弧插值。
先看平面上的两个2D2D2D向量v0v_0v0和v1v_1v1都是单位向量我们需要计算vtv_tvt它是沿着v0v_0v0到v1v_1v1弧的平滑插值。设www是v0v_0v0到vtv_tvt弧所截的角那么vtv_tvt就是v1v_1v1沿弧旋转twtwtw的结果。
需要考虑两点问题一是四元数qqq和−q-q−q代表同一方位但是作为slerp的参数时可能有不一样的结果是因为4D4D4D球面不是欧式空间的直接扩展而这种现象在2D3D2D 3D2D3D空间是不会发生的。解决方法是选择q0q_0q0和q1q_1q1的符号使得点乘q0⋅q1q_0\cdot q_1q0⋅q1的结果是非负。第二就是如果q0q_0q0和q1q_1q1非常接近sinθ\sin\thetasinθ会非常小这时除法会出现问题解决方法是此时采用线性插值。
在论文《从运动捕获数据中提取关键帧》也有介绍到这种四元数插值方法这里直接贴过来有兴趣去看看论文
若q1[w1,x1,y1,z1]q_1[w_1,x_1,y_1,z_1]q1[w1,x1,y1,z1]和q2[w2,x2,y2,z2]q_2[w_2,x_2,y_2,z_2]q2[w2,x2,y2,z2]为两个单位四元数它们之间的球面线性插值为 slerp(q1,q2;t)sin(1−t)θsinθq1sintθsinθq2slerp(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)θq1sinθsintθq2 其中θarccos(w1w2x1x2y1y2z1z2)\theta\arccos(w_1w_2x_1x_2y_1y_2z_1z_2)θarccos(w1w2x1x2y1y2z1z2)
直接撸代码
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 quaternionstheta acos(absD);sinTheta sin(theta);scale0 sin( ( 1.0 - t ) * theta) / sinTheta;scale1 sin( ( t * theta) ) / sinTheta;
end
if(d 0)scale1 -scale1;
endq3 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])
%将欧拉角转换为四元数
quatAjoint_euler2quat(skel,A);
quatBjoint_euler2quat(skel,B);
%执行四元数插值插20帧
internum20;
temp_quatzeros(31,4);%31个关节每个关节一个四元数
newMotionzeros(internum,62);%20帧每帧62维
for i1:internumti/internum;%对于角度采用四元数插值for j1:size(quatA,1)temp_quat(j,:)jointslerp(quatA(j,:),quatB(j,:),t); endnewMotion(i,:)joint_quat2euler(temp_quat);%对于位置采用线性插值posAA(1,1:3);posBB(1,1:3);newMotion(i,1:3)(1-t)*posAt*posB;
end
newMotion(find(isnan(newMotion)))0;
skelPlayDataA(skel,[A;newMotion;B])结果 后记
其实之前写过类似博客但是不是用markdown写的排版真的好丑我就把它们删掉写到此博客了。代码连接链接https://pan.baidu.com/s/1uLadyPL8yPlQWdPpLSWVrw 密码asph
代码也可以到我个人的CSDN上传空间去找或者微信公众号个人简介中的GitHub。此博客已同步更新至微信公众号
