%???Harris角點提取算法并精確至亞像素級?????????????????????????????

%t=input（‘z1.jpg‘‘s‘）;??%?提示輸入要處理的圖像
tic;?????????????????????????????????????????????????%計時開始
Image?=?imread（‘d:\2.jpg‘）;???????????????????????????????????%?讀取圖像
Image?=?rgb2gray（Image）;?????????????????????????????%?轉(zhuǎn)化為灰度圖像????
HdImage=Image;???????????????????????????????????????%?轉(zhuǎn)化后的灰度圖像
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用來求數(shù)組的平方
Iy2?=?Iy.^2;?????????????????????????????????????????%?.^2用來求數(shù)組的平方
Ixy?=?Ix.*Iy;????????????????????????????????????????%?數(shù)組相乘
h=?fspecial（‘gaussian‘[9?9]2）;?????????????????????%?產(chǎn)生9*9的高斯窗口，sigma=2，產(chǎn)生的窗口越大，得到的角點越少（7-11）
A?=?filter2（hIx2）;??????????????????????????????????%?用產(chǎn)生的高斯窗口處理Ix2得到A?
B?=?filter2（hIy2）;??????????????????????????????????%?用產(chǎn)生的高斯窗口處理Iy2得到B?
C?=?filter2（hIxy）;??????????????????????????????????%?用產(chǎn)生的高斯窗口處理Ixy得到C
height?=?size（HdImage1）;????????????????????????????%?計算圖像的第一維的元素數(shù)，即行數(shù)
width?=?size（HdImage2）;?????????????????????????????%?計算圖像的第二維的元素數(shù)，即列數(shù)
CRF?=?zeros（heightwidth）;???????????????????????????%?生成一個和圖像大小一致的全0的double型數(shù)組，用來保存角點響應函數(shù)值
CRFmax?=?0;??????????????????????????????????????????%?保存圖像中最大的角點響應函數(shù)值?
%M?=?[A（ij）?C（ij）;C（ij）?B（ij）];?????????????
%CRF（ij）?=?det（M）-0.05*（trace（M））^2;????????????????%?計算角點響應函數(shù)值，k的取值一般在0.04--0.06???
CRF=（A.*B?-?C.^2）?-?0.05*（A?+?B）.^2;?????????????????%代碼的優(yōu)化把for循環(huán)改為向量循環(huán)，k=0.05?
CRFmax=max（max（CRF））;????????????????????????????????%找到最大的角點響應函數(shù)（用來設置閾值時用）
l=ordfilt2（CRF7^2ones（7））;?????????????????????????%生成在7*7的窗口進行非最大值抑制（排序濾波器）
k=（l==CRF）&（CRF>0.001*CRFmax）;???????????????????????%設定閾值為0.01*CRFmax只有是局部最大值并且角點響應函數(shù)值大于閾值才是角點
[v?u]?=?find（k）;????????????????????????????????????%?找到角點是的位置，并保存
count?=?size（u1）;???????????????????????????????????%?用來記錄角點的個數(shù)
disp（‘檢測到的角點個數(shù)為：‘）??
disp（count）??????????????????????????????????????????%?輸出角點個數(shù)?
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）
???????????