sift--算法学习 2012 年5月草稿

这个是几个月前开始学习的东西，好久没有碰过了~~

现在来 自惭形秽一下。

这是队友给我的matlab的资料，网上流传的学习中有两份，一份是vedaldi的，一份是lowe的。

这里这份是仿v的，全部用matlab来实现，速度肯定比不上v的matlab与C的混编版，更别说纯c版

因为matlab的循环是非常费时间的，我有一次循环计算量大了差点死机。

废话不多说，我从do_demo1开始吧

1.读入一幅图

% Add subfolder path.
main; 
img1_dir = 'demo-data\';
img1_file = 'beaver11.bmp';

I1=imreadbw([img1_dir img1_file]) ; 

2.调整大小，我觉得这个是从运算量的角度来说
I1_rgb = imread([img1_dir img1_file]) ; 
I1=imresize(I1, [240 320]);

I1_rgb =imresize(I1_rgb, [240 320]);
I1=I1-min(I1(:)) ;                                     %min(I1) <0 我感觉是去除负值，我调试的时候发现有负值，为什么不知道
I1=I1/max(I1(:)) ;


%fprintf('CS5240 -- SIFT: Match image: Computing frames and descriptors.\n') ;

 %FRAMES(1:2,k)  center (X,Y) of the frame k,  k帧的中心元素坐标
%     FRAMES(3,k)    scale SIGMA of the frame k,    不知道这个中文怎么翻译，尺度空间的规模
%     FRAMES(4,k)    orientation THETA of the frame k.  方向

 DESCR is a DxK matrix stores one descriptor per columm (usually
%   D=128).  描述子

[frames1,descr1,gss1,dogss1] = do_sift( I1, 'Verbosity', 1, 'NumOctaves', 4, 'Threshold',  0.1/3/2 ) ; %0.04/3/2

figure(11) ; clf ; plotss(dogss1) ; colormap gray ;
drawnow ;


figure(2) ; clf ;
imshow(I1_rgb) ; axis image ;


hold on ;
h=plotsiftframe( frames1 ) ; set(h,'LineWidth',1,'Color','g') ;

function [frames,descriptors,scalespace,difofg]=sift(I,varargin)
%% file:       sift.m
% author:      Noemie Phulpin
% description: SIFT algorithm


warning off all;


[M,N,C] = size(I) ;

3.生成高斯尺度空间

O可以理解成组。S是不同的尺度。

% Lowe's choices
S=3 ;
omin= 0 ;
%O=floor(log2(min(M,N)))-omin-4 ; % Up to 16x16 images
O = 4;
这部分是算法里面有的
sigma0=1.6*2^(1/S) ;   
sigman=0.5 ;

这个门限值我觉得有点带经验色彩
thresh = 0.2 / S / 2 ; % 0.04 / S / 2 ;
r = 18 ;


NBP = 4 ;
NBO = 8 ;
magnif = 3.0 ;


% Parese input
compute_descriptor = 0 ;
discard_boundary_points = 1 ;
verb = 0 ;


% Arguments sanity check
if C > 1
  error('I should be a grayscale image') ;
end


frames = [] ;
descriptors = [] ;


% 
% Do the job
%
fprintf('---------------------------- CS5240 SIFT: Extract SIFT features from an image ------------------------------\n') ; tic ; 


fprintf('CS5240 -- SIFT: constructing scale space with DoG ...\n') ; tic ; 

这里是有关尺度空间的资料。


scalespace = do_gaussian(I,sigman,O,S,omin,-1,S+1,sigma0) ;

function L = do_gaussian(I,sigman,O,S,omin,smin,smax,sigma0)


%% file:        do_gaussian.m
 % author:      Noemie Phulpin
 % description: gaussian scale space of image I
 %%


if(nargin<7) 
   sigma0=1.6*k;
end

 omin是组的初始值，如果小于0就放大，大于0就缩小。
if omin<0
   for o=1:-omin
I=doubleSize(I);
   end
elseif omin>0
   for o=1:-omin
I=halveSize(I);
   end
end
[M,N] = size(I);                      %size of image


k = 2^(1/S);                          %scale space multiplicative step k
sigma0=1.6*k;                         %definition by Lowe
dsigma0 = sigma0*sqrt(1-1/k^2);       %scale step factor
sigmaN=0.5;                           %nominal smoothing of the image
so=-smin+1;                           %index offset  这是只是用来说明有多少层。

