点云滤波（二）——RMLS（鲁棒移动最小二乘）介绍

一、前言

在点云滤波（一）——MLS（移动最小二乘）/PSS介绍，我们讨论了MLS的基本原理，并给出了一个很直观的理解方式。虽然MLS可以显著的减少噪声，但我们也发现了其问题，即会将锐利边缘平滑掉。下图展示了我使用MLS在实际采集到数据上的处理结果，其中白色的点云是原始点云，红色的点云是经过MLS滤波之后的点云。

其原因在于MLS假设整个邻域处处光滑，只拟合了一个表面来表征邻域。但对于锐利边缘，这显然不合理。所以RMLS[1]对于这种锐利区域，采用了基于LMS（Least Medians of Squares）和前向搜索（forward-search paradigm）的鲁棒MLS方法来重建出分段光滑表面。下面的左图是直接使用MLS重建的结果，右图是使用RMLS重建的结果，可以明显看出RMLS能更好的保持锐利边缘。

我对这篇论文中的一些实现细节有一些疑惑，看了很久也想不明白，因此也没上手复现。网上好像也并没有复现这篇论文的工作。如果有大佬复现过或者有其他资源的话，请教教小弟！

二、关键技术

RMLS的关键点在于两个技术：最小中位数平方法（LMS）、前向搜索以及迭代重拟合（最后一个不是重点），本章对此进行介绍。

2.1 Least median of squares

2.1.1 目标函数

LMS是一种鲁棒回归方法，其崩溃点高达50%，也就是说它可以处理高达50%的异常数据点而不会导致估计结果出现巨大偏差。
$$argmi{n}_{\beta }median(f_{\beta }(x_i)-y_i{)}^2$$
但实际上论文中使用了更为鲁棒的$l_1$估计:
$$argmi{n}_{\beta }median|f_{\beta }(x_i)-y_i|$$

2.1.2 解法

文中使用随机采样的方法来对上式进行求解

1. 从邻域点中随机选择$k$个点，然后拟合出模型$\beta$
2. 然后计算剩余$N-k$个点的残差，取出中位数$r_{i,\beta }$
3. 重复计算$T$次，然后取出最小的$r_{i,\beta }$对应的$\beta$

2.1.3 局限性

尽管最小中值二乘法的崩溃点高达50%，但实际也没啥卵用。放在我们当前的表面重建任务中，尤其是遇到多个表面的问题，这个方法很容易就失效了。文中的4.1节参考了[2]中的方法，将邻域点的残差进行排序然后选出第k小的残差作为评价指标，而不是简单使用中值。
本人概率论学的一塌糊塌，因此参考的那篇文章我看不太懂，想不明白这么做的道理。并且一直以来我都觉得这种选某一个特定数值的方式过于暴力和抽象，并不是很敢用。

2.2 Forward search paradigm

2.2.1 前向与后向搜索

• 有了前向搜索，那自然要问“前向”二字体现在何处，是否有后向搜索。在我们当前面对的模型拟合问题中，前向与后向体现在数据点的选择方向上，前向是逐步添加数据；后向则是从所有的数据点开始，然后逐步删除数据。

• 论文中用了如图所示的线性回归例子来说明后向搜索存在的问题：
  
  按后向搜索的思想，首先会用所有点拟合出一条直线（图中黑线），然后计算剩余点的残差，逐步删除残差最大的点，于是图中红色的数据点就被当作离群点删除了，显然不太合适。

• 对于前向搜索，则会先从一个靠谱但却比较粗糙的拟合结果开始（这个初始解的获取是通过随机取点进行拟合，然后进行上面的LMS，迭代几次，选择指标最小的模型作为初始解），添加残差最小的点，然后对模型重新拟合（有点RANSAC的意思），接着基于新模型来观测一些指标（文中使用了residual-plot，还提到了$i$th Studentized residual, Cook’s distance 和 modified Cook distance），进而判断是否需要终止，如果不满足条件则继续添加点。示意图如下
  

• 残差与误差

  • 残差（Residual）是观测值与模型拟合值之间的差
  • 误差（Error）是观测值与真实值之间的差

  通常情况真实值都得不到，所以没办法衡量误差，只能计算残差。

2.2.2 终止条件

• 论文中指出可以通过观测几种统计量来判断前向搜索是否中止，文中举的例子是residual-plot。前向搜索会从全部邻域点$\mathcal{S}$（下图中的青色区域）中选出几个初始点$\mathcal{Q}$来拟合模型。拟合之后，计算邻域中每个点的残差，将残差最小的点添加到$\mathcal{Q}$中。不断重复上述步骤，就可以绘制出如下的residual-plot。其中纵轴代表残差的数值，横轴代表添加的点数，不同颜色的曲线代表$\mathcal{S}$中的不同点。
  
  观察residual-plot，可以总结出以下特点，在前向搜索的最开始阶段，点可以被分为两部分：高残差点和低残差点，可以发现图中存在一个明显的转折点，即高残差点残差变小，而低残差点残差变大，这个位置就代表离群点被加入到$\mathcal{Q}$中了，这时候就该停止搜索，保持上一次的估计结果作为最终结果。

