gusucode.com > KLT Tracking源码程序matlab > harris.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Harris角点提取算法 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [posc,posr1]=harris(X) %filename = 'Lena.jpg'; filename = 'Lena.jpg'; X = imread(filename); % 读取图像 % imshow(X); Info = imfinfo(filename); %获取图像相关信息 % if (Info.BitDepth > 8) % f = rgb2gray(X); % end f = X; %《基于特征点的图像配准与拼接技术研究》 %计算图像亮度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 k2 = 0.983; 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.001 * 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); posr1 = k2; % figure,imshow(ori_im); % hold on; % plot(posr, posc, 'r+');