三维物体追踪笔记(1)-基于边缘的三维物体追踪——理论、公式推导与实现
发布日期:2022-01-31 02:37:24 浏览次数:36 分类:技术文章

本文共 7154 字,大约阅读时间需要 23 分钟。

1. 基于边缘检测的三维跟踪建模

三维物体追踪是已知图像中某个物体在已知一系列空间三维点位置(或者是一个3D 模型面片集)的信息下,将这些点通过一个恰到好处的位姿(R,t)进行转换后投影到图像上。问题的求解目标是这个位姿(R,t),难点是并不知道这些三维点对应图像上的哪些像素点。

假如已知这个3D模型在图像对应物体具有清晰轮廓,特别是物体与背景之间有较好的灰度值区分度的(比如白色桌子上放一本MVG书),就可以利用基于边缘检测的三维物体跟踪。

令:x是重投影回去的轮廓(以下简称轮廓)上的离散点,x表示二维点。此处为方便起见,先命名x对应的三维点的坐标为 ( X , Y , Z ) (X,Y,Z) (X,Y,Z)

令:x在垂直于轮廓方向且向内的相邻点为 x p e x_{pe} xpe,平行于轮廓的相邻点为 x p a x_{pa} xpa

由于基于边缘检测的三维物体跟踪的最重要的特点是:垂直于边缘方向的梯度大,平行于边缘方向的梯度小,于是有:

E = 1 2 ∑ l i n e ∈ c o n t o u r ∑ x ∈ l i n e ( I ( x ) − I ( x p a ) ) 2 + ( I ( x ) − I ( x p e ) − 256 ) 2 E =\frac{1}{2} \sum_{line\in contour}\sum_{x \in line}(I(x)-I(x_{pa}))^2+(I(x) - I(x_{pe}) - 256)^2 E=21linecontourxline(I(x)I(xpa))2+(I(x)I(xpe)256)2

其中 ∣ I ( x ) − I ( x p a ) ∣ |I(x) - I(x_{pa})| I(x)I(xpa)表示梯度。

后者 ( 256 − I ( x ) − I ( x p e ) ) 2 (256-I(x) - I(x_{pe}))^2 (256I(x)I(xpe))2让候选点到边缘的梯度(即 I ( x ) − I ( x p e ) I(x) - I(x_{pe}) I(x)I(xpe))尽可能大(那样对应的E才会小)。灰度值最大值一般为256。

由于前者对结果的影响很小,因此我们考虑后者即可:

