1. 首先，我们不禁要问什么是harris角点？
对于角点，到目前为止还没有明确的数学定义。但是你可以认为角点就是极值点，即在某方面属性特别突出的点。一般的角点检测都是对有具体定义的、或者是能够具体检测出来的兴趣点的检测。这意味着兴趣点可以是角点，是在某些属性上强度最大或者最小的孤立点、线段的终点，或者是曲线上局部曲率最大的点。
通俗的来说，在一副图像中，我们可以认为角点是物体轮廓线的连接点（见图1），当拍摄视角变化的时候，这些特征点仍能很好地保持稳定的属性。
图1 corner
角点在保留图像图形重要特征的同时，可以有效地减少信息的数据量，使其信息的含量很高，有效地提高了计算的速度，有利于图像的可靠匹配，使得实时处理成为可能。它的各种应用，这里我就不再赘述了。
2. 如何检测出harris角点？
图2 角点检测的基本思想
角点检测最原始的想法就是取某个像素的一个邻域窗口，当这个窗口在各个方向上进行小范围移动时，观察窗口内平均的像素灰度值的变化（即，E(u,v)，Window-averaged change of intensity）。从上图可知，我们可以将一幅图像大致分为三个区域（‘flat’，‘edge’，‘corner’），这三个区域变化是不一样的。
其中，u,v是窗口在水平，竖直方向的偏移，
这里可以先简单复习一下泰勒级数展开的知识，因为马上就用到啦，
这是一维的情况，对于多元函数，也有类似的泰勒公式。
对I(x+u,y+v)进行二维泰勒级数展开，我们取一阶近似，有
图中蓝线圈出的部分我们称之为结构张量（structure tensor），用M表示。
讲到这里，先明确一点，我们的目的是什么？我们的目的是寻找这样的像素点，它使得我们的u，v无论怎样取值，E(u,v)都是变化比较大的，这个像素点就是我们要找的角点。不难观察，上式近似处理后的E(u,v)是一个二次型，而由下述定理可知，
令E(u,v)=常数，我们可用一个椭圆来描绘这一函数。
椭圆的长短轴是与结构张量M的两个特征值
相对应的量。通过判断的情况我们就可以区分出‘flat’，‘edge’，‘corner’这三种区域，因为最直观的印象：
corner：在水平、竖直两个方向上变化均较大的点，即Ix、Iy都较大；
edge ：仅在水平、或者仅在竖直方向有较大的点，即Ix和Iy只有其一较大 ；
flat ： 在水平、竖直方向的变化量均较小的点，即Ix、Iy都较小；
而结构张量M是由Ix，Iy构成，它的特征值正好可以反映Ix，Iy的情况，下面我以一种更容易理解的方式来讲述椭圆的物理意义。
这样是不是更清楚了呢^_^……，因此我们可以得出结论：
当然，大牛们并没有止步于此，这样通过判断两个变量的值来判断角点毕竟不是很方便。于是他们想出了一种更好的方法，对，就是定义角点响应函数R（corner response function），
针对三种区域，R的取值如何呢？
至此，我们就可以通过判断R的值来判断某个点是不是角点了。
角点：R为大数值整数
边缘：R为大数值负数
平坦区：绝对值R是小数值
3. harris角点检测算法步骤
1.利用Soble计算出XY方向的梯度值
2.计算出Ix^2,Iy^2,Ix*Iy
3.利用高斯函数对Ix^2,Iy^2,Ix*Iy进行滤波
4.计算局部特征结果矩阵M的特征值和响应函数C(i,j)=Det(M)-k(trace(M))^2 (0.04<=k<=0.06)
5.将计算出响应函数的值C进行非极大值抑制，滤除一些不是角点的点，同时要满足大于设定的阈值
下面放上源代码：
#include "opencv2/imgproc/imgproc.hpp" #include "opencv2/highgui/highgui.hpp" #include
#include
using namespace cv; using namespace std; /* RGB转换成灰度图像的一个常用公式是： Gray = R*0.299 + G*0.587 + B*0.114 */ //灰度转换函数* //第一个参数image输入的彩色RGB图像的引用； //第二个参数imageGray是转换后输出的灰度图像的引用； //* void ConvertRGB2GRAY(const Mat &image, Mat &imageGray); //Sobel卷积因子计算X、Y方向梯度和梯度方向角 //第一个参数imageSourc原始灰度图像； //第二个参数imageSobelX是X方向梯度图像； //第三个参数imageSobelY是Y方向梯度图像； //第四个参数pointDrection是梯度方向角数组指针 //* void SobelGradDirction(Mat &imageSource, Mat &imageSobelX, Mat &imageSobelY); //计算Sobel的X方向梯度幅值的平方* //第一个参数imageGradX是X方向梯度图像； //第二个参数SobelAmpXX是输出的X方向梯度图像的平方 //* void SobelXX(const Mat imageGradX, Mat_
&SobelAmpXX); //计算Sobel的Y方向梯度幅值的平方* //第一个参数imageGradY是Y方向梯度图像； //第二个参数SobelAmpXX是输出的Y方向梯度图像的平方 //* void SobelYY(const Mat imageGradY, Mat_
&SobelAmpYY); //计算Sobel的XY方向梯度幅值的乘积* //第一个参数imageGradX是X方向梯度图像； //第二个参数imageGradY是Y方向梯度图像； //第二个参数SobelAmpXY是输出的XY方向梯度图像 //* void SobelXY(const Mat imageGradX, const Mat imageGradY, Mat_
&SobelAmpXY); //计算一维高斯的权值数组* //第一个参数size是代表的卷积核的边长的大小 //第二个参数sigma表示的是sigma的大小 //* double *getOneGuassionArray(int size, double sigma); //高斯滤波函数的实现* //第一个参数srcImage是代表的输入的原图 //第二个参数dst表示的是输出的图 //第三个参数size表示的是卷积核的边长的大小 //* void MyGaussianBlur(Mat_
&srcImage, Mat_
&dst, int size); //计算局部特涨结果矩阵M的特征值和响应函数H = (A*B - C) - k*(A+B)^2 //M //A C //C B //Tr(M)=a+b=A+B //Det(M)=a*b=A*B-C^2 //计算输出响应函数的值得矩阵 // void harrisResponse(Mat_
&GaussXX, Mat_
&GaussYY, Mat_
&GaussXY, Mat_
&resultData,float k); //*非极大值抑制和满足阈值及某邻域内的局部极大值为角点 //第一个参数是响应函数的矩阵 //第二个参数是输入的灰度图像 //第三个参数表示的是输出的角点检测到的结果图 void LocalMaxValue(Mat_
&resultData, Mat &srcGray, Mat &ResultImage,int kSize); int main() { const Mat srcImage = imread("3.jpg"); if (!srcImage.data) { printf("could not load image...\n"); return -1; } imshow("srcImage", srcImage); Mat srcGray; ConvertRGB2GRAY(srcImage, srcGray); Mat imageSobelX; Mat imageSobelY; Mat resultImage; Mat_
imageSobelXX; Mat_
imageSobelYY; Mat_
imageSobelXY; Mat_
GaussianXX; Mat_
GaussianYY; Mat_
GaussianXY; Mat_
HarrisRespond; //计算Soble的XY梯度 SobelGradDirction(srcGray, imageSobelX, imageSobelY); //计算X方向的梯度的平方 SobelXX(imageSobelX, imageSobelXX); SobelYY(imageSobelY, imageSobelYY); SobelXY(imageSobelX, imageSobelY, imageSobelXY); //计算高斯模糊XX YY XY MyGaussianBlur(imageSobelXX, GaussianXX,3); MyGaussianBlur(imageSobelYY, GaussianYY, 3); MyGaussianBlur(imageSobelXY, GaussianXY, 3); harrisResponse(GaussianXX, GaussianYY, GaussianXY, HarrisRespond, 0.05); LocalMaxValue(HarrisRespond, srcGray, resultImage, 3); imshow("imageSobelX", imageSobelX); imshow("imageSobelY", imageSobelY); imshow("resultImage", resultImage); waitKey(0); return 0; } void ConvertRGB2GRAY(const Mat &image, Mat &imageGray) { if (!image.data || image.channels() != 3) { return; } //创建一张单通道的灰度图像 imageGray = Mat::zeros(image.size(), CV_8UC1); //取出存储图像像素的数组的指针 uchar *pointImage = image.data; uchar *pointImageGray = imageGray.data; //取出图像每行所占的字节数 size_t stepImage = image.step; size_t stepImageGray = imageGray.step; for (int i = 0; i < imageGray.rows; i++) { for (int j = 0; j < imageGray.cols; j++) { pointImageGray[i*stepImageGray + j] = (uchar)(0.114*pointImage[i*stepImage + 3 * j] + 0.587*pointImage[i*stepImage + 3 * j + 1] + 0.299*pointImage[i*stepImage + 3 * j + 2]); } } } //存储梯度膜长 void SobelGradDirction(Mat &imageSource, Mat &imageSobelX, Mat &imageSobelY) { imageSobelX = Mat::zeros(imageSource.size(), CV_32SC1); imageSobelY = Mat::zeros(imageSource.size(), CV_32SC1); //取出原图和X和Y梯度图的数组的首地址 uchar *P = imageSource.data; uchar *PX = imageSobelX.data; uchar *PY = imageSobelY.data; //取出每行所占据的字节数 int step = imageSource.step; int stepXY = imageSobelX.step; int index = 0;//梯度方向角的索引 for (int i = 1; i < imageSource.rows - 1; ++i) { for (int j = 1; j < imageSource.cols - 1; ++j) { //通过指针遍历图像上每一个像素 double gradY = P[(i + 1)*step + j - 1] + P[(i + 1)*step + j] * 2 + P[(i + 1)*step + j + 1] - P[(i - 1)*step + j - 1] - P[(i - 1)*step + j] * 2 - P[(i - 1)*step + j + 1]; PY[i*stepXY + j*(stepXY / step)] = abs(gradY); double gradX = P[(i - 1)*step + j + 1] + P[i*step + j + 1] * 2 + P[(i + 1)*step + j + 1] - P[(i - 1)*step + j - 1] - P[i*step + j - 1] * 2 - P[(i + 1)*step + j - 1]; PX[i*stepXY + j*(stepXY / step)] = abs(gradX); } } //将梯度数组转换成8位无符号整型 convertScaleAbs(imageSobelX, imageSobelX); convertScaleAbs(imageSobelY, imageSobelY); } void SobelXX(const Mat imageGradX, Mat_
&SobelAmpXX) { SobelAmpXX = Mat_
(imageGradX.size(), CV_32FC1); for (int i = 0; i < SobelAmpXX.rows; i++) { for (int j = 0; j < SobelAmpXX.cols; j++) { SobelAmpXX.at
(i, j) = imageGradX.at
(i, j)*imageGradX.at
(i, j); } } //convertScaleAbs(SobelAmpXX, SobelAmpXX); } void SobelYY(const Mat imageGradY, Mat_
&SobelAmpYY) { SobelAmpYY = Mat_
(imageGradY.size(), CV_32FC1); for (int i = 0; i < SobelAmpYY.rows; i++) { for (int j = 0; j < SobelAmpYY.cols; j++) { SobelAmpYY.at
(i, j) = imageGradY.at
(i, j)*imageGradY.at
(i, j); } } //convertScaleAbs(SobelAmpYY, SobelAmpYY); } void SobelXY(const Mat imageGradX, const Mat imageGradY, Mat_
&SobelAmpXY) { SobelAmpXY = Mat_
(imageGradX.size(), CV_32FC1); for (int i = 0; i < SobelAmpXY.rows; i++) { for (int j = 0; j < SobelAmpXY.cols; j++) { SobelAmpXY.at
(i, j) = imageGradX.at
(i, j)*imageGradY.at
(i, j); } } //convertScaleAbs(SobelAmpXY, SobelAmpXY); } //计算一维高斯的权值数组 double *getOneGuassionArray(int size, double sigma) { double sum = 0.0; //定义高斯核半径 int kerR = size / 2; //建立一个size大小的动态一维数组 double *arr = new double[size]; for (int i = 0; i < size; i++) { // 高斯函数前的常数可以不用计算，会在归一化的过程中给消去 arr[i] = exp(-((i - kerR)*(i - kerR)) / (2 * sigma*sigma)); sum += arr[i];//将所有的值进行相加 } //进行归一化 for (int i = 0; i < size; i++) { arr[i] /= sum; cout << arr[i] << endl; } return arr; } void MyGaussianBlur(Mat_
&srcImage, Mat_
&dst, int size) { CV_Assert(srcImage.channels() == 1 || srcImage.channels() == 3); // 只处理单通道或者三通道图像 int kerR = size / 2; dst = srcImage.clone(); int channels = dst.channels(); double* arr; arr = getOneGuassionArray(size, 1);//先求出高斯数组 //遍历图像 水平方向的卷积 for (int i = kerR; i < dst.rows - kerR; i++) { for (int j = kerR; j < dst.cols - kerR; j++) { float GuassionSum[3] = { 0 }; //滑窗搜索完成高斯核平滑 for (int k = -kerR; k <= kerR; k++) { if (channels == 1)//如果只是单通道 { GuassionSum[0] += arr[kerR + k] * dst.at
(i, j + k);//行不变，列变换，先做水平方向的卷积 } else if (channels == 3)//如果是三通道的情况 { Vec3f bgr = dst.at
(i, j + k); auto a = arr[kerR + k]; GuassionSum[0] += a*bgr[0]; GuassionSum[1] += a*bgr[1]; GuassionSum[2] += a*bgr[2]; } } for (int k = 0; k < channels; k++) { if (GuassionSum[k] < 0) GuassionSum[k] = 0; else if (GuassionSum[k] > 255) GuassionSum[k] = 255; } if (channels == 1) dst.at
(i, j) = static_cast
(GuassionSum[0]); else if (channels == 3) { Vec3f bgr = { static_cast
(GuassionSum[0]), static_cast
(GuassionSum[1]), static_cast
(GuassionSum[2]) }; dst.at
(i, j) = bgr; } } } //竖直方向 for (int i = kerR; i < dst.rows - kerR; i++) { for (int j = kerR; j < dst.cols - kerR; j++) { float GuassionSum[3] = { 0 }; //滑窗搜索完成高斯核平滑 for (int k = -kerR; k <= kerR; k++) { if (channels == 1)//如果只是单通道 { GuassionSum[0] += arr[kerR + k] * dst.at
(i + k, j);//行变，列不换，再做竖直方向的卷积 } else if (channels == 3)//如果是三通道的情况 { Vec3f bgr = dst.at
(i + k, j); auto a = arr[kerR + k]; GuassionSum[0] += a*bgr[0]; GuassionSum[1] += a*bgr[1]; GuassionSum[2] += a*bgr[2]; } } for (int k = 0; k < channels; k++) { if (GuassionSum[k] < 0) GuassionSum[k] = 0; else if (GuassionSum[k] > 255) GuassionSum[k] = 255; } if (channels == 1) dst.at
(i, j) = static_cast
(GuassionSum[0]); else if (channels == 3) { Vec3f bgr = { static_cast
(GuassionSum[0]), static_cast
(GuassionSum[1]), static_cast
(GuassionSum[2]) }; dst.at
(i, j) = bgr; } } } delete[] arr; } void harrisResponse(Mat_
&GaussXX, Mat_
&GaussYY, Mat_
&GaussXY, Mat_
&resultData,float k) { //创建一张响应函数输出的矩阵 resultData = Mat_
(GaussXX.size(), CV_32FC1); for (int i = 0; i < resultData.rows; i++) { for (int j = 0; j < resultData.cols; j++) { float a = GaussXX.at
(i, j); float b = GaussYY.at
(i, j); float c = GaussXY.at
(i, j); resultData.at
(i, j) = a*b - c*c - k*(a + b)*(a + b); } } } //非极大值抑制 void LocalMaxValue(Mat_
&resultData, Mat &srcGray, Mat &ResultImage, int kSize) { int r = kSize / 2; ResultImage = srcGray.clone(); for (int i = r; i < ResultImage.rows - r; i++) { for (int j = r; j < ResultImage.cols - r; j++) { if (resultData.at
(i, j) > resultData.at
(i - 1, j - 1) && resultData.at
(i, j) > resultData.at
(i - 1, j) && resultData.at
(i, j) > resultData.at
(i - 1, j - 1) && resultData.at
(i, j) > resultData.at
(i - 1, j + 1) && resultData.at
(i, j) > resultData.at
(i, j - 1) && resultData.at
(i, j) > resultData.at
(i, j + 1) && resultData.at
(i, j) > resultData.at
(i + 1, j - 1) && resultData.at
(i, j) > resultData.at
(i + 1, j) && resultData.at
(i, j) > resultData.at
(i + 1, j + 1)) { if ((int)resultData.at
(i, j) > 18000) { circle(ResultImage, Point(i, j), 5, Scalar(0,0,255), 2, 8, 0); } } } } }
发布者：全栈程序员-站长，转载请注明出处：https://javaforall.net/175987.html原文链接：https://javaforall.net
