Harris, Shi-Tomasi角点检测与亚像素角点检测笔记

角点检测

角点可以是两个边缘的角点；

角点是邻域内具有两个主方向的特征点；

角点检测
基本思想
Harris
Shi-Tomasi
Harris 和 Shi-Tomasi 的缺点
亚像素角点检测
- OpenCV接口
参考

基本思想

在各个方向上移动窗口，若灰度均发生了较大的变化，那么就认为在窗口内遇到了角点。（图3）
如果窗口在图像各个方向上移动时，窗口内图像的灰度没有发生变化，那么窗口内就不存在角点。（图1）
如果窗口在某一个方向移动时，窗口内图像的灰度发生了较大的变化，而在另一些方向上没有发生变化，那么，窗口内的图像可能就是一条直线的线段。（图2）

Harris

数学模型

上面的自相关函数\(c(x,y;\Delta x,\Delta y)\)表示窗口W平移\((\Delta x,\Delta y)\)后产生的灰度变化。

若对于任何\((\Delta x,\Delta y)\)，自相关函数c的值均较大，则说明该窗口内存在角点。

上面的M矩阵为x,y方向的梯度变化。

注意：此时还不能说M为梯度变化协方差矩阵，A，B，C不是方差和协方差。因此需要进行均值化处理：将各维随机变量（这里为\(I_x,I_y\)两维）减去其均值，使处理后的随机变量的均值为0。接下来就可以进行主成分分析。

假设下面的\((I_x,I_y)\)已经进行过均值化处理，则此时M即为梯度分布协方差矩阵。

针对平坦区域，边缘区域以及角点区域三种情形的梯度分布\((I_x,I_y)\)进行分析：

显然：

平坦区域的\((I_x,I_y)\)分布在原点附近，相应的其自相关函数c在\((\Delta x,\Delta y)\)相同时较小。
边缘区域的\((I_x,I_y)\)在某一个方向上较为分散，按照PCA来理解应该只有一个主方向。
角点区域的\((I_x,I_y)\)在x,y两个方向上分布均较为分散。

对协方差矩阵M进行对角化，则特征值即为主方向上的方差，对应的特征向量即为主方向（详见PCA原理）。

得出结论：

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

在实际判断角点时，无需计算矩阵M的特征值，使用下式近似计算角点响应值：

\[R=detM - \alpha (traceM)^2 \\ detM = \lambda_1\lambda_2 = AB-C^2 \\ traceM = \lambda_1 + \lambda_2 = A+B \]

\(\alpha\)为一常数，通常取0.04～0.06.

使用matlab绘制响应值函数图像如下：

则可以发现，该函数图像刚好满足协方差矩阵M特征值与角点的关系：

MATLAB代码如下：

x = linspace(0,10);
y = linspace(0,10);
[X,Y] = meshgrid(x,y);
Z = X.*Y - 0.05*(X+Y).^2;
contour(X,Y,Z,'ShowText','on');
xlabel("\lambda_1");ylabel("\lambda_2");

算法实现

性质

参数α对角点检测的影响：

增大α的值，将减小角点响应值R，减少被检测角点的数量；减小α的值，将增大角点响应值R，增加被检测角点的数量。
Harris角点检测对亮度和对比度的变化不敏感。（灰度不变性）
Harris 检测器具有旋转不变性，但不具有尺度不变性，也就是说尺度变化可能会导致角点变为边缘，如下图所示：

为保证角点的尺度不变性，提出了多尺度Harris角点检测。

OpenCV接口

void cornerHarris( 
    InputArray src, 
    OutputArray dst, 
    int blockSize,
    int ksize, 
    double k,
    int borderType=BORDER_DEFAULT );

src – 输入的单通道8-bit或浮点图像。
dst – 存储着Harris角点响应的图像矩阵，大小与输入图像大小相同，是一个浮点型矩阵。
blockSize – 邻域大小。
apertureSize – 扩展的微分算子大。
k – 响应公式中的，参数αα。
boderType – 边界处理的类型.

该接口很少使用，一般我们使用goodFeaturesToTrack函数：

void cv::goodFeaturesToTrack(
    cv::InputArray image, // 输入图像（CV_8UC1 CV_32FC1）
    cv::OutputArray corners, // 输出角点vector
    int maxCorners, // 最大角点数目
    double qualityLevel, // 质量水平系数（小于1.0的正数，一般在0.01-0.1之间）
    double minDistance, // 最小距离，小于此距离的点忽略
    cv::InputArray mask = noArray(), // mask=0的点忽略
    int blockSize = 3, // 使用的邻域数
    bool useHarrisDetector = false, // false ='Shi Tomasi metric'
    double k = 0.04 // Harris角点检测时使用
);

