Harris算子，用于检测特征点
特征值都比较大时，即窗口中含有角点
特征值一个较大，一个较小，窗口中含有边缘
特征值都比较小，窗口处在平坦区域
https://blog.csdn.net/yudingjun0611/article/details/7991601
讲明白了为什么窗口函数选择高斯函数而不是阶跃函数：https://blog.csdn.net/lwzkiller/article/details/54633670
讲明白了Harris角点检测的特点：
https://blog.csdn.net/woxincd/article/details/60754658
讲明白了M矩阵的由来公式推导：
https://blog.csdn.net/qq_36387683/article/details/80550078
讲明白了为什么矩阵M要进行对角化以及标准化：https://blog.csdn.net/NewThinker_wei/article/details/45603583
讲明白了Ix,Iy与M所画椭圆的长短轴之间的关系以及进行Harris检测的具体步骤：
https://blog.csdn.net/linqianbi/article/details/78930239
算法思路
利用偏导求导公式
Ix Iy分别为该点处的梯度。
过程
1、求x、y梯度
2、高斯加权（即w（x，y）
3、得到矩阵H
4、本应该求特征值，但是此处转换为响应函数
5、输出
代码：
C++
gradient.h文件
#pragma once
#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;
void gradient(Mat image);
gradient.cpp文件
#include "gradient.h"
void gradient(Mat image)
{
//cout << "测试成功"<<endl;
if (image.empty()) //如果为空
{
std::cout << "Could not open or find the image." << std::endl;
}
//int rows = image.rows - 1, cols = image.cols - 1;
Mat image_X, image_Y, image_XY,image_X2, image_Y2, image_GX2, image_GY2,image_GXY,M,M1,M2;
//Mat tr_x = (Mat_<int>(1, 5) << -2, -1, 0, 1, 2);
//Mat tr_y = (Mat_<char>(3, 1) << 1, 0, -1);
Mat tr_x = (Mat_<char>(1, 5) << -2,-1, 0, 1,2);
Mat tr_y = (Mat_<char>(5, 1) << -2, -1, 0, 1, 2);
//for (int i = 0; i <= image.rows - 1; i++)
//{
// for (int j = 2; j <= image.cols - 3; j++)
// {
// image_X.at<char>(i, j) = -2 * image.at<char>(i, j - 2) - 1 * image.at<char>(i, j - 1) + image.at<char>(i, j + 1) + 2 * image.at<char>(i, j + 1);
// }
//
//}
//cv::filter2D(image, image_X, CV_8S, tr_x, cv::Point(-1, -1), 0);
filter2D(image , image_X , CV_64F, tr_x);//x方向梯度,
image_X2 = image_X.mul(image_X);
//imshow("X方向梯度图", image_X2);
filter2D(image, image_Y, CV_64F, tr_y);//y方向梯度
image_Y2 = image_Y.mul(image_Y);
//imshow("Y方向梯度图", image_Y2);
//GaussianBlur
image_XY = image_X.mul(image_Y);
//imshow("XY方向梯度图", image_XY);
Mat gaussKernel = getGaussianKernel(3, 1);
filter2D(image_X2, image_GX2, CV_64F, gaussKernel);
filter2D(image_Y2, image_GY2, CV_64F, gaussKernel);
filter2D(image_XY, image_GXY, CV_64F, gaussKernel);
imshow("Y方向梯度图", image_GX2);
//imshow("Y方向梯度图", image_GY2);
double alpha = 0.04;
Mat cornerStrength(image_GX2.size(), image_GX2.type());
double Rmax = 0;
for (int i = 0; i < image.rows; i++)
{
for (int j = 0; j < image.cols; j++)
{
//hconcat(image_GX2.at<char>(i,j), image_GXY.at<char>(i, j), M1);
//hconcat( image_GXY.at<char>(i, j), image_GY2.at<char>(i, j), M2);
//vconcat(M1, M2, M);
double det_m = image_GX2.at<double>(i, j) * image_GY2.at<double>(i, j)- image_GXY.at<double>(i, j)* image_GXY.at<double>(i, j);
double trace_m = image_GX2.at<double>(i, j) + image_GY2.at<double>(i, j);
cornerStrength.at<double>(i, j) = det_m - alpha * trace_m * trace_m;
if (cornerStrength.at<double>(i, j) > Rmax)
Rmax = cornerStrength.at<double>(i, j);
}
}
//imshow("cornerStrength", cornerStrength);
//double maxStrength;
//minMaxLoc(cornerStrength, NULL, &maxStrength, NULL, NULL);//求取全局最大值，为下一步求解阈值做准备
//Mat dilated;
//Mat localMax, imgDst;
//默认3*3核膨胀，膨胀之后，除了局部最大值点和原来相同，其它非局部最大值点被
//3*3邻域内的最大值点取代
//dilate(cornerStrength, dilated, Mat());
//compare(cornerStrength, dilated, localMax, CMP_EQ);//3结果图像，类型是CV_8UC1，即单通道8位图，大小和src1和src2中最大的那个一样，比较结果为真的地方值为 255，否则为0
//cv::Mat cornerMap, cornerTh;
//double qualityLevel = 0.01;
//double threshold = qualityLevel * maxStrength;
//cv::threshold(cornerStrength, cornerTh,threshold, 255, cv::THRESH_BINARY);
// 转为8-bit图
//cornerStrength.convertTo(cornerMap, CV_8U);// 和局部最大值图与，剩下角点局部最大值图，即：完成非最大值抑制
//cv::bitwise_and(cornerStrength, localMax, cornerMap);
//imgDst = cornerMap.clone();
Mat localMax= Mat::zeros(image.rows, image.cols, CV_64F);
for (int i = 1; i < cornerStrength.rows-1; i++)
{
for (int j = 1; j < cornerStrength.cols-1; j++)
{
if (
cornerStrength.at<double>(i, j)>0.01* Rmax&&
cornerStrength.at<double>(i, j) > cornerStrength.at<double>(i - 1, j - 1)
&& cornerStrength.at<double>(i, j) > cornerStrength.at<double>(i - 1, j)
&& cornerStrength.at<double>(i, j) > cornerStrength.at<double>(i - 1, j + 1)
&& cornerStrength.at<double>(i, j) > cornerStrength.at<double>(i, j - 1)
&& cornerStrength.at<double>(i, j) > cornerStrength.at<double>(i, j + 1)
&& cornerStrength.at<double>(i, j) > cornerStrength.at<double>(i + 1, j - 1)
&& cornerStrength.at<double>(i, j) > cornerStrength.at<double>(i + 1, j)
&& cornerStrength.at<double>(i, j) > cornerStrength.at<double>(i + 1, j + 1)
)
localMax.at<double>(i, j) = 1;
}
}
for (int y = 0; y < localMax.rows; y++)
{
//const uchar* rowPtr = cornerMap.ptr<uchar>(y);
for (int x = 0; x < localMax.cols; x++)
{
// 非零点就是角点
if (localMax.at<double>(y, x) ==1 )
{
Point PT=Point(x, y);
circle(image, PT, 3, Scalar(255, 0, 0));
}
}
}
imshow("Harris", image);
}main.cpp
#include"gradient.h"
void main()
{
Mat M = imread("D:\\Users\\Ou\\Desktop\\角点测试图\\格网.png", IMREAD_GRAYSCALE);
//cv::cvtColor(M, M, cv::COLOR_BGR2GRAY);
gradient(M);
waitKey();
cv::waitKey(0);//等待用户需要多长时间毫秒,零意味着永远等待
}
Matlab
% clear all; clc ;tic;
% ori_im = imread("D:\Users\Ou\Desktop\角点测试图\格网.png");
% if(size(ori_im,3)==3)
% ori_im = rgb2gray(uint8(ori_im)); %转为灰度图像
% end
%
% fx = [-2 -1 0 1 2]; % x方向梯度算子(用于Harris角点提取算法)
% Ix = filter2(fx,ori_im); % x方向滤波
% fy = [-2;-1;0;1;2]; % y方向梯度算子(用于Harris角点提取算法)
% Iy = filter2(fy,ori_im); % y方向滤波
% Ix2 = Ix.^2;
% Iy2 = Iy.^2;
% imshow(Ix2)
% figure
% imshow(Iy2)
% Ixy = Ix.*Iy;
% imshow(Ixy)
% h= fspecial('gaussian',[7 7],2);
% Ix2 = filter2(h,Ix2);
% Iy2 = filter2(h,Iy2);
% Ixy = filter2(h,Ixy);
% imshow(Ix2)
% figure
% imshow(Iy2)
% figure
% imshow(Ixy)
% height = size(ori_im,1);
% width = size(ori_im,2);
% result = zeros(height,width);
% R = zeros(height,width);
% for i = 1:height
% for j = 1:width
% M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)]; % auto correlation matrix
% R(i,j) = det(M)-0.06*(trace(M))^2;
% end
% end
%
% cnt = 0;
% for i = 2:height-1
% for j = 2:width-1
% % 进行非极大抑制，窗口大小3*3
% if R(i,j) > R(i-1,j-1) && R(i,j) > R(i-1,j) && R(i,j) > R(i-1,j+1) && R(i,j) > R(i,j-1) && R(i,j) > R(i,j+1) && R(i,j) > R(i+1,j-1) && R(i,j) > R(i+1,j) && R(i,j) > R(i+1,j+1)
% result(i,j) = 1;
% cnt = cnt+1;
% end
% end
% end
% Rsort=zeros(cnt,1);
% [posr, posc] = find(result == 1);
% for i=1:cnt
% Rsort(i)=R(posr(i),posc(i));
% end
% [Rsort,ix]=sort(Rsort,1);
% Rsort=flipud(Rsort);
% ix=flipud(ix);
% ps=100;
% posr2=zeros(ps,1);
% posc2=zeros(ps,1);
% for i=1:ps
% posr2(i)=posr(ix(i));
% posc2(i)=posc(ix(i));
% end
% imshow(ori_im);
%
% hold on;
% plot(posc2,posr2,'g+');
% toc;
%MatLab角点检测程序harris。
ori_im2=rgb2gray(imread('D:\Users\Ou\Desktop\角点测试图\格网.png'));
%ori_im2=imresize(ori_im2',0.50,'bicubic'); %加上这句图就变成竖着的了
fx = [5 0 -5;8 0 -8;5 0 -5]; % % la gaucienne，ver axe x
Ix = filter2(fx,ori_im2); % la convolution vers axe x
fy = [5 8 5;0 0 0;-5 -8 -5]; % la gaucienne，ver axe y
Iy = filter2(fy,ori_im2); % la convolution vers axe y
Ix2 = Ix.^2;
Iy2 = Iy.^2;
Ixy = Ix.*Iy;
clear Ix;
clear Iy;
h= fspecial('gaussian',[3 3],2); % générer une fonction gaussienne，sigma=2
Ix2 = filter2(h,Ix2);
imshow(Ix2);
Iy2 = filter2(h,Iy2);
Ixy = filter2(h,Ixy);
height = size(ori_im2,1);
width = size(ori_im2,2);
result = zeros(height,width); % enregistrer la position du coin
R = zeros(height,width);
K=0.04;
Rmax = 0; % chercher la valeur maximale de R
for i = 1:height
for j = 1:width
M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];
R(i,j) = det(M)-K*(trace(M))^2; % % calcule R
if R(i,j) > Rmax
Rmax = R(i,j);
end;
end;
end;
cnt = 0;
for i = 2:height-1
for j = 2:width-1
% réduire des valuers minimales ，la taille de fenetre 3*3
if R(i,j) > 0.01*Rmax && R(i,j) > R(i-1,j-1) && R(i,j) > R(i-1,j) && R(i,j) > R(i-1,j+1) && R(i,j) > R(i,j-1) && R(i,j) > R(i,j+1) && R(i,j) > R(i+1,j-1) && R(i,j) > R(i+1,j) && R(i,j) > R(i+1,j+1)
result(i,j) = 1;
cnt = cnt+1;
end;
end;
end;
imshow(result);
[posr2, posc2] = find(result == 1);
cnt % compter des coins
figure
imshow(ori_im2);
hold on;
plot(posc2,posr2,'w+');
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
