角点检测（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⊗(−1 0 1)Iy=∂x∂I=I⊗(−1 0 1)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++:void cornerHarris(InputArray src, OutputArray dst, int blockSize, int ksize, double k, int borderType=BORDER_DEFAULT )

Name Disc.

Parameters src – Input single-channel 8-bit or floating-point image.
dst – Image to store the Harris detector responses. It has the type CV_32FC1 and the same size as src .
blockSize – Neighborhood size (see the details on [cornerEigenValsAndVecs()](https://docs.opencv.org/3.0-beta/modules/imgproc/doc/feature_detection.html?highlight=cornerharris#void cornerEigenValsAndVecs(InputArray src, OutputArray dst, int blockSize, int ksize, int borderType)) ).
ksize – Aperture parameter for the [Sobel()](https://docs.opencv.org/3.0-beta/modules/imgproc/doc/filtering.html#void Sobel(InputArray src, OutputArray dst, int ddepth, int dx, int dy, int ksize, double scale, double delta, int borderType)) operator.
k – Harris detector free parameter. See the formula below.
borderType – Pixel extrapolation method. See [borderInterpolate()](https://docs.opencv.org/3.0-beta/modules/core/doc/operations_on_arrays.html#int borderInterpolate(int p, int len, int borderType)) .

Name	Disc.
Parameters	src – Input single-channel 8-bit or floating-point image. dst – Image to store the Harris detector responses. It has the type `CV_32FC1` and the same size as `src` . blockSize – Neighborhood size (see the details on [`cornerEigenValsAndVecs()`](https://docs.opencv.org/3.0-beta/modules/imgproc/doc/feature_detection.html?highlight=cornerharris#void cornerEigenValsAndVecs(InputArray src, OutputArray dst, int blockSize, int ksize, int borderType)) ). ksize – Aperture parameter for the [`Sobel()`](https://docs.opencv.org/3.0-beta/modules/imgproc/doc/filtering.html#void Sobel(InputArray src, OutputArray dst, int ddepth, int dx, int dy, int ksize, double scale, double delta, int borderType)) operator. k – Harris detector free parameter. See the formula below. borderType – Pixel extrapolation method. See [`borderInterpolate()`](https://docs.opencv.org/3.0-beta/modules/core/doc/operations_on_arrays.html#int borderInterpolate(int p, int len, int borderType)) .

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++:**void goodFeaturesToTrack(InputArray image, OutputArray corners, int maxCorners, double qualityLevel, double minDistance, InputArray mask=noArray(), int blockSize=3, bool useHarrisDetector=false, double k=0.04 )

Name Disc.

Parameters: image – Input 8-bit or floating-point 32-bit, single-channel image.
eig_image – The parameter is ignored.
temp_image – The parameter is ignored.
corners – Output vector of detected corners.
maxCorners – Maximum number of corners to return. If there are more corners than are found, the strongest of them is returned.
qualityLevel – Parameter characterizing the minimal accepted quality of image corners. The parameter value is multiplied by the best corner quality measure, which is the minimal eigenvalue (see [cornerMinEigenVal()](https://docs.opencv.org/3.0-beta/modules/imgproc/doc/feature_detection.html?highlight=cornerharris#void cornerMinEigenVal(InputArray src, OutputArray dst, int blockSize, int ksize, int borderType)) ) or the Harris function response (see [cornerHarris()](https://docs.opencv.org/3.0-beta/modules/imgproc/doc/feature_detection.html?highlight=cornerharris#void cornerHarris(InputArray src, OutputArray dst, int blockSize, int ksize, double k, int borderType)) ). The corners with the quality measure less than the product are rejected. For example, if the best corner has the quality measure = 1500, and the qualityLevel=0.01 , then all the corners with the quality measure less than 15 are rejected.
minDistance – Minimum possible Euclidean distance between the returned corners.
mask – Optional region of interest. If the image is not empty (it needs to have the type CV_8UC1 and the same size as image ), it specifies the region in which the corners are detected.
blockSize – Size of an average block for computing a derivative covariation matrix over each pixel neighborhood. See [cornerEigenValsAndVecs()](https://docs.opencv.org/3.0-beta/modules/imgproc/doc/feature_detection.html?highlight=cornerharris#void cornerEigenValsAndVecs(InputArray src, OutputArray dst, int blockSize, int ksize, int borderType)) .
useHarrisDetector – Parameter indicating whether to use a Harris detector (see [cornerHarris()](https://docs.opencv.org/3.0-beta/modules/imgproc/doc/feature_detection.html?highlight=cornerharris#void cornerHarris(InputArray src, OutputArray dst, int blockSize, int ksize, double k, int borderType))) or [cornerMinEigenVal()](https://docs.opencv.org/3.0-beta/modules/imgproc/doc/feature_detection.html?highlight=cornerharris#void cornerMinEigenVal(InputArray src, OutputArray dst, int blockSize, int ksize, int borderType)).
k – Free parameter of the Harris detector.

Name	Disc.
Parameters:	image – Input 8-bit or floating-point 32-bit, single-channel image. eig_image – The parameter is ignored. temp_image – The parameter is ignored. corners – Output vector of detected corners. maxCorners – Maximum number of corners to return. If there are more corners than are found, the strongest of them is returned. qualityLevel – Parameter characterizing the minimal accepted quality of image corners. The parameter value is multiplied by the best corner quality measure, which is the minimal eigenvalue (see [`cornerMinEigenVal()`](https://docs.opencv.org/3.0-beta/modules/imgproc/doc/feature_detection.html?highlight=cornerharris#void cornerMinEigenVal(InputArray src, OutputArray dst, int blockSize, int ksize, int borderType)) ) or the Harris function response (see [`cornerHarris()`](https://docs.opencv.org/3.0-beta/modules/imgproc/doc/feature_detection.html?highlight=cornerharris#void cornerHarris(InputArray src, OutputArray dst, int blockSize, int ksize, double k, int borderType)) ). The corners with the quality measure less than the product are rejected. For example, if the best corner has the quality measure = 1500, and the `qualityLevel=0.01` , then all the corners with the quality measure less than 15 are rejected. minDistance – Minimum possible Euclidean distance between the returned corners. mask – Optional region of interest. If the image is not empty (it needs to have the type `CV_8UC1` and the same size as `image` ), it specifies the region in which the corners are detected. blockSize – Size of an average block for computing a derivative covariation matrix over each pixel neighborhood. See [`cornerEigenValsAndVecs()`](https://docs.opencv.org/3.0-beta/modules/imgproc/doc/feature_detection.html?highlight=cornerharris#void cornerEigenValsAndVecs(InputArray src, OutputArray dst, int blockSize, int ksize, int borderType)) . useHarrisDetector – Parameter indicating whether to use a Harris detector (see [`cornerHarris()`](https://docs.opencv.org/3.0-beta/modules/imgproc/doc/feature_detection.html?highlight=cornerharris#void cornerHarris(InputArray src, OutputArray dst, int blockSize, int ksize, double k, int borderType))) or [`cornerMinEigenVal()`](https://docs.opencv.org/3.0-beta/modules/imgproc/doc/feature_detection.html?highlight=cornerharris#void cornerMinEigenVal(InputArray src, OutputArray dst, int blockSize, int ksize, int borderType)). k – Free parameter of the Harris detector.

算法流程如下：

Function calculates the corner quality measure at every source image pixel using the cornerMinEigenVal() or cornerHarris() .
Function performs a non-maximum suppression (the local maximums in 3 x 3 neighborhood 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,y qualityMeasureMap(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 than maxDistance.

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角点检测原理分析