lsd:tracking

最新推荐文章于 2022-03-22 21:51:33 发布

可爱的小蚂蚁

最新推荐文章于 2022-03-22 21:51:33 发布

阅读量1.1k

点赞数 2

文章标签： lsd

本文链接：https://blog.youkuaiyun.com/u013004597/article/details/52303017

版权

svo_slam 专栏收录该内容

3 篇文章

订阅专栏

SE(3) Tracking new frame:
　　新图像与模板图像在映射后的光度残差表达式：

r p (x, ξ j i) = I i (x) - I j (ω (x, D i (x), ξ j i))

$\begin{equation} r_p(\mathbf{x,\xi_{ji}})=I_i(x)-I_j(\mathbf{\omega(x,D_i(x),\xi_{ji}))} \notag \end{equation}$
其中j是新图像帧，i是关键图像帧。forward additive的光度残差对迭代增量

Δξ $\Delta \mathbf{\xi}$ 求导得：

J Δ ξ = \nabla I | ω (x, (T) \partial π \partial p ' ∣ ∣ ∣ T p . \partial p ' \partial T ∣ ∣ ∣ Δ T T i . \partial T \partial Δ T ∣ ∣ ∣ I . \partial e x p ( Δ ξ ^ ) \partial Δ ξ ∣ ∣ ∣ 0 .

$\begin{equation} \begin{split} &\mathbf{J_{\Delta \xi}} =\nabla I|_{\omega(x,(\mathbf{T})} \left. \dfrac{\partial{\pi}}{\partial{p'}} \right |_{\mathbf{Tp}}. \left. \dfrac{\partial{p'}}{\partial{T}} \right |_{\mathbf{\Delta T T_i}}. \left. \dfrac{\partial{\mathbf{T}}}{\partial{\Delta T}} \right |_{\mathbf{I}}. \left. \dfrac{\partial{exp(\hat{\Delta \mathbf{\xi}})}}{\partial{\Delta \mathbf{\xi}}} \right |_{\mathbf{0}}. \\ \end{split} \notag \end{equation}$
其中：

∇I|ω(x,T)=(ΔIuΔIv) $\nabla I|_{\omega(x,\mathbf{T})}=(\Delta I_u \quad \Delta I_v)$ ：模板图像上像素坐标

x(u,v) $x(u,v)$ 在T进行warp后得到的新图像

I′ $I'$ 中像素坐标

x′(u′,v′) $x'(u',v')$ 的梯度.

∂π∂p′∣∣∣Tp.=⎛⎝⎜⎜fx1z′00fy1z′−fxx′z′2−fyy′z′2⎞⎠⎟⎟ $\left. \dfrac{\partial{\pi}}{\partial{p’}} \right |_{\mathbf{Tp}}.= \left( \begin{array}{ccc} f_x\dfrac{1}{z'}&0&-f_x\dfrac{x'}{z'^2}\\ 0&f_y\dfrac{1}{z'}&-f_y\dfrac{y'}{z'^2} \end{array} \right)$ :针孔相机模型求导。
这里令

Ti=⎛⎝⎜⎜⎜⎜r11r21r310r12r22r320r13r23r330txtytz1⎞⎠⎟⎟⎟⎟ $T_i= \left( \begin{array}{ccc|c} r_{11}&r_{12}&r_{13}&t_x\\ r_{21}&r_{22}&r_{23}&t_y\\ r_{31}&r_{32}&r_{33}&t_z\\ \hline 0&0&0&1 \end{array} \right)$ ,

p′=Tip $p’=T_ip$ , 所以有：

\partial p ' \partial T ∣ ∣ ∣ Δ T T i . \partial T \partial Δ T ∣ ∣ ∣ I . \partial e x p ( Δ ξ ^ ) \partial Δ ξ ∣ ∣ ∣ 0 . = ⎛ ⎝ ⎜ ⎜ x 00 0 x 0 00 x y 00 0 y 0 00 y z 00 0 z 0 00 z 100010001 ⎞ ⎠ ⎟ ⎟ = ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ r 11 00 r 12 00 r 13 00 t x 00 0 r 11 00 r 12 00 r 13 00 t x 0 00 r 11 00 r 12 00 r 13 00 t x r 21 00 r 22 00 r 23 00 t y 00 0 r 21 00 r 22 00 r 23 00 t y 0 00 r 21 00 r 22 00 r 23 00 t y r 31 00 r 32 00 r 33 00 t z 00 0 r 31 00 r 32 00 r 33 00 t z 0 00 r 31 00 r 32 00 r 33 00 t z 000000000100000000000010000000000001 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ = ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 000000000100000000000010000000000001 0000010 - 1 0000 00 - 1 000100000 010 - 1 00000000 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\begin{equation} \begin{split} \left. \dfrac{\partial{p'}}{\partial{T}} \right |_{\mathbf{\Delta T T_i}}.&= \left( \begin{array}{ccc|ccc|ccc|ccc} x&0&0&y&0&0&z&0&0&1&0&0\\ 0&x&0&0&y&0&0&z&0&0&1&0\\ 0&0&x&0&0&y&0&0&z&0&0&1 \end{array} \right) \\ \left. \dfrac{\partial{\mathbf{T}}}{\partial{\Delta T}} \right |_{\mathbf{I}}.&= \left( \begin{array}{ccc|ccc|ccc|ccc} r_{11}&0&0&r_{21}&0&0&r_{31}&0&0&0&0&0\\ 0&r_{11}&0&0&r_{21}&0&0&r_{31}&0&0&0&0\\ 0&0&r_{11}&0&0&r_{21}&0&0&r_{31}&0&0&0\\ r_{12}&0&0&r_{22}&0&0&r_{32}&0&0&0&0&0\\ 0&r_{12}&0&0&r_{22}&0&0&r_{32}&0&0&0&0\\ 0&0&r_{12}&0&0&r_{22}&0&0&r_{32}&0&0&0\\ r_{13}&0&0&r_{23}&0&0&r_{33}&0&0&0&0&0\\ 0&r_{13}&0&0&r_{23}&0&0&r_{33}&0&0&0&0\\ 0&0&r_{13}&0&0&r_{23}&0&0&r_{33}&0&0&0\\ t_x&0&0&t_y&0&0&t_z&0&0&1&0&0\\ 0&t_x&0&0&t_y&0&0&t_z&0&0&1&0\\ 0&0&t_x&0&0&t_y&0&0&t_z&0&0&1\\ \end{array} \right)\\ \left. \dfrac{\partial{exp(\hat{\Delta \mathbf{\xi}})}}{\partial{\Delta \mathbf{\xi}}} \right |_{\mathbf{0}}.&= \left( \begin{array}{ccc|ccc} 0&0&0&0&0&0\\ 0&0&0&0&0&1\\ 0&0&0&0&-1&0\\ 0&0&0&0&0&-1\\ 0&0&0&0&0&0\\ 0&0&0&1&0&0\\ 0&0&0&0&1&0\\ 0&0&0&-1&0&0\\ 0&0&0&0&0&0\\ 1&0&0&0&0&0\\ 0&1&0&0&0&0\\ 0&0&1&0&0&0\\ \end{array} \right) \end{split} \notag \end{equation}$
　　注意上面的表达式中在迭代过程中，

