关于Harris角点检测的实现
来源:互联网 发布:手机淘宝改中差评 编辑:程序博客网 时间:2024/06/05 23:02
harris角点检测算法非常有名,网上也有一堆源代码,这里自己贴下自己找到的一些资料,以备将来查看。
1.matlab中有内置的corner函数可以实现harris焦点的提取
2.网上也有很多博文中贴出了源代码,这里我们都做了测试,至于网上源代码中 的这段,研究了好久都没有研究出来,为啥这个高斯函数的一阶微分会是这个结果呢,不知道呀2333333333333
如果有大神知道,请劳烦告知一下,谢谢了
fx = [5 0 -5;8 0 -8;5 0 -5]; % 高斯函数一阶微分,x方向(用于改进的Harris角点提取算法)
3.自己参考大神们的源码自己也写了一个附在其中,其中也遇到一些问题,比如imfilter, conv2, filter2之间的关系,简单来讲conv2是卷积操作,filter2是相关操作,imfilter理论上都能实现以上两者的效果,在我们放出的代码中的test部分就是对matlab论坛中一个帖子的测试代码了。
% -----------------【harris corner detect】--------------% author : zhyh2010% date : 20150511% target : harris corner detect% --------------------------【end】----------------------function main % clc % close all % clear all %harris_matlab_inside;% test myharris_1;% harris_refer1;% harris_refer2;end% -----------------【matlab 内置corner函数进行角点检测】--------------% method : corner 只支持单通道% --------------------------【end】----------------------function harris_matlab_inside src = imread('1.jpg'); figure, imshow(src) title('original pic') src_corner = corner(rgb2gray(src), 'Harris', 1000); figure, imshow(src); hold on plot(src_corner(:, 1), src_corner(:, 2), 'b*') hold offend% -----------------【matlab conv2 filter2 imfilter 区别】--------------% 检验博文 http://www.ilovematlab.cn/thread-293710-1-1.html% --------------------------【end】----------------------function test A=[4 3 1 2;0 1 1 3;5 2 0 0] B=[1 2 3;0 -1 2;1 1 0] % kernel C = conv2(A, B, 'full') D = filter2(B, A, 'full') E_1 = imfilter(A, B, 0, 'full', 'conv') E_2 = imfilter(A, B, 0, 'full', 'corr')end% -----------------【myharris_1】--------------% 根据原理自行实现 harris 角点检测% --------------------------【end】----------------------function myharris_1 src = imread('1.jpg'); img = rgb2gray(src); h_x = [-2 -1 0 1 2];% Ix = imfilter(img, h_x, 0, 'full', 'corr');% % Ix = imfilter(img, h_x, 0, 'same', 'corr');% Ix = Ix(2:end-1, 2:end -1); Ix = filter2(h_x, img); h_y = [-2 -1 0 1 2]';% % Iy = imfilter(img, h_y, 0, 'same', 'corr');% % Iy = imfilter(img, h_y, 0, 'full', 'corr');% Iy = Iy(2:end-1, 2:end -1); Iy = filter2(h_y, img); Ix2 = Ix.^2; Iy2 = Iy.^2; Ixy = Ix.*Iy; kernel_size = 3; kernel_sigma = 2; h = fspecial('gaussian', kernel_size, kernel_sigma); Ix2 = filter2(h, Ix2); Iy2 = filter2(h, Iy2); Ixy = filter2(h, Ixy); pic_width = size(img, 2); pic_height = size(img, 1); alpha = 0.05; Judge_standard = zeros(size(img)); result = zeros(size(img)); for i = 1 : pic_height for j = 1 : pic_width M = [Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)]; Judge_standard(i, j) = det(M) - alpha * trace(M)^2; end end MaxJS = max(Judge_standard(:)); R = Judge_standard; threshold = 0.1 * MaxJS; for i = 2 : pic_height - 1 for j = 2 : pic_width - 1 if (Judge_standard(i, j) > threshold) ... && (bigger_around(Judge_standard, i, j)) result(i, j) = 1; end end end [posx, posy] = find(result == 1); fprintf('---------------find %d corners---------------\n', length(posx)); figure, imshow(src) hold on plot(posy, posx, 'b*') hold offend% -----------------【myharris_1 接口子函数】--------------% 本质 3 X 3 局部抑制操作% isbigger = bigger_around(A, i, j)% A : 存储角点评判标准的值% i,j : 当前坐标 % attention : 2 <= i <= size(A, 1) - 1% 2 <= j <= size(A, 2) - 1% --------------------------【end】----------------------function isbigger = bigger_around(A, i, j) assert(i >= 2 && i <= size(A, 1) - 1); assert(j >= 2 && j <= size(A, 2) - 1); cond = A(i, j) > A(i - 1, j - 1); cond = cond && (A(i, j) > A(i - 1, j)); cond = cond && (A(i, j) > A(i - 1, j + 1)); cond = cond && (A(i, j) > A(i, j - 1));% cond = cond && (A(i, j) > A(i, j)); cond = cond && (A(i, j) > A(i, j + 1)); cond = cond && (A(i, j) > A(i + 1, j - 1)); cond = cond && (A(i, j) > A(i + 1, j)); cond = cond && (A(i, j) > A(i + 1, j + 1)); isbigger = cond;end% -----------------【别人的harris算法实现1】--------------% 参考地址: http://blog.sina.com.cn/s/blog_6c73951701017bwi.html% --------------------------【end】----------------------function harris_refer1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Harris角点提取算法 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear; %filename = 'Lena.jpg'; filename='1.jpg'; 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', [3 3], 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+');end% -----------------【别人的harris算法实现2】--------------% 参考地址: http://blog.csdn.net/makenothing/article/details/12884331% --------------------------【end】----------------------function harris_refer2 %MatLab角点检测程序harris。 ori_im2=rgb2gray(imread('1.jpg')); %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); 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; [posr2, posc2] = find(result == 1); cnt % compter des coins figure imshow(ori_im2); hold on; plot(posc2,posr2,'w*'); end% -----------------【别人的harris算法实现3】--------------% 参考地址: http://blog.csdn.net/makenothing/article/details/12884331% --------------------------【end】----------------------function harris_refer3 %%%Prewitt Operator Corner Detection.m %%%时间优化--相邻像素用取差的方法求Harris角点 %% %clear; Image = imread('1.jpg'); % 读取图像 Image = im2uint8(rgb2gray(Image)); dx = [-1 0 1;-1 0 1;-1 0 1]; %dx:横向Prewitt差分模版 Ix2 = filter2(dx,Image).^2; Iy2 = filter2(dx',Image).^2; Ixy = filter2(dx,Image).*filter2(dx',Image); %生成 9*9高斯窗口。窗口越大,探测到的角点越少。 h= fspecial('gaussian',9,2); A = filter2(h,Ix2); % 用高斯窗口差分Ix2得到A B = filter2(h,Iy2); C = filter2(h,Ixy); nrow = size(Image,1); ncol = size(Image,2); Corner = zeros(nrow,ncol); %zeros用来产生一个全零矩阵,故矩阵Corner用来保存候选角点位置,初值全零,值为1的点是角点 %参数t:点(i,j)八邻域的“相似度”参数,只有中心点与邻域其他八个点的像素值之差在 %(-t,+t)之间,才确认它们为相似点,相似点不在候选角点之列 t=20; %我并没有全部检测图像每个点,而是除去了边界上boundary个像素,也就是从第8行第8列开始遍历。 %因为我们感兴趣的角点并不出现在边界上 %个人觉得这一部分是的主要目的是找出可能是角点的点,缩小范围,加快运算速度。 %具体思想是如果中心点(i,j)周围8个点中有7、8个点灰度值与之相似,那么该中心点应该处于平坦区域,不可能为角点, %如果中心点(i,j)周围只有1个点或者没有点与之相似,那么该中心点也不可能为角点。 boundary=8; for i=boundary:nrow-boundary+1 for j=boundary:ncol-boundary+1 nlike=0; %相似点个数 if Image(i-1,j-1)>Image(i,j)-t && Image(i-1,j-1)<Image(i,j)+t nlike=nlike+1; end if Image(i-1,j)>Image(i,j)-t && Image(i-1,j)<Image(i,j)+t nlike=nlike+1; end if Image(i-1,j+1)>Image(i,j)-t && Image(i-1,j+1)<Image(i,j)+t nlike=nlike+1; end if Image(i,j-1)>Image(i,j)-t && Image(i,j-1)<Image(i,j)+t nlike=nlike+1; end if Image(i,j+1)>Image(i,j)-t && Image(i,j+1)<Image(i,j)+t nlike=nlike+1; end if Image(i+1,j-1)>Image(i,j)-t && Image(i+1,j-1)<Image(i,j)+t nlike=nlike+1; end if Image(i+1,j)>Image(i,j)-t && Image(i+1,j)<Image(i,j)+t nlike=nlike+1; end if Image(i+1,j+1)>Image(i,j)-t && Image(i+1,j+1)<Image(i,j)+t nlike=nlike+1; end if nlike>=2 && nlike<=6 Corner(i,j)=1;%如果周围有2~6个相似点,那(i,j)就是角点 end; end; end; CRF = zeros(nrow,ncol); % CRF用来保存角点响应函数值,初值全零 CRFmax = 0; % 图像中角点响应函数的最大值,作阈值之用 k=0.05; % 计算CRF %工程上常用CRF(i,j) =det(M)/trace(M)计算CRF,那么此时应该将下面第105行的 %比例系数k设置大一些,k=0.1对采集的这几幅图像来说是一个比较合理的经验值 for i = boundary:nrow-boundary+1 for j = boundary:ncol-boundary+1 if Corner(i,j)==1 %只关注候选点 M = [A(i,j) C(i,j); C(i,j) B(i,j)]; CRF(i,j) = det(M)-k*(trace(M))^2; if CRF(i,j) > CRFmax CRFmax = CRF(i,j); end; end end; end; %CRFmax count = 0; % 用来记录角点的个数 t=0.01; % 下面通过一个3*3的窗口来判断当前位置是否为角点 for i = boundary:nrow-boundary+1 for j = boundary:ncol-boundary+1 if Corner(i,j)==1 %只关注候选点的八邻域 if CRF(i,j) > t*CRFmax && CRF(i,j) >CRF(i-1,j-1) ......%?????为什么要CRF(i,j) > t*CRFmax啊?求大神告知 && CRF(i,j) > CRF(i-1,j) && CRF(i,j) > CRF(i-1,j+1) ...... && CRF(i,j) > CRF(i,j-1) && CRF(i,j) > CRF(i,j+1) ...... && CRF(i,j) > CRF(i+1,j-1) && CRF(i,j) > CRF(i+1,j)...... && CRF(i,j) > CRF(i+1,j+1) count=count+1;%这个是角点,count加1 else % 如果当前位置(i,j)不是角点,则在Corner(i,j)中删除对该候选角点的记录 Corner(i,j) = 0; end; end; end; end; % disp('角点个数'); % disp(count) figure,imshow(Image); % display Intensity Image hold on; % toc(t1) for i=boundary:nrow-boundary+1 for j=boundary:ncol-boundary+1 column_ave=0; row_ave=0; k=0; if Corner(i,j)==1 for x=i-3:i+3 %7*7邻域 for y=j-3:j+3 if Corner(x,y)==1 % 用算数平均数作为角点坐标,如果改用几何平均数求点的平均坐标,对角点的提取意义不大 row_ave=row_ave+x; column_ave=column_ave+y; k=k+1; end end end end if k>0 %周围不止一个角点 plot( column_ave/k,row_ave/k ,'g.'); end end; end; end% http://blog.csdn.net/jinshengtao/article/details/17717731
参考文献
1.imfilter 与conv2 和 filter2 之间的区别
2.harris源码1
3.harris源码2
0 0
- 关于Harris角点检测的实现
- 优化的Harris角点检测算法的MATLAB实现
- vc实现harris角点检测算法
- harris角点检测算法实现VC
- OpenCV2.4.4实现Harris角点检测
- Harris角点检测及实现
- OpenCV2.4.4实现Harris角点检测
- Python 实现Harris 角点检测
- Harris角点检测原理及实现
- Harris角点检测实现详解
- Harris 角点检测
- Harris角点检测
- Harris角点检测
- Harris角点检测
- Harris角点检测
- harris角点检测
- Harris角点检测
- Harris角点检测
- 嵌入式WIFI与普通WIFI的区别
- POJ 2226 Muddy Fields
- wifi网络接入原理(上)——扫描Scanning
- 八大排序算法
- Matlab中的静态变量和全局变量
- 关于Harris角点检测的实现
- wifi网络接入原理(中)——认证Authentication
- c#中的ref、out、params参数
- POJ 3186 Treats for the Cows
- 安卓博客
- base64加密代码
- 日常养胃三大误区
- seo属于哪个行业
- wifi网络接入原理(下)——关联Association