%???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