E = 1 2 ∑ l i n e ∈ c o n t o u r ∑ x i ∈ l i n e ( ( I ( x i ) − I ( x i p e ) − 256 ) 2 E =\frac{1}{2} \sum_{line\in contour}\sum_{x_i \in line}((I(x_i) - I(x_i{pe}) - 256)^2 E=21linecontourxiline((I(xi)I(xipe)256)2

**注:**在实验室常用做法中使用的是根据该marker是黑色还是白色而决定使用:

E i = 1 2 ( I ( x ) − I ( x p p e ) + 256 ) 2 ( 白 色 ) E_i = \frac{1}{2}(I(x)-I(x_{p_{pe}}) + 256) ^2(白色) Ei=21(I(x)I(xppe)+256)2

还是 E i = 1 2 ( I ( x ) − I ( x p e ) − 256 ) 2 ( 黑 色 ) E_i = \frac{1}{2}(I(x)-I(x_{pe}) - 256)^2(黑色) Ei=21(I(x)I(xpe)256)2,以便让像素差更接近与0.(我们以下假设是黑色)

2. 基于模型求导

2.1 链式分解

f i = ( I ( x i ) − I ( x i p e ) − 256 ) f_i = (I(x_i) - I(x_i{pe}) - 256) fi=(I(xi)I(xipe)256),则 E = 1 2 f T f = 1 2 [ f 1 f 2 … f n ] [ f 1 f 2 … f n ] E=\frac {1}{2}f^Tf= \frac{1}{2} \begin{bmatrix} f_1& f_2&…&f_n\end{bmatrix}\begin{bmatrix} f_1\\ f_2\\…\\f_n\end{bmatrix} E=21fTf=21[f1f2fn]f1f2fn

E的最小化问题,转换为对无数个 E i E_i Ei最小化问题求解。
已知:
f i f_i fi对李代数空间上的位姿 ϕ \phi ϕ求导:
∂ f i ∂ ϕ = J i = ∂ I ( x p e ) ∂ ϕ − ∂ I ( x ) ∂ ϕ \frac{\partial f_i}{\partial\phi} = J_i = \frac{ \partial I(x_{pe})}{\partial \phi} - \frac{ \partial I(x)}{\partial \phi} ϕfi=Ji=ϕI(xpe)ϕI(x)

那么 f f f向量整体对李代数空间上的位姿 ϕ \phi ϕ求导得到该问题的雅克比矩阵

J = [ J 1 J 2 … J n ] J = \begin{bmatrix} J_1\\ J_2\\…\\J_n\end{bmatrix} J=J1J2Jn
那么根据高斯牛顿推理可以得到E的梯度
g = J T f = [ J 1 J 2 … J n ] [ f 1 f 2 … f n ] = ∑ ( J i f i ) g = J^Tf =\begin{bmatrix} J_1& J_2&…&J_n\end{bmatrix}\begin{bmatrix} f_1\\ f_2\\…\\f_n\end{bmatrix} = \sum (J_if_i) g=JTf=[J1J2Jn]f1f2fn=(Jifi)
和E对应的海塞矩阵
H = J T J = [ J 1 J 2 … J n ] [ J 1 J 2 … J n ] = ∑ J i T J i H = J^TJ =\begin{bmatrix} J_1& J_2&…&J_n\end{bmatrix}\begin{bmatrix} J_1\\ J_2\\…\\J_n\end{bmatrix}= \sum J_i^TJ_i H=JTJ=[J1J2Jn]J1J2Jn=JiTJi

2.2 子模块计算:

(3)式的这两项的计算分别如下:

I. x处像素值对位姿求导

∂ I ( x ) ∂ ϕ = ∂ I ∂ x ∂ x ∂ ϕ \frac{ \partial I(x)}{\partial \phi} = \frac{\partial I}{\partial x} \frac{\partial x}{\partial \phi} ϕI(x)=xIϕx

$\frac{\partial I}{\partial x} $表征此处的图像梯度,由sobel求导可得;

∂ x ∂ ϕ \frac{\partial x}{\partial \phi} ϕx表示三维点对李代数的求导:

J x = ∂ x ∂ ϕ = [ − f x X Y Z 2 f x ( 1 + X 2 Z 2 ) − f x Y Z f x Z 0 − f x X Z 2 − f x ( 1 + Y 2 Z 2 ) f y X Y Z 2 f y X Z 0 f y Z − f y Y Z 2   ] J_{x} = \frac{\partial x}{\partial \phi} = \begin{bmatrix} -\frac{f_xXY}{Z^2} & f_x(1+\frac{X^2}{Z^2}) & -\frac{f_xY}{Z} & \frac{f_x}{Z} & 0 & -\frac{f_xX}{Z^2}\\ - f_x(1+\frac{Y^2}{Z^2}) &\frac{f_yXY}{Z^2} & \frac{f_yX}{Z} &0 & \frac{f_y}{Z} & -\frac{f_yY}{Z^2} \\\ \end{bmatrix} Jx=ϕx=Z2fxXYfx(1+Z2Y2) fx(1+Z2X2)Z2fyXYZfxYZfyXZfx00ZfyZ2fxXZ2fyY

其中 f x , f y f_x,f_y fx,fy是内参数。

**II. x p e x_{pe} xpe处像素值对位姿求导 **

$J_{dxpe} = \frac{ \partial I(x_{pe})}{\partial \phi} $

因为 x p e = x + d x p e x_{pe} = x + d_{x_{pe}} xpe=x+dxpe,而II中已经求出 ∂ I ( x ) ∂ ϕ \frac{ \partial I(x)}{\partial \phi} ϕI(x),所以只用对 d x p e dx_{pe} dxpe求关于 ϕ \phi ϕ的求导结果 J d x p e J_{dxpe} Jdxpe

对于某一条边,设其两端顶点为 x 1 , x 2 x_1,x_2 x1,x2,令

J x 1 = ∂ I ( x 1 ) ∂ ϕ , J x 2 = ∂ I ( x 2 ) ∂ ϕ J_{x1}=\frac{ \partial I(x_1)}{\partial \phi} ,J_{x2} = \frac{ \partial I(x_2)}{\partial \phi} Jx1=ϕI(x1),Jx2=ϕI(x2)

J d x p e = ∂ I ( d x p e ) ∂ ϕ = ∂ I ( R 2 x 2 ( − π / 2 ) d x ) ∂ ϕ = ∂ ( R 2 x 2 ( − π / 2 ) d x ) x 1 − x 2 ∣ x 1 − x 2 ∣ ) ∂ ϕ J_{dxpe} = \frac{ \partial I(dx_{pe})}{\partial \phi} = \frac{ \partial I(R_{2x2}(-\pi /2)dx)}{\partial \phi} =\frac{\partial ({R_{2x2}(-\pi /2)dx) \frac{x1-x2}{|x1-x2|}})}{\partial{\phi}} Jdxpe=ϕI(dxpe)=ϕI(R2x2(π/2)dx)=ϕ(R2x2(π/2)dx)x1x2x1x2)

