嗯.....本来打算一两天写一个简单的算法,没想到Harris进行的这么顺利,很快就完成了,那今天就写两个吧。承接上一篇的Forstner算子,继续完成Harris算子,参考资料为:
张祖勋《数字摄影测量学》
c++代码如下:
#include <iostream>
#include<opencv2/opencv.hpp>
using namespace std;
using namespace cv;
vector<Point2f>featurePoints;
int k = 2;//5*5的窗口
int ksize = 5;//5*5的极值检测窗口
void Harris(Mat src)
{
Mat I = Mat::zeros(src.rows, src.cols, CV_64FC1);//存储各像素的兴趣值
Mat Ixx = Mat::zeros(src.rows, src.cols, CV_64FC1);//M左上角(高斯滤波前)
Mat Iyy = Mat::zeros(src.rows, src.cols, CV_64FC1);//M右下角
Mat Ixy = Mat::zeros(src.rows, src.cols, CV_64FC1);//M副对角线
double GaussKernel[5][5] = //定义高斯卷积核
{
0.00390625, 0.015625, 0.0234375, 0.015625, 0.00390625,
0.015625, 0.0625, 0.09375, 0.0625, 0.015625,
0.0234375, 0.09375, 0.140625, 0.09375, 0.0234375,
0.015625, 0.0625, 0.09375, 0.0625, 0.015625,
0.00390625, 0.015625, 0.0234375, 0.015625, 0.00390625,
};
//1.一阶差分
for (int r = k; r < src.rows - k; r++)
{
for (int c = k; c < src.cols - k; c++)
{
double gx = 0, gy = 0;
for (int i = -k; i <= k - 1; i++)
{
gx += src.at<uchar>(r, c + i) - src.at<uchar>(r, c + i + 1);//x方向梯度
gy += src.at<uchar>(r + i, c) - src.at<uchar>(r + i + 1, c);//y方向梯度
}
Ixx.at<double>(r, c) = gx * gx;
Ixy.at<double>(r, c) = gx * gy;
Iyy.at<double>(r, c) = gy * gy;
}
}
for (int r = k; r < src.rows - 5; r++)
{
for (int c = k; c < src.cols - 5; c++)
{
double M11 = 0, M12 = 0, M22 = 0;
for (int j = 0; j < 5; j++)
{
for (int i = 0; i < 5; i++)
{
//2.对梯度值进行高斯滤波
M11 += Ixx.at<double>(r + j, c + i) * GaussKernel[j][i];//高斯滤波后的M值
M12 += Ixy.at<double>(r + j, c + i) * GaussKernel[j][i];
M22 += Iyy.at<double>(r + j, c + i) * GaussKernel[j][i];
}
}
double detM = M11 * M22 - M12 * M12;
double trM = M11 + M22;
double k = 0.04;
//3.计算兴趣值
I.at<double>(r, c) = detM - k * trM * trM;//计算兴趣值
}
}
//cout << "I:" << I << endl;
//4.选取极值点
for (int r = ksize; r < src.rows - ksize; r += ksize)
{
for (int c = ksize; c < src.cols - ksize; c += ksize)
{
double max = 0;
int Flag = 0;
Point2f point;
for (int j = -ksize; j <= ksize; j++)
{
for (int i = -ksize; i <= ksize; i++)
{
double value = I.at<double>(r + j, c + i);
//局部非最大值抑制
if (value > max)
{
max = value;
point.x = c + i;
point.y = r + j;
Flag = 1;
}
}
}
if (Flag == 1)
{
featurePoints.push_back(point);
}
}
}
cout << "特征点个数:" << featurePoints.size() << endl;
}
void drawPoint(Mat img)
{
for (int i = 0; i < featurePoints.size(); i++)
{
int centerx = featurePoints[i].x;
int centery = featurePoints[i].y;
circle(img, Point(centerx, centery), 4, Scalar(255, 0, 0), 1);
}
}
int main()
{
Mat src = imread("C:\\Users\\Administrator\\Desktop\\images\\标定板.png", IMREAD_GRAYSCALE);
Harris(src);
Mat img = imread("C:\\Users\\Administrator\\Desktop\\images\\标定板.png", IMREAD_COLOR);
drawPoint(img);
imwrite("D:\\Software\\VS2019\\source\\repos\\摄影测量与三维重建\\3.Harris点特征提取算子\\result\\Harris点特征提取算子.jpg", img);
waitKey(0);
return 0;
}
这个算法比较简单,也可能有前边两个做铺垫,这个明显写起来很顺手,也没遇到太大的阻碍,原图和结果如下: