MATLAB实现图像的SIFT特征提取源码

function [frames,descriptors,scalespace,difofg]=do_sift(I,varargin)
 %% file:       sift.m
 % author:      Noemie Phulpin
 % description: SIFT algorithmwarning off all;
[M,N,C] = size(I) ;
% 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 应该是灰度图') ;
 endframes = [] ;
 descriptors = [] ;% 
 %   开始工作
 %
 fprintf('---------------------------- 开始 SIFT: 从图像中提取SIFT特征 ------------------------------\n') ; tic ; fprintf('SIFT: 用DoG构造尺度空间 ...\n') ; tic ; 
scalespace = do_gaussian(I,sigmaN,O,S,omin,-1,S+1,sigma0) ;
fprintf('                高斯尺度空间计时: (%.3f s)\n',toc) ; tic ; 
difofg = do_diffofg(scalespace) ;
fprintf('                构建相减尺度空间: (%.3f s)\n',toc) ;
for o=1:scalespace.O
     
     
     fprintf('CS5240 -- SIFT: 计算 “组”  %d\n', o-1+omin) ;
                 tic ;
     
   %  DOG octave 的局部极值检测
     oframes1 = do_localmax(  difofg.octave{o}, 0.8*thresh, difofg.smin  ) ;
     oframes2 = do_localmax( -difofg.octave{o}, 0.8*thresh, difofg.smin  ) ;
     oframes = [oframes1 ,oframes2 ] ; 
     
     
     fprintf('CS5240 -- SIFT: 初始化关键点 # %d.  \n', ...
       size(oframes, 2)) ;
     fprintf('                用时 (%.3f s)\n', ...
        toc) ;
     tic ;
     
     if size(oframes, 2) == 0
         continue;
     end
     
   % 移除靠近边界的关键点
     rad = magnif * scalespace.sigma0 * 2.^(oframes(3,:)/scalespace.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) ;%把不是靠近边界点的极值点重新放入oframes中
         
     fprintf('CS5240 -- SIFT: 移除靠近边界关键点后 # %d \n', size(oframes,2)) ;
       tic ;
         
   % 精简局部, 阈值强度 和移除边缘关键点
        oframes = do_extrefine(...
          oframes, ...
          difofg.octave{o}, ...
          difofg.smin, ...
          thresh, ...
          r) ;
    
        fprintf('CS5240 -- SIFT:  移除低对比度和边缘上关键点后 # %d \n',size(oframes,2)) ;
     fprintf('                Time (%.3f s)\n',  toc) ;
     tic ;
   
     fprintf('CS5240 -- SIFT: 计算特征点方向\n');  % 计算方向
     oframes = do_orientation(...
         oframes, ...
         scalespace.octave{o}, ...
         scalespace.S, ...
         scalespace.smin, ...
         scalespace.sigma0 ) ;
     fprintf('                用时: (%.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:计算方向后的特征点 # %d  \n', size(frames,2)) ;
   % Descriptors
     
     fprintf('CS5240 -- SIFT: 计算 描述子...\n') ;
     tic ;
         
     sh = do_descriptor(scalespace.octave{o}, ...
                     oframes, ...
                     scalespace.sigma0, ...
                     scalespace.S, ...
                     scalespace.smin, ...
                     'Magnif', magnif, ...
                     'NumSpatialBins', NBP, ...
                     'NumOrientBins', NBO) ;
     
     descriptors = [descriptors, sh] ;%每一组计算描述子向量后补充到descriptors数组中
     
     fprintf('                用时: (%.3f s)\n\n\n',toc) ; 
     
       
     
 end fprintf('CS5240 -- SIFT: 关键点总数: %d \n\n\n', size(frames,2)) ;