basalt中的相对位姿和绝对位姿导数
在作者的这部分代码是相对难以和代码对应的,使用相对位姿的导数目的:为了节约很多计算量,假设两帧之间有很多特征点,vins等的做法是特征点残差直接对绝对位姿求导数,basalt是残差对相对位姿求导数×相对位姿对绝对位姿求导数,这样在计算过程中虽然左半部分残差对相对位姿求导不一样但是相对位姿对绝对位姿求导数是一样的因此这部分就节省了很多的计算量。这部分代码要注意的要点:
- 作者使用 S O 3 + t r a n s SO3 + trans SO3+trans的形式表述位姿
- 作者采用了巧妙的坐标准换,使得式子更加相对
s o 3 + t r a n s so3+trans so3+trans和 s e 3 se3 se3导数的差别
basalt中定义的位姿的加法,可以看出其是先进行平移加再进行旋转的左乘,但更新方式和 S E ( 3 ) SE(3) SE(3)是有区别的。做这个步骤的操作目的:在推倒公式的时候SE形式具有很好的公式和库的支持,so+trans的形式就需要自己分开写表达式再写,麻烦!因此可以对SE整体求导再转换回来。
inline static void incPose(const Vec6& inc, SE3& T) {
T.translation() += inc.template head<3>();// 位移更新
T.so3() = SO3::exp(inc.template tail<3>()) * T.so3();// 位姿进行左扰动更新
}
S E ( 3 )增量 SE(3)增量 SE(3)增量 δ ξ 1 = [ δ ρ 1 δ ϕ 1 ] \delta \xi_1 = \begin{bmatrix} \delta \rho_1 \\ \delta \phi_1 \\ \end{bmatrix} δξ1=[δρ1δϕ1] S O ( 3 ) + t r a n s 增量 SO(3)+trans增量 SO(3)+trans增量 δ ξ 2 = [ δ ρ 2 δ ϕ 2 ] \delta \xi_2 = \begin{bmatrix} \delta \rho_2 \\ \delta \phi_2 \\ \end{bmatrix} δξ2=[δρ2δϕ2]
S E ( 3 ) 更新 SE(3)更新 SE(3)更新 E x p ( δ ξ 1 ) T = [ ( I + [ δ ϕ 1 ] × ) δ ρ 1 0 I ] [ ( R t 0 I ] Exp(\delta \xi_1)T = \begin{bmatrix} (I+[\delta \phi_1]^{×}) & \delta\rho_1 \\\ 0 & I\\ \end{bmatrix} \begin{bmatrix} (R & t \\ 0 & I\\ \end{bmatrix} Exp(δξ1)T=[(I+[δϕ1]×) 0δρ1I][(R0tI] 式子1
S O ( 3 ) + t r a n s 更新 SO(3)+trans更新 SO(3)+trans更新 E x p ( δ ξ 2 ) T = [ ( I + [ δ ϕ 2 ] × ) R t + δ t 2 0 I ] Exp(\delta \xi_2)T = \begin{bmatrix} (I+[\delta \phi_2]^{×})R & t + \delta t_2 \\\ 0 & I\\ \end{bmatrix} Exp(δξ2)T=[(I+[δϕ2]×)R 0t+δt2I] 式子2
将式子1和式子2进行对比得到
δ ξ 1 = [ δ ρ 1 δ ϕ 1 ] = [ I [ t ] × 0 I ] [ δ ρ 2 δ ϕ 2 ] = [ I [ t ] × 0 I ] δ ξ 2 \delta \xi_1 = \begin{bmatrix} \delta \rho_1 \\ \delta \phi_1 \\ \end{bmatrix} = \begin{bmatrix} I & [t]^{×} \\ 0 & I \\ \end{bmatrix} \begin{bmatrix} \delta \rho_2 \\ \delta \phi_2 \\ \end{bmatrix} =\begin{bmatrix} I & [t]^{×} \\ 0 & I \\ \end{bmatrix} \delta \xi_2 δξ1=[δρ1δϕ1]=[I0[t]×I][δρ2δϕ2]=[I0[t]×I]δξ2
因此得到的一个关系是: ∂ δ ξ 1 ∂ δ ξ 2 = [ I [ t ] × 0 I ] \frac{\partial \delta \xi_1}{ \partial \delta \xi_2} = \begin{bmatrix} I & [t]^{×} \\ 0 & I \\ \end{bmatrix} ∂δξ2∂δξ1=[I0[t]×I]
一个变换
其实相对位姿对位姿的导数在slam14的位姿图优化部分已经介绍了,但是和代码还是有一点对应不上的地方,对应不上的就是一个小转换。把 T i w T_i^w Tiw和 T j w T_j^w Tjw都乘以一个 T ′ = [ I − t w i w 0 1 ] T' = \begin{bmatrix} I & -t_{wi}^w \\ 0 & 1 \end{bmatrix} T′=[I0−twiw1]的变换 其作用相当于转换到在i位置新建一个和世界坐标系姿态一样但原点在i位置处的坐标系中。感性的认识是在这个坐标系下的旋转和平移增量如果是定义在解耦的so3+trans的形式下就是在世界坐标系下的增量
感性认识理解:原来basalt定义的残差是 e r r = l n [ T j i ( ( T j w ) T ∗ T i w ) ] err = ln[ T^i_j ((T^w_j)^T * T^w_i) ] err=ln[Tji((Tjw)T∗Tiw)] 加上这个变换后 e r r n e w = l n [ T j i ( ( T ′ T j w ) T ∗ T ′ T i w ) ] = l n [ T j i ( ( T j w ) T T ′ T T ′ T i w ) ] = l n [ T j i ( ( T j w ) T ∗ T i w ) ] err_{new} = ln[ T^i_j ((T' T^w_j)^T * T' T^w_i) ]=ln[T_j^i((T_j^w)^T T'^T T'T_i^w)]=ln[ T^i_j ((T^w_j)^T * T^w_i) ] errnew=ln[Tji((T′Tjw)T∗T′Tiw)]=ln[Tji((Tjw)TT′TT′Tiw)]=ln[Tji((Tjw)T∗Tiw)] 即这两个误差是一样的,但此时再对其进行求雅克比的时候则对 T ′ T x w T'T_x^w T′Txw一起求导得到的都是,再带入se3对so3+trans的导数关系,这样确实对i的导数trans部分和j的导数在Adj上的区别了。用人话来说就是,变换之前世界坐标系下的导数和当前转换后世界坐标系下的导数及残差应该是一样的,因此选择转换后的坐标系为基准,对转换后的坐标内求相对对绝对位姿的导数是和代码的表述是一致的。
证明。首先对于导数表达式子中对se3的导数只是正负号有差别,带入 ∂ s e 3 ∂ s o 3 + t r a n s = [ I [ t ] × 0 I ] \frac{\partial se3}{\partial so3+trans} = \begin{bmatrix}I & [t]^{×} \\ 0 & I \\ \end{bmatrix} ∂so3+trans∂se3=[I0[t]×I] 到雅克比矩阵中进行简写为(…为 J r J_r Jr的那一部分) . . A d j ( T i T ) [ I [ t w ] × 0 I ] = . . . [ R i T [ − R i T ∗ t i w ] × 0 R i T ] [ I [ t w ] × 0 I ] ..Adj(Ti^T) \begin{bmatrix} I & [t^w]^{×} \\ 0 & I\end{bmatrix} = ...\begin{bmatrix} R^T_i & [-R^T_i*t_i^w]^{×} \\ 0 & R^T_i \end{bmatrix} \begin{bmatrix} I & [t^w]^{×} \\ 0 & I\end{bmatrix} ..Adj(TiT)[I0[tw]×I]=...[RiT0[−RiT∗tiw]×RiT][I0[tw]×I] 对式子进行打开得到:
[ R i T R i T [ t w ] × + ( − R i T t i w ) × R i T 0 R i T ] = [ R i T R i T [ t w ] × ∗ R i ∗ R i T + ( − R i T t i w ) × R i T 0 R i T ] \begin{bmatrix} R_i^T & R_i^T[t^w]^{×} + (-R_i^Tt_i^w)^{×}R_i^T \\ 0 & R_i^T\end{bmatrix} = \begin{bmatrix} R_i^T & R_i^T[t^w]^{×}*R_i *R_i^T + (-R_i^Tt_i^w)^{×}R_i^T \\ 0 & R_i^T\end{bmatrix} [RiT0RiT[tw]×+(−RiTtiw)×RiTRiT]=[RiT0RiT[tw]×∗Ri∗RiT+(−RiTtiw)×RiTRiT]
= [ R i T ( R i T t w ) × ∗ R i T + ( − R i T t i w ) × R i T 0 R i T ] = \begin{bmatrix} R_i^T & (R_i^Tt^w)^{×} *R_i^T + (-R_i^Tt_i^w)^{×}R_i^T \\ 0 & R_i^T\end{bmatrix} =[RiT0(RiTtw)×∗RiT+(−RiTtiw)×RiTRiT]
= [ R i T ( R i T t w − R i T t i w ) × ∗ R i T 0 R i T ] = \begin{bmatrix} R_i^T & (R_i^Tt^w-R_i^Tt_i^w)^{×} *R_i^T \\ 0 & R_i^T\end{bmatrix} =[RiT0(RiTtw−RiTtiw)×∗RiTRiT]
= [ R i T ( R i T ( t w − t i w ) ) × ∗ R i T 0 R i T ] = \begin{bmatrix} R_i^T & (R_i^T(t^w-t_i^w))^{×} *R_i^T \\ 0 & R_i^T\end{bmatrix} =[RiT0(RiT(tw−tiw))×∗RiTRiT]
当对 i i i求导数的时候, t w = t i w t^w = t_i^w tw=tiw此时矩阵为 = [ R i T 0 0 R i T ] = \begin{bmatrix} R_i^T & 0 \\ 0 & R_i^T\end{bmatrix} =[RiT00RiT]
当对 j j j求导数的时候, t w = t j w t_w = t_j^w tw=tjw此时矩阵为 = [ R i T ( R i T ( t i j w ) ) × ∗ R i T 0 R i T ] = [ R i T ( t i j i ) × ∗ R i T 0 R i T ] = \begin{bmatrix} R_i^T & (R_i^T(t_{ij}^w))^{×} *R_i^T \\ 0 & R_i^T\end{bmatrix} = \begin{bmatrix} R_i^T & (t_{ij}^i)^{×} *R_i^T \\ 0 & R_i^T\end{bmatrix} =[RiT0(RiT(tijw))×∗RiTRiT]=[RiT0(tiji)×∗RiTRiT]
总结
导数为:
∂ e r r T i w = J r − 1 ( r e s ) ∗ [ R i T 0 0 R i T ] \frac{\partial err}{T^w_i} = J_r^{-1}(res) * \begin{bmatrix} R_i^T & 0 \\ 0 & R_i^T\end{bmatrix} Tiw∂err=Jr−1(res)∗[RiT00RiT]
∂ e r r T j w = J r − 1 ( r e s ) ∗ [ R i T ( t i j i ) × ∗ R i T 0 R i T ] \frac{\partial err}{T^w_j} = J_r^{-1}(res) * \begin{bmatrix} R_i^T & (t_{ij}^i)^{×} *R_i^T \\ 0 & R_i^T\end{bmatrix} Tjw∂err=Jr−1(res)∗[RiT0(tiji)×∗RiTRiT]
对应代码为:
inline Sophus::Vector6d relPoseError(
const Sophus::SE3d& T_i_j, const Sophus::SE3d& T_w_i,
const Sophus::SE3d& T_w_j, Sophus::Matrix6d* d_res_d_T_w_i = nullptr,
Sophus::Matrix6d* d_res_d_T_w_j = nullptr) {
Sophus::SE3d T_j_i = T_w_j.inverse() * T_w_i;// 相对位姿
Sophus::Vector6d res = Sophus::se3_logd(T_i_j * T_j_i);// 计算残差
if (d_res_d_T_w_i || d_res_d_T_w_j) {// 如果需要计算jacobian
Sophus::Matrix6d J;
Sophus::rightJacobianInvSE3Decoupled(res, J);// 算jr的逆
Eigen::Matrix3d R = T_w_i.so3().inverse().matrix();
Sophus::Matrix6d Adj;
Adj.setZero();
Adj.topLeftCorner<3, 3>() = R;
Adj.bottomRightCorner<3, 3>() = R;
if (d_res_d_T_w_i) {// 计算左边的导数
*d_res_d_T_w_i = J * Adj;
}
if (d_res_d_T_w_j) {// 如果计算残差对右边的导数
// 在右上角填充
Adj.topRightCorner<3, 3>() = Sophus::SO3d::hat(T_j_i.inverse().translation()) * R;
*d_res_d_T_w_j = -J * Adj;
}
}
return res;
}