李群资料:
https://en.wikipedia.org/wiki/3D_rotation_group
http://www.ethaneade.com/latex2html/lie_groups/lie_groups.html
http://www.ethaneade.com/latex2html/lie_groups/lie_groups.html
https://github.com/ankurhanda/gvnn
https://en.wikipedia.org/wiki/3D_rotation_group
扰动法:
因此扰动法的导数可以理解为对0处的导数,容易求解,
此时0处的导数容易计算:
只是最终求解出以后:
而不是简单相加,这就是Dynamic Fusion里面的都是在0处进行Taylor展开。
最近在做SMPL相关的工作,SMPL中每个关节的旋转角度采用李代数来表达,其中李代数到李群的转换与罗德里格斯变换相同,可以表达为。而其中微分模型中,SO(3)的微分写作:
具体证明参考BCH公式: https://en.wikipedia.org/wiki/Baker-Campbell-Hausdorff_formula
其中为左雅各比矩阵,
所以可以根据泰勒展开可以得到SO3矩阵对于每个李代数中分量的导数:
其中
其中
根据实际运算中我们也可以通过微小扰动法求得,
根据微小扰动方法得的结果通常与导数法有一定差别,是由于步长选择和数值计算截断误差造成的,由于两种误差来源互相冲突,因此步长的选择有一个合适的范围,通常我们认为解析法计算的结果较数值扰动法为准确,在这种情况下,可以采用Matlab进行实验求出最佳的步长范围:
图1 数值计算误差与step
计算结果证明步长大约等于EPS的10的十次方倍到十一次方倍可以得到最小计算误差,大概为10^-10左右。 在Matlab中EPS通常为2.2204e-16 。
附录1: 罗德里格斯的matlab实现
%{
% 罗德里格斯变换 % 输入:1*3的罗德里格斯向量 % 输出3*3矩阵 为旋转矩阵
% 参考文献:
% 新浪微博:http://blog.sina.com.cn/s/blog_5fb3f125010100hp.html
% CSDN空间:http://blog.csdn.net/bugrunner/article/details/7359100
% 维基百科:http://en.wikipedia.org/wiki/Rodrigues'_rotation_formula
%}
function Rotation_matrix=Rodrigues2Rotation(rod)
r = rod;
theta=norm(r);
% r = r/theta ;
%%
% rx=[ 0 -r(3) r(2) ;
% r(3) 0 -r(1) ;
% -r(2) r(1) 0 ]; %%
% Rotation_matrix = eye(3,3)+sin(theta)/theta*rx+(1-cos(theta))/theta^2*rx*rx;
one_r = r./theta;
Rotation_matrix=cos(theta)*eye(3,3)+(1-cos(theta))*(one_r*one_r')+sin(theta)*[0 -one_r(3) one_r(2);one_r(3) 0 -one_r(1); -one_r(2) one_r(1) 0] ;
% 下面这种办法也可以:
%{
rx=[0 -r(3) r(2);r(3) 0 -r(1); -r(2) r(1) 0];
R2=eye(3,3)+sin(theta)/theta*rx+(1-cos(theta))/theta^2*rx*rx;
(R-R')/2-sin(theta)*[0 -one_r(3) one_r(2);one_r(3) 0 -one_r(1);-one_r(2) one_r(1) 0 ]
Rx=(R-R')/2;
%}
附录2: 数值导数计算方法:
function matrix = diff_rod(vec , i , delta_rank)
delta = eps*10^delta_rank ;
vec1 = vec ;
'vec2 = vec ;
vec1(i) = vec(i) - delta ;
vec2(i) = vec(i) + delta ;
matrix = (Rodrigues2Rotation(vec2) - Rodrigues2Rotation(vec1) ) / (delta*2.0) ; %
end
附录3:解析法计算SO(3)导数
function matrix = deriviive( vec , i )
matrix_l = Jaccobi_l(vec) ;
vec_l = matrix_l(:,i) ; %%%
matrix = skew_matrix(vec_l) * Rodrigues2Rotation( vec ) ; %%
function skew_mat = skew_matrix(vec_n)
skew_mat=[ 0 -vec_n(3) vec_n(2) ;
vec_n(3) 0 -vec_n(1) ;
-vec_n(2) vec_n(1) 0 ];
function J_l = Jaccobi_l(vec)
theta = norm(vec);
vec_n = vec / theta ;
rx=[ 0 -vec_n(3) vec_n(2) ;
vec_n(3) 0 -vec_n(1) ;
-vec_n(2) vec_n(1) 0 ];
J_l = sin(theta)/theta*eye(3)+( 1-sin(theta)/theta )*(vec_n*vec_n') + (1-cos(theta))/theta * rx ;
% J_l = inv(J_l) ;
function matrix = i_matrix(i)
if i ==1
matrix=[ 0 0 0 ;
0 0 -1 ;
0 1 0];
else
if i==2
matrix=[ 0 0 1;
0 0 0;
-1 0 0];
else
matrix=[ 0 -1 0 ;
1 0 0 ;
0 0 0] ;
end
end
其中写作如下: