视觉SLAM总结——视觉特征子综述

视觉特征子综述

这篇文章是一篇2D特征子和3D特征子的综述，每个特征子的介绍都是以性质——算法简述——实验的构架介绍的，性质主要是简单总结这个特征子的特点，算法简述是参考了其他博客和文章的一个总结，实验是基于PCL和其他博客的一些Demo实现的一些简单Demo（这里特地感谢下我的同门提供了部分代码），实验的代码在我的Github上，本人水平有限，如果文中有任何错误，欢迎大家指正。

需要明确的一个概念是特征描述子是由关键点(keypoints)+特征描述子(feature)两部分构成，下面所介绍的算法中，一部分是同时介绍了关键点和特征描述子的算法，比如SIFT，另一部分是更加强调关键点或者特征描述子的提取算法，比如Harris，ISS等，偏重是不一样的。

第一部分：2D特征子

1. Harris

性质
算法的核心是利用局部窗口在图像上进行移动判断灰度发生较大的变化，其性质主要是
（1）阈值决定检测点数量，
（2）Harris角点检测算子对亮度和对比读变换不敏感，
（3）Harris角点检测具有旋转不变性（通过特征值判定的优势，可以理解为椭圆的长短轴与方向无关），
（4）Harris角点检测不具有尺度不变性。

算法简述

如上图所示：算法描述像素点平移后的窗口内的灰度变换是通过自相关函数给出，其中 $w(x,y)$为窗口加权函数。
$$c(x,y;\Delta x,\Delta y)=\sum _{(u,v)\in W(x,y)}w(u,v)(I(u,v)–I(u+\Delta x,v+\Delta y){)}^2$$通过一阶近似并简化为 $$c(x,y;\Delta x,\Delta y)\approx \sum _w(I_x(u,v)\Delta x+I_y(u,v)\Delta y{)}^2=[\Delta x,\Delta y]M(x,y)\left[\begin{array}{c}\Delta x\\ \Delta y\end{array}\right]$$其中 $$M(x,y)=\sum _w\left[\begin{array}{cc}{I}_x(x,y{)}^2& I_x(x,y)I_y(x,y)\\ I_x(x,y)I_y(x,y)& I_y(x,y{)}^2\end{array}\right]$$而上述方程可以表示成一个椭圆， $M(x,y)$的两个特征值决定了椭圆的形状，如下

根据特征值可以将角点的判定分为三种情况：
（1）图像中的直线。一个特征值大，另一个特征值小， $\lambda 1\gg \lambda 2$或 $\lambda 2\gg \lambda 1$。自相关函数值在某一方向上大，在其他方向上小。
（2）图像中的平面。两个特征值都小，且近似相等；自相关函数数值在各个方向上都小。
（3）图像中的角点。两个特征值都大，且近似相等，自相关函数在所有方向都增大。
根据以上结论设计了一个响应函数如下，进行角点的阈值判定 $$R=det\mathbf{M}-\alpha (trace\mathbf{M}{)}^2$$最后在邻域内进行非极大值抑制即获得图像角点

实验
实验结果如下图所示：

图中，绿色点为初始提取的的Harris角点，红色空心点为通过亚像素处理后的角点，仔细观察后，红色空心点确实位置确实更加准确。

2. SIFT（Scale-invariant feature transform）

性质
算法的核心是高斯差分以及邻域内的直方图统计，其性质主要如下：
（1）SIFT特征是图像的局部特征，对旋转、尺度缩放、亮度变化保持不变性，对视角变化、仿射变换、噪声也保持一定程度的稳定性（其中尺度不变性主要是通过尺度空间保证的，旋转不变性主要是通过主方向保证的）
（2）信息量丰富（128维特征），适用于在海量特征数据库中进行匹配
（3）实时性不高，对边缘光滑的目标无法准确提取特征点

算法简述
（1）第一步：检测尺度空间的极值点： 首先在使用高斯金字塔构建尺度空间时，主要分成两部分，对图像做降采样（获得Octave），以及对图像做不同尺度的高斯模糊（获得Interval）。通过同一Octave内相邻的Interval作差即可高斯差分图像，如下图所示：

公式为 $$D(x,y,\sigma )=(G(x,y,k\sigma )-G(x,y,\sigma ))\ast I(x,y)=L(x,y,k\sigma )-L(x,y,\sigma )$$其中 $$G(x,y,\sigma )=\frac{1}{2\pi {\sigma }^2}{e}^{-\frac{x^2+y^2}{2{\sigma }^2}}$$这里有一点值得注意的是，在降采样时，高斯金字塔中一Octive的底层图像是由前一组图像的倒数第3张图像(S+1层)隔点采样得到。这样也保证了尺度变化的连续性。如下图所示取（ $k=1/3$，同时注意下图中DoG空间和高斯空间的区别)

关键点的检测见识检测DOG空间内同一个Octive领域26个像素的极值点，如下图所示，获得极值点后会通过（1）非极大值抑制（2）海森矩阵去除边界点进行不稳定点的去除。

（2）第二步：精确定位特征点（亚像素）
SIFT是采用插值法进行亚像素的定位，对DoG函数进行泰勒展开有$$D(X)=D+\frac{\partial D^T}{\partial X}\Delta X+\frac{1}{2}\Delta X^T\frac{{\partial }^{2}D}{\partial X^2}\Delta X$$$D(x)$相当于获得了一个对于$\Delta X$相关的函数，对$\Delta X$求导等于零即可获得极值点的偏移量$$\hat{X}=-(\frac{{\partial }^{2}D}{\partial X^2}{)}^{-1}\frac{\partial D}{\partial X}$$通过迭代最后获得的精确关键点为：$$D(\hat{X})=D+\frac{1}{2}\frac{\partial D^T}{\partial X}\hat{X}$$
（3）第三步：设定主方向
对于通过上述方法获得的关键点，在其邻近高斯金字塔图像的$3\sigma$邻域窗口内计算所有像素点梯度幅值和方向，公式如下$$m(x,y)=\sqrt{[L(x+1,y)-L(x-1,y){]}^2+[L(x,y+1)-L(x,y-1){]}^2}$$$$\theta (x,y)=ta{n}^{-1}[L(x,y+1)-L(x,y-1)]/[L(x+1,y)-L(x-1,y)]$$模值$m(x,y)$按$\sigma =1.5{\sigma }_{oct}$、邻域窗口为$3\sigma =3\times 1.5\sigma oct$的高斯分布加权,在完成关键点的梯度计算后，使用直方图统计邻域内像素的梯度和方向，梯度直方图将梯度方向$(0,360^{\circ })$分为$36$柱(bins)，如下图所示，直方图的峰值所在的方向代表了该关键点的主方向。

（4）第四步：生成描述子
在经过上述流程后，检测到的每个关键点有三个信息： 位置、尺度以及方向，接下来要做的就是抽象出一组特征向量描述每个关键点。SIFT特征描述子是关键点 邻域梯度经过高斯加权后统计结果的一种表示。通过对关键点周围图像区域分块，计算块内的梯度直方图，得到表示局部特征点信息的特征向量。
 特征描述子与特征点所在的尺度空间有关，因此对梯度的求取应该在特征点对应的高斯图像上进行。将关键点附近的邻域划分为 $d\times d(d=4)$个子区域，每个子区域的大小与关键点方向分配时相同，即边长为 $3{\sigma }^{oct}$。
下一步特别重要，保证了特征的旋转不变形，就是 旋转坐标轴至关键点的主方向，旋转后采样点的新坐标为 $$\left[\begin{array}{c}x’\\ y’\end{array}\right]=\left(\begin{array}{cc}cos\theta & -sin\theta \\ sin\theta & cos\theta \end{array}\right)\left[\begin{array}{c}x\\ y\end{array}\right],x,y\in [-radius,radius]$$

如上图所示，原本是按照左图进行采样，经过主方向旋转之后按照右图方向所示进行采样。将旋转后的采样点坐标分配到对应的子区域，计算影响子区域的采样点的梯度和方向，分配到8个方向上。

得到128维特征向量后，为了去除光照变化的影响，需要对向量进行归一化处理。非线性光照变化仍可能导致梯度幅值的较大变化，但对梯度方向影响较小。因此对于超过阈值0.2的梯度幅值设为0.2，然后再进行一次归一化。最后将特征向量按照对应高斯金字塔的尺度大小排序。至此，SIFT特征描述子形成。
参考:
图像特征之SIFT特征匹配

实验
实验结果如下图所示：

如图所示，可以看出SIFT的匹配是比较准确的，几乎所有的匹配线都是平行线。

3. SURF（Speed up robust feature）

性质
SURF算法是对SIFT算法加强版，同时加速的具有鲁棒性的特征，其同样具有尺度、亮度和旋转不变性。SURF的核心在于采用了harr特征以及积分图像integral image的概念，这大大加快了程序的运行时间。

算法简述
SURF的算法流程和SIFT几乎是一致的，只是中间某些模块使用了更加鲁棒高效的算法
（1）尺度空间的极值检测
首先要讨论下SURF的核心积分图像的概念，首先对于一个图像$I(x,y)$,其Hessen矩阵如下：$$H(I(x,y))=\left[\begin{array}{cc}\frac{{\partial }^{2}I}{\partial x^2}& \frac{{\partial }^{2}I}{\partial x\partial y}\\ \frac{{\partial }^{2}I}{\partial x\partial y}& \frac{{\partial }^{2}I}{\partial y^2}\end{array}\right]$$$H$矩阵的判别式是$$Det(H)=\frac{{\partial }^{2}I}{\partial x^2}\ast \frac{{\partial }^{2}I}{\partial y^2}-\frac{{\partial }^{2}I}{\partial x\partial y}\ast \frac{{\partial }^{2}I}{\partial x\partial y}$$在构建Hessian矩阵前需要对图像进行高斯滤波，经过滤波后的Hessian矩阵表达式为：$$H(x,y,\sigma )=\left[\begin{array}{cc}{L}_{xx}(x,y,\sigma )& L_{xy}(x,y,\sigma )\\ L_{xy}(x,y,\sigma )& L_{yy}(x,y,\sigma )\end{array}\right]$$其中$(x,y)$为像素位置,$L(x,y,\sigma )=G(\sigma )\ast I(x,y)$,代表着图像的高斯尺度空间，是由图像和不同的高斯卷积得到。在离散数学图像中，一阶导数是相邻像素的灰度差$$L_x=L(x+1,y)-L(x,y)$$二阶导数是对一阶导数的再次求导$$L_{xx}=[L(x+1,y)-L(x,y)]-[L(x,y)-L(x-1,y)]$$$$=L(x+1,y)+L(x-1,y)-2L(x,y)$$反过来看Hessian矩阵的判别式，其实就是当前点对水平方向二阶偏导数乘以垂直方向二阶偏导数再减去当前水平、垂直二阶偏导的二次方：$$Det(H)=L_{xx}\ast L_{yy}-L_{xy}\ast L_{xy}$$通过这种方法可以为图像中每个像素计算出其H行列式的决定值，并用这个值来判别图像局部特征点。Hesson矩阵判别式中的$L(x,y)$是原始图像的高斯卷积，由于高斯核服从正太分布，从中心点往外，系数越来越小，为了提高运算速度，SURF算法使用了盒式滤波器来替代高斯滤波器$L$，所以在$L_{x}y$上乘了一个加权系数0.9，目的是为了平衡因使用盒式滤波器近似所带来的误差，则H矩阵判别式可表示为：$$Det(H)=L_{xx}\ast L_{yy}-(0.9\ast L_{xy}{)}^2$$盒式滤波器和高斯滤波器的示意图如下：

那么为什么盒式滤波器可以提高运算速度呢？这就涉及到积分图的使用， 盒式滤波器对图像的滤波转化成计算图像上不同区域间像素的加减运算问题，这正是积分图的强项，只需要简单积分查找积分图就可以完成。
所谓积分图就是只遍历一次图像就可以求出图像中所有区域像素和的快速算法，大大的提高了图像特征值计算的效率。积 分图主要的思想是将图像从起点开始到各个点所形成的矩形区域像素之和作为一个数组的元素保存在内存中，当要计算某个区域的像素和时可以直接索引数组的元素，不用重新计算这个区域的像素和，从而加快了计算（这有个相应的称呼，叫做 动态规划算法）。积分图能够在多种尺度下，使用相同的时间（常数时间）来计算不同的特征，因此大大提高了检测速度。

然后就是构建尺度空间，SURF的尺度空间构建和SIFT是不同的，如下图所示：

SIFT算法 下一组图像的长宽均是上一组的一半，同一组不同层图像之间尺寸一样，但是所使用的尺度空间因子(高斯模糊系数σ)逐渐增大；而在SURF算法中， 不同组间图像的尺寸都是一致的，不同的是不同组间使用的盒式滤波器的模板尺寸逐渐增大，同一组不同层图像使用相同尺寸的滤波器，但是滤波器的尺度空间因子逐渐增大, SURF的尺度空间构建方式省去了下采样的过程，因此提升了处理速度

（2）特征点过滤并进行精确定位。
SURF特征点的定位过程和SIFT算法一致，将经过Hessian矩阵处理的每个像素点(即获得每个像素点Hessian矩阵的判别式值)与其图像域（相同大小的图像）和尺度域（相邻的尺度空间）的所有相邻点进行比较**（注意这里就不再是高斯差分图像），当其大于（或者小于）所有相邻点时，该点就是极值点。初步定位出特征点后，再经过滤除能量比较弱的关键点以及错误定位的关键点，筛选出最终的稳定的特征点。

（3）设定主方向
SURF算法中，采用的是统计特征点圆形邻域内的Harr小波特征，即在特征点的圆形邻域内，统计60度扇形内所有点的水平、垂直Harr小波特征**总和，然后扇形以0.2弧度大小的间隔进行旋转并再次统计该区域内Harr小波特征值之后，最后将值最大的那个扇形的方向作为该特征点的主方向。该过程示意图如下

Harr特征分为三类：边缘特征、线性特征、中心特征和对角线特征，组合成特征模板，如下图所示。特征模板内有白色和黑色两种矩形，并定义该模板的特征值为白色矩形像素和减去黑色矩形像素和。Harr特征值反映了图像的灰度变化情况。

（4）生成描述子
SURF算法中，也是提取特征点周围4×4个矩形区域块，但是所取得矩形区域方向是沿着特征点的主方向。每个子区域统计25个像素点水平方向和垂直方向的Haar小波特征，这里的水平和垂直方向都是相对主方向而言的。该Harr小波特征为水平方向值之和、垂直方向值之和、水平方向值绝对值之和以及垂直方向绝对之和4个方向。该过程示意图如下：

把这4个值作为每个子块区域的特征向量，所以一共有4×4×4=64维向量作为SURF特征的描述子，比SIFT特征的描述子减少了一半。这里还有个区别可能是，SIFT中仅仅是对坐标轴进行了旋转，采样点随之进行了旋转，但是计算采样点的梯度方向是还是按照原来的坐标系进行描述，而SURF中的Harr小波模块是随着坐标轴旋转到了新的坐标系下面进行特征的计算。
参考
Surf特征提取分析
SURF特征提取原理详细分析及opencv API调用

实验
实验结果如下图所示：

如图所示，SURF的匹配效果也是比较好的

4. ORB（Oriented FAST and rotated BRIEF）

性质
ORB特征是将FAST特征点的检测方法与BRIEF特征描述子结合起来，并在它们原来的基础上做了改进与优化。ORB的最大特点就是速度快，且ORB的特征描述子具备旋转和尺度不变性。

算法流程
（1）oFast关键点提取
ORB算法的特征提取是由FAST算法改进的，这里成为oFAST（FASTKeypoint Orientation）。也就是说，在使用FAST提取出特征点之后，给其定义一个特征点方向，以此来实现特征点的旋转不变形。FAST算法是公认的最快的特征点提取方法。FAST算法提取的特征点非常接近角点类型。oFAST算法如下：

步骤一：粗提取。该步能够提取大量的特征点，但是有很大一部分的特征点的质量不高。下面介绍提取方法。从图像中选取一点P，如上图。我们判断该点是不是特征点的方法是，以P为圆心画一个半径为3pixel的圆。圆周上如果有连续n个像素点的灰度值比P点的灰度值大或者小，则认为P为特征点。一般n设置为12。为了加快特征点的提取，快速排出非特征点，首先检测1、9、5、13位置上的灰度值，如果P是特征点，那么这四个位置上有3个或3个以上的的像素值都大于或者小于P点的灰度值。如果不满足，则直接排出此点。

步骤二：机器学习的方法筛选最优特征点。简单来说就是使用ID3算法训练一个决策树，将特征点圆周上的16个像素输入决策树中，以此来筛选出最优的FAST特征点。

步骤三：非极大值抑制去除局部较密集特征点。使用非极大值抑制算法去除临近位置多个特征点的问题。为每一个特征点计算出其响应大小。计算方式是特征点P和其周围16个特征点偏差的绝对值和。在比较临近的特征点中，保留响应值较大的特征点，删除其余的特征点。

步骤四：特征点的尺度不变形。建立金字塔，来实现特征点的多尺度不变性。设置一个比例因子scaleFactor（opencv默认为1.2）和金字塔的层数nlevels（pencv默认为8）。将原图像按比例因子缩小成nlevels幅图像。nlevels幅不同比例的图像提取特征点总和作为这幅图像的oFAST特征点。

步骤五：实现特征点的旋转不变性。ORB算法提出使用矩（moment）法来确定FAST特征点的方向。也就是说通过矩来计算特征点以r为半径范围内的质心，特征点坐标到质心形成一个向量作为该特征点的方向。矩定义如下：$$m_{pq}=\sum _{x,y\in r}{x}^p{y}^{q}I(x,y)$$其中，$I(x,y)$为图像表达式。该矩的质心为：$$C=(\frac{m_{10}}{m_{00}},\frac{m_{01}}{m_{00}})$$假设关键点坐标为$O$，则向来的角度即为该特征点的方向。计算公式如下$$\Theta =arctan(\frac{m_{10}}{m_{00}}/\frac{m_{01}}{m_{00}})=arctan(m_{01}/m_{10})$$

（2）rBRIEF特征描述子
BRIEF算法计算出来的是一个二进制串的特征描述符。它是在一个特征点的邻域内，选择n对像素点$p_i、q_i（i=1,2,\dots ,n）$。然后比较每个点对的灰度值的大小。如果I(pi)> I(qi)，则生成二进制串中的1，否则为0。所有的点对都进行比较，则生成长度为n的二进制串。一般n取128、256或512，opencv默认为256。另外，值得注意的是为了增加假设原始的BRIEF算法在特征点$S\times S$（一般S取31）邻域内选取n对点集。特征描述符的抗噪性，算法首先需要对图像进行高斯平滑处理。在ORB算法中，在这个地方进行了改进，在使用高斯函数进行平滑后，又用了其他操作，使其更加的具有抗噪性。

如上图所示分别是五种不同的点对选取方式，分别是（1）在图像块内平均采样；（2）p和q都符合(0,S2/25)的高斯分布；（3）p符合(0,S2/25)的高斯分布，而q符合(0,S2/100)的高斯分布；（4）在空间量化极坐标下的离散位置随机采样；（5）把p固定为(0,0)，q在周围平均采样。通过实验比较第二种匹配效果最好。

rBRIEF就是在原始的算法上进行了旋转不变性的改进，即在使用oFast算法计算出的特征点中包括了特征点的方向角度。假设原始的BRIEF算法在特征点SxS（一般S取31）邻域内选取n对点集。$$D=((x_1,y_1),(x_2,y_2)(x_3,y_3)...(x_{2n},y_{2n}))$$经过旋转角度$\theta$旋转后，得到新的点对$$D_{\theta }=R_{\theta }D$$在新的点集位置上比较点对的大小形成二进制串的描述符。这里需要注意的是，在使用oFast算法是在不同的尺度上提取的特征点。因此，在使用BRIEF特征描述时，要将图像转换到相应的尺度图像上，然后在尺度图像上的特征点处取$S\times S$邻域，然后选择点对并旋转，得到二进制串描述符。

rBRIEF同时增强了描述符的可区分性（不相关性），ORB论文中没有使用5种方法中的任意一种，而是使用统计学习的方法来重新选择点对集合。

首先建立300k个特征点测试集。对于测试集中的每个点，考虑其31x31邻域。这里不同于原始BRIEF算法的地方是，这里在对图像进行高斯平滑之后，使用邻域中的某个点的5x5邻域灰度平均值来代替某个点对的值，进而比较点对的大小。这样特征值更加具备抗噪性。另外可以使用积分图像加快求取5x5邻域灰度平均值的速度。

从上面可知，在31x31的邻域内共有(31-5+1)x(31-5+1)=729个这样的子窗口，那么取点对的方法共有M=265356种，我们就要在这Ｍ种方法中选取256种取法，选择的原则是这256种取法之间的相关性最小。怎么选取呢？

（1）在300k特征点的每个31x31邻域内按M种方法取点对，比较点对大小，形成一个300kxM的二进制矩阵Q。矩阵的每一列代表300k个点按某种取法得到的二进制数。
（2）对Q矩阵的每一列求取平均值，按照平均值到0.5的距离大小重新对Q矩阵的列向量排序，形成矩阵T。
（3）将T的第一列向量放到R中。
（4）取T的下一列向量和R中的所有列向量计算相关性，如果相关系数小于设定的阈值，则将T中的该列向量移至R中。
（5）按照（4）的方式不断进行操作，直到R中的向量数量为256。

通过这种方法就选取了这256种取点对的方法。这就是rBRIEF算法。
参考：
ORB特征提取详解

实验
实验结果如下图所示：

如图所示，ORB算法虽然快，但是匹配效果确实差了很多。

5. 一部分思考总结

在前文所介绍的四种2D特征子中，Harris是最为简单的，仅仅具备旋转不变形而不具备尺度不变形。而SIFT、SURF和ORB均具备旋转和尺度不变性，旋转不变性通常是通过求得关键点的主方向获得的，而尺度不变形是通过高斯金字塔获得的。另外，旋转和尺度不变性应该都是针对特征描述子来说的，关键点的提取在不同尺度就应该是不同的关键点，而对于旋转不变性来说，关键点都是通过邻域的阈值特性确定的，所以即使图片旋转了也不会影响到关键点的提取。

第二部分：3D特征子

在接下来介绍的八种3D特征子中，前四种主要是关于关键点的提取算法，后四种主要是关于特征描述子的提取算法，前后是可以相互结合的。例如你可以提取3D SIFT的关键点然后用PFH对关键点进行特征描述，然后再进行匹配。

1. 3D Harris

性质
3D Harris与2D Harris具备同样的性质，他们都不具备尺度不变性，是一种比较简单基础的特征点检测算法。

算法简述
算法的推到过程和2D Harris是一致的，这里主要来分析一下区别，再2D Harris中主要是对图像灰度的协方差矩阵进行特征值求解$$M(x,y)=\sum _w\left[\begin{array}{cc}{I}_x(x,y{)}^2& I_x(x,y)I_y(x,y)\\ I_x(x,y)I_y(x,y)& I_y(x,y{)}^2\end{array}\right]$$如果两个特征值都大，证明图像在两个主方向上灰度变化程度都大，因此是角点。

在点云中，没有灰度的概念，因此假定点云密度是均匀的，采用小正方体内的点云数量变化程度作为角点判定标准，类似2D Harris的推到过程可以获得如下梯度协方差矩阵。$$M(x,y)=\sum _w\left[\begin{array}{ccc}{I}_x(x,y{)}^2& I_x(x,y)I_y(x,y)& I_x(x,y)I_z(x,y)\\ I_x(x,y)I_y(x,y)& I_y(x,y{)}^2& I_y(x,y)I_z(x,y)\\ I_x(x,y)I_z(x,y)& I_y(x,y)I_z(x,y)& I_z(x,y{)}^2\end{array}\right]$$同样对特征值进行判定，某个特征值大则证明再该特征值对应的主方向上密度梯度变换明显，如果三个特征值都大，则证明该点是角点，当然，在实际操作中3D Harris中角点的定义应该情况更加复杂，最后的判定方式需要阅读源码才能知道。

实验
实验结果如下图所示：

如图所示，根据Harris主要是根据点云变化的梯度情况来判定是否为角点，因此Harris主要分布在变换较为剧烈的位置，例如耳朵尖尖。

2. 3D SIFT（3D Scale-invariant feature transform）

性质
3D SIFT同2D SIFT一样具有尺度和旋转不变性，对视角变化，仿射变换和噪声等保持一定程度的稳定性。

算法简述
3D SIFT同2D SIFT的算法流程基本一致，只是再2D图像和3D点云的处理细节上稍有些不同，比如2D高斯金字塔换成3D高斯金字塔，如下图所示：

再比如，再2D SIFT算法中特征描述子的窗口大小是16×16， 特征向量的维度是4×4×8=128维，而在3D SIFT中窗口大小维16×16×16，特征向量维度是4×4×4×8×8=4096维，算法流程上在此就不再赘述了。

实验
实验结果如下图所示：

SIFT角点也是主要分布在变化相对剧烈的部分，SIFT计算时间相对也较长。

3. ISS（Intrinsic shape signatures）

性质
通过构建关键点的局部坐标系来实现特征的描述，进而提取关键点，思路比较简单。

算法简述
（1）通过协方差矩阵建立模型
其思路其实就是通过协方差矩阵分析某点周围点云的分布情况，步骤如下：
1）把$p_i$和周围点$p_j$的坐标相减：本质上这生成了许多从$p_i->p_j$的向量，理想情况下pi的法线应该是垂直于这些向量的
2）构建这些向量的协方差矩阵，矩阵的特征值所对应的特征向量就是这个邻域点云的三个垂直的主方向（其中一个是法向量）

（2）利用特征值之间的关系提取特征点
这里只定性的解释，如果上述协方差的三个特征值大小相差不大，意味着三个方向的点云分布情况相似，桌角这样的角点就可能产生这样的分布。而三个特征值某一个特征值相对较大，就意味着该方向上点云分布相对密集，就有可能是边界。

这里我们注意到ISS算法和Harris算法中都出现了一个协方差矩阵，那这两个协方差矩阵一样吗？当然不一样，Harris中的协方差矩阵描述的点云变化梯度的相关性，其最大特征值对应的是点云密度变化最大的方向，而ISS中的协方差矩阵描述的是点云分布的相关性，其最大特征值对应的是点云密度分布最大的方向。

实验
实验结果如下图所示

如图所示，ISS算法所提取的角点分布是相对均匀的。

4. NARF（Normal aligned radial feature）

性质
NARF 特征点提取算法的两个目标是：1). 提取出的特征点应位于表面稳定的区域以确保对法线的估计具有较强的鲁邦性，同时在其邻域点的法线应具有明显的变化。2). 提取出的点应尽可能的位于物体的边缘处，从而使得在不同的视点下都能提取到大部分相同的特征点。

算法简述
（1）边缘提取
其步骤主要如下：
1）使用启发式方法寻找不跨越边缘的相邻点的典型3D距离
首先根据点云生成Range Image，对于Range Image中的每一个点 $p_i$ 选择一个2D方形窗口内，计算方形窗口内的所有邻域点距离 $p_i$ 的距离，并对距离进行排序。$d_0$ 表示与 $p_i$ 距离最近的点，显然是 $p_i$ 自己。 $d_{s^2}$是与$p_i$最远的点，这就有可能是跨越边缘的点了。 选择一个$d_m$，作为与$p_i$同平面但距离最远的点。也就是说，如果$d_0-d_{s^2}$是一个连续递增的数列，那么$d_m$可以取平均值。如果这个数列存在某个阶跃跳动（可能会形成类似阶跃信号）那么则发生阶跃的地方应该是有边缘存在，不妨取阶跃点为$d_m$(距离较小的按个阶跃点）原文并未如此表述此段落，原文取s=5, m=9 作为m点的一个合理估计。

2）使用 1) 得到的距离信息，计算该点是边缘点的可信度(score)
对点 $p_i$ ，计算四个分数值，分别表示点的上、左、右、下方存在边缘的可能(下面仅以右边方向为例，其他方向计算方法相似)：
首先计算右边点的平均距离$$p_{right}=\frac{1}{m_p}\sum _{i=1}^{m_p}{p}_{x+i,y}$$然后计算3D距离$$d_{right}=||p_{x,y}-p_{right}||$$计算$S_{right}$$$S_{right}=max(0,1-\frac{\delta }{d_{right}})$$$S_{right}$越大，表明点$p_i$右边距离变化越明显，存在边缘的概率越大，最后再对$S_{right}$图进行平滑操作

3）如果是边缘点，将该点分为上object borders , shadow borders , veil points 3类点中的一类
根据 $p_x,y$​ 与 $p_{right}$​ 的3D距离的比较

4）非极大抑制以获得精确的边缘点位置
如果 $s_{right}$ 大于阈值(论文中取 0.8) ，并且是 $(p_{x-1},y,p_{x+1},y)$邻域内的极大值，则认为 $p_x,y$​ 为 obstacle border

（2）关键点提取
关键点提取的指导原则如下：
1）必须考虑边缘和表面结构信息; 2）必须是能在不同视角下都能被检测到的点 3）必须位于稳定的区域，从而可以得到的稳定的法向量

关键点提取的指导原则如下：
1） 对 range image 中的每一个点，计算一个表示其邻域表面变化的权重分数和这种变化的一个主方向(dominant direction)。那么对于边缘点，可以取其 权重 为1 ， 主方向为边缘方向。对于其他点，取权重为 $(1-lamda{)}^3$ ， 方向为法向在平面 p上的投影。 其中，曲率$lamda$和法向由去 $p_i$点及其周边与之距离小于 $2\delta$的点进行PCA主成分分析平面获得， $p$垂直于 $p_i$与原点连线。
2）对每个点，根据邻域点的主方向计算一个特征值，该值表达了两个性质 a）邻域内点的主方向之间的差异；b）该点所处的表面的局部变化量(表面的稳定程度) 。

3）对特征图进行平滑并进行非极大抑制，找到最终的关键点。

（3）特征描述子提取
特征描述子可以参看下面博客，因为其设计设计也比较繁琐，所以不在此赘述。

参考
3D点云——NARF特征

实验
实验结果如下图所示：

上图是根据点云图生成的RangeImage图

如图所示，关键点的分布是区域的，原因是NARF的实现边缘提取部分是根据RangeImage获得的，因此造成了关键点区域分布的情况，所以NARF应该更加使用于实际相机坐标系的使用。

5. PFH（Point feature histogram）

性质
PFH 的目的是将一个点的邻域均值曲率几何特性编码到多维的直方图中，这样的高维数据提供了非常具有信息量的特征表达。对点云不仅在六维变换（旋转及平移各三维）下具有不变性，还能很好地适应不同程度的采样或不同程度的噪声影响

