<span style="font-size:10px; font-family: Arial, Helvetica, sans-serif; background-color: rgb(255, 255, 255); font-weight: normal;">之前在斯坦福的機器視覺演算法與應用課程上學了一些東西，並用matlab程式設計實現，沒來得及整理，現在把它整理一下，出來的效果可能不太完善，這有待後續的研究與改進。</span>

一.特徵檢測（提取）

        基於特徵的影象配準方法是影象配準中最常見的方法之一。它不是直接利用影象畫素值，二十通過畫素值匯出的符號特徵（如特徵點、特徵線、特徵區域）來實現影象配準，因此可以克服利用灰度資訊進行影象配準的缺點，主要體現在以下三個方面：（1）利用特徵點而不是圖

像灰度資訊，大大減少了在匹配過程中的計算量；（2）特徵點的匹配度量值相對位置變化比較敏感，可以提高匹配的精度；（3）特徵點的提取過程可以減少噪聲的影響，對灰度變化、影象形變以及遮擋等都有較好的適應能力。

       一類重要的點特徵：角點（corner points），其定義主要有以下：

1.  區域性視窗沿各方向移動，灰度均產生明顯變化的點      
2. 影象局部曲率突變的點
3. 典型的角點檢測演算法：Harris角點檢測、CSS角點檢測
4. Harris角點檢測基本思想
   從影象區域性的小視窗觀察影象特徵，角點定義：視窗向任意方向的移動都導致影象灰度的明顯變化（如下圖）

Harris檢測：數學表達

將影象視窗平移[u,v]產生灰度變化E(u,v)

由泰勒展開，得：

利用角點響應函式：
判斷特徵點是否為角點的依據：R只與M值有關，R為大數值正數時特徵點為角點，R為大數值負數時為邊緣，R為小數值時為平坦區，如下圖：
尋找R位於一定閾值之上的區域性最大值，去除偽角點。 Harris角點檢測流程：
1.通過高斯函式的導數對原始影象進行卷積計算；影象在水平方向和垂直方向的導數Ix和Iy；
2.計算對應這些梯度外積即（Ix2 、Iy2、IxIy）三個影象如下圖：
4.使用高斯函式對以上影象進行卷積濾波；
3.使用前面的公式計算角點響應函式R值；
5.對計算到的角點影象進行區域性極大值抑制。
以下是我的特徵檢測程式碼，由於程式設計能力有限，可能過程中會出現程式碼不規範以及不完善的情況，望見諒！

function [x, y, confidence, scale, orientation] = get_interest_points(image,feature_width)

feature_width = 16;
fx = [-1,0,1;-2,0,2;-1,0,1];
fy = [-1,-2,-1;0,0,0;1,2,1];
Ix = imfilter(image,fx);
Iy = imfilter(image,fy);
Ix2 = Ix.* Ix;
%subplot(2,2,1);
%imshow(Ix2);
Ixy = Ix.* Iy;
% subplot(2,2,2);
% imshow(Ixy);
Iy2 = Iy.* Iy;
% subplot(2,2,3);
% imshow(Iy2);
h = fspecial('gaussian',25,2);
% subplot(2,2,4);
% imshow(h);
Ix2 = imfilter(Ix2,h);
Ixy = imfilter(Ixy,h);
Iy2 = imfilter(Iy2,h);

%尋找最大R值
Rmax = 0;
l = 0.05;
[img_height,img_width] = size(image);
M = zeros(img_height,img_width);
for i = 1:img_height
    for j = 1:img_width
        A = [Ix2(i,j),Ixy(i,j);Ixy(i,j),Iy2(i,j)];
        M(i,j) = det(A) -  l*(trace(A))^2;
        if M(i,j) > Rmax
            Rmax = M(i,j);
        end
    end
end


%區域性非極大值抑制
k = 0.01;
cnt = 0;   %記錄角點數目
result = zeros(img_height,img_width);
for i=2:img_height-1
    for j =2:img_width-1
        if M(i,j)>k*Rmax && M(i,j)>M(i-1,j-1) && M(i,j)>M(i-1,j)&&M(i,j)>M(i-1,j+1)&&M(i,j)>M(i,j-1)&&M(i,j)>M(i,j+1)&&M(i,j)>M(i+1,j-1)&&M(i,j)>M(i+1,j)&&M(i,j)>M(i+1,j+1)
            result(i,j) = 1;
            cnt = cnt + 1;
        end
    end
end

[y,x] = find(result==1);
x_indexleft = find(x <= feature_width/2 - 1);   %尋找x方向上左邊邊緣角點
%去除邊緣角點
x(x_indexleft) = [];
y(x_indexleft) = [];
y_indexup = find(y <= feature_width/2 - 1);    %尋找y方向上部邊緣角點
x(y_indexup) = [];
y(y_indexup) = [];
x_indexright = find(x > img_width - 8);       %尋找x方向右邊的邊緣角點
x(x_indexright) = [];
y(x_indexright) = [];
y_indexdown = find(y > img_height - 8);       %尋找y方向下部的邊緣角點
x(y_indexdown) = [];
y(y_indexdown) = [];
%x(find(x<feature_width/2 - 1)) = [];
%y(find(y<feature_width/2 - 1)) = [];
disp(length(x));    %顯示角點數目
figure;
imshow(image);
hold on;
plot(x,y,'ro');
hold off;
end
 

