Matlab圖像處理學習筆記（八）：用廣義霍夫變換篩選sift特征點

經過幾天的學習研究，終於完成了廣義霍夫變換（Generalised Hough transform）對特征點的篩選。此法不僅僅針對sift特征點，surf，Harris等特征點均可適用。

這幾天我發現關於廣義霍夫變換的資料少之又少，不過經過仔細研讀各方的資料，我對Generalised Hough transform也有了一點認識，如果有理解的不對的地方，還請指正。

在介紹Generalised Hough transform之前，先簡要介紹一下霍夫變換。我們都知道，用霍夫變換可以檢測直線，圓，橢圓等，或者說只要該形狀是可解析的，都可以用霍夫變換來進行識別。我認為，霍夫變換的關鍵點有兩點：

1、投票機制的引入。

2、參數空間的轉換。

舉個例子，好比一條線段，在x-y坐標系下，一條直線的特征並不是那么明顯的，但在轉換到p，theta參數空間后，一條線段的特征就很明顯，再加上投票機制，如果一個累加單元中有較大值，則可以判定存在一條線段。圓也一樣。但直線、圓、橢圓等這些形狀其實都是可解析的，也就是可以用一個方程式來表達，如果現在有一個不規則形狀，他就沒辦法了。因此，D.H. Ballard在1981年將霍夫變換進行推廣，提出了廣義霍夫變換（Generalised Hough transform）。好了，下面進入本文的關鍵，廣義霍夫變換（Generalised Hough transform）。

本文的主要參考資料如下：

%referrence：
% David G. Lowe，Distinctive Image Features from Scale-Invariant Keypoints
% David G. Lowe，Object Recognition from Local Scale-Invariant Features
% D.H. Ballard, "Generalizing the Hough Transform to Detect Arbitrary Shapes"
% wiki：http://en.wikipedia.org/wiki/Generalised_Hough_transform#cite_note-1

轉載請注明出處：http://blog.csdn.net/u010278305

廣義霍夫變換的目的是要在一個具體場景中尋找一個不規則已知物體的位置，我們把這個不規則已知物體叫做模板，這個已知物體的模板由一系列點組成。通過參數空間轉換，我們把目標轉換為尋找一個轉換矩陣，通過該矩陣，可以使得模板與場景中模板出現的位置相匹配。只要我們確定了這個轉換矩陣的參數，那么模板在場景中的位置也就確定了。

那么這個轉換矩陣的參數如何確定呢？這就要通過投票機制，如果一個參數獲得的投票最多，那么就認為這個參數是正確的。

一般情況下，要進行Generalised Hough transform，先要選取一個參考點，然后將普通的坐標轉換為(r，, ɸ)這種形式，r為邊界點到選取參考點的長度，ɸ為參考點到邊界點的連線與邊界點切線的夾角。如下圖所示：

之后建立一個R-table，列出邊界點上所有點的角度與其長度的對應關系，這樣便能完整的表述我們的模板物體。R-table的建立過程如下：

R-table如下：

在之后的步驟中，我們可以通過角度來尋找匹配，再用計算出的值填充累加單元。但由於我們已經尋得各特征點的匹配關系，在實施的過程中，我們跳過這一步。直接進入下一步。

該算法的具體實施過程如下：

1、選取模板的（0，0）為參考點，代表模板的位置。模板的每一個特征點與參考點相減。（由於參考點選為（0，0），故各特征點不變。

2、初始化累加表.A[xcmin. . . xcmax][ycmin. . . ycmax][qmin . . .qmax][smin . . . smax]。其中q代表旋轉角度，s代表縮放尺度。在A的大小選取上，我們按David G. Lowe給出的建議，We use broad bin sizes of 30 degrees for orientation, a factor of 2 for scale, and 0.25 times the maximum projected training image dimension (using the predicted scale) for location.

3、遍歷場景中的每一點，記x',y'為模板中特征點的坐標。則可以用下式推出場景中模板的坐標，並對累加單元進行累加。

4、找到所有累加單元的最大值，並找到其索引，則其對應的（xc，yc）為模板中（0，0）點在場景中對應的坐標。

5、找到沒有為最大累加單元做出貢獻的點，並認為它是奇異點，進行排除。

6、存儲剩下的點。

OK，Generalised Hough transform過程到此為止。之后，我按照David G. Lowe在Distinctive Image Features from Scale-Invariant Keypoints一文中給出的說明，對最終的匹配點進行仿射變換，求出仿射變換的參數。模板到場景坐標的轉換矩陣為：

我們將上式記為Ax=b；那么我們的目標就是求解x，x的求解方式如下：

到此，我們完成了所有工作。

下面給出matlab源代碼：