算法简述
图中显示了一个中心点 $(p_q)$ 的 PFH 计算影响区域，其中该点是 $(p_q)$是三维空间中半径为$r$的球体中心点，即红色点所示。该中心点与其所有在半径为$r$的球体内部的$k$个邻域点两两互相连接组成一个网络，最终的 PFH 几何特征将通过计算邻域内这些所有点对之间的变化关系而得到的直方图。因此，计算每个点的 PFH 的计算复杂度为 $O(k^2)$。

为了计算两点之间的相对差异，给定两点 $p_1$ 和 $p_2$ ，且各自的法向分别为 $n_1$ 和 $n_2$，并在其中一个点上定义了一个局部的坐标系，三个单位向量 $u、v、w$的建立规则如下： $$u=n_1$$ $$v=u\times {\frac{p_2-p_1}{\parallel p_2-p_1\parallel }}_2$$ $$w=u\times v$$其中， ${\parallel p_2-p_1\parallel }_2$的含义是两点 $p_1$ 和 $p_2$ 之间的欧式距离。然后利用上述的 $uvw$ 坐标系， $n_1$ 和 $n_2$的差异可用三个角度 $(\alpha ,\varphi ,\theta )$ 来表达，其中的创建规则如下： $$a=v\cdot n_2$$ $$u\cdot \frac{p_2-p_1}{\parallel p_2-p_1{\parallel }_2}$$ $$\theta =\arctan (w\cdot n_2,u\cdot n_2)$$两点的空间坐标系如下图所示。

随后，对于 $k$ 邻域内的所有点对，计算这个三元组 $⟨\alpha ,\varphi ,\theta ⟩$。可见，这种表达方式仅仅用了3个参数，就能涵盖两个点原本信息中的12个参数（每个点的位置和法线各需要3个参数表达空间信息）。最后，将三元组 $⟨\alpha ,\varphi ,\theta ⟩$ 放入到直方图的各个子区间之中，以形成最终的 PFH 特征表达。

计算完中心点领域内 $n$ 个点之间的所有三元组，一共会得到 $C_n^2$ 个三元组。通常情况下，三元组中的每个特征值会被分成 $b$ 等分，所以该三元组会形成一个 $b^3$ 维的直方图，每一个维代表了其中某个值的某一个范围区间。然后，再去统计出现在各个子区间的频率即可。在实际项目中计算 PFH 时，我们的设置 $b=5$，即把每个特征值分为5等分，因此 PFH 是一个125维的向量。

参考：
几何特征系列：Point Feature Histogram（点特征直方图）

实验
实验结果如下图所示：

如图所示，红色兔子是模型，绿色兔子是场景，绿色兔子上红色的点则是红色兔子通过PFH匹配后对其的点。可以看出，对其效果是较好的，但是并没有精确对准，因此只能算作粗提取，进行位姿估计时还需要进行ICP的精估计。

6. SHOT（Signature of histogram of orientation）

性质
建立了一个独特的局部坐标系和一种新的特征描述子，具有较好的鲁棒性和描述性，在它的一些扩展方法中实现了几何和颜色纹理特征的组合

算法简述
（1）建立局部坐标系
首先获取关键点的K-近邻点，采用局部坐标系来作为几何信息的参考原点。
局部坐标系的Z轴确定方式：对K-近邻点拟合平面，平面法向量作为局部坐标系的Z轴，平面法向量即求协方差矩阵最小特征值对应的特征向量（同上面ISS算法中描述的一致）；此种解法对数据噪声有很强的鲁棒性，关键点在于要对数据去中心化处理，将坐标原点移动到数据重心。根据特征点P到重心Oi形成的向量与法向量的点乘来判断法向量正负。SHOT在计算Z轴的时候，加了距离权重。协方差矩阵如下：$$M=\frac{1}{{\sum }_{i:d_i<R}(R-d_i)}\sum _{i:d_i<R}(R-d_i)(P_i-P)(Pi-P{)}^T$$局部坐标系X轴确定方式：计算法向量n与每个近邻的点乘，找夹角最大的的P（即cos最小），然后计算P到过特征点F点平面的投影P′。特征点F与投影P′组成的方向为X轴方向。
局部坐标系Y轴确定方式：x轴和Z轴叉乘得到。

这和ISS算法中描述的局部坐标系的建立方式确实不一样，但是好处在哪里还是有点不明觉厉…

（2）进行直方图编码
借鉴2D图像的SIFT特征，邻域内法向量整体的直方图特征鲁棒性较低，因此采用局部直方图形式进行编码。SHOT将整个邻域划分为32小块，纬度2分，半径2分，经度8分，如下图所示。

计算每个近邻点相对于法向量的夹角（特征值），对每个小块内的特征值构建直方图（11级），累积直方图采用插值的形式，分别是角度差值，纬度方向插值、经度方向插值和半径方向插值。
起始就是定义一个规则，判断当前特征值落在哪个区间里面，每个曲面的特征值按照下图方式计算。最后得到的特征向量为32x11 = 352维。

（3）将描述子在纹理上进行扩展
由于最初提出的SHOT描述子是不针对纹理和颜色信息的，单纯的几何信息编码影响识别的准确率，因此后来SHOT加入了颜色信息。SHOT中的颜色编码主要是LAB空间(CIELAB计算)中的颜色；直方图31级，因此颜色描述子长度32x31 = 992;添加颜色信息的SHOT描述子共1344维。