∂p′∂T∣∣∣ΔTTi. $\left. \dfrac{\partial{p'}}{\partial{T}} \right |_{\mathbf{\Delta T T_i}}.$ 和

∂exp(Δξ^)∂Δξ∣∣∣0. $\left. \dfrac{\partial{exp(\hat{\Delta \mathbf{\xi}})}}{\partial{\Delta \mathbf{\xi}}} \right |_{\mathbf{0}}.$ 是始终不变的。但是

∇I|ω(x,T) $\nabla I|_{\omega(x,\mathbf{T})}$ 、

∂π∂p′∣∣∣Tp. $\left. \dfrac{\partial{\pi}}{\partial{p'}} \right |_{\mathbf{Tp}}.$ 和

∂T∂ΔT∣∣∣I. $\left. \dfrac{\partial{\mathbf{T}}}{\partial{\Delta T}} \right |_{\mathbf{I}}.$ 每次迭代都是不一样的。由

Ti $T_i$ 的表达式可知：

⎧ ⎩ ⎨ x' = r 11 x + r 12 y + r 13 z + t x y' = r 21 x + r 22 y + r 23 z + t y z' = r 31 x + r 32 y + r 33 z + t z

$\begin{equation} \begin{cases} x'=r_{11}x+r_{12}y+r_{13}z+t_x\\ y'=r_{21}x+r_{22}y+r_{23}z+t_y\\ z'=r_{31}x+r_{32}y+r_{33}z+t_z \end{cases} \notag\\ \end{equation}$
所以：

∂p′∂T∣∣∣ΔTTi.∂T∂ΔT∣∣∣I.= ⎛⎝⎜⎜x′000x′000x′y′000y′000y′z′000z′000z′100010001⎞⎠⎟⎟ $\left. \dfrac{\partial{p’}}{\partial{T}} \right |_{\mathbf{\Delta T T_i}}. \left. \dfrac{\partial{\mathbf{T}}}{\partial{\Delta T}} \right |_{\mathbf{I}}.=\ \left( \begin{array}{ccc|ccc|ccc|ccc} x'&0&0&y'&0&0&z'&0&0&1&0&0\\ 0&x'&0&0&y'&0&0&z'&0&0&1&0\\ 0&0&x'&0&0&y'&0&0&z'&0&0&1 \end{array} \right)$

exp(Δξ^)=(e[ω]×0Vv1) $exp(\hat{\Delta \xi})= \left( \begin{array}{cc} e^{[\omega]_{\times}}&\mathbf{V}v\\ 0&1 \end{array} \right)$
其中：

[ω]×=⎛⎝⎜0ω3−ω2−ω30ω1ω2−ω10⎞⎠⎟;v=⎛⎝⎜v1v2v3⎞⎠⎟ $[\omega]_{\times}= \left( \begin{array}{ccc} 0&-\omega_3&\omega_2\\ \omega_3&0&-\omega_1\\ -\omega_2&\omega_1&0 \end{array} \right); \mathbf{v}= \left( \begin{array}{c} v_1\\ v_2\\ v_3 \end{array} \right)$

e [ω] \times V = I + sin ( | | ω | | ) | | ω | | [ω] \times + 1 - cos ( | | ω | | ) | | ω | | 2 [ω] 2 \times = I + 1 - cos ( | | ω | | ) | | ω | | 2 [ω] \times + | | ω | | - sin ( | | ω | | ) | | ω | | 3 [ω] 2 \times

$\begin{equation} \begin{split} e^{[\mathbf{\omega}]_\times}&= \mathbf{I}+\dfrac{\sin(||\mathbf{\omega}||)}{||\mathbf{\omega}||}[\mathbf{\omega}]_\times+\dfrac{1-\cos(||\mathbf{\omega}||)}{||\mathbf{\omega}||^2}[\mathbf{\omega}]^2_\times \\ \mathbf{V}&=\mathbf{I}+\dfrac{1-\cos(||\mathbf{\omega}||)}{||\mathbf{\omega}||^2}[\mathbf{\omega}]_\times+\dfrac{||\mathbf{\omega}||-\sin(||\mathbf{\omega}||)}{||\mathbf{\omega}||^3}[\mathbf{\omega}]^2_\times \end{split} \notag \end{equation}$
上式在

||ω||=0 $||\mathbf{\omega}||=0$ 附近，对

sin(||ω||) $\sin(||\mathbf{\omega}||)$ 和

cos(||ω||) $\cos(||\mathbf{\omega}||)$ 进行一阶泰勒展开得：

e [ω] \times V \approx I + 1 2 [ω] \times + 1 6 [ω] 2 \times \approx I + [ω] \times + 1 2 [ω] 2 \times

$\begin{equation} \begin{split} e^{[\mathbf{\omega}]_\times}& \approx \mathbf{I}+\dfrac{1}{2}[\mathbf{\omega}]_\times+\dfrac{1}{6}[\mathbf{\omega}]^2_\times \\ \mathbf{V}& \approx \mathbf{I}+[\mathbf{\omega}]_\times+\dfrac{1}{2}[\mathbf{\omega}]^2_\times \end{split} \notag \end{equation}$
所以矩阵求导得：
这里写图片描述

其中

e1=[1,0,0]T $\mathbf{e_1=[1,0,0]^T}$ ,

e2=[0,1,0]T $\mathbf{e_2=[0,1,0]^T}$ ,

e1=[0,0,1]T $\mathbf{e_1=[0,0,1]^T}$ .

J Δ ξ = = = (\nabla I u \nabla I v) T ⎛ ⎝ ⎜ ⎜ f x 1 z ' 0 0 f y 1 z ' - f x x ' z ' 2 - f y y ' z ' 2 ⎞ ⎠ ⎟ ⎟ ⎛ ⎝ ⎜ ⎜ x' 00 0 x' 0 00 x' y' 00 0 y' 0 00 y' z' 00 0 z' 0 00 z' 100010001 ⎞ ⎠ ⎟ ⎟ ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ 0 3 \times 3 0 3 \times 3 0 3 \times 3 I 3 \times 3 - [e 1] \times - [e 2] \times - [e 2] \times 0 3 \times 3 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ (\nabla I u \nabla I v) ⎛ ⎝ ⎜ ⎜ f x 1 z ' 0 0 f y 1 z ' - f x x ' z ' 2 - f y y ' z ' 2 ⎞ ⎠ ⎟ ⎟ ⎛ ⎝ ⎜ ⎜ 100010001 0 - z y' z' 0 - x' - y' x' 0 ⎞ ⎠ ⎟ ⎟ 1 z ' (\nabla I u f x \nabla I v f y) T ⎛ ⎝ ⎜ ⎜ ⎜ 1001 - x ' z ' - y ' z ' - x ' y ' z ' - (z' + y ' 2 z ') z' + x ' 2 z ' x ' y ' z ' - y' x' ⎞ ⎠ ⎟ ⎟ ⎟

$\begin{equation} \begin{split} \mathbf{J_{\Delta \xi}} &=\left( \begin{array}{c} \nabla I_u \\ \nabla I_v \end{array} \right)^\mathbb{T} \left( \begin{array}{ccc} f_x\dfrac{1}{z'}&0&-f_x\dfrac{x'}{z'^2}\\ 0&f_y\dfrac{1}{z'}&-f_y\dfrac{y'}{z'^2} \end{array} \right) \left( \begin{array}{ccc|ccc|ccc|ccc} x'&0&0&y'&0&0&z'&0&0&1&0&0\\ 0&x'&0&0&y'&0&0&z'&0&0&1&0\\ 0&0&x'&0&0&y'&0&0&z'&0&0&1 \end{array} \right) \left( \begin{array}{c|c} \mathbf{0_{3\times 3}}&-[\mathbf{e_1}]_\times \\ \mathbf{0_{3\times 3}}&-[\mathbf{e_2}]_\times \\ \mathbf{0_{3\times 3}}&-[\mathbf{e_2}]_\times \\ \hline \mathbf{I_{3\times 3}}&\mathbf{0_{3\times 3}} \end{array} \right)\\ =&(\nabla I_u \quad \nabla I_v) \left( \begin{array}{ccc} f_x\dfrac{1}{z'}&0&-f_x\dfrac{x'}{z'^2}\\ 0&f_y\dfrac{1}{z'}&-f_y\dfrac{y'}{z'^2} \end{array} \right) \left( \begin{array}{ccc|ccc} 1&0&0&0&z'&-y'\\ 0&1&0&-z&0&x'\\ 0&0&1&y'&-x'&0 \end{array} \right)\\ =&\dfrac{1}{z'} \left( \begin{array}{c} \nabla I_u f_x\\ \nabla I_v f_y \end{array} \right)^\mathbb{T} \left( \begin{array}{ccc|ccc} 1&0&-\dfrac{x'}{z'}&-\dfrac{x'y'}{z'}&z'+\dfrac{x'^2}{z'}&-y'\\ 0&1&-\dfrac{y'}{z'}&-(z'+\dfrac{y'^2}{z'})&\dfrac{x'y'}{z'}&x' \end{array} \right)\\ \end{split} \notag \end{equation}$
　　注意：forward additive的本质是将第一幅图像向第二幅图像warp后，比较warp得到的图像和第二幅图像之间的残差。在迭代过程中,warp后图像的像素坐标

(u′,v′) $(u',v')$ ,点坐标

(x′,y′,z′) $(x',y',z')$ 需要每次都重新求。
　　光度残差协方差的来源分别为：模板图像光度协方差、新图像光度协方差、逆深度估计引起的模板光度协方差。根据残差协方差传播的方程：

ΣI≈JΣDTJT $\mathbf{\Sigma_I \approx J \Sigma_{D_T} J^T}$ ,\和函数:

I(ω(x,DT(X),T))=I(π(Tπ−1(x,1DT(x)))) $I(\omega(x,D_T(X),\mathbf{T}))=I(\pi(\mathbf{T}\pi^{-1}(x,\dfrac{1}{D_T(x)})))$ ，所以光度对逆深度的求导为：

\partial r p ( T ) \partial D T = \nabla I | ω (x, T i) \partial π \partial p ' ∣ ∣ ∣ T i p . \partial p ' \partial p ∣ ∣ ∣ p . \partial π - 1 \partial Z ∣ ∣ ∣ 1 D T ( x ) . \partial 1 x \partial x ∣ ∣ ∣ ∣ ∣ D T (x) .

$\begin{equation} \dfrac{\partial{\mathbf{r_p(T)}}}{\partial{D_T}}=\nabla I|_{\omega(x,\mathbf{T_i})} \left. \dfrac{\partial{\pi}}{\partial{p'}} \right |_{\mathbf{T_ip}}. \left. \dfrac{\partial{p'}}{\partial{p}} \right |_{\mathbf{p}}. \left. \dfrac{\partial{\pi^{-1}}}{\partial{Z}} \right |_{\dfrac{1}{D_T(x)}}. \left. \dfrac{\partial{\dfrac{1}{x}}}{\partial{x}} \right |_{D_T(x)}. \notag \end{equation}$
其中：
这里写图片描述

因为

{ππ−1:(u,v)=Kp:p=K−1(u,v) $\begin{cases}\pi&:(u,v)=Kp \\ \pi^{-1}&:p=K^{-1}(u,v)\end{cases}$ ,所以有

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪x=z(u−cx)fxy=z(u−cy)fyz=z $\begin{cases}x=\dfrac{z(u-c_x)}{f_x} \\ y=\dfrac{z(u-c_y)}{f_y}\\z=z\end{cases}$ ,
因而：
这里写图片描述

该导数化简为：

\partial r p ( T ) \partial D T = = = = \nabla I | ω (x, T i) \partial π \partial p ' ∣ ∣ ∣ T i p . \partial p ' \partial p ∣ ∣ ∣ p . \partial π - 1 \partial Z ∣ ∣ ∣ 1 D T ( x ) . \partial 1 x \partial x ∣ ∣ ∣ ∣ ∣ D T (x) . (\nabla I u \nabla I v) ⎛ ⎝ ⎜ ⎜ f x 1 z ' 0 0 f y 1 z ' - f x x ' z ' 2 - f y y ' z ' 2 ⎞ ⎠ ⎟ ⎟ ⎛ ⎝ ⎜ r 11 r 21 r 31 r 12 r 22 r 32 r 13 r 23 r 33 ⎞ ⎠ ⎟ ⎛ ⎝ ⎜ x y z ⎞ ⎠ ⎟ D T (x) - 1 D T ( x ) 2 (\nabla I u \nabla I v) T ⎛ ⎝ ⎜ ⎜ f x 1 z ' 0 0 f y 1 z ' - f x x ' z ' 2 - f y y ' z ' 2 ⎞ ⎠ ⎟ ⎟ ⎛ ⎝ ⎜ x' - t x y' - t y z' - t z ⎞ ⎠ ⎟ 1 D T ( x ) 1 D T ( x ) z ' 2 (\nabla I u f x (t x z' - t z x') + \nabla I v f y (t y z' - t z y'))

$\begin{equation} \begin{split} \dfrac{\partial{\mathbf{r_p(T)}}}{\partial{D_T}}=&\nabla I|_{\omega(x,\mathbf{T_i})} \left. \dfrac{\partial{\pi}}{\partial{p'}} \right |_{\mathbf{T_ip}}. \left. \dfrac{\partial{p'}}{\partial{p}} \right |_{\mathbf{p}}. \left. \dfrac{\partial{\pi^{-1}}}{\partial{Z}} \right |_{\dfrac{1}{D_T(x)}}. \left. \dfrac{\partial{\dfrac{1}{x}}}{\partial{x}} \right |_{D_T(x)}.\\ =&(\nabla I_u \quad \nabla I_v) \left( \begin{array}{ccc} f_x\dfrac{1}{z'}&0&-f_x\dfrac{x'}{z'^2}\\ 0&f_y\dfrac{1}{z'}&-f_y\dfrac{y'}{z'^2} \end{array} \right) \left( \begin{array}{ccc} r_{11}&r_{12}&r_{13}\\ r_{21}&r_{22}&r_{23}\\ r_{31}&r_{32}&r_{33} \end{array} \right) \left( \begin{array}{c} x\\ y\\ z \end{array} \right)D_T(x) \dfrac{-1}{D_T(x)^2}\\ =& \left( \begin{array}{c} \nabla I_u\\ \nabla I_v \end{array} \right)^T \left( \begin{array}{ccc} f_x\dfrac{1}{z'}&0&-f_x\dfrac{x'}{z'^2}\\ 0&f_y\dfrac{1}{z'}&-f_y\dfrac{y'}{z'^2} \end{array} \right) \left( \begin{array}{c} x'-tx\\ y'-ty\\ z'-tz \end{array} \right)\dfrac{1}{D_T(x)}\\ =& \dfrac{1}{D_T(x)z'^2}(\nabla I_uf_x(t_xz'-t_zx')+\nabla I_vf_y(t_yz'-t_zy')) \end{split} \notag \end{equation}$
综上，光度残差的协方差计算公式是：

σ 2 r p (T) = = σ 2 I + σ 2 J + (\partial r p ( T ) \partial D T) 2 σ 2 D T 2 σ 2 I + (\partial r p ( T ) \partial D T) 2 V D T

$\begin{equation} \begin{split} \sigma^2_{r_p(\mathbf{T})}=&\sigma^2_I+\sigma^2_J+(\dfrac{\partial{\mathbf{r_p(T)}}}{\partial{D_T}})^2 \sigma^2_{D_T} \\ =&2\sigma^2_I+(\dfrac{\partial{\mathbf{r_p(T)}}}{\partial{D_T}})^2 V_{D_T} \end{split} \notag \end{equation}$
Sim(3) Tracking keyframe:
　　当新选定一帧图像作为关键帧时，就会依据以前关键帧对新选择的关键帧上某些像素位置的深度给出估计。有了深度估计后，考虑tracking时会考虑两幅图像的光度匹配残差和新图像深度估计与warp后深度的深度残差。使用这两个残差去驱动求解、修正这两个关键帧之间的位姿变换。\
　　首先看下残差代价函数：

E (T) : = \sum x \in Ω D T [α x ( T ) σ 2 r p , x ( T ) r 2 p, x (T) + α x ( T ) σ 2 r d , x ( T ) r 2 d, x (T)]

$\begin{equation} E(\mathbf{T}):=\sum\limits_{\mathbf{x}\in\Omega D_T}[ \dfrac{\alpha_\mathbf{x}(\mathbf{T})}{\sigma^2_{r_{p,\mathbf{x}}(\mathbf{T})}}r^2_{p,\mathbf{x}}(\mathbf{T})+ \dfrac{\alpha_\mathbf{x}(\mathbf{T})}{\sigma^2_{r_{d,\mathbf{x}}(\mathbf{T})}}r^2_{d,\mathbf{x}}(\mathbf{T})] \notag \end{equation}$
其中：

r p, x (T) r d, x (T) α x σ 2 r p, x (T) σ 2 r d, x (T) : = I T (x) - I I ([x'] 1, 2) : = [x'] 3 - D I ([x'] 1, 2) : = ρ δ (r 2 p , x ( T ) σ 2 r p , x ( T ) + r 2 d , x ( T ) σ 2 r d , x ( T )) : = 2 σ 2 I + (\partial r p , x ( T ) \partial D T ( x ) ∣ ∣ ∣ D T (x) .) 2 V D T (x) : = (\partial r d , x ( T ) \partial D I ∣ ∣ ∣ D I ([x'] 1, 2) .) 2 V I ([x'] 1, 2) + (\partial r d , x ( T ) \partial D T ( x ) ∣ ∣ ∣ D T (x) .) 2 V D T (x)

$\begin{equation} \begin{split} r_{p,\mathbf{x}}(\mathbf{T})&:=I_T(\mathbf{x})-I_I([\mathbf{x'}]_{1,2})\\ r_{d,\mathbf{x}}(\mathbf{T})&:=[\mathbf{x'}]_3-D_I([\mathbf{x'}]_{1,2})\\ \alpha_{\mathbf{x}}&:=\rho_{\delta}(\dfrac{r^2_{p,\mathbf{x}}(\mathbf{T})}{\sigma^2_{r_{p,\mathbf{x}}(\mathbf{T})}}+\dfrac{r^2_{d,\mathbf{x}}(\mathbf{T})}{\sigma^2_{r_{d,\mathbf{x}}(\mathbf{T})}})\\ \sigma^2_{r_{p,\mathbf{x}}(\mathbf{T})}&:=2\sigma^2_I+ (\left. \dfrac{\partial{\mathbf{r_{p,\mathbf{x}}(T)}}}{\partial{D_T(\mathbf{x})}}\right |_{D_T(\mathbf{x})}.)^2 V_{D_T}(\mathbf{x})\\ \sigma^2_{r_{d,\mathbf{x}}(\mathbf{T})}&:= (\left. \dfrac{\partial{\mathbf{r_{d,\mathbf{x}}(T)}}}{\partial{D_I}}\right |_{D_I([\mathbf{x'}]_{1,2})}.)^2 V_I(\mathbf{[x']_{1,2}})\\ &\quad +(\left. \dfrac{\partial{\mathbf{r_{d,\mathbf{x}}(T)}}}{\partial{D_T(\mathbf{x})}}\right |_{D_T(\mathbf{x})}.)^2 V_{D_T}(\mathbf{x}) \end{split} \notag \end{equation}$
其中

x′=ωs(x,DT(x),T) $\mathbf{x'}=\omega_s(\mathbf{x},D_T(\mathbf{x}),\mathbf{T})$ 表示变换和映射过后的点，前两个坐标是其在新图像中的像素坐标，第三个是其深度值。

ρδ $\rho_\delta$ 是huber加权的权重函数。\
　　下面来看下 光度残差函数对于变换增量的导数：

J p, x = \nabla I | ω (x, T) \partial π \partial p ' ∣ ∣ ∣ T p . \partial p ' \partial T ∣ ∣ ∣ Δ T T i . \partial T \partial Δ T ∣ ∣ ∣ I . \partial e x p ( Δ ξ ^ ) \partial Δ ξ ∣ ∣ ∣ 0 .

$\begin{equation} \begin{split} &\mathbf{J_{p, \mathbf{x}}} =\nabla I|_{\omega(x,\mathbf{T})} \left. \dfrac{\partial{\pi}}{\partial{p'}} \right |_{\mathbf{Tp}}. \left. \dfrac{\partial{p'}}{\partial{T}} \right |_{\mathbf{\Delta T T_i}}. \left. \dfrac{\partial{\mathbf{T}}}{\partial{\Delta T}} \right |_{\mathbf{I}}. \left. \dfrac{\partial{exp(\hat{\Delta \mathbf{\xi}})}}{\partial{\Delta \mathbf{\xi}}} \right |_{\mathbf{0}}. \\ \end{split} \notag \end{equation}$
上式与SE(3)中等式唯一的不同是

ξ $\mathbf{\xi}$ 有7个变量，增加了变量s,其他表达式不变，但：

\partial e x p ( Δ ξ ^ ) \partial Δ ξ ∣ ∣ ∣ 0 . = ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 000000000100000000000010000000000001 0000010 - 1 0000 00 - 1 000100000 010 - 1 00000000 100010001000 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\begin{equation} \left. \dfrac{\partial{exp(\hat{\Delta \mathbf{\xi}})}}{\partial{\Delta \mathbf{\xi}}} \right |_{\mathbf{0}}.= \left( \begin{array}{ccc|ccc|c} 0&0&0&0&0&0&1\\ 0&0&0&0&0&1&0\\ 0&0&0&0&-1&0&0\\ 0&0&0&0&0&-1&0\\ 0&0&0&0&0&0&1\\ 0&0&0&1&0&0&0\\ 0&0&0&0&1&0&0\\ 0&0&0&-1&0&0&0\\ 0&0&0&0&0&0&1\\ 1&0&0&0&0&0&0\\ 0&1&0&0&0&0&0\\ 0&0&1&0&0&0&0\\ \end{array} \right) \notag \end{equation}$
与前面的计算过程一样，经过化简后得到：

J p, x = 1 z ' (\nabla I u f x \nabla I v f y) T ⎛ ⎝ ⎜ ⎜ ⎜ 1001 - x ' z ' - y ' z ' - x ' y ' z ' - (z' + y ' 2 z ') z' + x ' 2 z ' x ' y ' z ' - y' x' 00 ⎞ ⎠ ⎟ ⎟ ⎟

$\begin{equation} \begin{split} \mathbf{J_{p, \mathbf{x}}}= \dfrac{1}{z'} \left( \begin{array}{c} \nabla I_u f_x\\ \nabla I_v f_y \end{array} \right)^\mathbb{T} \left( \begin{array}{ccc|ccc|c} 1&0&-\dfrac{x'}{z'}&-\dfrac{x'y'}{z'}&z'+\dfrac{x'^2}{z'}&-y'&0\\ 0&1&-\dfrac{y'}{z'}&-(z'+\dfrac{y'^2}{z'})&\dfrac{x'y'}{z'}&x'&0 \end{array} \right) \end{split} \notag \end{equation}$
　　 光度残差函数度对于逆深度的导数与SE(3)一样：

\partial r p , x ( T ) \partial D T = 1 D T ( x ) z ' 2 (\nabla I u f x (t x z' - t z x') + \nabla I v f y (t y z' - t z y'))

$\begin{equation} \dfrac{\partial{r_{p,\mathbf{x}}(T)}}{\partial{D_T}}= \dfrac{1}{D_T(\mathbf{x})z'^2}(\nabla I_uf_x(t_xz'-t_zx')+\nabla I_vf_y(t_yz'-t_zy'))\notag \end{equation}$
　　 深度残差对于增量位姿变换的导数：

J d, x = (\partial π z \partial p ' ∣ ∣ ∣ T p . - \nabla D | ω (x, T i) \partial π \partial p ' ∣ ∣ ∣ T p .) \partial p ' \partial T ∣ ∣ ∣ Δ T T i . \partial T \partial Δ T ∣ ∣ ∣ I . \partial e x p ( Δ ξ ^ ) \partial Δ ξ ∣ ∣ ∣ 0 .

$\begin{equation} \begin{split} \mathbf{J_{d, \mathbf{x}}} =& (\left. \dfrac{\partial{\pi_z}}{\partial{p'}} \right |_{\mathbf{Tp}}. -\nabla D|_{\omega(\mathbf{x,T_i})}\left. \dfrac{\partial{\pi}}{\partial{p'}} \right |_{\mathbf{Tp}}. ) \left. \dfrac{\partial{p'}}{\partial{T}} \right |_{\mathbf{\Delta T T_i}}. \left. \dfrac{\partial{\mathbf{T}}}{\partial{\Delta T}} \right |_{\mathbf{I}}. \left. \dfrac{\partial{exp(\hat{\Delta \mathbf{\xi}})}}{\partial{\Delta \mathbf{\xi}}} \right |_{\mathbf{0}}. \\ \end{split} \notag \end{equation}$
其中：\

∂πz∂p′∣∣∣Tp.=∂1z′∂p′∣∣∣∣∣Tp.=(00−1z′2) $\left. \dfrac{\partial{\pi_z}}{\partial{p’}} \right |_{\mathbf{Tp}}. =\left. \dfrac{\partial{\dfrac{1}{z’}}}{\partial{p’}} \right |_{\mathbf{Tp}}. =\left( \begin{array}{ccc} 0&0&-\dfrac{1}{z'^2} \end{array} \right)$ \
基于性能原因，我们忽略

∇D|ω(x,Ti)∂π∂p′∣∣∣Tp. $\nabla D|_{\omega(\mathbf{x,T_i})}\left. \dfrac{\partial{\pi}}{\partial{p'}} \right |_{\mathbf{Tp}}.$
计算化简得：

J d, x = (00 - 1 z ' 2 - y ' z ' 2 x ' z ' 2 0 - 1 z ')

$\begin{equation} \mathbf{J_{d, \mathbf{x}}}= \left( \begin{array}{ccccccc} 0&0&-\dfrac{1}{z'^2}&-\dfrac{y'}{z'^2}&\dfrac{x'}{z'^2}&0&-\dfrac{1}{z'} \end{array} \right) \notag \end{equation}$
　　 深度残差对于新关键帧I逆深度的导数：

\partial r d , x ( T ) \partial D I ∣ ∣ ∣ D I ([x'] 1, 2) . = - 1

$\begin{equation} \left . \dfrac{\partial{r_{d,\mathbf{x}}(\mathbf{T})}}{\partial{D_I}} \right |_{D_I([\mathbf{x'}]_{1,2})}.=-1 \notag \end{equation}$
　　 深度残差对于模板关键帧T逆深度的导数：

\partial r d , x ( T ) \partial D T ∣ ∣ ∣ D T (x) = \partial π z \partial p ' ∣ ∣ ∣ T i p . \partial p ' \partial p ∣ ∣ ∣ p . \partial π - 1 \partial Z ∣ ∣ ∣ 1 D T ( x ) . \partial 1 D \partial D ∣ ∣ ∣ ∣ ∣ D T (x) = ⎛ ⎝ ⎜ ⎜ ⎜ 00 - 1 z ' 2 ⎞ ⎠ ⎟ ⎟ ⎟ T ⎛ ⎝ ⎜ r 11 r 21 r 31 r 12 r 22 r 32 r 13 r 23 r 33 ⎞ ⎠ ⎟ ⎛ ⎝ ⎜ x y z ⎞ ⎠ ⎟ D T (x) - 1 D T ( x ) 2 = (00 - 1 z ' 2) ⎛ ⎝ ⎜ x' - t x y' - t y z' - t z ⎞ ⎠ ⎟ 1 D T ( x ) = z ' - t z D T ( x ) z ' 2

$\begin{equation} \begin{split} \left . \dfrac{\partial{r_{d,\mathbf{x}}(\mathbf{T})}}{\partial{D_T}} \right |_{D_T(\mathbf{x})}&= \left. \dfrac{\partial{\pi_z}}{\partial{p'}} \right |_{\mathbf{T_ip}}. \left. \dfrac{\partial{p'}}{\partial{p}} \right |_{\mathbf{p}}. \left. \dfrac{\partial{\pi^{-1}}}{\partial{Z}} \right |_{\dfrac{1}{D_T(\mathbf{x})}}. \left. \dfrac{\partial{\dfrac{1}{D}}}{\partial{D}} \right |_{D_T(\mathbf{x})}\\ &= \left( \begin{array}{c} 0\\ 0\\ -\dfrac{1}{z'^2} \end{array} \right)^\mathbf{T} \left( \begin{array}{ccc} r_{11}&r_{12}&r_{13}\\ r_{21}&r_{22}&r_{23}\\ r_{31}&r_{32}&r_{33} \end{array} \right) \left( \begin{array}{c} x\\ y\\ z \end{array} \right)D_T(x) \dfrac{-1}{D_T(\mathbf{x})^2}\\ &= \left( \begin{array}{ccc} 0&0&-\dfrac{1}{z'^2} \end{array} \right) \left( \begin{array}{c} x'-t_x\\ y'-t_y\\ z'-t_z \end{array} \right)\dfrac{1}{D_T(\mathbf{x})}\\ &=\dfrac{z'-t_z}{D_T(\mathbf{x})z'^2} \end{split} \notag \end{equation}$
Storage Jacobian Parameter:
　　对于

sim(3) $sim(3)$ 的求解过程中的存储参数在这里介绍下，方便看懂作者的源代码。由于

Jp,x $\mathbf{J}_{p,\mathbf{x}}$ 在s上的导数始终为0，所以该雅可比仍然采用SE(3)的6参数格式。由于

Jd,x $\mathbf{J}_{d,\mathbf{x}}$ 只在4个位置上有数据，所以采用4参数存储格式。
　　由上一部分可知：
　　这里写图片描述

所以

Δ ξ = - [(J p, x J d, x) (J p, x J d, x)] - 1 (J p, x J d, x) (r p, x r d, x) = - [J T p, x J p, x + J T d, x J d, x] - 1 [J p, x r p, x + J d, x r d, x] = - [A p, x + A d, x] - 1 [b p, x + b d, x] = - A - 1 b

$\begin{equation} \begin{split} \Delta \mathbf{\xi}&=- [\left( \begin{array}{cc} \mathbf{J}_{p,\mathbf{x}}&\mathbf{J}_{d,\mathbf{x}} \end{array} \right) \left( \begin{array}{c} \mathbf{J}_{p,\mathbf{x}}\\ \mathbf{J}_{d,\mathbf{x}} \end{array} \right)]^{-1} \left( \begin{array}{cc} \mathbf{J}_{p,\mathbf{x}}&\mathbf{J}_{d,\mathbf{x}} \end{array} \right) \left( \begin{array}{c} r_{p,\mathbf{x}}\\ r_{d,\mathbf{x}} \end{array} \right) \\ &=-[\mathbf{J}^{\mathbf{T}}_{p,\mathbf{x}}\mathbf{J}_{p,\mathbf{x}}+ \mathbf{J}^{\mathbf{T}}_{d,\mathbf{x}}\mathbf{J}_{d,\mathbf{x}}]^{-1} [\mathbf{J}_{p,\mathbf{x}}r_{p,\mathbf{x}}+\mathbf{J}_{d,\mathbf{x}}r_{d,\mathbf{x}}]\\ &=-[\mathbf{A}_{p,\mathbf{x}}+\mathbf{A}_{d,\mathbf{x}}]^{-1} [\mathbf{b}_{p,\mathbf{x}}+\mathbf{b}_{d,\mathbf{x}}]\\ &=-\mathbf{A^{-1}b} \end{split} \notag \end{equation}$
其中

A b = A p, x + A d, x = ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ A p, 11 A p, 21 A p, 31 A p, 41 A p, 51 A p, 61 0 A p, 12 A p, 22 A p, 32 A p, 42 A p, 52 A p, 62 0 A p, 13 A p, 23 A p, 33 A p, 43 A p, 53 A p, 63 0 A p, 14 A p, 24 A p, 34 A p, 44 A p, 54 A p, 64 0 A p, 15 A p, 25 A p, 35 A p, 45 A p, 55 A p, 65 0 A p, 16 A p, 26 A p, 36 A p, 46 A p, 56 A p, 66 0 0000000 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ + ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 00000000000000 00 A d, 11 A d, 21 A d, 31 0 A d, 41 00 A d, 12 A d, 22 A d, 32 0 A d, 42 00 A d, 13 A d, 23 A d, 33 0 A d, 43 0000000 00 A d, 14 A d, 24 A d, 34 0 A d, 44 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ = b p, x + b d, x = ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ b p 1 b p 2 b p 3 b p 4 b p 5 b p 6 0 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ + ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 00 b d 1 b d 2 b d 3 0 b d 4 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\begin{equation} \begin{split} \mathbf{A}&=\mathbf{A}_{p,\mathbf{x}}+\mathbf{A}_{d,\mathbf{x}}\\ &= \left( \begin{array}{cccccc|c} \mathbf{A}_{p,11}&\mathbf{A}_{p,12}&\mathbf{A}_{p,13}&\mathbf{A}_{p,14}&\mathbf{A}_{p,15}&\mathbf{A}_{p,16}&0\\ \mathbf{A}_{p,21}&\mathbf{A}_{p,22}&\mathbf{A}_{p,23}&\mathbf{A}_{p,24}&\mathbf{A}_{p,25}&\mathbf{A}_{p,26}&0\\ \mathbf{A}_{p,31}&\mathbf{A}_{p,32}&\mathbf{A}_{p,33}&\mathbf{A}_{p,34}&\mathbf{A}_{p,35}&\mathbf{A}_{p,36}&0\\ \mathbf{A}_{p,41}&\mathbf{A}_{p,42}&\mathbf{A}_{p,43}&\mathbf{A}_{p,44}&\mathbf{A}_{p,45}&\mathbf{A}_{p,46}&0\\ \mathbf{A}_{p,51}&\mathbf{A}_{p,52}&\mathbf{A}_{p,53}&\mathbf{A}_{p,54}&\mathbf{A}_{p,55}&\mathbf{A}_{p,56}&0\\ \mathbf{A}_{p,61}&\mathbf{A}_{p,62}&\mathbf{A}_{p,63}&\mathbf{A}_{p,64}&\mathbf{A}_{p,65}&\mathbf{A}_{p,66}&0\\ \hline 0&0&0&0&0&0&0\\ \end{array} \right) + \left( \begin{array}{cc|ccc|cc} 0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0\\\hline 0&0&\mathbf{A}_{d,11}&\mathbf{A}_{d,12}&\mathbf{A}_{d,13}&0&\mathbf{A}_{d,14}\\ 0&0&\mathbf{A}_{d,21}&\mathbf{A}_{d,22}&\mathbf{A}_{d,23}&0&\mathbf{A}_{d,24}\\ 0&0&\mathbf{A}_{d,31}&\mathbf{A}_{d,32}&\mathbf{A}_{d,33}&0&\mathbf{A}_{d,34}\\\hline 0&0&0&0&0&0&0\\ 0&0&\mathbf{A}_{d,41}&\mathbf{A}_{d,42}&\mathbf{A}_{d,43}&0&\mathbf{A}_{d,44} \end{array} \right)\\ \mathbf{b}&=\mathbf{b}_{p,\mathbf{x}}+\mathbf{b}_{d,\mathbf{x}}\\ &= \left( \begin{array}{c} \mathbf{b}_{p1}\\ \mathbf{b}_{p2}\\ \mathbf{b}_{p3}\\ \mathbf{b}_{p4}\\ \mathbf{b}_{p5}\\ \mathbf{b}_{p6}\\ \hline 0\\ \end{array} \right)+ \left( \begin{array}{c} 0\\ 0\\ \mathbf{b}_{d1}\\ \mathbf{b}_{d2}\\ \mathbf{b}_{d3}\\ 0\\ \hline \mathbf{b}_{d4}\\ \end{array} \right) \end{split} \notag \end{equation}$
简单起见，这里没有考虑权重，但代码里是加了权重的。
Huber Weight Function:
　　huber函数的形式

αx:=ρδ(r2p,x(T)σ2rp,x(T)+r2d,x(T)σ2rd,x(T)) $\alpha_{\mathbf{x}}:=\rho_{\delta}(\dfrac{r^2_{p,\mathbf{x}}(\mathbf{T})}{\sigma^2_{r_{p,\mathbf{x}}(\mathbf{T})}}+\dfrac{r^2_{d,\mathbf{x}}(\mathbf{T})}{\sigma^2_{r_{d,\mathbf{x}}(\mathbf{T})}})$
　　虽然写成这样，但真正的变量是：

x=|rp,x(T)σrp,x(T)|+|rd,x(T)σrd,x(T)| $\mathbf{x}=|\dfrac{r_{p,\mathbf{x}}(\mathbf{T})}{\sigma_{r_{p,\mathbf{x}}(\mathbf{T})}}|+|\dfrac{r_{d,\mathbf{x}}(\mathbf{T})}{\sigma_{r_{d,\mathbf{x}}(\mathbf{T})}}|$ .
所以权重函数是：

α x = ω (x) = ⎧ ⎩ ⎨ 1 k | x | | x | < = k | x | > = k

$\begin{equation} \alpha_{\mathbf{x}}=\omega(\mathbf{x})= \begin{cases} 1 \quad &|\mathbf{x}|<=k\\ \dfrac{k}{|\mathbf{x}|} \quad &|\mathbf{x}|>=k \end{cases} \notag \end{equation}$
Affine Lighting Correction:
　　考虑光照变化的场景，建立了光照仿射模型:

r p, x (ξ) = a I T (x) + b - I I (x')

$\begin{equation} r_{p,\mathbf{x}}(\mathbf{\xi})=aI_{T}(\mathbf{x})+b-I_{I}(\mathbf{x'}) \notag \end{equation}$
　　这里a,b可以消除场景光照变化引起的共性光度残差，而

ξ $\mathbf{\xi}$ 可以消除由于几何运动引起的异性光度残差。引入a,b消除共性残差后，会使

ξ $\mathbf{\xi}$ 对异性残差更为敏感，在较少的迭代次数内就能求出

ξ $\mathbf{\xi}$ 的真值。
　　考虑a,b引起的共性残差的迭代最小化公式是：

E a, b (a, b) = \sum x \in Ω D T ρ a, b (a I T (x) + b - I I (x'))

$\begin{equation} E_{a,b}(a,b)=\sum\limits_{\mathbf{x}\in\Omega_{D_T}}\rho_{a,b}(aI_{T}(\mathbf{x})+b-I_{I}(\mathbf{x'})) \notag \end{equation}$
其中

ρa,b(r)=min{δmax,r2} $\rho_{a,b}(r)=min\{\delta_{max},r^2\}$ .
关于a,b的值作者给出的封闭解是这样的：

⎧ ⎩ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ a * b * = \sum x \in Ω L I T ( x ) I I ( x ' ) \sum x \in Ω L I I ( x ' ) I I ( x ' ) = 1 | Ω L | \sum i (I T (x') - a * I I (x))

$\begin{equation} \begin{cases} a^*&=\dfrac{\sum_{\mathbf{x\in\Omega_L}}I_T(\mathbf{x})I_I(\mathbf{x'})}{\sum_{\mathbf{x\in\Omega_L}}I_I(\mathbf{x'})I_I(\mathbf{x'})}\\ b^*&=\dfrac{1}{|\Omega_L|}\sum\limits_{i}(I_T(\mathbf{x'})-a^*I_I(\mathbf{x})) \end{cases} \notag \end{equation}$

ΩL:={x∈ΩT|ρa,b(aIT(x)+b−II(x′))<δmax} $\Omega_L:=\{\mathbf{x}\in\Omega_T|\rho_{a,b}(aI_{T}(\mathbf{x})+b-I_{I}(\mathbf{x'}))<\delta_{max}\}$
上面关于a,b应该是有错误的.
　　下面给出两种可能正确的解：
　　一、根据线性回归函数性质可知：

⎧ ⎩ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ a * b * = \sum x \in Ω L I T ( x ) I I ( x ' ) - n I T ( x ¯ ¯ ¯ ¯ ¯ ¯ ¯ ) I I ( x ' ¯ ¯ ¯ ¯ ¯ ¯ ¯ \sum x \in Ω L I T ( x ' ) I T ( x ' ) - n I T ( x ' ) ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ 2 = 1 | Ω L | \sum i (I I (x') - a * I T (x))

$\begin{equation} \begin{cases} a^*&=\dfrac{\sum_{\mathbf{x\in\Omega_L}}I_T(\mathbf{x})I_I(\mathbf{x'})-n\overline{I_T(\mathbf{x}})\overline{I_I(\mathbf{x'}}}{\sum_{\mathbf{x\in\Omega_L}}I_T(\mathbf{x'})I_T(\mathbf{x'})-n\overline{I_T(\mathbf{x'})}^2}\\ b^*&=\dfrac{1}{|\Omega_L|}\sum\limits_{i}(I_I(\mathbf{x'})-a^*I_T(\mathbf{x})) \end{cases} \notag \end{equation}$
　　二、依据源代码，提出a,b封闭解的一种可能假设。(仍然存在不足，仅供参考。)
给定一组样本

xi $\mathbf{x}_i$ 和每个样本的权重

pi $p_i$ ,做如下简记：

⎧ ⎩ ⎨ s x s w s x x = \sum x i \cdot p i = \sum p i = \sum x 2 i \cdot p i

$\begin{equation} \begin{cases} sx&=\sum\mathbf{x}_i\cdot p_i\\ sw&=\sum p_i\\ sxx&=\sum\mathbf{x}^2_i \cdot p_i \end{cases}\notag \end{equation}$
则这组样本的平均值是：

E(x)=∑xi⋅pi∑pi=sxsw $E(\mathbf{x})=\dfrac{\sum\mathbf{x}_i\cdot p_i}{\sum p_i}=\dfrac{sx}{sw}$
方差为：

D (x) = E [(x - E (x)) 2] = \sum ( x i - E ( x ) ) 2 \cdot p i \sum p i = \sum x 2 i \cdot p i - 2 E ( x ) \sum x i \cdot p i + E 2 ( x ) \sum p i s w = s x x - 2 E ( x ) s x + E 2 ( x ) s w s w = s x x s w - 2 (s x s w) 2 + (s x s w) 2 = s x x s w - (s x s w) 2

$\begin{equation} \begin{split} D(\mathbf{x})&=E[(\mathbf{x}-E(\mathbf{x}))^2]\\ &=\dfrac{\sum(\mathbf{x}_i-E(\mathbf{x}))^2\cdot p_i}{\sum p_i}\\ &=\dfrac{\sum\mathbf{x}^2_i\cdot p_i-2E(\mathbf{x})\sum\mathbf{x}_i\cdot p_i+E^2(\mathbf{x})\sum p_i}{sw}\\ &=\dfrac{sxx-2E(\mathbf{x})sx+E^2(\mathbf{x})sw}{sw}\\ &=\dfrac{sxx}{sw}-2(\dfrac{sx}{sw})^2+(\dfrac{sx}{sw})^2\\ &=\dfrac{sxx}{sw}-(\dfrac{sx}{sw})^2 \end{split} \notag \end{equation}$
理想情况下有：

II(x′)=aIT(x)+b $I_{I}(\mathbf{x'})=aI_{T}(\mathbf{x})+b$
根据上面的推导，令

D (I I (x')) D (I T (x)) = s y y s w - (s y s w) 2 = s x x s w - (s x s w) 2

$\begin{equation} \begin{split} D(I_{I}(\mathbf{x'}))&=\dfrac{syy}{sw}-(\dfrac{sy}{sw})^2\\ D(I_{T}(\mathbf{x}))&=\dfrac{sxx}{sw}-(\dfrac{sx}{sw})^2 \end{split} \notag \end{equation}$
由

D(II(x′))=a2⋅D(IT(x)) $D(I_{I}(\mathbf{x'}))=a^2\cdot D(I_{T}(\mathbf{x}))$ 得：

⎧ ⎩ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ a * b * = D ( I I ( x ' ) ) D ( I T ( x ) ) - - - - - - - - \sqrt = s y y s w - ( s y s w ) 2 s x x s w - ( s x s w ) 2 - - - - - - - - - - - -  ⎷     = s y y - ( s y ) 2 s w s x x - ( s x ) 2 s w - - - - - - - - - - -  ⎷      = s y s w - a * \cdot s x s w = s y - a * \cdot s x s w

$\begin{equation} \begin{cases} a^*&=\sqrt{\dfrac{D(I_{I}(\mathbf{x'}))}{D(I_{T}(\mathbf{x}))}} =\sqrt{\dfrac{\dfrac{syy}{sw}-(\dfrac{sy}{sw})^2}{\dfrac{sxx}{sw}-(\dfrac{sx}{sw})^2}} =\sqrt{\dfrac{syy-\dfrac{(sy)^2}{sw}}{sxx-\dfrac{(sx)^2}{sw}}}\\ b^*&=\dfrac{sy}{sw}-a^*\cdot \dfrac{sx}{sw} =\dfrac{sy-a^*\cdot sx}{sw} \end{cases} \notag \end{equation}$