这里是生成一个结构体，我是怎么理解的。
%scale space structure
L.O = O;
L.S = S;
L.sigma0 = sigma0;
L.omin = omin;
L.smin = smin;
L.smax = smax;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%First Octave
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%initilize the octave with S sub-levels
L.octave{1} = zeros(M,N,smax-smin+1);  %我不记得这里我改过没有，好像就是加上so


%initilize the first sub-level
sig=sqrt( (sigma0*k^smin)^2 - (sigmaN/2^omin)^2 );
%b=smooth2(I,sig) ;
%[N1,M1]=size(b)
%b(1:4,1:4)
%c=imsmooth(I,sig) ;
%[N2,M2]=size(c)
%c(1:4,1:4)
L.octave{1}(:,:,1) = smooth(I,sig);
smooth也是作者自己写的，根据sig的值建立滤波算子，进行滤波
function J = smooth(I,s)
%% file:        smooth.m
% author:      Noemie Phulpin
% description: smoothing image
%filter 
h=fspecial('gaussian',ceil(4*s),s);
%convolution
J=imfilter(I,h);
return;

%other sub-levels 
for s=smin+1:smax
    dsigma = k^s * dsigma0;
    L.octave{1}(:,:,s+so) = smooth( squeeze(L.octave{1}(:,:,s-1+so)) ,dsigma);
end
squeeze这里可否去掉呢，我去掉了运行demo-1是没有问题的



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Folowing Octaves
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We need to initialize the first level of octave (o,smin) from
  % the closest possible level of the previous octave. A level (o,s)
  % in this octave corrsponds to the level (o-1,s+S) in the previous
  % octave. In particular, the level (o,smin) correspnds to
  % (o-1,smin+S). However (o-1,smin+S) might not be among the levels
  % (o-1,smin), ..., (o-1,smax) that we have previously computed.
  % The closest pick is
  %
  %                       /  smin+S    if smin+S <= smax
  % (o-1,sbest) , sbest = |
  %                       \  smax      if smin+S > smax
  %
  % The amount of extra smoothing we need to apply is then given by
  %
  %  ( sigma0 2^o 2^(smin/S) )^2 - ( sigma0 2^o 2^(sbest/S - 1) )^2
  %
  % As usual, we divide by 2^o to cancel out the effect of the
  % downsampling and we get
  %
  %  ( sigma 0 k^smin )^2 - ( sigma0 2^o k^(sbest - S) )^2
  %

%convert all octaves
for o=2:O
    
    sbest = min(smin+S,smax);
    TMP = halvesize( squeeze(L.octave{o-1}(:,:,sbest+so)) );
    sigma_next = sigma0*k^smin;
    sigma_prev = sigma0*k^(sbest-S);

这里如果smin+S《smax，则sbest=smin+S，那么sigma_prev=sigma_next

否则，说明smin+S》smax，则要调整TMp的值，不然就是上一组的最后一层直接下采样
    
    if (sigma_next>sigma_prev)
       sig=sqrt(sigma_next^2-sigma_prev^2);
       TMP= smooth( TMP,sig);
    end
    
    [M,N] = size(TMP);
    L.octave{o} = zeros(M,N,smax-smin+1); 
    L.octave{o}(:,:,1) = TMP;
    
    %other sub-levels 

这里我把v的注释加上了，其实就是提取了下公因式

 % Here we go from (omin,s-1) to (omin,s). The extra smoothing
  % standard deviation is
  %
  %  (sigma0 2^omin 2^(s/S) )^2 - (simga0 2^omin 2^(s/S-1/S) )^2
  %
  % Aftred dividing by 2^omin (to take into account the fact
  % that the image has been pre-scaled omin octaves), the
  % standard deviation of the smoothing kernel is
  %
  %   dsigma = sigma0 k^s sqrt(1-1/k^2)
  %

    for s=smin+1:smax
        dsigma = k^s * dsigma0;
        L.octave{o}(:,:,s+so) = smooth( squeeze(L.octave{o}(:,:,s-1+so)) ,dsigma);
    end


end







下面的这个图象缩小放大的方法是最简单的了。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Auxiliary functions


function J = halvesize(I)
J=I(1:2:end,1:2:end);


function J = doubleSize(I)
[M,N]=size(I) ;
J = zeros(2*M,2*N) ;
J(1:2:end,1:2:end) = I ;
J(2:2:end-1,2:2:end-1) = ...
0.25*I(1:end-1,1:end-1) + ...
0.25*I(2:end,1:end-1) + ...
0.25*I(1:end-1,2:end) + ...
0.25*I(2:end,2:end) ;
J(2:2:end-1,1:2:end) = ...
0.5*I(1:end-1,:) + ...
    0.5*I(2:end,:) ;
