点云滤波(三)——EAR（边缘感知重采样）介绍_edge-aware point set resampling

一、前言

这是点云滤波系列的第三篇学习笔记，在之前我们学习了非常经典的MLS算法（点云滤波（一）——MLS（移动最小二乘）/PSS介绍），以及为解决MLS不具备保持边缘问题所提出的变体RMLS（点云滤波（二）——RMLS（鲁棒移动最小二乘）介绍）。RMLS由于没有找到相关的开源代码，并且我对方法也存在一些疑惑，害怕复现出来依托答辩，所以就没有花时间去复现。但是寻找能够保持边缘的边缘的点云平滑算法依然是我目前面对的问题，因此我会继续找一些相关的论文进行学习并形成笔记。这次阅读的文章是2013年发表于TOG上的论文Edge-aware point set resampling（EAR），作者也将算法实现发布在了CGAL库中。

二、论文摘要+引言+相关研究

看了这么多论文，感觉写摘要还是有公式可以套用的，我现在觉得大部分摘要都会包含以下几部分：

1. 一句“废话”介绍主要背景
   准确的法线估计对于表面重建和point set rendering很重要，但是法线的获取很容易受到噪声的干扰，并且锐利边缘处法线的准确计算更是一个难题
2. 现有方法存在的问题，凸显本文研究价值
   指出即使是著名的双边滤波在边缘处也会有不理想的结果
3. 指明本文的主要贡献
   提出了一个可以感知边缘并有效处理噪声和离群点的重采样方法
4. 简明扼要的指出方法的核心思想
   不直接计算边缘处的法线，而是通过平坦处可靠的法线来渐进地对边缘处进行重采样
5. 基于xxx数据集证明了方法的xxx特性（鲁棒性/实时性/高精度）
   EAR算法可以输出不受噪声影响的法线，以及被很好保持的锐利边缘
6. 方法具有实际应用价值
   EAR算法可以提升感知边缘类重建方法的以及点集渲染类技术的效果

引言部分是对摘要的展开，这篇文章的组织性很强，摘要中的每一点都在引言部分展开。

在related work部分，本文也强调了MLS、LOP，这些方法都假设下垫面是光滑的，因此无法实现边缘保持。同RMLS一样，EAR也要将下垫面处理成分段光滑表面。

需要注意EAR和MLS是重采样算法，APSS/RIMLS/Cocone/Poisson是更下游的重建算法，其中APSS/Cocone/Poisson是edge-oblivious方法，RIMLS是edge-aware方法。

三、EAR算法

3.1 Overview

EAR的核心主要有两部分：用于在远离边缘处采样的edge-aware projection operator以及用于填补边缘区域的bilateral projection operator

• Resampling away from edges
  
  由上图可以看出，经过这一步，会输出具有可靠法线的定向点，并且在锐利边缘处会留下间隙，流程如下:

  1. 使用PCA来估计法线方向；
  2. 调整法线方向使得其连续 ；
  3. 对法线应用双边平滑，以粗略的揭示边缘位置；
  4. 使用局部各向异性投影来重新采样，这样可以在不处理边缘数据的情况下可靠的估计法线。

• Edge-aware upsampling
  
  上采样过后，锐利边缘处的空隙被重新填补

  1. 先对上一步得出来的有向点做滤波；
  2. 计算待插入点的基础位置、投影方向以及投影距离；
  3. 迭代进行上采样过程，直到将边缘处的空隙填补完毕。

3.2 Resampling away from edges

这一环节的输入是包含缺陷的无序点云 Q = { q j } j ∈ J Q=\{q_j\}_{j\in J} Q={qj​}j∈J​，输出的是平滑后的远离边缘的有向点云 S = { s i } i ∈ I = { ( p i , n i ) } S=\{s_i\}_{i\in I}=\{(p_i,n_i)\} S={si​}i∈I​={(pi​,ni​)}达到的效果如下：

实现对噪声的滤除，同时将点推离边缘，并在边缘留下空隙。其中有一个重要参数${\sigma }_p$，数值越大边缘处的空隙越大，在噪声更多的情况下可以获得更好的分离效果。接下来看看这种推离是怎么实现的。

1. 使用PCA来获得法线的初值
2. 用双边滤波对法线进行调整
   $$n_i\leftarrow \frac{{\sum }_{s_i^{\prime }\in N_{s_i}}\theta (\Vert p_i-p_{i^{\prime }}\Vert )\psi (n_i,n_{i^{\prime }})n_{i^{\prime }}}{{\sum }_{s_i^{\prime }\in N_{s_i}}\theta (\Vert p_i-p_{i^{\prime }}\Vert )\psi (n_i,n_{i^{\prime }})}$$
   其中 N s i = { s i ′ ∣ s i ′ ∈ S ∧ ∥ p i − p i ′ ∥ < σ p } N_{s_i} = \{s_{i'} | s_{i'} \in S \land \|p_i - p_{i'}\| < \sigma_p\} Nsi​​={si′​∣si′​∈S∧∥pi​−pi′​∥<σp​}，$\theta (r)=e^{-r^2/{\sigma }_p^2},\psi (n_i,n_{i^{\prime }})=e^{-{\left(\frac{1-n_i^T{n}_{i^{\prime }}}{1-\cos ({\sigma }_n)}\right)}^2}$，$\theta$和$\psi$分别是欧式距离权重与法线差异权重，文中设定${\sigma }_n=15°$，$\psi$的图像如下,${\sigma }_n$控制权重衰减的位置。
   
3. 对LOP算子进行优化，获得可以感知边缘的能力，来在远离边缘处重采样
   先简单介绍下LOP算子，他的一个重要特点是无需参数化，回想MLS，需要对每个查询点的近邻拟合一个局部平面，这个过程即为参数化。但是对于一些复杂结构或者噪声混乱的情况，获取局部平面很困难。下面的第一幅图展示的是因为噪声复杂导致局部拟合平面失败，第二幅图中两个圆柱中间也不能够成功参数化。而LOP在应对更混乱的噪声以及复杂结构时候表现良好。
   
   LOP是基于优化的方法目标函数如下（为避免混乱，我将原来的符号修改成了EAR中的），其中$P^k$是第$k$次迭代时候的投影点，$P^0$是从$Q$中随机选取的子集。
   P k + 1 = G ( P k ) = arg min ⁡ P = { p i } { E 1 ( P k , Q ) + E 2 ( P k , Q ) } , E 1 ( P k , Q ) = ∑ i ∈ I ∑ j ∈ J ∥ p i k − q j ∥ θ ( ∥ p i k − q j ∥ ) , E 2 ( P k , Q ) = ∑ i ∈ I λ i ∑ i ∈ I \ { i } η ( ∥ p i k − p i ′ k ∥ ) θ ( ∥ p i k − p i ′ k ∥ ) P^{k+1}=G(P^k) = \argmin_{P=\{p_i\}} \{ E_1(P^k, Q) + E_2(P^k, Q) \},\\ E_1(P^k, Q) = \sum_{i \in I} \sum_{j \in J} \|p_i^k - q_j\| \theta(\|p_i^k - q_j\|),\\ E_2(P^k, Q) = \sum_{i\in I} \lambda_{i} \sum_{i \in I \backslash \{i\}} \eta(\|p_i^k - p_{i'}^k\|) \theta(\|p_i^k - p_{i'}^k\|) Pk+1=G(Pk)=P={pi​}argmin​{E1​(Pk,Q)+E2​(Pk,Q)},E1​(Pk,Q)=i∈I∑​j∈J∑​∥pik​−qj​∥θ(∥pik​−qj​∥),E2​(Pk,Q)=i∈I∑​λi​i∈I\{i}∑​η(∥pik​−pi′k​∥)θ(∥pik​−pi′k​∥)

• $E_1$来源于几何中位数（L1-median），到$Q$中n个点距离和最小的点$p$被称作$Q$的几何中位数，$p={\mathrm{argmin}}_x{\sum }_{j\in J}||q_j-x||$。$E_1$不同的地方在于需要找到一个点集$P$而不只是一个点 $p$，添加一个距离权重项$\theta (r)=e^{-r^2/(h/4{)}^2}$之后，相当于$p_i$只会看到他附近的一个小邻域中的点， 可以认为$p_i$是其邻域点的几何中位数。

• $E_2$的作用是为了使得$P$中的点分布均匀$\eta (r)=1/(3r^3)$是一个斥力函数（这个地方原文好像写错了，论文中笔误好像挺多的。。。。）会惩罚距离过近的投影点

  但是尽管LOP对于原始数据的几何重建效果很好，但他仍然是一个各向同性的操作子，因为距离权重$\theta (r)$不考虑锐利的几何特征。因此EAR做出的改进是考虑法线防线来推断出边缘的大致位置。具体做法是使用$\varphi (n_i,p_i-q_j)=e^{-(n_i^T(p_i-q_j){)}^2/{\sigma }_p^2}$代替$E_1$中的距离权重项。所以完整形式如下：
  P k + 1 = G ( P k ) = arg min ⁡ P = { p i } { E 1 ( P k , Q ) + E 2 ( P k , Q ) } , E 1 ( P k , Q ) = ∑ i ∈ I ∑ j ∈ J ∥ p i k − q j ∥ ϕ ( n i , p i − q j ) , E 2 ( P k , Q ) = ∑ i ∈ I λ i ∑ i ∈ I \ { i } η ( ∥ p i k − p i ′ k ∥ ) θ ( ∥ p i k − p i ′ k ∥ ) P^{k+1}=G(P^k) = \argmin_{P=\{p_i\}} \{ E_1(P^k, Q) + E_2(P^k, Q) \},\\ E_1(P^k, Q) = \sum_{i \in I} \sum_{j \in J} \|p_i^k - q_j\| \phi(n_i,p_i-q_j),\\ E_2(P^k, Q) = \sum_{i\in I} \lambda_{i} \sum_{i \in I \backslash \{i\}} \eta(\|p_i^k - p_{i'}^k\|) \theta(\|p_i^k - p_{i'}^k\|) Pk+1=G(Pk)=P={pi​}argmin​{E1​(Pk,Q)+E2​(Pk,Q)},E1​(Pk,Q)=i∈I∑​j∈J∑​∥pik​−qj​∥ϕ(ni​,pi​−qj​),E2​(Pk,Q)=i∈I∑​λi​i∈I\{i}∑​η(∥pik​−pi′k​∥)θ(∥pik​−pi′k​∥)

  $e^{-(n_i^T(p_i-q_j){)}^2/{\sigma }_p^2}$可以用来衡量$p_i$与$q_j$之间的几何差异，如下图所示$\varphi (n_p,p-q_1)=1$而$\varphi (n_p,p-q_2)$权重较小，所以最后$P$的选择会倾向于远离锐利特征的位置（但有个问题，这样操作的话不是就失去了局部性么）
  
  这里的$\eta (r)=-r$，用于平衡$E_1$和$E_2$的${\lambda }_i$的取值取决于$\mu$，LOP指出$\mu$反应了在几何精确度($\mu \in [0.1,0.25]$)和均匀分布$\mu \in [0.3,0.45]$之间的权衡，EAR中选择的是$\mu$=0.45。另一个重要参数是$h=d_{bb}/\sqrt{m}$，其中$d_{bb}$是bounding box的对角线长度，$m$是输入点云的个数，${\sigma }_p=h$。

3.3 Edge-Preserving Upsampling

上一步我们获得了边缘处留有空隙的点云，在这一步中需要将空隙填补上。每次都会插入一个新点$(p_k,n_k)$，这个点要满足三个目标：1）$p_k$位于下垫面上；2）$n_k$垂直于$p_k$处的下垫面；3）插入点要均匀的分布。直接求解需要在5D解空间中进行搜索，这非常困难，尤其是在锐利特征附近。本文的做法分三步：找到一个near-sparsest的插入基点$b_k$（目标3），优化出投影距离$d_k$（目标1），计算投影方向$n_k$（目标2），最终的投影点即为$p_k=b_k+d_k{n}_k$。

• $b_k$的选取
  这一部分参考了MLS的上采样，MLS会在局部找到一个最大外接球，用数学语言表示如下：$C(b)=mi{n}_{s_{i^{\prime }}\in {\mathcal{N}}_{s_i}}D(b,s_{i^{\prime }})$，这里的$D(b,s_{i^{\prime }})$是欧式距离。然后在其球心添加新点，但用到当前情况下会出现问题，左图所示的就是最大外接球，由于$b_k$有可能会被推出球外，所以很有可能每次找到的都是同一个点。本文的改进方式是替换欧式距离为正交距离$D(b,s_{i^{\prime }})=\Vert b-p_{i^{\prime }}-n_{i^{\prime }}^T(b-p_{i^{\prime }})n_{i^{\prime }}\Vert$，用红色向量标示出来了。
  
  
  对于$b_k$的求解，文中使用的方式是找到有着largest clearance的$s_i$与邻域点连线的中点。如何选择下一个$s_i$，本文设计了一个优先函数$P(s_i)$:
  $$P(s_i)=\underset{s_{i^{\prime }}\in {\mathcal{N}}_{s_i}}{\max }(2-n_i^T{n}_{i^{\prime }}{)}^{\rho }C(\frac{p_i+p_{i^{\prime }}}{2})$$
  当$\rho =0$时候，优先级只取决于最大的midpoint clearance。当$\rho$大于零时，靠近边缘的点因为法线变化较大而有更高的优先级。
  
• 确定投影距离$d_k$
  选好$b_k$之后，给定法线方向$n$，$d_k$的确定方式如下：
  $$d_k(b_k,n)=\frac{{\sum }_{s_i\in N_{b_k}}(n^T(b_k-p_i))\theta (\Vert b_k-p_i\Vert )\psi (n,n_i)}{{\sum }_{s_i\in N_{b_k}}\theta (\Vert b_k-p_i\Vert )\psi (n,n_i)}$$
  其实就是将$b_k$邻域点到$n$的正交距离做了一个双边加权
  
• 确定投影方向$n_k$
  一个好的投影方向应该满足以下的要求，1）投影方向符合邻域的法线分布；2）在这个投影方向上的移动距离要尽可能小
  为了提高效率，搜索会在生成$b_k$的两个端点$s_i,s_j$的法线邻域附近进行，使用投影距离来选择在在哪个邻域搜索 l = arg min ⁡ l ∈ { i , j } d k ( b k , n l ) l=\argmin_{l\in\{i,j\}}d_k(b_k,n_l) l=argminl∈{i,j}​dk​(bk​,nl​)，从下图可以一目了然的看出这种做法的合理性。
  
  然后选择最符合已有法线分布的方向，符合程度靠$f(b_k,n)={\sum }_{s_i\in N_{b_k}}\Vert n-n_i{\Vert }^2\theta (\Vert b_k-p_i\Vert )\psi (n,n_i)$来量化，选择使得$f$最小的$n$。