角点检测（Harris & Shi-Tomas）的原理及OpenCV API 的应用
这篇博客的内容主要时对参考中多篇博客的总结。
1. 角点
在现实世界中，角点对应于物体的拐角，道路的十字路口、丁字路口等。关于角点的具体描述可以有几种：
1). 一阶导数（即灰度的梯度）的局部最大所对应的像素点；
2). 两条及两条以上边缘的交点；
3). 图像中梯度值和梯度方向的变化速率都很高的点；
4).角点处的一阶导数最大，二阶导数为零，指示物体边缘变化不连续的方向。
从求解思路上可以分为两种，一种基于边缘，一种基于图像灰度。前者往往需要对图像边缘进行编码，这在很大程度上依赖于图像的分割与边缘提取，具有相当大的难度和计算量，且一旦待检测目标局部发生变化，很可能导致操作的失败。早期主要有 Rosenfeld 和 Freeman 等人的方法，后期有 CSS 等方法，参考下图。
基于图像灰度的方法通过计算点的曲率及梯度来检测角点，避免了第一类方法存在的缺陷，此类方法主要有Moravec 算子、Forstner 算子、Harris 算子、SUSAN 算子等。
这篇文章主要介绍 Harris 角点检测的算法原理。
2. Harris角点
2.1 基本原理
Harris 角点检测是一种直接基于灰度图的角点提取算法，稳定性高，尤其对L型角点（也就是直角）检测精度高。缺点也是明显的，就是运算速度慢。
人眼对角点的识别通常是在一个局部的小区域或小窗口完成的。
如果在各个方向上移动这个特定的小窗口，窗口内区域的灰度发生了较大的变化，那么就认为在窗口内遇到了角点；如果这个特定的窗口在图像各个方向上移动时，窗口内图像的灰度没有发生变化，那么窗口内就不存在角点；如果窗口在某一个方向移动时，窗口内图像的灰度发生了较大的变化，而在另一些方向上没有发生变化，那么，窗口内的图像可能就是一条直线的线段。
对于图像 ，当在点(x,y)(x,y)(x,y)处平移(Δx,Δy)(Δx,Δy)(Δx,Δy)后的自相似性，可以通过自相关函数给出：
c(x,y;Δx,Δy)=∑(u,v)∈W(x,y)w(u,v)(I(u,v)–I(u+Δx,v+Δy))2c(x,y;Δx,Δy)=\sum_{(u,v)∈W(x,y)}w(u,v)(I(u,v)–I(u+Δx,v+Δy))^2 c(x,y;Δx,Δy)=(u,v)∈W(x,y)∑w(u,v)(I(u,v)–I(u+Δx,v+Δy))2
其中，w(x,y)w(x,y)w(x,y)是以点(x,y)(x,y)(x,y)为中心的窗口，w(u,v)w(u,v)w(u,v)为加权函数，它既可是常数，也可以是高斯加权函数。
根据泰勒展开，对图像I(x,y)I(x,y)I(x,y)在平移(Δx,Δy)(Δx,Δy)(Δx,Δy)后进行一阶近似：
I(u+Δx,v+Δy)=I(u,v)+Ix(u,v)Δx+Iy(u,v)Δy+O(Δx2,Δy2)≈I(u,v)+Ix(u,v)Δx+Iy(u,v)Δy\begin{aligned} I(u+Δx,v+Δy) =&I(u,v)+I_x(u,v)Δx+I_y(u,v)Δy+O(Δx2,Δy2)\\ ≈&I(u,v)+I_x(u,v)Δx+I_y(u,v)Δy \end{aligned} I(u+Δx,v+Δy)=≈I(u,v)+Ix(u,v)Δx+Iy(u,v)Δy+O(Δx2,Δy2)I(u,v)+Ix(u,v)Δx+Iy(u,v)Δy
其中，Ix,IyI_x,I_yIx,Iy是图像I(x,y)I(x,y)I(x,y)的偏导数，这样的话，自相关函数则可以简化为：
c(x,y;Δx,Δy)≈∑w(Ix(u,v)Δx+Iy(u,v)Δy)2=[Δx,Δy]M(x,y)[ΔxΔy]\begin{aligned} c(x,y;Δx,Δy)≈&∑_w(I_x(u,v)Δx+I_y(u,v)Δy)^2\\ = &[Δx,Δy]M(x,y) \left[ \begin{matrix} Δx \\ Δy\end{matrix} \right] \end{aligned} c(x,y;Δx,Δy)≈=w∑(Ix(u,v)Δx+Iy(u,v)Δy)2[Δx,Δy]M(x,y)[ΔxΔy]
其中
M(x,y)=∑w[Ix(x,y)2Ix(x,y)Iy(x,y)Ix(x,y)Iy(x,y)Iy(x,y)2]=[∑wIx(x,y)2∑wIx(x,y)Iy(x,y)∑wIx(x,y)Iy(x,y)∑wIy(x,y)2]]=[ACCB]\begin{aligned} M(x,y)= &\sum_w \left[ \begin{matrix} 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{matrix} \right] \\ =& \left[ \begin{matrix} \sum_w I_x(x,y)^2 & \sum_wI_x(x,y)I_y(x,y) \\ \sum_w I_x(x,y)I_y(x,y) &\sum_w I_y(x,y)^2] \end{matrix} \right] \\ =& \left[ \begin{matrix} A &C \\ C &B \end{matrix} \right] \end{aligned} M(x,y)===w∑[Ix(x,y)2Ix(x,y)Iy(x,y)Ix(x,y)Iy(x,y)Iy(x,y)2][∑wIx(x,y)2∑wIx(x,y)Iy(x,y)∑wIx(x,y)Iy(x,y)∑wIy(x,y)2]][ACCB]
也就是说图像 I(x,y)I(x,y)I(x,y)在点(x,y)(x,y)(x,y)处平移(Δx,Δy)(Δx,Δy)(Δx,Δy)后的自相关函数可以近似为二项函数：
c(x,y;Δx,Δy)≈AΔx2+2CΔxΔy+BΔy2c(x,y;Δx,Δy)≈AΔx^2+2CΔxΔy+BΔy^2 c(x,y;Δx,Δy)≈AΔx2+2CΔxΔy+BΔy2
其中
A=∑wI2xB=∑wI2yC=∑wIxIyA=∑_wI^2x\\ B=∑_wI^2y\\ C=∑_wI_xI_y A=w∑I2xB=w∑I2yC=w∑IxIy
二次项函数本质上就是一个椭圆函数。椭圆的扁率和尺寸是由M(x,y)M(x,y)M(x,y)的特征值λ1、λ2λ1、λ2λ1、λ2决定的，椭圆的方向是由M(x,y)M(x,y)M(x,y)的特征矢量决定的，如下图所示，椭圆方程为：
[Δx,Δy]M(x,y)[ΔxΔy]=1[Δx,Δy]M(x,y) \left[ \begin{matrix} Δx \\ Δy \end{matrix} \right] = 1 [Δx,Δy]M(x,y)[ΔxΔy]=1
椭圆函数特征值与图像中的角点、直线（边缘）和平面之间的关系如下图所示。共可分为三种情况：
图像中的直线。一个特征值大，另一个特征值小，λ1≫λ2λ1≫λ2λ1≫λ2 或者λ2≫λ1λ2≫λ1λ2≫λ1。自相关函数值在某一方向上大，在其他方向上小。
图像中的平面。两个特征值都小，且近似相等；自相关函数数值在各个方向上都小。
图像中的角点。两个特征值都大，且近似相等，自相关函数在所有方向都增大。
根据二次项函数特征值的计算公式，我们可以求M(x,y)M(x,y)M(x,y)矩阵的特征值。但是Harris给出的角点求取方法并不需要计算具体的特征值，而是计算一个角点响应值RRR来判断角点。RRR的计算公式为：
R=detM−α(traceM)2R=detM−α(traceM)^2 R=detM−α(traceM)2
式中，detMdetMdetM为矩阵M=[A,BB,C]M = \left[ \begin{matrix} A, \space B \\ B, \space C \end{matrix} \right]M=[A,BB,C]的行列式，traceMtraceMtraceM为矩阵的直迹；α\alphaα为常数，取值范围为0.04~0.06。
**只要在某一点(x,y)(x,y)(x,y)处，RRR超过某一设定的阈值,即认为该点为角点。**显然增大ααα的值，将减小角点响应值RRR，降低角点检测的灵敏性，减少被检测角点的数量；减小ααα值，将增大角点响应值RRR，增加角点检测的灵敏性，增加被检测角点的数量。
事实上，特征值隐含在detMdetMdetM和traceMtraceMtraceM中，因为
detM=λ1λ2=AC−B2traceM=λ2+λ2=A+CdetM=λ_1λ_2=AC−B^2 \\ traceM=λ_2+λ_2=A+C detM=λ1λ2=AC−B2traceM=λ2+λ2=A+C
2.2 算法流程
Harris算子的计算步骤:
计算图像I(x,y)I(x,y)I(x,y)在XXX和YYY两个方向的梯度Ix、IyIx、IyI_x、I_yI_x、I_yIx、IyIx、Iy，其中
Ix=∂I∂x=I⊗(−101)Iy=∂I∂x=I⊗(−101)TIx=\frac{∂I}{∂x}=I⊗(−1\space 0 \space 1) \\ Iy=\frac{∂I}{∂x}=I⊗(−1 \space0 \space1)^T Ix=∂x∂I=I⊗(−101)Iy=∂x∂I=I⊗(−101)T
计算图像两个方向梯度的乘积。
I2x=Ix⋅IxI2y=Iy⋅IyIxy=Ix⋅IyI^2x=I_x⋅I_x\\ I^2y=I_y⋅I_y\\ I_{xy}=I_x⋅I_y I2x=Ix⋅IxI2y=Iy⋅IyIxy=Ix⋅Iy
使用高斯函数对Ix2、Iy2I^2_x、I^2_yIx2、Iy2和IxyI_{xy}Ixy进行高斯加权（取σ=1σ=1σ=1），生成矩阵MMM的元素A、B、CA、B、CA、B、C。
A=g(Ix2)=Ix2⊗wC=g(Iy2)=Iy2⊗wB=g(Ix,y)=Ixy⊗wA=g(I^2_x)=I^2_x⊗w\\ C=g(I^2_y)=I^2_y⊗w\\ B=g(I_{x,y})=I_{xy}⊗w A=g(Ix2)=Ix2⊗wC=g(Iy2)=Iy2⊗wB=g(Ix,y)=Ixy⊗w
计算每个像素的 Harris 响应值 RRR，并对小于某一阈值ttt的RRR置为零。
R={R:detM−α(traceM)2<t}R=\{R:detM−α(traceM)^2<t\} R={R:detM−α(traceM)2<t}
在3×33×33×3或5×55×55×5的邻域内进行非最大值抑制，局部最大值点即为图像中的角点。
2.3 Opencv API
opencv 已经提供了相关的 API，具体可以参考：
1).cornerHarris()
2). cornerHarris demo
C++:voidcornerHarris(InputArraysrc, OutputArraydst, intblockSize, intksize, doublek, intborderType=BORDER_DEFAULT )
The function runs the Harris corner detector on the image. Similarly to cornerMinEigenVal()
and cornerEigenValsAndVecs() , for each pixel (x,y)(x,y)(x,y) it calculates a 2×22 \times 22×2 gradient covariance matrix M(x,y)M^{(x,y)}M(x,y) over a blockSize×blockSizesblockSize \times blockSizesblockSize×blockSizes neighborhood. Then, it computes the following characteristic:
dst(x,y)=detM(x,y)−k⋅(trMx,y)2dst(x,y) = detM^{(x,y)} - k\cdot(trM^{x,y})^2 dst(x,y)=detM(x,y)−k⋅(trMx,y)2
Corners in the image can be found as the local maxima of this response map.
参数说明
blockSize就是窗口函数W(x, y)大小；
Ksize-sobel算子核的大小，可取1，3，5。。。（用sobel函数来得到Ix， Iy）；
K表示计算角度响应时候的参数大小，默认在0.04~0.06。
关于Harris算法的源代码可以参考Harris源代码分析。
2.4 API 调用实例
#include <opencv2/opencv.hpp> #include "opencv2/highgui/highgui.hpp" #include "opencv2/imgproc/imgproc.hpp" using namespace cv;using namespace std;Mat g_srcImage, g_srcImage1, g_grayImage;int thresh = 30; //当前阈值 int max_thresh = 175; //最大阈值 void on_CornerHarris(int, void*);//回调函数 int main(int argc, char** argv){g_srcImage = imread("lol19.jpg", 1);if (!g_srcImage.data){printf("读取图片错误！ \n");return -1;}imshow("原始图", g_srcImage);g_srcImage1 = g_srcImage.clone();//存留一张灰度图 cvtColor(g_srcImage1, g_grayImage, CV_BGR2GRAY);//创建窗口和滚动条 namedWindow("角点检测", CV_WINDOW_AUTOSIZE);createTrackbar("阈值: ", "角点检测", &thresh, max_thresh, on_CornerHarris);//调用一次回调函数，进行初始化 on_CornerHarris(0, 0);waitKey(0);return(0);}void on_CornerHarris(int, void*){Mat dstImage;//目标图 Mat normImage;//归一化后的图 Mat scaledImage;//线性变换后的八位无符号整型的图 //置零当前需要显示的两幅图，即清除上一次调用此函数时他们的值 dstImage = Mat::zeros(g_srcImage.size(), CV_32FC1);g_srcImage1 = g_srcImage.clone();//进行角点检测 //第三个参数表示邻域大小，第四个参数表示Sobel算子孔径大小，第五个参数表示Harris参数cornerHarris(g_grayImage, dstImage, 2, 3, 0.04, BORDER_DEFAULT);// 归一化与转换 normalize(dstImage, normImage, 0, 255, NORM_MINMAX, CV_32FC1, Mat());convertScaleAbs(normImage, scaledImage);//将归一化后的图线性变换成8位无符号整型 // 将检测到的，且符合阈值条件的角点绘制出来 for (int j = 0; j < normImage.rows; j++){for (int i = 0; i < normImage.cols; i++){//Mat::at<float>(j,i)获取像素值，并与阈值比较if ((int)normImage.at<float>(j, i) > thresh + 80){circle(g_srcImage1, Point(i, j), 5, Scalar(10, 10, 255), 2, 8, 0);circle(scaledImage, Point(i, j), 5, Scalar(0, 10, 255), 2, 8, 0);}}}imshow("角点检测", g_srcImage1);imshow("角点检测2", scaledImage);}
3. Shi-Tomas 角点检测
3.1 Shi-Tomas 原理
Shi-Tomas 角点检测是基于 Harris 角点检测的改进版本。上面已经提到 Harris 的角点响应值 RRR 是和 kkk 值相关的，可是 kkk 是一个经验值，如何设置一个最佳值就成了应用时的难点。
Shi-Tomas 发现，角点的稳定性其实和矩阵 MMM 的较小特征值有关，于是直接用较小的那个特征值作为分数。这样就不用调整kkk值了。
所以 Shi-Tomas 将分数公式改为如下形式：
R=min(λ1,λ2)R = min(\lambda_1, \lambda_2 ) R=min(λ1,λ2)
3.2 OpenCV API 及算法流程
**C++:**voidgoodFeaturesToTrack(InputArrayimage, OutputArraycorners, intmaxCorners, doublequalityLevel, doubleminDistance, InputArraymask=noArray(), intblockSize=3, booluseHarrisDetector=false, doublek=0.04 )
算法流程如下：
Function calculates the corner quality measure at every source image pixel using thecornerMinEigenVal()orcornerHarris().
Function performs a non-maximum suppression (the local maximums in3 x 3neighborhood are retained).
The corners with the minimal eigenvalue less than qulityLevel⋅maxx,yqualityMeasureMap(x,y)qulityLevel \cdot max_{x,y} \space qualityMeasureMap(x,y)qulityLevel⋅maxx,yqualityMeasureMap(x,y) are rejected.
The remaining corners are sorted by the quality measure in the descending order.
Function throws away each corner for which there is a stronger corner at a distance less thanmaxDistance.
3.3 API 应用
#include "opencv2/highgui/highgui.hpp"#include "opencv2/imgproc/imgproc.hpp"#include <iostream>#include <stdio.h>#include <stdlib.h>using namespace cv;using namespace std;Mat src, src_gray;int maxCorners = 23;int maxTrackbar = 100;RNG rng(12345); //RNG：random number generator，随机数产生器char* source_window = "Image";void goodFeaturesToTrack_Demo(int, void*);int main(){//转化为灰度图src = imread("E:\\image\\timg.jpg", 1);cvtColor(src, src_gray, CV_BGR2GRAY);namedWindow(source_window, CV_WINDOW_AUTOSIZE);//创建trackbarcreateTrackbar("MaxCorners:", source_window, &maxCorners, maxTrackbar, goodFeaturesToTrack_Demo);imshow(source_window, src);goodFeaturesToTrack_Demo(0, 0);waitKey(0);return(0);}void goodFeaturesToTrack_Demo(int, void*){if (maxCorners < 1) { maxCorners = 1; }//初始化 Shi-Tomasi algorithm的一些参数vector<Point2f> corners;double qualityLevel = 0.01;double minDistance = 10;int blockSize = 3;bool useHarrisDetector = false;double k = 0.04;//给原图做一次备份Mat copy;copy = src.clone();// 角点检测goodFeaturesToTrack(src_gray, corners, maxCorners, qualityLevel, minDistance, Mat(), blockSize, useHarrisDetector, k);//画出检测到的角点cout << "** Number of corners detected: " << corners.size() << endl;int r = 4;for (int i = 0; i < corners.size(); i++){circle(copy, corners[i], r, Scalar(rng.uniform(0, 255), rng.uniform(0, 255),rng.uniform(0, 255)), -1, 8, 0);}namedWindow(source_window, CV_WINDOW_AUTOSIZE);imshow(source_window, copy);}
参考
1.Harris角点
2.OpenCV角点检测：Harris角点及Shi-Tomasi角点检测
3.Opencv学习笔记（五）Harris角点检测
4.Harris角点检测原理分析