%function:
%       用廣義霍夫變換（Generalised Hough transform）篩選sift得到的特征點。
%       並用最小平方法（The least-squares solution）求取仿射變換的參數
%注意：
%      由於matlab沒有自帶sift特征提取，sift特征提取調用了該算法作者提供的底層調用。
%referrence：
%      David G. Lowe，Distinctive Image Features from Scale-Invariant Keypoints
%      David G. Lowe，Object Recognition from Local Scale-Invariant Features
%      D.H. Ballard, "Generalizing the Hough Transform to Detect Arbitrary Shapes"
%      wiki：http://en.wikipedia.org/wiki/Generalised_Hough_transform#cite_note-1
%date:2015-1-15
%author:chenyanan
%轉載請注明出處：http://blog.csdn.net/u010278305

%清空變量，讀取圖像
clear;close all
fprintf('/******************************\n**It''s writed by chenyn2014.\n******************************/\n');

%讀取模板
rgb2=imread ('head.jpg');
gray2=rgb2gray(rgb2);
imwrite(gray2,'scene2.jpg');

%讀取場景
rgb1=imread ('girl.jpg');
gray1=rgb2gray(rgb1);

%對場景進行仿射變換，以驗證sift特征的尺度不變性、旋轉不變性
scale_value=0.7;
scale_tform=[scale_value 0 0;0 scale_value 0;0 0 1];
rotate_value=2*pi/30;
rotate_tform=[cos(rotate_value) sin(rotate_value) 0;-sin(rotate_value) cos(rotate_value) 0;0 0 1];
shift_x=0;  shift_y=0;
shift_tform=[1 0 0; 0 1 0;shift_x shift_y 1];
tform = affine2d(scale_tform*rotate_tform*shift_tform);
gray1 = imwarp(gray1,tform);
imwrite(gray1,'scene1.jpg');

image1='scene1.jpg';
image2='scene2.jpg';

% Find SIFT keypoints for each image
%keypoint location (row, column, scale, orientation)
[im1, des1, loc1] = sift(image1);
[im2, des2, loc2] = sift(image2);

% For efficiency in Matlab, it is cheaper to compute dot products between
%  unit vectors rather than Euclidean distances.  Note that the ratio of 
%  angles (acos of dot products of unit vectors) is a close approximation
%  to the ratio of Euclidean distances for small angles.
%
% distRatio: Only keep matches in which the ratio of vector angles from the
%   nearest to second nearest neighbor is less than distRatio.
distRatio = 0.6;   

% For each descriptor in the first image, select its match to second image.
des2t = des2';                          % Precompute matrix transpose
match_point=1;
for i = 1 : size(des1,1)
   dotprods = des1(i,:) * des2t;        % Computes vector of dot products
   [vals,indx] = sort(acos(dotprods));  % Take inverse cosine and sort results

   % Check if nearest neighbor has angle less than distRatio times 2nd.
   if (vals(1) < distRatio * vals(2))
      %將有效的匹配點存進數組，方便后邊處理
      matchedPoints1(match_point,:)=[loc1(i,1) loc1(i,2)];
      matchedPoints2(match_point,:)=[loc2(indx(1),1) loc2(indx(1),2)];
      match_point=match_point+1;
   end
end

% Create a new image showing the two images side by side.
im3 = appendimages(im1,im2);

%繪制剔除奇異點之前的結果
% Show a figure with lines joining the accepted matches.
figure('name','original','Position', [100 100 size(im3,2) size(im3,1)]);
colormap('gray');
imagesc(im3);title('original');
hold on;
cols1 = size(im1,2);
for i = 1: size(matchedPoints1,1)
    line([matchedPoints1(i,2) matchedPoints2(i,2)+cols1], ...
         [matchedPoints1(i,1) matchedPoints2(i,1)], 'Color', 'c');
end
hold off;
num = size(matchedPoints1,1);
fprintf('Found %d matches.\n', num);


