Harris角点

1. 特征定义

   • 轮廓之间角点
   • 视角变化，仍然能被识别
   • 梯度变化大

2. 提取思想

   使用一个固定窗口在图像上进行任意方向上的滑动，比较滑动前后窗口中像素灰度变化程度，若任意方向滑动都有较大灰度变化，则认为该窗口内存在角点。这里自定义的变量是滑动的步长与滑动窗口的大小。

3. 推导计算

   窗口滑动前后的像素灰度变化：

   $E(u,v)=\sum\limits_{(x,y)€W}w(x,y)[I(x+u,y+v)-I(x,y)]^2$

   参数分析：

   • (x,y)是窗口对应的像素坐标位置；
   • (u,v)是窗口的偏移量；
   • I(x,y)是像素坐标位置(x,y)的图像灰度值；
   • w(x,y)是窗口函数，相当于不同类型的权重算子，高斯or平均；能表示窗口不同的贡献度；

   化简

   根据泰勒公式一阶展开

   $f(x+u,y+v)≈f(x,y)+uf_x(x,y)+vf_y(x,y)$

   则

   $E(u,v)\\=\sum\limits_{(x,y)€W}w(x,y)[I(x+u,y+v)-I(x,y)]^2\\ ≈\sum\limits_{(x,y)€W}w(x,y)[I(x,y)+uI_x+vI_y-I(x,y)]^2\\ =\sum\limits_{(x,y)€W}w(x,y)[u^2I_x^2+2uvI_xI_y+v^2I_y^2]\\ =\sum\limits_{(x,y)€W}w(x,y)\begin{bmatrix}u & v\end{bmatrix}\begin{bmatrix} I_x^2 & I_xIy \\ I_xI_y & I_y^2\end{bmatrix}\begin{bmatrix}u \\ v\end{bmatrix}\\ =\begin{bmatrix}u & v\end{bmatrix}(\sum\limits_{(x,y)€W}w(x,y)\begin{bmatrix} I_x^2 & I_xIy \\   I_xI_y & I_y^2\end{bmatrix})\begin{bmatrix}u \\ v\end{bmatrix}$

   化简后：

   $E(u,v)=\begin{bmatrix}u & v\end{bmatrix}M\begin{bmatrix}u \\ v\end{bmatrix}$

   $M=\sum\limits_{(x,y)€W}w(x,y)\begin{bmatrix} I_x^2 & I_xIy \\   I_xI_y & I_y^2\end{bmatrix}$,$I_x,I_y$分别是窗口内像素点(x,y)在x方向和y方向上的梯度值。

4. 但事实并非如此

   如上就是判断角点的理论方法，而实际上，是通过对窗口内的每个像素在xy方向上的梯度统计分析做判断。

   • 特征值都比较大时，即窗口中含有角点；
   • 特征值一个较大，一个较小，窗口中含有边缘；
   • 特征值都比较小，窗口处在平坦区域；

   对每一个窗口计算得到一个分数R，根据R的大小来判定窗口内是否存在harris特征角。

   $R=det(M)-k(trace(M))^2$

   $det(M)=λ_1λ_2$

   $trace(M)=λ_1+λ_2$

   R取决于M的特征值，对于角点|R|很大，平坦的区域|R|很小，边缘的|R|为负值，同时设定R的阈值进行角点判断。

5. 总流程：

   • 计算图像在xy方向的梯度

     $I_x=\frac{∂I}{∂x}=I×\begin{bmatrix} -1 & 0 & 1 \end{bmatrix}\\ I_y=\frac{∂I}{∂y}=I×\begin{bmatrix} -1 & 0 & 1 \end{bmatrix}^T$

   • 计算梯度乘积

     $I_x^2=I_xI_x\\ I_y^2=I_yI_y\\ I_xI_y=I_x*I_y$

   • 使用窗口函数对梯度乘积进行加权，计算M矩阵

     $M=\begin{bmatrix} A& C \\C  & B\end{bmatrix}=\sum\limits_{(x,y)€W}w(x,y)\begin{bmatrix} I_x^2 & I_xIy \\   I_xI_y & I_y^2\end{bmatrix}$

     $A=\sum\limits_{(x,y)€W}g(I_x^2)=\sum\limits_{(x,y)€W}I_x^2w(x,y)\\ B=\sum\limits_{(x,y)€W}g(I_y^2)=\sum\limits_{(x,y)€W}I_y^2w(x,y)\\ C=\sum\limits_{(x,y)€W}g(I_xI_y)=\sum\limits_{(x,y)€W}I_xI_y*w(x,y)$

   • 计算每个像素点的Harris响应值R

     $R=det(M)-k(trace(M))^2$

   • 过滤

     $R=\{R:det(M)-k(trace(M))^2 > t \}$

6. Harris3D思想

   Harris3D 确定角点的方法是方块体内点云的数量变化。

   在Harris2D里，我们使用了图像梯度构成的协方差矩阵。

   想象一下，如果在 点云中存在一点p，在p上建立一个局部坐标系：z方向是法线方向，x,y方向和z垂直，在p上建立一个小正方体。

   • 如果小正方体沿z方向移动，那小正方体里的点云数量应该不变
   • 如果小正方体位于边缘上，则沿边缘移动，点云数量几乎不变，沿垂直边缘方向移动，点云数量改
   • 如果小正方体位于角点上，则有两个方向都会大幅改变点云数量

   由二维的M矩阵想到用法向量x,y,z构成协方差矩阵，那么它应该是一个对称矩阵，该矩阵的特征向量有一个方向是法线方向，另外两个方向和法线垂直，那么直接用协方差矩阵替换掉图像里的M矩阵，就得到了点云的Harris算法。其中，半径r可以用来控制角点的规模，r小，则对应的角点越尖锐（对噪声更敏感），r大，则可能在平缓的区域也检测出角点。

   在一阶信息量的情况下描述了两个相邻像素的关系。显然这个思想可以轻易的移植到点云上来。

   2D中的M矩阵也有相应的3D形式，2D使用梯度，拓展到3D中则使用法向量（包含法线和方向两个信息）。

   $A=N_x^2=N_xN_x\\ B=N_y^2=N_yN_y\\ C=N_z^2=N_zN_z\\ D=N_xN_y=N_xN_y\\ E=N_xN_z=N_xN_z\\ F=N_yN_z=N_yN_z\\$

   $M= \begin{bmatrix} A& F& E \\F &B &D \\ E&D&C\end{bmatrix} =\sum\limits_{(x,y,z)€W}w(x,y,z) \begin{bmatrix} N_x^2 & N_yN_z & N_xN_z\\ N_yN_z & N_y^2 & N_xN_y\\ N_xN_z & N_xN_y & N_z^2\\ \end{bmatrix}$

   • 计算每个像素点的Harris响应值R

     $R=det(M)-k(trace(M))^2$

     $det(M)=ABC+2DEF-AD^2-BE^2-CF^2\\ trace(M)=A+B+C$

   在PCL源码中，也能看到是此设计思路，只摘取部分

   template <typename PointInT, typename PointOutT, typename NormalT> void
   pcl::HarrisKeypoint3D<PointInT, PointOutT, NormalT>::calculateNormalCovar (const std::vector<int>& neighbors, float* coefficients) const
   {
   	unsigned count = 0;
     // indices        0   1   2   3   4   5   6   7
     // coefficients: xx  xy  xz  ??  yx  yy  yz  ??
   	// ......省略1万字符
   	memset (coefficients, 0, sizeof (float) * 8);
     for (std::vector<int>::const_iterator iIt = neighbors.begin(); iIt != neighbors.end(); ++iIt)
     {
       if (pcl_isfinite (normals_->points[*iIt].normal_x))
       {
         coefficients[0] += normals_->points[*iIt].normal_x * normals_->points[*iIt].normal_x;
         coefficients[1] += normals_->points[*iIt].normal_x * normals_->points[*iIt].normal_y;
         coefficients[2] += normals_->points[*iIt].normal_x * normals_->points[*iIt].normal_z;
   
         coefficients[5] += normals_->points[*iIt].normal_y * normals_->points[*iIt].normal_y;
         coefficients[6] += normals_->points[*iIt].normal_y * normals_->points[*iIt].normal_z;
         coefficients[7] += normals_->points[*iIt].normal_z * normals_->points[*iIt].normal_z;
   
         ++count;
       }
     }
   	// ......省略1万字符
   }
   
   template <typename PointInT, typename PointOutT, typename NormalT> void
   pcl::HarrisKeypoint3D<PointInT, PointOutT, NormalT>::responseHarris (PointCloudOut &output) const
   {
   // ......省略1万字符
   	for (int pIdx = 0; pIdx < static_cast<int> (input_->size ()); ++pIdx)
     {
       const PointInT& pointIn = input_->points [pIdx];
       output [pIdx].intensity = 0.0; //std::numeric_limits<float>::quiet_NaN ();
       if (isFinite (pointIn))
       {
         std::vector<int> nn_indices;
         std::vector<float> nn_dists;
         tree_->radiusSearch (pointIn, search_radius_, nn_indices, nn_dists);
         calculateNormalCovar (nn_indices, covar);
   
         float trace = covar [0] + covar [5] + covar [7];
         if (trace != 0)
         {
           float det = covar [0] * covar [5] * covar [7] + 2.0f * covar [1] * covar [2] * covar [6]
                     - covar [2] * covar [2] * covar [5]
                     - covar [1] * covar [1] * covar [7]
                     - covar [6] * covar [6] * covar [0];
   
           output [pIdx].intensity = 0.04f + det - 0.04f * trace * trace;
         }
       }
       output [pIdx].x = pointIn.x;
       output [pIdx].y = pointIn.y;
       output [pIdx].z = pointIn.z;
     }
   // ......省略1万字符
   }
   

7. 代码样例

   • Harris2DDIY，按逻辑一步步来，用的是https://www.cnblogs.com/zyly/p/9508131.html的程序

     • Code

   • Harris2D 用opencv接口实现

     • Code

   • Harris3D 用pcl接口实现

     • Code

   • 效果图

     

     

   本节参考来源于以下网友，十分感谢

   • 第十一节、Harris角点检测原理(附源码)
   • 协方差wiki百科解释
   • 机器人三维视觉 pcl 点云库，云关键点 /Keypoints，Harris3D、 NARF