其中, R 2 × 2 ( − π / 2 ) d x ) R_{2\times2}(-\pi /2)dx) R2×2(π/2)dx)是因为 x p e x_{pe} xpe的变化相当于x的变化旋转90度。

逆时针还是顺时针旋转90度?取决于x点的位置。我们可以通过x与投影到像素平面的几何中心的连线的斜率得到它是顺时针还是逆时针。

我们规定所有点投影下来,所有点先后顺序按照逆时针连接,也就是连线总是: x 1 − > x 2 x_1->x_2 x1>x2

则当 x 2 x_2 x2的斜率大于 x 1 x_1 x1时,则令 z = x 1 − x 2 z = x_1 - x_2 z=x1x2

J d x p e = R 2 x 2 ( − π / 2 ) ⋅ ∂ Z ∣ Z ∣ ∂ Z ∂ ( x 1 − x 2 ) ∂ ϕ J_{dxpe} =R_{2x2}(-\pi /2) \cdot \frac{\partial \frac{ Z}{|Z|}}{\partial Z} \frac{\partial (x1-x2)}{\partial \phi} Jdxpe=R2x2(π/2)ZZZϕ(x1x2)

其中,旋转矩阵以及 ∂ ( x 1 − x 2 ) ∂ ϕ = J x 1 − J x 2 \frac{\partial (x1-x2)}{\partial \phi} =J_{x1} - J_{x2} ϕ(x1x2)=Jx1Jx2容易求得,只需要专心求第二项。

f = ( f x , f y ) T = ∂ Z ∣ Z ∣ f =(f_x,f_y)^T= \frac{\partial Z}{|Z|} f=(fx,fy)T=ZZ表示一个分布在单位圆上的点。

如果Z为n维,则f表示落在n维“球”面上的点。

令fxx为 f x 对 Z x f_x对Z_x fxZx的求导,以此类推fxy,fyy(易证fyx = fxy);同时令 n = 1 ( x 1 2 + x 2 2 + . . . + x n 2 ) n=\frac{1}{\sqrt{(x_1^2+x_2^2+...+x_n^2)}} n=(x12+x22+...+xn2) 1

由于Z_x具有共同形式 x i ( x 1 2 + x 2 2 + . . . + x n 2 \frac{x_i}{\sqrt{(x_1^2+x_2^2+...+x_n^2}} (x12+x22+...+xn2 xi:对于这种类型的求导:

f i i = 1 ∑ x i 2 + x i ( − 1 2 ) ∑ x i 2 3 ⋅ 2 x i = n − n 3 x i 2 fii = \frac{1}{\sqrt{\sum_{x_i^2}}}+ x_i (-\frac{1}{2}) \sqrt{

{\sum_{x_i^2}}}^3 \cdot 2x_i=n-n^3x_i^2 fii=xi2 1+xi(21)xi2 32xi=nn3xi2

f i j = x i ( − 1 2 ) ∑ x i 2 3 ⋅ 2 x j = − n 3 ∗ x i ∗ x j fij = x_i (-\frac{1}{2}) \sqrt{

{\sum_{x_i^2}}}^3 \cdot 2x_j = -n^3*x_i*x_j fij=xi(21)xi2 32xj=n3xixj

所以,对于之前的两维的 Z = [ d x , d y ] T Z = [dx, dy]^T Z=[dx,dy]T,有:

∂ Z ∣ Z ∣ ∂ Z = [ n − n 3 d x 2 − n 3 d x d y − n 3 d x d y n − n 3 d y 2 ] \frac{ \partial \frac{ Z}{|Z|}}{\partial Z} = \begin{bmatrix} n-n^3d_x^2 & -n^3dxdy \\ -n^3dxdy & n-n^3dy^2 \end{bmatrix} ZZZ=[nn3dx2n3dxdyn3dxdynn3dy2]

所以: P = R 2 × 2 ( − π / 2 ) ⋅ ∂ Z ∣ Z ∣ ∂ Z = [ 0 1 − 1 0 ] ⋅ ∂ Z ∣ Z ∣ ∂ Z = [ − n 3 d x d y − ( n − n 3 d y 2 ) n − n 3 d x 2 n 3 d x d y ] P = R_{2\times 2}(-\pi /2) \cdot \frac{\partial \frac{ Z}{|Z|}}{\partial Z} = \begin{bmatrix} 0 & 1\\ -1 &0 \end{bmatrix} \cdot \frac{ \partial \frac{ Z}{|Z|}}{\partial Z} = \begin{bmatrix} -n^3dxdy & -(n-n^3dy^2) \\ n-n^3d_x^2 & n^3dxdy \end{bmatrix} P=R2×2(π/2)ZZZ=[0110]ZZZ=[n3dxdynn3dx2(nn3dy2)n3dxdy]

而由于前边我们已经解释了:

J d x p e = P ⋅ ( J x 1 − J x 2 ) J_{dxpe} = P \cdot (J_{x1} - J_{x2}) Jdxpe=P(Jx1Jx2)

这时就可以得到e对$\phi $的求导了:

J x p e = J x + J d x p e J_{xpe} = J_x + J_{dxpe} Jxpe=Jx+Jdxpe

3. 工程实现需要留意的地方

3.1 实现优化算法

LM调参对追踪效果其实是有比较大的影响,比如LM的初始 λ \lambda λ值(一般用海塞矩阵的对角元素的最大值$\tau 来 作 为 初 始 值 来作为初始值 \lambda_0 $),对结果影响比较大;另外一个需要调整的参数是收敛条件,比如当前error和上一个error的差相对于上一个error的比值。

转载地址:https://blog.csdn.net/Kevin_cc98/article/details/79582906 如侵犯您的版权,请留言回复原文章的地址,我们会给您删除此文章,给您带来不便请您谅解!

上一篇:论文阅读-位姿估计-SE3-Nets Learning Rigid Body Motion using Deep Neural Networks
下一篇:双目标定(三)标定流程(含矫正)

发表评论

最新留言

做的很好,不错不错
[***.243.131.199]2024年04月24日 14时27分30秒

关于作者

    喝酒易醉,品茶养心,人生如梦,品茶悟道,何以解忧?唯有杜康!
-- 愿君每日到此一游!

推荐文章