第一个参数是输入图像（8位或32位单通道图）。

第二个参数是检测到的所有角点，类型为vector或数组，由实际给定的参数类型而定。如果是vector，那么它应该是一个包含cv::Point2f的vector对象；如果类型是cv::Mat,那么它的每一行对应一个角点，点的x、y位置分别是两列。

第三个参数用于限定检测到的点数的最大值。

第四个参数表示检测到的角点的质量水平（通常是0.10到0.01之间的数值，不能大于1.0）。

第五个参数用于区分相邻两个角点的最小距离（小于这个距离得点将进行合并）。

第六个参数是mask，如果指定，它的维度必须和输入图像一致，且在mask值为0处不进行角点检测。

第七个参数是blockSize，表示在计算角点时参与运算的区域大小，常用值为3，但是如果图像的分辨率较高则可以考虑使用较大一点的值。

第八个参数用于指定角点检测的方法，如果是true则使用Harris角点检测，false则使用Shi Tomasi算法。

第九个参数是在使用Harris算法时使用，最好使用默认值0.04。

Shi-Tomasi

Harris的改进，考虑到Harris算法的稳定性与经验值k有关，k不好设定最佳值。

Shi-Tomasi 发现，角点的稳定性其实和矩阵 M 的较小特征值有关，于是直接用较小的那个特征值作为分数。这样就不用调整k值了。

所以 Shi-Tomasi 将分数公式改为如下形式：

\[R = min(\lambda_1,\lambda_2) \]

Harris 和 Shi-Tomasi 的缺点

Harris 和 Shi-Tomasi 都是基于梯度计算的角点检测方法，Shi-Tomasi 的效果要好一些。

基于梯度的检测方法有一些缺点: 计算复杂度高，图像中的噪声可以阻碍梯度计算。

想要提高检测速度的话，可以考虑基于模板的方法：FAST角点检测算法。

亚像素角点检测

对检测到的角点作进一步的优化计算，可使角点的精度达到亚像素（subPixel）级别。

上图中q为真实角点，p为角点附近任意一点。角点附近必然存在灰度不同区域，则p点或位于灰度均匀区域，或位于灰度不同区域的分界线上。

若p在灰度均匀区域内部，则如左图所示：p点梯度为0；若p在两个灰度不同区域的边缘，如右图所示，则p点的梯度方向垂直边缘方向，而边缘方向与向量q-p方向一致。因此，q-p向量与p点的梯度向量点积运算结果始终为0。在初始角点附近邻域内，可以找到多个不同的点p_i，计算其梯度方向\(\nabla I(p_i)\)，则可以得到方程：

\[\nabla I(p_i) \cdot (p_i - q)= 0 \]

多个方程联立，使用最小二乘法即可求解出点q。

然后以求解出的点q为邻域中心，继续迭代，即可获得满足精度的亚像素角点q。

关于求解的详细公式推导和思想，可以参考亚像素角点的求法

OpenCV接口

void cv::cornerSubPix(
    cv::InputArray image, // 输入图像
    cv::InputOutputArray corners, // 角点（既作为输入也作为输出）
    cv::Size winSize, // 区域大小为 NXN; N=(winSize*2+1)
    cv::Size zeroZone, // 类似于winSize，但是总具有较小的范围，Size(-1,-1)表示忽略
    cv::TermCriteria criteria // 停止优化的标准
);

第一个参数是输入图像，和cv::goodFeaturesToTrack()中的输入图像是同一个图像。

第二个参数是检测到的角点，即是输入也是输出。

第三个参数是Size类型，表示搜索窗口的半径。若winSize=Size(5，5)，那么就表示使用（52+1）x（52+1）=11*11大小的搜索窗口。

第四个参数为Size类型，表示死区的一半尺寸。而死区为不对搜索区的中央位置做求和运算的区域，用来避免自相关矩阵出现的某些可能的奇异性。值为（-1，-1）表示没有死区。

第五个参数用于表示计算亚像素时停止迭代的标准，可选的值有cv::TermCriteria::MAX_ITER 、cv::TermCriteria::EPS（可以是两者其一，或两者均选），前者表示迭代次数达到了最大次数时停止，后者表示角点位置变化的最小值已经达到最小时停止迭代。二者均使用cv::TermCriteria()构造函数进行指定