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

時間 2019-12-15

標籤 matlab 圖像處理學習筆記廣義變換篩選 sift 特徵欄目 MATLAB 简体版

原文原文鏈接

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

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

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

一、投票機制的引入。less

二、參數空間的轉換。ide

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

本文的主要參考資料以下：idea

%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-1spa

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

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

那麼這個轉換矩陣的參數如何肯定呢？這就要經過投票機制，若是一個參數得到的投票最多，那麼就認爲這個參數是正確的。

通常狀況下，要進行Generalised Hough transform，先要選取一個參考點，而後將普通的座標轉換爲(r，, ɸ)這種形式，r爲邊界點到選取參考點的長度，ɸ爲參考點到邊界點的連線與邊界點切線的夾角。以下圖所示：

以後創建一個R-table，列出邊界點上全部點的角度與其長度的對應關係，這樣便能完整的表述咱們的模板物體。R-table的創建過程以下：

R-table以下：

在以後的步驟中，咱們能夠經過角度來尋找匹配，再用計算出的值填充累加單元。但因爲咱們已經尋得各特徵點的匹配關係，在實施的過程當中，咱們跳過這一步。直接進入下一步。

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

一、選取模板的（0，0）爲參考點，表明模板的位置。模板的每個特徵點與參考點相減。（因爲參考點選爲（0，0），故各特徵點不變。

二、初始化累加表.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.

三、遍歷場景中的每一點，記x',y'爲模板中特徵點的座標。則能夠用下式推出場景中模板的座標，並對累加單元進行累加。

四、找到全部累加單元的最大值，並找到其索引，則其對應的（xc，yc）爲模板中（0，0）點在場景中對應的座標。

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

六、存儲剩下的點。

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