計算機視覺之（一）利用Harris檢測子進行角點特徵檢測（含matlab原始碼）
好的點特徵如圖1所示，在角點處，影象上的一個小視窗沿任何方向移動，視窗內的影象都有明顯的灰度變化，因此角點是好的點特徵。
圖1
那麼，我們要怎樣識別角點呢？由圖2可以看出，角點具有類似圖2中的第一列所展示的性質：
圖2
假設視窗W發生位置偏移(u,v)；比較偏移前後視窗中每一個畫素點的灰度變化值；使用灰度誤差平方和來構造一個誤差函式E(u,v)，其中的視窗函式是用來濾波的。
H稱為自相關矩陣, λmax和λmin是自相關矩陣的特徵值。其中E（u,v）是一個二次型函式，二次型函式的本質就是一個橢圓，橢圓的扁率和尺寸是由H的特徵λmax和λmin值λmax和λmin決定的，橢圓的方向由H的特徵向量決定。圖3-3表示的是點線面與橢圓的關係，其中（d）圖中的橢圓跑哪去了，我也不知道。因為圖是從掃描版的PDF中擷取的，有點模糊。
下面這幅圖是王永明等在《影象區域性不變性特徵與描述》一書中闡述橢圓與點線面的關係用的。第一行是不同典型影象的灰度梯度分佈圖，第二行是對梯度分別的橢圓擬合。不過有一點沒有明白，如果按照哈里斯的理論線性邊緣時，R不應該為負值的嗎？那王永明的這幅圖中的R代表什麼含義啊到底？為什麼線性邊緣時，R=0.3328？搞不懂啊搞不懂！！！
正如圖3所示，只有當兩者都比較大，並且大小相當時對應點才為角點，兩者都非常小時為平坦區域；一大一小時為邊界。
圖3
在1988年，哈里斯在其論文《A combined corner and edge detector》裡給出了更有效的角點響應函式
R為正值時，檢測到的是角點；R為負時檢測到的是邊；R很小時檢測到的是平坦區域，由此也就有了更便於計算的數學公式。哈老先生對此做的圖如下：
下面是對Harris演算法的Matlab實現，原始碼參考自http://blog.csdn.net/aflyeagle/article/details/5116799
個人對一些地方做了修改說明。
執行環境windows8.1+Matlab R2013b
%%%Prewitt Operator Corner Detection.m
%%%時間優化--相鄰畫素用取差的方法求Harris角點
%%
clear;
Image = imread('884.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;
採用國際會議中心的圖片進行實驗，效果如下：
原圖：
標記效果圖：
為了驗證Harris的旋轉不變性，將原圖旋轉後實驗效果如下：