J(1:2:end,2:2:end-1) = ...
0.5*I(:,1:end-1) + ...
    0.5*I(:,2:end) ;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf('                Time for Gaussian scale space construction: %.3f s\n',toc) ; tic ; 

4.生成的尺度空间求差。


difofg = do_diffofg(scalespace) ;

function D = do_diffofg(L)


%% file:        do_diffofg.m
 % author:      Noemie Phulpin
 % description: substraction of consecutive levels of the scale space SS.
 %%




D.smin = L.smin;
D.smax = L.smax-1;
D.omin =L.omin;
D.O = L.O;
D.S = L.S;
D.sigma0 = L.sigma0;


for o=1:D.O
    
    [M,N,S] = size(L.octave{o});
    D.octave{o} = zeros(M,N,S-1);
    
    for s=1:S-1
        D.octave{o}(:,:,s) = L.octave{o}(:,:,s+1) -  L.octave{o}(:,:,s);   
    end;
    
end;

fprintf('                Time for Differential scale space construction: %.3f s\n',toc) ;


for o=1:scalespace.O
    
    
fprintf('CS5240 -- SIFT: computing octave %d\n', o-1+omin) ;  %相当于从omin开始
                tic ;

  % Local maxima of the DOG octave
    oframes1 = do_localmax(  difofg.octave{o}, 0.8*thresh, difofg.smin  ) ;
oframes = [oframes1 , do_localmax( - difofg.octave{o}, 0.8*thresh, difofg.smin)] ; 
    
    
    fprintf('CS5240 -- SIFT: initial keypoints # %d.  \n', ...
      size(oframes, 2)) ;
    fprintf('                Time (%.3f s)\n', ...
       toc) ;
    tic ;

    if size(oframes, 2) == 0
        continue;
    end
    
  % Remove points too close to the boundary
    rad = magnif * scalespace.sigma0 * 2.^(oframes(3,:)/scalespace.S) * NBP / 2 ;
    %magnif=3;
    %sigma0*2.^(s/S)*nbp/2
    sel=find(...
      oframes(1,:)-rad >= 1                     & ...
      oframes(1,:)+rad <= size(scalespace.octave{o},2) & ...      %列数
      oframes(2,:)-rad >= 1                     & ...
      oframes(2,:)+rad <= size(scalespace.octave{o},1)     ) ;     %行数
    oframes=oframes(:,sel) ;

fprintf('CS5240 -- SIFT: keypoints # %d after discarding from boundary\n', size(oframes,2)) ;
      tic ;

  % Refine the location, threshold strength and remove points on edges
    oframes = do_extrefine(...
  oframes, ...
  difofg.octave{o}, ...
  difofg.smin, ...
  thresh, ...
  r) ;
   
    fprintf('CS5240 -- SIFT: keypoints # %d after discarding from low constrast and edges\n',size(oframes,2)) ;
    fprintf('                Time (%.3f s)\n',  toc) ;
    tic ;
  
    fprintf('CS5240 -- SIFT: compute orientations of keypoints\n');


  % Compute the orientations
oframes = do_orientation(...
oframes, ...
scalespace.octave{o}, ...
scalespace.S, ...
scalespace.smin, ...
scalespace.sigma0 ) ;
fprintf('                time: (%.3f s)\n', toc);tic;

  % Store frames
x     = 2^(o-1+scalespace.omin) * oframes(1,:) ;
y     = 2^(o-1+scalespace.omin) * oframes(2,:) ;
sigma = 2^(o-1+scalespace.omin) * scalespace.sigma0 * 2.^(oframes(3,:)/scalespace.S) ;
frames = [frames, [x(:)' ; y(:)' ; sigma(:)' ; oframes(4,:)] ] ;


fprintf('CS5240 -- SIFT: keypoints # %d after orientation computation \n', size(frames,2)) ;
  % Descriptors

    fprintf('CS5240 -- SIFT: compute descriptors...\n') ;
    tic ;

sh = do_descriptor(scalespace.octave{o}, ...
                    oframes, ...
                    scalespace.sigma0, ...
                    scalespace.S, ...
                    scalespace.smin, ...
                    'Magnif', magnif, ...
                    'NumSpatialBins', NBP, ...
                    'NumOrientBins', NBO) ;
    
    descriptors = [descriptors, sh] ;
    
    fprintf('                time: (%.3f s)\n\n\n',toc) ; 
    
      
    
end 


fprintf('CS5240 -- SIFT: total keypoints: %d \n\n\n', size(frames,2)) ;