opencv3/C++ Harris角点、Shi-Tomasi角点&亚像素角点

角点检测在图像匹配、目标识别、目标跟踪、运动估计与三维重建等CV领域起着非常重要的作用。

角点定义

关于角点的定义有以下几种：

1、角点是两条及两条以上的边缘的交点；

2、角点处的一阶导数最大，二阶导数为零；

3、角点是一阶导数（即灰度梯度）的局部最大对应的像素点；

4、角点指示了物体边缘变化不连续的方向；

5、角点指图像梯度值和梯度方向的变化速率都很高的点；

Harris角点

Harris角点检测原理：

当一个窗口在图像上移动，在平滑区域上，窗口在各个方向没有变化（如图a），在边缘上，窗口在边缘方向上没有变化（如图b），在角点处，窗口在各个方向都有变化。



用I(x,y)表示图像灰度，w(x,y)表示图像窗口，用E(u,v)表示将窗口平移[u,v]造成的图像灰度的平均变化(自相关函数)，则：

E(u,v)=∑x,y
1ddc3
w(x,y)[I(x+u,y+v)−I(x,y)]2 E(u,v)=∑x,yw(x,y)[I(x+u,y+v)−I(x,y)]2

进行一次泰勒多项式展开：

I(x+u,y+v)=I(x,y)+Ixu+Iyv+O(u2,v2) I(x+u,y+v)=I(x,y)+Ixu+Iyv+O(u2,v2)

则：

E(u,v)=∑x,yw(x,y)[Ixu+Iyv+O(u2,v2)]2 E(u,v)=∑x,yw(x,y)[Ixu+Iyv+O(u2,v2)]2

E(u,v)=Ax2+By2+2Cxy E(u,v)=Ax2+By2+2Cxy

其中：

A=I2x⨂h(x,y) A=Ix2⨂h(x,y)

B=I2y⨂h(x,y) B=Iy2⨂h(x,y)

C=IxI2y⨂h(x,y) C=IxIy2⨂h(x,y)

h(x,y)是一个高斯平滑滤波函数；

则对于微小移动量[u,v]：

E(u,v)≅[u,v]M[u v] E(u,v)≅[u,v]M[u v]

其中：

M=∑x,yw(x,y)[I2xIxIyIxIyI2y] M=∑x,yw(x,y)[Ix2IxIyIxIyIy2]

导数使用Sobel算子计算。

M(u,v)=[ACCB] M(u,v)=[ACCB]

M为自相关函数E(x,y)的近似Hessian矩阵(M为2*2矩阵)。

设 λ1、λ2 λ1、λ2为M的特征值，定义角点相应函数R为：

R=λ1λ2−k(λ1+λ2)2 R=λ1λ2−k(λ1+λ2)2

为避免求解M的特征值，使用矩阵的迹与行列式进行进行替换；

矩阵的迹 tr(M)=A+B tr(M)=A+B

矩阵的行列式 det(M)=AB−C2 det(M)=AB−C2

则函数R为：

R=det(M)−k(tr(M))2 R=det(M)−k(tr(M))2

R值越大表明该点越是角点。当R大于零且较大时对应角点，若R绝对值较小时对应平坦区域，若R较小但小于零则对应边缘。

R值与点类型的关系：



Harris角点检测算法就是对角点响应函数R进行阈值处理：R > threshold。

OpenCV3角点检测

自定义角点检测

cornerEigenValsAndVecs()计算特征值和特征向量

函数cornerEigenValsAndVecs()计算角点检测的图像块的特征值和特征向量。

对于每个像素p，函数cornerEigenValsAndVecs考虑一个 blockSize×blockSize blockSize×blockSize邻域 S（p） S（p）。它计算邻域上的导数的协方差矩阵为：

M = [∑S(p)(dI/dx)2∑S(p)dI/dxdI/dy∑S(p)dI/dxdI/dy∑S(p)(dI/dy)2][∑S(p)(dI/dx)2∑S(p)dI/dxdI/dy∑S(p)dI/dxdI/dy∑S(p)(dI/dy)2]

其中导数使用Sobel算子计算。

之后，它找到 M M的特征向量和特征值，并将它们存储在目标图像中作为 (λ1，λ2，x1，y1，x2，y2) (λ1，λ2，x1，y1，x2，y2)（输出dst类型为CV_32FC（6）类型）

λ1，λ2 λ1，λ2是 M M的非排序特征值

x1，y1 x1，y1是对应于 λ1 λ1的特征向量

x2，y2 x2，y2是对应于 λ2 λ2的特征向量

