Harris角点检测学习总结
人们普遍认为角点是二维图像亮度变化剧烈的点或图像边缘曲线上曲率极大值的点。这些点在保留图像图形重要特征的同时,可以有效地减少信息的数据量,使其信息的含量很高,有效地提高了计算的速度,有利于图像的可靠匹配,使得实时处理成为可能。其在三维场景重建、运动估计、目标跟踪、目标识别、图像配准与匹配等计算机视觉领域起着非常重要的作用。
角点的检测主要有两类基于图像边缘的方法和基于图像灰度的方法。前者很大程度上依赖于图像的分割和边缘提取，一旦待检测目标发生局部变化，很可能导致操作失败，因此该类方法使用范围较小；后者有很多方法，包括Harris算子，Moravec算子，Susan算子等等。
参考英文文献：《A COMBINED CORNER AND EDGE DETECTOR》，1988，Chris Harris & Mike Stephens.
Harris角点检测的实现步骤：
（1）计算图像I(x,y)在x和y两个方向的梯度Ix、Iy.
(2)计算图像两个方向的梯度的乘积.
（3）使用高斯函数对偏导进行高斯加权，生成矩阵M.
(4)计算每个像元的Harris响应值R，并对小于某一阈值的R置为零.
(5)在3*3或5*5的领域内进行非极大值抑制，局部极大值点即为图像中的角点。
相关MATLAB代码有两个版本
function [posr, posc]=Harris1(in_image,a) % 功能：检测图像的Harris角点 % 输入：in_image－待检测的RGB图像数组 % a－角点参数响应，取值范围为：0.04～0.06 % 输出：posr－所检测出角点的行坐标向量 % posc－所检测出角点的列坐标向量 % 将RGB图像转化成灰度图像 in_image=rgb2gray(\'in_image\'); % unit8型转化为双精度double64型 ori_im=double(in_image); %%%%%%计算图像在x、y 两个方向的梯度%%%%%% % x方向梯度算子模板 fx = [-1 0 1]; % x方向滤波 Ix = filter2(fx,ori_im); % y方向梯度算子 fy = [-1;0;1]; % y方向滤波 Iy = filter2(fy,ori_im); %%%%%%计算两个方向的梯度乘积%%%%%% Ix2 = Ix.^2; Iy2 = Iy.^2; Ixy = Ix.*Iy; %%%%%%使用高斯函数对梯度乘积进行加权%%%%%% % 产生7*7的高斯窗函数，sigma=2 h= fspecial(\'gaussian\',[7 7],2); Ix2 = filter2(h,Ix2); Iy2 = filter2(h,Iy2); Ixy = filter2(h,Ixy); %%%%%%计算每个像元的Harris响应值%%%%%% [height,width]=size(ori_im); R = zeros(height,width); % 像素(i,j)处的Harris响应值 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)-a*(trace(M))^2; end end %%%%%%去掉小于阈值的Harris响应值%%%%%% Rmax=max(max(R)); % 阈值 t=0.01* Rmax; for i = 1:height for j = 1:width if R(i,j)<t R(i,j) = 0; end end end %%%%%%进行3×3邻域非极大值抑制%%%%%% % 进行非极大抑制，窗口大小3*3 corner_peaks=imregionalmax(R); countnum=sum(sum(corner_peaks)); %%%%%%显示所提取的Harris角点%%%%%% % posr是用于存放行坐标的向量 [posr, posc] = find(corner_peaks== 1); % posc是用于存放列坐标的向量 figure imshow(in_image) hold on for i = 1 : length(posr) plot(posc(i),posr(i),\'r+\'); end
这个函数程序在运行过程中始终有错误，有待我进一步修改。
可以调用C=cornermetric(I,\’Harris\’)来检测图像的Harris角点特征。其中，I为输入的灰度图像矩阵；C为角点量度矩阵，用来探测I中的角点信息，并与I同尺寸，C的值越大表示图像I中的像素越有可能是一个角点。
相应代码如下：
%%%%%%读入图像并显示%%%%%%
I=imread(\'hua.jpg\');
%将RGB图像转换成灰度图像
I=rgb2gray(I);
subplot(1,3,1);
imshow(I);
title(\'输入图像\');
%%%%%%生成角点度量矩阵并进行调整%%%%%%
%生成角点度量矩阵
C = cornermetric(I,\'Harris\');
C_adjusted = imadjust(C);
subplot(1,3,2);
imshow(C_adjusted);
title(\'角点矩阵\');
%%%%%%寻找并显示Harris角点%%%%%%
corner_peaks = imregionalmax(C);
corner_idx = find(corner_peaks == true);
[r g b] = deal(I);
r(corner_idx) = 255;
g(corner_idx) = 255;
b(corner_idx) = 0;
RGB = cat(3,r,g,b);
subplot(1,3,3);
imshow(RGB);
title(\'检测出的Harris角点\');
下面是第二个MATLAB版本程序。能够运行成功。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Harris角点提取算法 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
%filename = \'Lena.jpg\';
filename=\'xx.png\';
X = imread(filename); % 读取图像
% imshow(X);
Info = imfinfo(filename); %获取图像相关信息
if (Info.BitDepth > 8)
f = rgb2gray(X);
end
%《基于特征点的图像配准与拼接技术研究》
%计算图像亮度f(x,y)在点(x,y)处的梯度-----------------------------------------------
% fx = [5 0 -5;8 0 -8;5 0 -5]; % 高斯函数一阶微分，x方向(用于改进的Harris角点提取算法)
ori_im = double(f) / 255; %unit8转化为64为双精度double64
fx = [-2 -1 0 1 2]; % x方向梯度算子(用于Harris角点提取算法)
Ix = filter2(fx, ori_im); % x方向滤波
% fy = [5 8 5;0 0 0;-5 -8 -5]; % 高斯函数一阶微分，y方向(用于改进的Harris角点提取算法)
fy = [-2; -1; 0; 1; 2]; % y方向梯度算子(用于Harris角点提取算法)
Iy = filter2(fy, ori_im); % y方向滤波
%构造自相关矩阵---------------------------------------------------------------
Ix2 = Ix .^ 2;
Iy2 = Iy .^ 2;
Ixy = Ix .* Iy;
clear Ix;
clear Iy;
h= fspecial(\'gaussian\', [7 7], 2); % 产生7*7的高斯窗函数，sigma=2
Ix2 = filter2(h,Ix2);
Iy2 = filter2(h,Iy2);
Ixy = filter2(h,Ixy);
%提取特征点---------------------------------------------------------------
height = size(ori_im, 1);
width = size(ori_im, 2);
result = zeros(height, width); % 纪录角点位置，角点处值为1
R = zeros(height, width);
Rmax = 0; % 图像中最大的R值
k = 0.06; %k为常系数，经验取值范围为0.04~0.06
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) - k * (trace(M)) ^ 2; % 计算R
if R(i,j) > Rmax
Rmax = R(i, j);
end;
end;
end;
%T = 0.01 * Rmax;%固定阈值，当R(i, j) > T时，则被判定为候选角点
T = 0.1 * Rmax;%固定阈值，当R(i, j) > T时，则被判定为候选角点
%在计算完各点的值后，进行局部非极大值抑制-------------------------------------
cnt = 0;
for i = 2 : height-1
for j = 2 : width-1
% 进行非极大抑制，窗口大小3*3
if (R(i, j) > T && 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;
i = 1;
for j = 1 : height
for k = 1 : width
if result(j, k) == 1;
corners1(i, 1) = j;
corners1(i, 2) = k;
i = i + 1;
end;
end;
end;
[posc, posr] = find(result == 1);
figure,imshow(ori_im);
hold on;
plot(posr, posc, \'r+\');
C语言程序代码如下
// My_corner.cpp : 定义控制台应用程序的入口点。 // #include "stdafx.h" #include<iostream> #include <opencv2/core/core.hpp> #include <opencv2/highgui/highgui.hpp> #include <opencv2/imgproc/imgproc.hpp> #include <GL/glut.h> using namespace std; using namespace cv; void m_filter(double *src,double *&dst,int height,int width,double *mask,int mW,int mH) { double a; for (int i = 0;i < height;i++) { for (int j = 0;j < width;j++) { a = 0.0; //去掉靠近边界的行列，无法滤波，超出范围 if (i > int(mH/2) && i < height - int(mH/2) && j > int(mW) && j < width - int(mW/2)) { for (int m = 0;m < mH;m++) { for(int n = 0;n < mW;n++) { a += src[(i+m-int(mH/2))*width+(j+n-int(mW))]*mask[m*mW+n]; } } } dst[i*width+j] = a; } } } int _tmain(int argc, _TCHAR* argv[]) { Mat img_src = imread("test1.jpg"); Mat img_g,corner; double *src,*img_x,*img_y,*img_x2,*img_y2,*img_xy; Size size = img_src.size(); img_g = Mat(size,CV_8UC1); corner = Mat(size,CV_8UC1); src = new double[size.width*size.height]; img_x = new double[size.width*size.height]; img_y = new double[size.width*size.height]; img_x2 = new double[size.width*size.height]; img_y2 = new double[size.width*size.height]; img_xy = new double[size.width*size.height]; cvtColor(img_src,img_g,CV_BGR2GRAY); for (int i=0;i<size.height;i++) { for (int j=0;j<size.width;j++) { src[i*size.width+j]=(double)img_g.data[i*size.width+j]/255; } } //定义水平方向差分算子并求img_x double dx[25] = {-2,-1,0,1,2,-2,-1,0,1,2,-2,-1,0,1,2,-2,-1,0,1,2,-2,-1,0,1,2}; m_filter(src,img_x,size.height,size.width,dx,5,5); //定义垂直方向差分算子并求img_y double dy[25] = {-2,-2,-2,-2,-2,-1,-1,-1,-1,-1,0,0,0,0,0,1,1,1,1,1,2,2,2,2,2}; m_filter(src,img_y,size.height,size.width,dy,5,5); FILE *fp; //fp=fopen("test.txt","w+"); //计算img_x2、img_y2、img_xy for (int i = 0;i < size.height;i++) { for (int j = 0;j < size.width;j++) { img_x2[i*size.width+j] = img_x[i*size.width+j] * img_x[i*size.width+j]; img_y2[i*size.width+j] = img_y[i*size.width+j] * img_y[i*size.width+j]; img_xy[i*size.width+j] = img_x[i*size.width+j] * img_y[i*size.width+j]; //fprintf(fp,"%f,%f,%f,%f,%f",img_x[i*size.width+j],img_y[i*size.width+j],img_x2[i*size.width+j],img_y2[i*size.width+j],img_xy[i*size.width+j]); //fprintf(fp,"\n"); } } //fclose(fp); //对img_x2、img_y2、img_xy进行高斯平滑,本例中使用5×5的高斯模板 // 计算模版参数 int gausswidth = 7; double sigma = 2; double *g = new double[gausswidth*gausswidth]; for(int i = 0;i < gausswidth;i++) for(int j = 0;j < gausswidth;j++) { g[i*gausswidth+j] = exp(-((i-int(gausswidth/2))*(i-int(gausswidth/2))+(j-int(gausswidth/2))*(j-int(gausswidth/2)))/(2*sigma)); } //归一化：使模板参数之和为1 double total=0; for(int i = 0;i < gausswidth*gausswidth;i++) total += g[i]; for(int i = 0;i < gausswidth;i++) for(int j = 0;j < gausswidth;j++){ g[i*gausswidth+j] /= total; } //进行高斯平滑 m_filter(img_x2,img_x2,size.height,size.width,g,gausswidth,gausswidth); m_filter(img_y2,img_y2,size.height,size.width,g,gausswidth,gausswidth); m_filter(img_xy,img_xy,size.height,size.width,g,gausswidth,gausswidth); //计算角点量R，R(i,j) = det(M) - k * (trace(M)) ^ 2 //Tr(M)=img_x2+img_y2 ,det(M)=img_x2*img_y2-img_xy^2 double *R = new double[size.width*size.height]; double k = 0.06; //k为常系数，经验取值范围为0.04~0.06 double Rmax=0; for(int i = 0; i < size.height; i++) { for(int j = 0; j < size.width; j++) { R[i*size.width+j] = img_x2[i*size.width+j] * img_y2[i*size.width+j] - img_xy[i*size.width+j] * img_xy[i*size.width+j] - k * (img_x2[i*size.width+j] + img_y2[i*size.width+j]) * (img_x2[i*size.width+j] + img_y2[i*size.width+j]); if (R[i*size.width+j] > Rmax) { Rmax = R[i*size.width+j]; } //printf("%f\n",R[i*size.width+j]); } } double max; double T = 0.1 * Rmax;//固定阈值，当R(i, j) > T时，则被判定为候选角点 double *mx = new double[size.width*size.height]; //在计算完各点的值后，进行局部非极大值抑制 for (int i=1;i<size.height-1;i++) { for (int j=1;j<size.width-1;j++) { //写不下了，分两行。。。进行非极大抑制，窗口大小3*3 if (R[i*size.width+j]> T && R[i*size.width+j] > R[(i-1)*size.width+j-1] && R[i*size.width+j] > R[(i-1)*size.width+j] && R[i*size.width+j] >R[(i-1)*size.width+j+1] && R[i*size.width+j] >R[i*size.width+j-1]) { if (R[i*size.width+j]>R[i*size.width+j+1]&&R[i*size.width+j]>R[(i+1)*size.width+j-1]&&R[i*size.width+j]>R[(i+1)*size.width+j]&&R[i*size.width+j]>R[(i+1)*size.width+j+1]) { corner.data[i*size.width+j]=1; } } } } //作图 Scalar sca(0,255,0); for(int i = 0; i < size.height; i++) { for(int j = 0; j < size.width; j++) { if (corner.data[i*size.width+j] == 1) { Point p(j,i); circle(img_src,p,20,sca); } } } imshow("find corner :)",img_src); waitKey(0); free(g);free(R);free(mx); return 0; }
版权声明：本文为natalie原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。