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