函数cornerEigenValsAndVecs()的输出可用于鲁棒的边缘或角点检测。

cornerEigenValsAndVecs()参数说明：

void cornerEigenValsAndVecs(
InputArray src, //输入单通道8位或浮点图像。
OutputArray dst,//与src具有相同的大小(CV_32FC（6）类型)
int blockSize, //邻域大小
int ksize, //Sobel算子的孔径参数
int borderType = BORDER_DEFAULT //像素外插方法
);

利用cornerEigenValsAndVecs()自定义角点检测

利用cornerEigenValsAndVecs()自定义角点检测示例：

#include<opencv2/opencv.hpp>
#include<math.h>
using namespace cv;

Mat src, gray, dst, harrisRspImg;
double harrisMinRsp;
double harrisMaxRsp;
int qualityLevel = 30;
int maxCount = 100;
void cornerTrack(int, void*);

int main()
{
src = imread("E:/image/image/bdb.jpg");
if (src.empty())
{
printf("can not load image \n");
return -1;
}
namedWindow("input", WINDOW_AUTOSIZE);
imshow("input", src);
cvtColor(src, gray, CV_BGR2GRAY);

float k = 0.04;
dst = Mat::zeros(src.size(), CV_32FC(6));
harrisRspImg = Mat::zeros(src.size(), CV_32FC1);
//计算角点检测的图像块的特征值和特征向量
cornerEigenValsAndVecs(gray, dst, 3, 3, 4);
for (int r = 0; r < dst.rows; r++)
{
for (int c = 0; c < dst.cols; c++)
{
double lambda1 = dst.at<Vec6f>(r,c)[0];
double lambda2 = dst.at<Vec6f>(r,c)[1];
harrisRspImg.at<float>(r, c) = lambda1*lambda2 - k*pow((lambda1+lambda2),2);

}
}
minMaxLoc(harrisRspImg, &harrisMinRsp, &harrisMinRsp, 0, 0, Mat());

namedWindow("output", WINDOW_AUTOSIZE);
createTrackbar("QualityValue", "output", &qualityLevel, maxCount, cornerTrack);
cornerTrack(0, 0);
waitKey(0);
return 0;
}

void cornerTrack(int, void*)
{
if (qualityLevel < 10)
{
qualityLevel = 10;
}
Mat showImage = src.clone();
float t = harrisMinRsp + ((((double)qualityLevel)/maxCount)*(harrisMaxRsp - harrisMinRsp));
for (int r = 0; r < src.rows; r++)
{
for (int c = 0; c < src.cols; c++)
{
float value = harrisRspImg.at<float>(r, c);
if (value > t)
{
circle(showImage, Point(c, r), 2, Scalar(0,255,255), 2, 8, 0);
}
}
}
imshow("output", showImage);
}

亚像素角点检测