• 上图是Noise-free区域，实际上在Noisy区域也有一样的特点
  

• 论文作者绘制这幅图其实只是为了说明当离群点进入$\mathcal{Q}$时候发生的变化。实际操作时并没有寻找residual-plot的拐点，作者先在整个点云中标记一个已知的平滑区域，然后拟合出该区域的多项式函数，接着测量出最大的残差值作为$r_t$（maximal tolerated residual）。将$r_t$与每次$\mathcal{Q}$中的最大残差值$r_i$比较，如果小于，则认为离群点进入了$\mathcal{Q}$。如下图是将每次$r_i$记录下来形成的曲线图，橙色线代表离群点进入的时候。
  

2.3 迭代重拟合

对于锐利边缘区域，需要拟合出分段的多项式函数。所以当$\mathcal{Q}$终止时，从$\mathcal{S}$中移除这部分点$\mathcal{S}=\mathcal{S}/\mathcal{Q}$，然后重新对$\mathcal{S}$进行拟合，直到$\mathcal{S}$为空。

三、RMLS流程

RMLS主要分两部分，一部分是要为点$x$的邻域$\mathcal{S}$拟合出一个或多个多项式，多项式的数量取决于$\mathcal{S}$中平滑区域的个数；另一部分是拓展了MLS中的投影算子，使得$x$可以准确的投影到合适的表面上。

3.1 local classification

• 实际就是用了LMS和前向搜索对$\mathcal{S}$中的点进行分类，区分出属于不同表面上的点。
• 其实有了上边的讲解，论文中第四章的算法流程读起来就没有太大障碍了，流程如下：
  

但还是要对其中的一些细节进行说明：

1. $p$是模型参数的个数（自由度），对于直线需要两个参数，平面则需要三个参数，在LMS随机选点数$k$可以设为$p$，太小的$k$值没有充分利用已有信息；太大的$k$值又会引入噪声并增加迭代次数；
2. 如果想要为3D空间的数据点拟合多项式表面，需要在参数域下进行，这个其实在MLS中有提到，当时会先拟合一个平面，然后将邻域点投影到这个平面上，获得他们的2D坐标，这个就是图中的Parameter Domain；
3. 2.1.2中写到，为了求出满足最小化LMS的$\beta$，我们需要在点集中随机选出几个点，作者指出可以利用几何先验的知识来优化选点的过程，进而加速整个进程并且提升稳定性。他的出发点是这样，因为表面一定是由点和他的邻域组成的，所以选候选点的时候要集中在一个邻域里。这个思想被应用在了两个地方，一个时LMS的初始化选点上；另一个是前向搜索添加候选点中，也就是图中的ImmediateNeighbors。

3.2 Projection operator

与MLS不同，RMLS为$x$的邻域拟合出了分段平滑表面，所以$x$往哪里投影就是一个问题。作者这部分写的并没有很详细，接下来的某些内容只是我自己的理解，欢迎大家多多交流。

投影流程：

1. 判断邻域是否平滑：如果平滑则直接按照MLS的方式来对$x$进行投影，如果不平滑则按本文的方式来迭代重拟合出多个表面。判断是否平滑，论文中的方法是对整个邻域做一次拟合，然后计算最大残差$r_{max}$，如果$r_{max}>r_t$，则认为该区域不平滑（直接计算曲率可能会更简单准确）。
2. 对于多个表面的情况如何投影：
   论文中有这样一句话我并不是很懂，“When we identify two or more subsets, the surface is defined as the intersection of the smooth surfaces defined by the different neighborhoods.”。我猜测这里的intersection指的是图中黑色实线部分
   
   将$x$点向平滑表面上做投影，如果投影点都位于实线区域上，则选取距离最近的投影点，图中的$x_1$；如果一个投影点位于实线上，另一个位于虚线上，则保留实线上的投影点；最后一种情况，如果两个点都落在虚线上，那就投影到表面的交线上，这部分作者参考了[3]中的方法；下面是这三个投影方式对应区域的示意图。
   
   红色区域具体投影方式如下：
   
   (b)图中的$q_1$和$q_2$是红色区域点到两个MLS表面上的投影点，计算在这两点上MLS表面的切平面，找到两平面交线上使得$\Vert r-q_1\Vert -\Vert r-q_2\Vert$最小的点$r$，也就是图上的绿色点。将$r$点向两个表面再做投影，重复上述步骤，直到$q_1$和$q_2$之间的距离小于设定的阈值，迭代结束。

四、参考文献

[1]Fleishman S, Cohen-Or D, Silva C T. Robust moving least-squares fitting with sharp features[J]. ACM transactions on graphics (TOG), 2005, 24(3): 544-552.

[2]Miller J V, Stewart C V. MUSE: Robust surface fitting using unbiased scale estimates[C]//Proceedings CVPR IEEE Computer Society Conference on Computer Vision and Pattern Recognition. IEEE, 1996: 300-306.

[3]Pauly M, Keiser R, Kobbelt L P, et al. Shape modeling with point-sampled geometry[J]. ACM Transactions on Graphics (TOG), 2003, 22(3): 641-650.