Harris算子实现（C++）

---

一、Harris算子原理

Harris算子受信号处理中自相关函数的启发，给出与自相关函数相联系的矩阵M。M阵的特征值是自相关函数的一阶曲率，如果两个曲率值都高，那么就认为该点是角点特征。

① 首先确定一个n×n大小的影像窗口，对窗口内的每一个像素点进行一阶差分运算，求得在x，y方向的梯度；

② 对梯度值进行高斯滤波，高斯卷积模板的值取0.3~0.9；

③ 选取局部极值点，在窗口内取最大值。局部极值点的数目往往很多，也可以根据特征点提取数目的要求，对所有的极值点排序，根据要求选出兴趣值最大的若干个点作为最后的结果。

二、C++实现

1. 步骤

1. 导入图像
2. 计算图像在XY方向上的梯度
3. 计算方向梯度矩阵的元素
4. 利用高斯函数进行滤波
5. 计算方向梯度矩阵的迹和行列式
6. 每个像素的Harris响应值R，并对小于某一阈值t的R置为零
7. 进行非最大值抑制，局部最大值点即为图像中的角点
8. 导出图像

2. 实现代码

#include "gdal_priv.h"
#include "cpl_conv.h"
#include <iostream>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/opencv.hpp>
#include<cmath>
#include <string>
#include <vector>
using namespace std;
using namespace cv;

//main(int argc, char *argv[ ], char **env)才是UNIX和Linux中的标准写法。
int main(int argc, char** argv)
{
	void Harris(Mat image);
	//将所有需要的库函数等，导入进来并改变位数为X64
	//包括VC++目录，C/C++常规、连接器常规等，不同的库路径不同，引入的位置也不同
	string filename = "G:/myself/Major/MathPhotograph/Experiment/Picture02.jpg";//斜杠的方向不能反	//argc: 整数, 用来统计你运行程序时送给main函数的命令行参数的个数
	//* argv[ ]: 指针数组，用来存放指向字符串参数的指针，每一个元素指向一个参数
	//cout << filename << std::endl;	//测试语句
	if (argc > 1)	//如果参数比1大的话，就把第一个参数赋值给filename
	{
		filename = argv[1];	//将参数赋值给filename
		//cout <<"Hey! I'm here."<< filename << std::endl;  //测试语句
	}
	Mat image = imread(filename, IMREAD_GRAYSCALE);	//灰色样式读取图像到Mat中IMREAD_COLOR//IMREAD_GRAYSCALE//IMREAD_UNCHANGED
	if (image.empty())	//如果为空
	{
		std::cout << "Could not open or find the image." << std::endl;
		return -1;
	}
	//namedWindow("Display window", WINDOW_AUTOSIZE);	//建立大小固定的窗体显示图像,可删除
	//imshow("Display window", image);
	std::cout << "图像的类型编号为：" << image.type() << endl;	//网上有类型对应的编号
	std::cout << "图像的通道数量为：" << image.channels() << endl;

	//计算各个像元的兴趣值
	Harris(image);
	cv::waitKey(0);//等待用户需要多长时间毫秒,零意味着永远等待
	return 0;
}


//https://www.cnblogs.com/ronny/p/4009425.html
void Harris(Mat image)
{
	int rows = image.rows - 1, cols = image.cols - 1;//行列数
	int row, col, i, j;
	int Ix, Iy;

	Mat M(rows + 1, cols + 1, CV_32FC3, Scalar::all(0));	//构建M矩阵:三个通道的int类型的
	//cout << M.at<Vec3f>(1,2) << endl;
	int sideAdd, sideDes;
	for (row = 0; row <= rows; row++)
	{
		for (col = 0; col <= cols; col++)
		{
			//step1: 计算方向梯度
			sideAdd = 1, sideDes = 1;
			if (row == 0 || col == 0)	sideDes = 0;
			if (row == rows || col == cols)		sideAdd = 0;
			Ix = (image.at<uchar>(row + sideAdd, col) - image.at<uchar>(row - sideDes, col));
			Iy = (image.at<uchar>(row, col + sideAdd) - image.at<uchar>(row, col - sideDes));
			//cout <<row<<" "<<col<<" "<< Ix <<" "<< Iy << endl;	//测试语句

			//step2: 计算两个方向梯度乘积
			M.at<Vec3f>(row, col)[0] = pow(Ix, 2);
			M.at<Vec3f>(row, col)[1] = Ix*Iy;
			M.at<Vec3f>(row, col)[2] = pow(Iy, 2);
			//cout << row << " " << col << endl;
		}
	}
	//step3: 利用高斯函数进行高斯加权
	Mat M_Gauss(rows + 1, cols + 1, CV_32FC3, Scalar::all(0));
	cv::GaussianBlur(M, M_Gauss, cv::Size(5, 5), 3, 3);	//3为sigma
	//cout << M_Gauss.type() << endl;	//测试语句

	//step4: 每个像素的Harris响应值R，并对小于某一阈值t的R置为零 
	Mat R(rows + 1, cols + 1, CV_32FC1, Scalar::all(0));
	float t = 1000000, temp = 0;	//阈值大小
	float detM = 0, traceM = 0;		//行列式和迹
	for (row = 0; row <= rows; row++)
	{
		for (col = 0; col <= cols; col++)
		{
			//计算行列式和迹
			detM = M_Gauss.at<Vec3f>(row, col)[0] * M_Gauss.at<Vec3f>(row, col)[2] - pow(M_Gauss.at<Vec3f>(row, col)[1], 2);
			traceM = M_Gauss.at<Vec3f>(row, col)[0] + M_Gauss.at<Vec3f>(row, col)[2];
			//计算图像每个像元的响应值，α为经常常数，取值范围为0.04~0.06
			temp = detM - 0.05*pow(traceM, 2);		//过渡一下
			if (temp > t)	//如果大于规定的阈值
				R.at<float>(row, col) = temp;
		}
	}

	//step5：在3×3或5×5的邻域内进行非最大值抑制，局部最大值点即为图像中的角点。
	float winMax = 0;
	int maxR, maxC;		//窗口内最大值的行列值
	int windowR, windowC;	//窗口内的移动位置
	for (row = 0; row <= rows - 7; row += 7)	//窗口大小为.。
	{
		for (col = 0; col <= cols - 7; col += 7)
		{
			maxR = 0, maxC = 0;
			winMax = 0;
			for (windowR = 0; windowR < 7; windowR++)
			{
				for (windowC = 0; windowC < 7; windowC++)
				{
					if (winMax < R.at<float>(row + windowR, col + windowC))
					{
						winMax = R.at<float>(row + windowR, col + windowC);
						maxR = row + windowR;
						maxC = col + windowC;
					}
				}
			}
			//画图
			if (maxR > 0 && maxC > 0)
			{
				Point pt = Point(maxC, maxR);
				int pointColor = (0, 0, 255);
				circle(image, pt, 3, pointColor);
				std::cout << maxR << "   " << maxC << ":   " << winMax << endl;
			}
		}
	}
	imshow("Harris算子", image);
	cv::waitKey(0);
}

原理比较麻烦，实现比较容易。

三、处理结果