函数cornerSubPix(）通过迭代找到角点或径向鞍点精确的亚像素位置。

函数cornerSubPix(）参数说明：

void cornerSubPix(
InputArray image, //输入图像
InputOutputArray corners,//输入角点的初始坐标和为输出的精确坐标
Size winSize, //搜索窗口边长的一半
Size zeroZone,//搜索区域中间的死区大小的一半，(-1,-1)表示没有这样的大小。
TermCriteria criteria //终止角点优化迭代的条件
);

winSize为搜索窗口边长的一半，如果winSize = Size（5,5），则使用 5∗2+1×5∗2+1=11×11 5∗2+1×5∗2+1=11×11大小的搜索窗口。

实例：

#include<opencv2/opencv.hpp>
using namespace cv;

Mat src, gray, dst;
int maxCorners = 10;
void SubPixels(int, void*);
int main()
{

src = imread("E:/image/image/bdb.jpg");
if (src.empty())
{
printf("can not load image \n");
return -1;
}
namedWindow("input", WINDOW_AUTOSIZE);
imshow("input", src);
cvtColor(src, gray, CV_BGR2GRAY);
namedWindow("output", WINDOW_AUTOSIZE);
createTrackbar("corners", "output", &maxCorners, 200, SubPixels);

waitKey(0);

return 0;
}

void SubPixels(int, void*)
{
if (maxCorners < 5)
{
maxCorners = 5;
}
std::vector<Point2f> corners;
double qualityLevel = 0.01;
//获取角点
goodFeaturesToTrack(gray, corners, maxCorners, qualityLevel, 10, Mat(), 3, false, 0.04);
//清除控制台显示
system("cls");
printf("number of corners: %d \n", corners.size());
Mat result = src.clone();
for (int i = 0; i < corners.size(); i++)
{
circle(result, corners[i], 2, Scalar(0,255,0),2,8,0);
}
TermCriteria tc = TermCriteria(TermCriteria::EPS + TermCriteria::MAX_ITER, 40, 0.001);
//计算并打印出亚像素角点
cornerSubPix(gray, corners, Size(5,5), Size(-1,-1), tc);
for (int i = 0; i < corners.size(); i++)
{
printf("%d point(x,y) = (%f, %f) \n", i+1,corners[i].x,corners[i].y);
}
imshow("output", result);
}

Harris角点检测

R=det(M)−k(tr(M))2 R=det(M)−k(tr(M))2

参数k默认0.04~0.06

cornerHarris()函数参数说明

void cornerHarris(
InputArray src, //输入图像（8位单通道图像或浮点图像）
OutputArray dst, //用于存储harris探测器响应的图像（大小与src相同，类型为CV_32FC1）
int blockSize,//邻域大小
int ksize, //Sobel算子的孔径参数
double k,//计算角度响应时候的参数大小（默认0.04~0.06）
int borderType = BORDER_DEFAULT //像素外插方法
);

Harris角点检测示例：

#include<opencv2/opencv.hpp>
using namespace cv;

Mat src , gray;
int threSet = 130;
void harris(int, void*);
int main()
{

src = imread("E:/image/image/bdb.jpg");
if (src.empty())
{
printf("can not load image \n");
return -1;
}
namedWindow("input", WINDOW_AUTOSIZE);
imshow("input", src);
cvtColor(src, gray, CV_BGR2GRAY);
namedWindow("output", WINDOW_AUTOSIZE);

createTrackbar("threshold:", "output", &threSet, 255, harris);
harris(0,0);
waitKey(0);
return 0;
}

void harris(int, void*)
{
Mat dst, normdst;
dst = Mat::zeros(gray.size(), CV_32FC1);
//Haar角点检测，结果存放到dst
cornerHarris(gray, dst, 2, 3, 0.04, BORDER_DEFAULT);
//规范数组的值范围
normalize(dst, dst, 0, 255, NORM_MINMAX, CV_32FC1, Mat());
//缩放，计算绝对值，并将结果转换为8位
convertScaleAbs(dst, dst);
Mat result = src.clone();
//画出角点检测结果
for (int r = 0; r < result.rows; r++)
{
uchar* currentRow = dst.ptr(r);
for (int c = 0; c < result.cols; c++)
{
int value = (int)*currentRow;
if (value > threSet)
{
circle(result, Point(c,r), 2, Scalar(0,255,0), 2, 8, 0);
}
currentRow++;
}

}
imshow("output", result);
}

Shi-Tomasi角点检测

与Harris角点检测的不同在于在使用矩阵特征值 λ1λ2 λ1λ2计算角度响应的时候，Shi-Tomasi角点检测时候计算角点响应时使用的公式为：

R=min(λ1,λ2) R=min(λ1,λ2)

goodFeaturesToTrack()函数参数说明

函数功能：查找图像或指定图像区域中最突出的角点.。

goodFeaturesToTrack()角点检测参数说明：

void goodFeaturesToTrack(
InputArray image, //输入单通道8位或32位浮点型图像。
OutputArray corners,//输出检测到的角点
int maxCorners, //要返回的最大角点数
double qualityLevel, //图像角点的品质因子
double minDistance,//角点之间的最小距离（删除该范围内更强的角点）
InputArray mask = noArray(), //感兴趣区
int blockSize = 3,//计算协方差矩阵时的窗口大小
bool useHarrisDetector = false, //是否使用Harris角点检测（默认计算shi-tomasi角点）
double k = 0.04 //Harris角点检测需要的k值
);

Shi-Tomasi角点检测示例：

#include<opencv2/opencv.hpp>
using namespace cv;

void trackBar(int, void*);
int thre = 0;
Mat src, dst;
int main()
{
src = imread("E:/image/bdb.jpg");
if (src.empty())
{
printf("can not load image \n");
return -1;
}
namedWindow("input");
namedWindow("output");
imshow("input", src);
cvtColor(src, dst, COLOR_BGR2GRAY);

destroyWindow("output");
namedWindow("output");
createTrackbar("threshold:","output",&thre,250,trackBar);
waitKey();
return 0;
}

void trackBar(int, void*)
{
std::vector<Point2f> corners;
goodFeaturesToTrack(dst,corners, thre, 0.01, 10, Mat());
for (int i = 0; i < corners.size(); i++)
{
circle(src, corners[i], 2, Scalar(0,255,255), 2);
}
imshow("output", src);
}