Harris角点检测，及其Matlab和OpenCV实现

1、Harris角点检测算法实现步骤
（1）计算图像$I(x,y)$在$X$和$Y$两个方向的梯度$I_x,I_y$；$$I_x=\frac{\partial I}{\partial x}=I\otimes \left[\begin{array}{ccc}-1& 0& 1\end{array}\right]$$$$I_y=\frac{\partial I}{\partial y}=I\otimes {\left[\begin{array}{ccc}-1& 0& 1\end{array}\right]}^T$$（2）计算图像两个方向梯度的乘积；
$$I_x^2=I_x{I}_x,I_y^2=I_y{I}_y,I_{xy}=I_x{I}_y$$（3）使用高斯函数对$I_x^2,I_y^2,I_{xy}$进行高斯加权（取$\sigma =1$），生成矩阵$M$的元素$A,B,C$；
$$A=g(I_x^2)=I_x^2\otimes w$$$$B=g(I_y^2)=I_y^2\otimes w$$$$C=g(I_{xy})=I_{xy}\otimes w$$（5）在$3\times 3$或$5\times 5$的邻域内进行非最大值抑制，局部最大值点记为图像中的角点。
2、Matlab实现

close all
clear
clc
%% 0、读取图像并进行灰度处理
img=imread('lena.png');
[m,n,c]=size(img);
if c>1
   img_gray=rgb2gray(img);
else
   img_gray=img;
end
%% 1、计算两个方向的梯度Ix,Iy
temp=zeros(m+2,n+2);
Ix=temp;
Iy=temp;
temp(2:m+1,2:n+1)=img_gray;
Ix(2:m+1,2:n+1)=temp(2:m+1,3:n+2)-temp(2:m+1,1:n);
Iy(2:m+1,2:n+1)=temp(3:m+2,2:n+1)-temp(1:m,2:n+1);
%% 2、计算图像两个方向的梯度乘积Ix2,Iy2,Ixy
Ix2=Ix(2:m+1,2:n+1).*Ix(2:m+1,2:n+1);
Iy2=Iy(2:m+1,2:n+1).*Iy(2:m+1,2:n+1);
Ixy=Ix(2:m+1,2:n+1).*Iy(2:m+1,2:n+1);
%% 3、使用高斯函数对Ix2,Iy2,Ixy进行高斯加权（sigma=1），生成矩阵M的元素A,B,C
h=fspecial('gaussian',[5,5],1);
A=filter2(h,Ix2);
B=filter2(h,Iy2);
C=filter2(h,Ixy);
%% 4、计算每个像素的Harris响应值R，并对小于阈值t的R置零
R=zeros(m,n);
k=0.06;
Rmax=0;
for i=1:m
   for j=1:n
      M=[A(i,j),C(i,j);C(i,j),B(i,j)];
      R(i,j)=det(M)-k*(trace(M))^2;
      if R(i,j)>Rmax
         Rmax=R(i,j); 
      end
   end
end
t=0.02*Rmax;
%% 5、非极大值抑制
s=5;    % 邻域大小
corner=zeros(m,n);
for i=(s+1)/2:m-(s+1)/2+1
    for j=(s+1)/2:n-(s+1)/2+1
       if R(i,j)==max(R(i-(s+1)/2+1:i+(s+1)/2-1,j-(s+1)/2+1:j+(s+1)/2-1),[],'all') && R(i,j)>t
          corner(i,j)=1;
       else
          corner(i,j)=0;
       end
    end
end
%% 角点显示
[x,y]=find(corner>0);   % 角点坐标
figure
imshow(img)
hold on
plot(y,x,'r+','MarkerSize',10);
hold off

结果如图所示：

3、OpenCV实现

// Harris角点检测
#include <iostream>
#include <opencv2\opencv.hpp>
using namespace std;
using namespace cv;

int main()
{
	Mat img = imread("lena.png");
	Mat img_gray;
	if (img.channels() > 1)
		cvtColor(img, img_gray, COLOR_BGR2GRAY);
	else
		img_gray = img;
	Size mn = size(img_gray);
	int n = mn.height, m = mn.width;
	float k = 0.06, Rmax = 0, t;
	img_gray.convertTo(img_gray, CV_32FC1);

	// 1、计算两个方向的梯度Ix,Iy
	Mat temp = Mat::zeros(mn + Size(2, 2), CV_32FC1);
	Mat R = Mat::zeros(mn, CV_32FC1), corner = R.clone();
	Mat Ix, Iy, Ix2, Iy2, Ixy, A, B, C, M;
	temp.copyTo(Ix);
	temp.copyTo(Iy);
	img_gray.copyTo(temp(Rect(1, 1, m, n)));
	Ix(Rect(1, 1, m, n)) = temp(Rect(2, 1, m, n)) - temp(Rect(0, 1, m, n));
	Iy(Rect(1, 1, m, n)) = temp(Rect(1, 2, m, n)) - temp(Rect(1, 0, m, n));

	// 2、计算图像两个方向的梯度乘积Ix2,Iy2,Ixy
	Ix2 = Ix(Rect(1, 1, m, n)).mul(Ix(Rect(1, 1, m, n)));
	Iy2 = Iy(Rect(1, 1, m, n)).mul(Iy(Rect(1, 1, m, n)));
	Ixy = Ix(Rect(1, 1, m, n)).mul(Iy(Rect(1, 1, m, n)));
	
	// 3、使用高斯函数对Ix2,Iy2,Ixy进行高斯加权（sigma=1），生成矩阵M的元素A,B,C
	GaussianBlur(Ix2, A, Size(5, 5), 1, 0);
	GaussianBlur(Iy2, B, Size(5, 5), 1, 0);
	GaussianBlur(Ixy, C, Size(5, 5), 1, 0);

	// 4、计算每个像素的Harris响应值R，并对小于阈值t的R置零
	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < m; j++)
		{
			M = (Mat_<float>(2, 2) << A.at<float>(i, j), C.at<float>(i, j), C.at<float>(i, j), B.at<float>(i, j));
			R.at<float>(i, j) = determinant(M) - k * (trace(M).val[0]) * (trace(M).val[0]);
			if (R.at<float>(i, j) > Rmax)
				Rmax = R.at<float>(i, j);
		}
	}
	t = 0.02*Rmax;

	// 5、 非极大值抑制
	int s = 5;	// 邻域大小
	for (int i = (s + 1) / 2 - 1; i < n - (s + 1) / 2 + 1; i++)
	{
		for (int j = (s + 1) / 2 - 1; j < m - (s + 1) / 2 + 1; j++)
		{
			double Rm;
			minMaxIdx(R(Rect(j - (s - 1) / 2, i - (s - 1) / 2, s, s)), NULL, &Rm, NULL, NULL);
			if (fabsf(R.at<float>(i, j) - (float)Rm) < 1e-10 && (R.at<float>(i, j) > t))
			{
				corner.at<float>(i, j) = 255;
				drawMarker(img, Point(j, i), Scalar(0, 0, 255), MARKER_CROSS, 10, 1);
			}
			else
				corner.at<float>(i, j) = 0;
		}
	}
	imshow("img", img);
	waitKey(0);
	return 0;
}

结果如图所示：

注意：（1）Matlab中图像坐标系$x$轴竖直向下，$y$轴水平向右，OpenCV中图像坐标系$x$轴水平向右，$y$轴水平向右竖直向下；（2）Matlab和OpenCV中的Harris角点检测算法有些许出入，“猜想”是因为：在高斯加权时，Matlab和OpenCV所用的高斯核不一样。