目录
- 0.友情链接
- 1. 论文核心思想
- 1.1. 点云特征匹配
- 1.2. 两个校验
- 1.3. 鲁棒函数与BR对偶
- 1.4.1. Black-Rangarjan Duality (BR对偶性)
- 1.4.2.Derivation of Φρ\Phi_\rhoΦρ
- 1.4.3.E(T,L)E(\bm{T},L)E(T,L)的优化求解方法
- 4.写在后面
0.友情链接
- FGR基本介绍 CSDN博客:Fast Global Registration
- FGR的github链接:
- FGR的论文下载地址
- FGR复现代码(open3d)
- FPFH论文原文
1. 论文核心思想
1.1. 点云特征匹配
\qquadFast Global Registration(FGR)是ECCV 2016的一篇关于点云全局配准的论文。所谓点云配准就是将点云A通过六自由度刚性变换矩阵T,将其变换到点云B,达到场景合并的目的,常常用于点云场景重建、SLAM、多传感器融合领域。
\qquadFGR的基本理论建立在特征匹配的基础上的,FGR使用的特征是比较快速的Hand-crafted的FPFH(Fast Point Feature Histogram,快速点云特征直方图,论文见友情链接)。FPFH相比其他learning-based的特征Recall比较低,但是优点是不需要GPU,运算速度快。
\qquad对于源点云PPP和模板点云QQQ,首先使用FPFH提取每个点的特征组成F(P)F(P)F(P)和F(Q)F(Q)F(Q),然而通过特征距离的最近邻对二者进行匹配(这点和ICP是类似的),FPFH原文的特征距离法则略为复杂,而FGR实现代码中则是采用是更简单的向量二范数计算特征距离,感兴趣的可以看下FPFH的原文。
\qquad对于F(Q)F(Q)F(Q)中的每个特征f(qi),(qi∈Q)f(q_i),(q_i\in Q)f(qi),(qi∈Q),找出F(P)F(P)F(P)中最接近f(qi)f(q_i)f(qi)的f(pj)f(p_j)f(pj),即
argminpj∈P∥f(pj)−f(qi)∥22\arg\min_{p_j\in P}{\Vert f(p_j)-f(q_i)\Vert_2^2} argpj∈Pmin∥f(pj)−f(qi)∥22
\qquad以上即是FGR特征匹配的基本思想了。需要知悉的是,整个点云配准过程,FPFH的匹配只会操作一次,这便是FGR相比其他Robust Point Cloud Registration快的原因,然而这个操作其实也占用了FGR大部分的运行时间。
1.2. 两个校验
\qquad接下来进入FGR的核心部分——outlier rejection。由于特征匹配存在较多的outlier (如下图红色连线所示)。对于correspondence-based 方法而言,这些outlier会对最后的SVD过程造成干扰,因而FGR的目的就是滤除这些红色的错误匹配。
本图摘自论文(Graduated Non-Convexity for Robust Spatial Perception: From Non-Minimal Solvers to Global Outlier Rejection)
\qquadFGR的第一校验是双向校验(dual validation),在大多数特征匹配方法里还是比较常见的。记1.1节得到的匹配好的特征为K1={(pt1,q1),(pt2,q2),...,(ptn,qn)}K_1=\lbrace(p_{t_1},q_1),(p_{t_2},q_2),...,(p_{t_n},q_n)\rbraceK1={(pt1,q1),(pt2,q2),...,(ptn,qn)},(n应该等于点云QQQ的点数),对于K1K_1K1中的每个特征对(pti,qi)(p_{t_i},q_i)(pti,qi),寻找
qi∗=argminqk∈Q∥pti−qk∥22q_i^*=\arg\min_{q_k\in Q}{\Vert p_{t_i}-q_k \Vert_2^2} qi∗=argqk∈Qmin∥pti−qk∥22
若qi∗≠qiq_i^*\neq q_iqi∗=qi,则排除这个点云对。
\qquad第一校验筛选完之后的点云对记为K2K_2K2。
\qquad第二校验是距离差校验,选取K2K_2K2中任意三个点云对(p1,p2,p3)(p_1,p_2,p_3)(p1,p2,p3)和(q1,q2,q3)(q_1,q_2,q_3)(q1,q2,q3),它们需要满足如下的不等式:
∀i≠j,0.9<∥pi−pj∥∥qi−qj∥<10.9\forall i\neq j,\quad 0.9<\frac{\Vert p_i-p_j \Vert}{\Vert q_i-q_j \Vert}<\frac{1}{0.9}∀i=j,0.9<∥qi−qj∥∥pi−pj∥<0.91
若不满足,则排除(p1,q1),(p2,q2),(p3,q3)(p_1,q_1),(p_2,q_2),(p_3,q_3)(p1,q1),(p2,q2),(p3,q3)三对点。以此得到K3K_3K3。
\qquad如果是理想情况,应该满足∥pi−pj∥∥qi−qj∥=1\frac{\Vert p_i-p_j \Vert}{\Vert q_i-q_j \Vert}=1∥qi−qj∥∥pi−pj∥=1,FGR通过一个距离范围排除那些明显不符合距离检验的点,但需要注意的是,这个检验并不是始终有效的。TEASER那篇论文就证明,当outlier ratio特别大时,距离检验也失效了,原因可能是K2K_2K2中每三组对应点中都包含outlier,导致选不出几组正确的点。
1.3. 鲁棒函数与BR对偶
\qquad这可能是FGR论文最难懂也是最核心的部分。在完成了一次匹配和两轮校验之后,FGR需要通过Gradually Non-Convexity(GNC,渐进非凸过程)的方式完成outlier rejection。经典ICP在已完成特征匹配后的目标函数是
argminT∈SE(3)∑(p,q)∈K∥q−Tp∥\arg\min_{T\in SE(3)}{\sum_{(p,q)\in K}}{\Vert q-\textbf{T}p\Vert} argT∈SE(3)min(p,q)∈K∑∥q−Tp∥
其中KKK是matched feature的集合,在FGR中指的是K3K_3K3。
\qquad然而最小二乘的函数对于outlier非常敏感,即使经过两轮校验,仍然会存在大量的outlier。
\qquadFGR的目标函数是鲁棒函数Geman-McClure,即ρ=μx2μ+x2\rho=\frac{\mu x^2}{\mu+x^2}ρ=μ+x2μx2,即
argminT∈SE(3)∑(p,q)∈Kρ(∥q−Tp∥)\arg\min_{T\in SE(3)}{\sum_{(p,q)\in K}}{\rho(\Vert q-\textbf{T}p\Vert)} argT∈SE(3)min(p,q)∈K∑ρ(∥q−Tp∥)
对于鲁棒函数的使用,FGR原文是这样解释的:
使用适当的鲁棒惩罚函数是至关重要的,因为目标函数中很多匹配约束都是伪约束。为了实现高计算效率,我们不希望在优化过程中采样、验证、修剪或重新计算关联。一个精心选择的估计器ρ将自动执行验证和剪枝,而不施加额外的计算成本。
说的直白一些,鲁棒函数可以在不改变匹配对集合KKK的情况下,使得最优解T∗T^*T∗自然地满足内点优先匹配的原则。其中的原因其实是由鲁棒函数非凸性造成的。对于outlier具有较高的、较为统一的损失,对于inlier具有较低的损失。
\qquadFGR使用的Geman-McClure是一个经典的GNC函数,随着μ\muμ的减小,越来越接近于非凸函数,而当μ→∞\mu \rightarrow \inftyμ→∞时,等价于x2x^2x2,即理想凸函数。为什么要使用GNC而不直接采用一个鲁棒函数(如Huber Loss)优化呢?因为一开始得到的T可能是个非常错误的值,就会造成∑(p,q)∈Kρ(∥q−Tp∥)\sum_{(p,q)\in K}{\rho(\Vert q-\textbf{T}p\Vert)}(p,q)∈K∑ρ(∥q−Tp∥)也是个错误的函数,因此需要逐步降低函数ρ\rhoρ的凸性逐步逼近正确的T。
\qquad承接上文,FGR为了不增加额外的计算量(其实主要是不重新计算correspondence),采用了BR对偶的方式完成鲁棒函数的优化求解。
\qquad其做法就是增加额外的线性过程L\mathcal{L}L,构造目标函数E(T,L)E(\textbf{T},\mathcal{L})E(T,L),并通过对L\mathcal{L}L和T\textbf{T}T交替优化的方式完成的推导过程及优化求解。相关截图如下:
这一套操作对于非专业人士确实非常不友好,文中的解释也只是a prior(一个先验知识),实际上它牵涉到BR对偶性的知识,想完全看懂FGR的读者务必得往下看。
1.4.1. Black-Rangarjan Duality (BR对偶性)
原文【PDF】链接
\qquad这是构造鲁棒函数和线性过程关联性的非常有名的方法,用于解决early-vision(例如平面重建、三维重建、图像恢复)滤除outlier的问题。这篇20世纪论文的引用量高达900+,足以证明它完全没有过时。
\qquad原文的理论对于非数学专业的同学可能非常难懂,本文只讲解它当中最重要的部分——定理和鲁棒函数的构造方法。定理如下:
图片源于:https://arxiv.org/pdf/1909.08605?ref=https://githubhelp.com
\qquad这条定理说的是给定一个鲁棒函数ρ\rhoρ,定义ϕ(z)=ρ(z)\phi(z)=\rho(\sqrt{z})ϕ(z)=ρ(z),若ϕ(z)\phi(z)ϕ(z)满足一定的条件,则可以将原鲁棒函数优化的问题转化为定理中的weighted process的等价形式。而这个新的二元函数的优化过程详见1.4.3节。
本文并不会证明这条定理,但需要读者知悉以下几点:
- ϕ(⋅)\phi(\cdot)ϕ(⋅)和鲁棒函数ρ(⋅)\rho(\cdot)ρ(⋅)的关系
- 等价命题的形式(ϕ(⋅)\phi(\cdot)ϕ(⋅)需满足的条件)
- 等价命题Φρ(⋅)\Phi_\rho(\cdot)Φρ(⋅)的含义(关于wiw_iwi的惩罚项)
最难理解的是Φρ\Phi_\rhoΦρ的构造过程,下面将会给出FGR中选取Φρ\Phi_\rhoΦρ(详见1.4节)的推导过程。
1.4.2.Derivation of Φρ\Phi_\rhoΦρ
Black-Rangarjan 对偶性原文给出寻找的close-form的Φρ\Phi_\rhoΦρ的算法步骤如下:
实际上,大部分论文都将τ\tauτ默认为1处理。FGR中给出的形式为
这里的lp,ql_{\text{p},\text{q}}lp,q对应原文中的zzz,ρ(x)=μx2μ+x2\rho(x)=\frac{\mu x^2}{\mu+x^2}ρ(x)=μ+x2μx2
根据ϕ(⋅)\phi(\cdot)ϕ(⋅)与ρ(⋅)\rho(\cdot)ρ(⋅)的关系
ϕ(z)=ρ(x2)=μzμ+z(z≥0)\phi(z)=\rho(\sqrt{x^2})=\frac{\mu z}{\mu + z} \quad (z\geq 0)ϕ(z)=ρ(x2)=μ+zμz(z≥0)
ϕ(z)\phi(z)ϕ(z)需要满足Black-Rangarajan对偶性的三个条件,即
{limx←0ϕ′(x)=limx←0μ2(μ+x)2=1limx←∞ϕ′(x)=limx←∞μ2(μ+x)2=0ϕ′′(x)=−2μ2(μ+x)3<0\begin{cases} & \lim_{x \leftarrow 0}\phi'(x)=\lim_{x\leftarrow 0}\frac{\mu^2}{(\mu+x)^2}=1 \\[2ex] & \lim_{x \leftarrow \infin}\phi'(x)=\lim_{x \leftarrow \infin}\frac{\mu^2}{(\mu+x)^2}=0 \\[2ex] & \phi''(x)=-\frac{2\mu^2}{(\mu+x)^3}<0 \end{cases} ⎩⎨⎧limx←0ϕ′(x)=limx←0(μ+x)2μ2=1limx←∞ϕ′(x)=limx←∞(μ+x)2μ2=0ϕ′′(x)=−(μ+x)32μ2<0
根据上述推导,显然ϕ(⋅)\phi(\cdot)ϕ(⋅)都满足了,接下来需要根据Black-Rangarajan原文中提供的方法解出Φρ\Phi_\rhoΦρ。
首先,ϕ(w)=μwμ+w\phi(w)=\frac{\mu w}{\mu + w}ϕ(w)=μ+wμw, 令z=ϕ′(w)=μ2(μ+w)2z=\phi'(w)=\frac{\mu^2}{(\mu+w)^2}z=ϕ′(w)=(μ+w)2μ2,则w=(ϕ′)−1(z)=μ/z−μw=(\phi')^{-1}(z)=\mu/\sqrt{z}-\muw=(ϕ′)−1(z)=μ/z−μ
Φρ(z)=ϕ(w)−zw=μwμ+w−zw=μ(μz−μ)μ+μz−μ−z(μz−μ)=μ(1−z−z+z)=μ(z−1)2\begin{aligned} \Phi_\rho(z) =&\phi(w)-zw\\ =&\frac{\mu w}{\mu+w}-zw\\ =&\frac{\mu(\frac{\mu}{\sqrt{z}}-\mu)}{\mu+\frac{\mu}{\sqrt{z}}-\mu}-z(\frac{\mu}{\sqrt{z}}-\mu) \\ =&\mu (1-\sqrt{z}-\sqrt{z}+z)=\mu(\sqrt{z}-1)^2 \end{aligned}Φρ(z)====ϕ(w)−zwμ+wμw−zwμ+zμ−μμ(zμ−μ)−z(zμ−μ)μ(1−z−z+z)=μ(z−1)2
即FGR原文中的鲁棒函数公式。
1.4.3.E(T,L)E(\bm{T},L)E(T,L)的优化求解方法
利用BR对偶性求解离群值过程的方法是交替优化(Alternating optimization),即对x\bm{x}x和wiw_iwi进行交替优化,具体方法截图如下。
需要说明的是,第二步wiw_iwi的优化一般是存在闭解的,针对FGR的目标函数求wiw_iwi的偏导可得:
∂E(T,L)∂wi=∥p−Tq∥2+μlp,q−1lp,q\frac{\partial{E(\bm{T},L)}}{\partial w_i}=\Vert\text{p}-\bm{T}\text{q}\Vert^2+\mu\frac{\sqrt{l_{p,q}}-1}{\sqrt{l_{p,q}}}∂wi∂E(T,L)=∥p−Tq∥2+μlp,qlp,q−1
导函数在(0,+∞)(0,+\infin)(0,+∞)上单调且存在唯一的零点,即lp,ql_{p,q}lp,q的闭合形式最优解:
lp,q=(μμ+∥p−Tq∥2)2l_{p,q}=\left(\frac{\mu}{\mu+\Vert \text{p}-\bm{T}\text{q}\Vert^2 } \right) ^2lp,q=(μ+∥p−Tq∥2μ)2
第一步对于T\bm{T}T的优化,FGR采用的是一种线性近似的方式,具体方式如下:
这个为什么属于线性近似呢?可以将这个矩阵拆成两个部分看:
M=[RT01]M=\begin{bmatrix} R &T\\ 0 &1 \end{bmatrix} M=[R0T1]
平移部分就等于它本身,旋转部分根据李群的计算法则,设
Rx=[0−γβγ0−α−βα0]R_x=\begin{bmatrix} 0 &-\gamma &\beta \\ \gamma &0 &-\alpha \\ -\beta &\alpha & 0 \end{bmatrix} Rx=⎣⎡0γ−β−γ0αβ−α0⎦⎤
R=exp(Rx)≈R0+R11!+...=[1−γβγ1−α−βα1]R=exp(R_x)\approx R^0+\frac{R^1}{1!}+...=\begin{bmatrix} 1 &-\gamma &\beta \\ \gamma &1 &-\alpha \\ -\beta &\alpha & 1 \end{bmatrix} R=exp(Rx)≈R0+1!R1+...=⎣⎡1γ−β−γ1αβ−α1⎦⎤
因此这个近似在(α,β,γ)(\alpha,\beta,\gamma)(α,β,γ)较小的时候还是可以成立的。
其中TkT^kTk是最后一次迭代的TTT,而ξ=(α,β,γ,a,b,c)\xi=(\alpha,\beta,\gamma,a,b,c)ξ=(α,β,γ,a,b,c)实际上是通过高斯-牛顿法这种数值优化方法求解的,求解方法就是上图中的(8)。
4.写在后面
\qquadFGR的全部讲解到此就结束了,实际上,FGR在实现过程逐步地减小μ\muμ以实现Graduate Non-Convexity (GNC)。将BR对偶性与GNC相结合,可以解决非凸鲁棒函数的全局最优解的问题,同时还是最优性可证的。本文只是相关理论基础的一个补充,欢迎在阅读相关文献的基础上在本文评论区留言交流学习心得。