二、特徵描述 

 在檢測到特徵（關鍵點）之後，我們必須匹配它們，也就是說，我們必須確定哪些特徵來自於不同影象中的對應位置。物體識別的核心問題是將同一目標在不同時間、不同解析度、不同光照、不同位姿情況下所成的影象相匹配。而為了進行匹配，我們首先要合理的表示影象。

SIFT（Scale invariant feature transform）特徵通過計算檢測到的關鍵點周圍16x16視窗內每一個畫素的梯度得到。在這裡我只是簡單的實現類似於SIFT特徵描述子的特徵描述方法，即我通過每4x4的四分之一象限，通過將加權梯度值加到直方圖八個方向區間中的一個，計算出一個梯度方向直方圖，因此在每一個特徵點都會形成一個128維的非負值形成了一個原始版本的SIFT描述子向量如下圖，並且將其歸一化以減少對比度和增益的影響，最後為了使描述子對其他各種光度變化魯棒，再將這些值以0.2截尾，然後再歸一化到單位長度。

function [features] = get_features(image, x, y, feature_width)
y_gradient = imfilter(image,[-1,0,1]');
features = zeros(length(x),128);
for k = 1:length(x)                         %window of 16x16 of each feature 
    x_subwindow = x(k)-7:x(k)+8;   
    y_subwindow = y(k)-7:y(k)+8;
    x_subwindow_gradient = x_gradient(y_subwindow,x_subwindow);
    y_subwindow_gradient = y_gradient(y_subwindow,x_subwindow);
    angles = atan2(y_subwindow_gradient,x_subwindow_gradient); 
   %for i = 1:16
       %for j = 1:16
          % if angles(i,j)<0
              % angles(i,j) = angles(i,j) + 2*pi;
           %end
       %end
   %end
   angles(angles<0) = angles(angles<0) + 2*pi;
   %angles_bin = [0 pi/4 pi/2 3*pi/4 pi 5*pi/4 3*pi/2 7*pi/4 2*pi];
   angles_binranges = 0:pi/4:2*pi;
   B = zeros(1,128);
   for i = 1:feature_width/4:feature_width
       for j = 1:feature_width/4:feature_width
          %A = reshape(angles(i:i+3,j:j+3)',1,16);
          subwindow = angles(i:i+3,j:j+3);
          angles_bincounts = histc(subwindow(:),angles_binranges);
          begin = 1 + 32*(floor(i/4)) + 8*(floor(j/4));  %tab the orientation of each 4x4 window
          B(begin:begin+7) = angles_bincounts(1:8);    %get the eight orientation of each window
          %figure;
          %bar(angles_binranges,angles_bincounts,'histc');
       end
   end
   %disp(length(B/norm(B)));
   features(k,:) = B/norm(B);
end
%cut the end
power = 0.8;
features = features.^power;
end

三、特徵匹配

一旦我們從兩幅或者多幅影象中提取到特徵及其描述子後，下一步就是要在這些影象之間建立一些初始特徵之間的匹配。
匹配策略一：對前面提取到的兩幅影象的128維特徵描述子向量做歐式距離度量，最簡單的一個策略就是先設定一個閾值（最大距離），然後返回在這個閾值範圍之內的另外一個影象中的所有匹配。
匹配策略二：做最近鄰匹配，即比較最近鄰距離和次近鄰距離的比值，即最近鄰比率（NNDR）。 匹配策略一的缺點是，如果閾值設得太高，就會產生誤報，也就是說會出現不正確的匹配。如果閾值設得太低，就會產生很多“漏報”，也就是說，很多正確的匹配被丟失。 固定閾值，最近鄰和最近鄰比率匹配。在固定閾值（虛線圓）下，描述子DA未能與DB匹配，DD錯誤地與DC和DE匹配。如果我們選擇最近鄰，DA和DE匹配。使用最近鄰比率（NNDR），小的NNDR（d1/d2）正確地將DA和DB匹配，大的NNDR（d1'/d2'）正確地拒絕DD與DC、DE的匹配。 程式碼如下：

function [matches, confidences] = match_features(features1, features2)
ratio = 0.8;
num_matches = 0;
num_features1 = size(features1,1);
num_features2 = size(features2,1);
matches = zeros(num_features1,2);
confidences = zeros(num_features1);
for i = 1:num_features1
    distances = dist(features1(i,:),features2');
    [dists,index] = sort(distances);
    nndr = dists(1)/dists(2); 
    if nndr < ratio
      %num_matches = num_matches + 1;
      matches(i,:) = [i,index(1)];
      confidences(i) = 1 - nndr;
%     else 
%         confidences(i) = [];
%         matches(i,:) = [];
    end
end

%keep only those matched
matches_index = find(confidences>0);
matches = matches(matches_index,:);
confidences = confidences(matches_index);

% Sort the matches so that the most confident onces are at the top of the
% list. You should probably not delete this, so that the evaluation
% functions can be run on the top matches easily.
[confidences, ind] = sort(confidences, 'descend');
matches = matches(ind,:);

出來的效果如下：