%Initialize the Accumulator table（初始化累加表A）
%由於我們已知點之間的對應關系，故不再構造R-table;
%（個人認為R-table的作用就是通過中心點與邊緣點的連線r與該點切線夾角尋找對應關系，故我們沒必要再構建）
%David G. Lowe建議，We use broad bin sizes of 30 degrees for orientation, a factor
%of 2 for scale, and 0.25 times the maximum projected training image dimension (using the
%predicted scale) for location.
A=zeros(4,4,12,5);
%For each edge point (x, y)（遍歷每個場景中的特征點）
for i=1:size(matchedPoints1,1)
    for sita_index=1:12
        sita=2*pi/12*(sita_index-1);
        for scale_index=1:5
            scale=(2^(scale_index-1))*0.25;
            %Using the gradient angle ɸ, retrieve from the R-table all the (α, r) values indexed under ɸ.
            %用R-table中的夾角ɸ尋找與模板中哪幾個點匹配，此處我們已知對應關系，直接應用即可
            %matchedPoints2(i,1) :row  y
            %matchedPoints2(i,2) :col  x
            %用該遍歷時刻的尺度、旋轉變換尋找在場景中目標的坐標
            xc=matchedPoints1(i,2)-(matchedPoints2(i,2)*cos(sita)-matchedPoints2(i,1)*sin(sita))*scale;
            yc=matchedPoints1(i,1)-(matchedPoints2(i,1)*sin(sita)+matchedPoints2(i,2)*cos(sita))*scale;
            x=0;y=0;
            if(xc<=0.25*size(im1,2))
                x=1;
            elseif (xc<=0.5*size(im1,2))
                x=2;
            elseif (xc<=0.75*size(im1,2))
                x=3;
            elseif (xc<=size(im1,2))
                x=4;
            end
                    
            if(yc<=0.25*size(im1,1))
                y=1;
            elseif (yc<=0.5*size(im1,1))
                y=2;
            elseif (yc<=0.75*size(im1,1))
                y=3;
            elseif (yc<=size(im1,1))
                y=4;
            end
            
            %對累加表進行更新。
            if(x>=1&&x<=4&&y>=1&&y<=4&&sita_index>=1&&sita_index<=12&&scale_index>=1&&scale_index<=5)
                A(x,y,sita_index,scale_index)=A(x,y,sita_index,scale_index)+1;
                %存儲每個點各種旋轉角度與縮放系數下的locations 、旋轉角度sita、縮放尺度scale
                Apoint(i,sita_index,scale_index,:)=[x y sita_index scale_index];
            end
        end
    end
end

%搜尋累加表中的最大值，並認為其對應參數是正確的參數
maxA=0;
for x=1:4
    for y=1:4
        tmpA=reshape(A(x,y,:,:),12,5);
        tmp=max(max(tmpA));
        if(tmp>maxA)
            maxA=tmp;
            locate=[x y];
            [sita ,scale]=find(tmpA==tmp);
            sita_scale=[sita(1) scale(1)];
        end
    end
end

%剔除不屬於累加表中最大值所在bin的點
iner_point=1;
%For each edge point (x, y)
for i=1:size(matchedPoints1,1)
    for sita_index=1:12
        for scale_index=1:5
            x=Apoint(i,sita_index,scale_index,1);
            y=Apoint(i,sita_index,scale_index,2);
            sita_tmp=Apoint(i,sita_index,scale_index,3);
            scale_tmp=Apoint(i,sita_index,scale_index,4);
            if(x==locate(1)&&y==locate(2)&&sita_tmp==sita_scale(1)&&scale_tmp==sita_scale(2))
                point1(iner_point,:)=matchedPoints1(i,:);
                point2(iner_point,:)=matchedPoints2(i,:);
                iner_point=iner_point+1;
            end 
        end
    end
end

%用最小平方法求得仿射變換的參數
%The least-squares solution for the parameters x can be determined by solving the corresponding
%normal equations,
for i = 1: size(point2,1)
   AA(2*(i-1)+1,:)=[point2(i,2) point2(i,1) 0 0 1 0];
   AA(2*(i-1)+2,:)=[0 0 point2(i,2) point2(i,1) 0 1];
   b(2*(i-1)+1)=point1(i,2);
   b(2*(i-1)+2)=point1(i,1);
end
tmp=AA'*AA;
tmp=tmp^-1;
tmp=tmp*AA';
affine=tmp*b';
m1=affine(1);
m2=affine(2);
m3=affine(3);
m4=affine(4);
tx=affine(5);
ty=affine(6);

%用放射變換的結果求出場景中與模板中坐標（x,y）對應的坐標（u，v）
x=120;y=200;
u=m1*x+m2*y+tx;
v=m3*x+m4*y+ty;

%繪制出剔除奇異點之后的結果
% Show a figure with lines joining the accepted matches.
figure('name','result','Position', [100 100 size(im3,2) size(im3,1)]);
colormap('gray');
imagesc(im3);title('Generalised Hough transform');
hold on;
cols1 = size(im1,2);
for i = 1: size(point2,1)
    line([point1(i,2) point2(i,2)+cols1], ...
         [point1(i,1) point2(i,1)], 'Color', 'c');
end
%繪制模板與場景的對應點並連線
plot(x+cols1,y,'r*');
plot(u,v, 'r*');
line([x+cols1 u], ...
    [y v], 'Color', 'r','LineWidth',3);
hold off;
num = size(point1,1);
fprintf('Found %d matches.\n', num);

運行效果如下：

剔除奇異點之前：

剔除奇異點之后：

可以看到，仍有少許的錯誤匹配點，David G. Lowe提出可以用剛剛得到的變換矩陣繼續篩選，在此不再給出。

本文的源圖片可以在之前的博客中找到，之前版本的源碼包可在我的資源中下載。

轉載請注明出處：http://blog.csdn.net/u010278305