harris角点检测并精确到亚像素级

自己用matlab编写的Harris角点检测的源代码，并根据附近点的关系把角点位置精确到亚像素级，内附详细说明文档。

%   Harris角点提取算法并精确至亚像素级                             

t=input（‘请输入你要处理的图像全名（含扩展名）:‘‘s‘）;  % 提示输入要处理的图像
tic;                                                 %计时开始
Image = imread（t）;                                   % 读取图像
Image = rgb2gray（Image）;                             % 转化为灰度图像    
HdImage=Image;                                       % 转化后的灰度图像
fx = [-1 0 1;-1 0 1;-1 0 1];                         % x方向的Prewitt算子，用于对x方向滤波 
Ix = filter2（fxHdImage）;                            % 对x方向滤波 
fy = [-1 -1 -1;0 0 0;1 1 1];                         % y方向的Prewitt算子，用于对y方向滤波 
Iy = filter2（fyHdImage）;                            % 对y方向滤波 
Ix2 = Ix.^2;                                         % .^2用来求数组的平方
Iy2 = Iy.^2;                                         % .^2用来求数组的平方
Ixy = Ix.*Iy;                                        % 数组相乘
h= fspecial（‘gaussian‘[9 9]2）;                     % 产生9*9的高斯窗口，sigma=2，产生的窗口越大，得到的角点越少（7-11）
A = filter2（hIx2）;                                  % 用产生的高斯窗口处理Ix2得到A 
B = filter2（hIy2）;                                  % 用产生的高斯窗口处理Iy2得到B 
C = filter2（hIxy）;                                  % 用产生的高斯窗口处理Ixy得到C
height = size（HdImage1）;                            % 计算图像的第一维的元素数，即行数
width = size（HdImage2）;                             % 计算图像的第二维的元素数，即列数
CRF = zeros（heightwidth）;                           % 生成一个和图像大小一致的全0的double型数组，用来保存角点响应函数值
CRFmax = 0;                                          % 保存图像中最大的角点响应函数值 
%M = [A（ij） C（ij）;C（ij） B（ij）];             
%CRF（ij） = det（M）-0.05*（trace（M））^2;                % 计算角点响应函数值，k的取值一般在0.04--0.06   
CRF=（A.*B - C.^2） - 0.05*（A + B）.^2;                 %代码的优化把for循环改为向量循环，k=0.05 
CRFmax=max（max（CRF））;                                %找到最大的角点响应函数（用来设置阈值时用）
l=ordfilt2（CRF7^2ones（7））;                         %生成在7*7的窗口进行非最大值抑制（排序滤波器）
k=（l==CRF）&（CRF>0.001*CRFmax）;                       %设定阈值为0.01*CRFmax只有是局部最大值并且角点响应函数值大于阈值才是角点
[v u] = find（k）;                                    % 找到角点是的位置，并保存
count = size（u1）;                                   % 用来记录角点的个数
disp（‘检测到的角点个数为：‘）  
disp（count）                                          % 输出角点个数 
figureimshow（HdImage）;                              % 显示灰度图像
hold on; 
plot（uv‘r.‘）;                                      % 在灰度图像上用红色‘.‘标出角点的位置

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   以下程序把角点精确到亚像素级，所用窗口为3*3+5*5                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

u1=ones（count1）;                                    %保存改进后的角点坐标
v1=ones（count1）;                                    %保存改进后的角点坐标
HdImage=double（HdImage）; 
for i=1:count
    H=zeros（12）;                                    %保存梯度
    m=zeros（12）;
    n=0;
    if（（u（i）+3<=width）&（v（i）+3<=height）&（u（i）>3）&（v（i）>3））
        for p=（v（i）-1）:（v（i）+1）
            for q=（u（i）-1）:（u（i）+1）
        

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     目录           0  2010-11-26 11:45  Harris角点检测\
     文件        4303  2010-01-01 09:26  Harris角点检测\Harris.m
     文件      272229  2009-12-04 15:19  Harris角点检测\z1.jpg
     文件      266546  2009-12-04 15:19  Harris角点检测\z2.jpg
     文件      278441  2009-12-04 15:20  Harris角点检测\z3.jpg
     文件      275108  2009-12-04 15:20  Harris角点检测\z4.jpg
     文件      281826  2009-12-04 15:20  Harris角点检测\z5.jpg
     文件      264229  2009-12-04 15:21  Harris角点检测\z6.jpg
     文件      140800  2010-01-01 09:30  Harris角点检测\魏青 0710710 Harris角点检测.doc