参考：
3D特征描述子SHOT详解

实验
实验结果如下图所示：

如果所示是两只原地旋转45度位姿态不同的兔子的匹配图，通过图可能感受不出来什么，但是通过打印的数值结果来看，匹配结果还是较为准确的。

7. PPF（Point pair feature）

性质
对于缺乏表面纹理信息，或局部曲率变化很小，或点云本身就非常稀疏的物体，采用局部特征描述子很难有效的提取到匹配对。所以就有了所谓基于Point Pair 的特征，该特征使用了一些全局信息来进行匹配，最终的位姿估计结果并不会陷入局部最小值，但是计算量极大。
值得注意的是，这种方法不再涉及到关键点的概念。

算法简述
算法的思想很简单：
（1）ppf 特征为$[||d|{|}_2,<d,n_1>,<d,n_2>,<n_1,n_2>]$,其中，$d$两点连线向量，$n_1,n_2$分别为两点的法向量，该特征由一个长度和三个夹角构成；

（2）针对目标模型，在两两点之间构造点对特征F，如果有 $N$个点，那么就有 $N\times N$个特征（说明此算法是 $O(N^2)$的）， $N\times N$个特征形成特征集F_Set（哈希表结构）；
（3）在场景中任意取1定点a,再任意取1动点b，构造ppf特征，并从F_set中寻找对应的，那么理想情况下，如果找到了完全匹配的特征，则可获得点云匹配的结果；
（4）此算法是一种Hough投票算法，将位姿估计结果离散化，每次匹配都能得到一个旋转角度，如果m个b都投票给了某一旋转角度则可认为匹配成功。

PPF算法时间复杂度较高，因此如果能对点云进行预分割再进行匹配的话可以大大降低时间复杂度。

实验
下图是PPFMap的demo的实现效果，

因为我的电脑上没有GPU，这个demo需要CUDA支持，因此我暂时跑不了，由图可以知道PPF是可以从环境中匹配到目标模型的。

8. 4PCS（4-Points Congruent Sets）

性质
该方法适用于重叠区域较小或者重叠区域发生较大变化场景点云配准，无需对输入数据进行预滤波和去噪，算法能够快速准确的完成点云配准。这种方法也是不涉及关键点的概念。

算法简述
（1）在基准点云$P$中寻找满足要求的长基线（基线的确定与输入参数中overlap有很大关系，overlap越大，基线选择越长，长基线能够保证匹配的鲁棒性，且匹配数量较少）共面四点基，如图用$B=(a,b,c,d)$表示

（2）提取共面四点基的拓扑信息，按下式计算四点基间的两个比例因子$r_1、r_2$，两比例因子在点云旋转和平移变化中具有仿射不变性。$$r_1=\frac{||a-e||}{||a-b||}$$$$r_2=\frac{||c-e||}{||c-d||}$$（3）按下式计算$q1,q2\in Q$四种可能存在的交点位置，进而计算所有$Q$中长基线点对交点坐标，比较交点坐标确定匹配集合，$e_i\approx e_j$表示寻找到对应的一致全等四点，$i,j$分别表示第$i$个和第$j$个$Q$中长基线点对。$$e_1=q_1+r_1(q_2-q_1)$$$$e_2=q_1+r_2(q_2-q_1)$$（4）上图中$B=(a,b,c,d)$的全等四点对为$C=(q_1,q_3,q_4,q_5)$。寻找点云中所有$P$的共面四点集合记为$E=(B_1,B_2,...,B_m)$，$m$为$P$中四点集总数，重复上述步骤可得到全等四点集合$D=(C_1,C_2,...,C_n)$，$n$为全等四点集合总数。

在集合$D=(C_1,C_2,...,C_n)$寻找最优全等四点对，4PCS使用LCP策略寻找最优全等四点匹配，即计算全等四点旋转和平移变化参数（这里的求解过程就相当于是ICP求解问题，最少需要三对匹配点，最后是通过SVD恢复R和t），将四点转化应用到全局点云转化，记录全局配准中包含最大的一致匹配记为最优匹配，至此完成4PCS算法局部粗配准工作。
参考
4PCS点云粗配准算法介绍
4-Points Congruent Sets for Robust Pairwise Surface Registration——4PCS阅读笔记

实验
4PCS算法应该是已经集成到PCL里面了，但是能够查到的参考Demo比较，从实验室同门那里了解到目前跑4PCS相关的代码主要还是跑的源码包，相关的实验配置可以参考这个博客
安装测试点云配准4PCS算法

9. 一部分思考总结

上述讨论的8种只是3D特征子中的一部分，由于时间和篇幅的原因还有很多没有涉及到，比如PFH有加速版的FPFH，4PCS有Super-4PCS， Generalized-4PCS， V-4PCS等升级版本，不同的方法就是一种不同的设计思想。另外一个感受是，3D特征子相较于2D特征子更为复杂，但是很多3D特征子有借鉴与2D特征子的思想，例如PPF描述就有点像ORB中的BRIEF描述子的3D版。不管是2D还是3D特征子都是SLAM知识框架里面最基础的一部分，之后再慢慢地继续挖掘

此外，对SLAM算法感兴趣的同学可以看考我的博客SLAM算法总结——经典SLAM算法框架总结