www.pudn.com > MatlabSIFTcode_1.rar > find_features.m

 
%///////////////////////////////////////////////////////////////////////////////////////////// 
% 
% find_features - scale space feature detector based upon difference of gaussian filters. 
%                 selects features based upon their maximum response in scale space 
% 
% Usage:  maxima = find_features(pyr, img, thresh, radius, radius2, min_sep, edgeratio, disp_flag, img_flag) 
% 
% Parameters:   
%            pyr :      cell array of filtered image pyramid (built with build_pyramid) 
%            img :      original image (only used for visualization) 
%            thresh :   threshold value for maxima search (minimum filter response considered) 
%            radius :   radius for maxima comparison within current scale 
%            radius2:   radius for maxima comparison between neighboring scales 
%            disp_flag: 1- display each scale level on separate figure.  0 - display nothing 
%            img_flag: 1 - display filter responses. 0 - display original images. 
% 
% Returns: 
% 
%            maxima  - cell array of nX2 matrices of row,column coordinates of selected points on each scale level 
% 
% Author:  
% Scott Ettinger 
% scott.m.ettinger@intel.com 
% 
% May 2002 
%///////////////////////////////////////////////////////////////////////////////////////////// 
 
function maxima = find_features(pyr, img, scl, thresh, radius, radius2, disp_flag, img_flag) 
  %                pts = find_features(pyr,img,scl,thresh,radius,radius2,disp_flag,1); 
 
 
levels = size(pyr); 
levels = levels(2); 
 
mcolor = [ 0 1 0;   %color array for display of features at different scales 
           0 1 0; 
           1 0 0; 
           .2 .5 0; 
           0 0 1; 
           1 0 1; 
           0 1 1; 
           1 .5 0 
           .5 1 0 
          0 1 .5 
           .5 1 .5]; 
 
 
[himg,wimg] = size(img);                                            %get size of images 
[h,w] = size(pyr{2}); 
 
     
for i=2:levels-1 
 
    [h,w] = size(pyr{i}); 
    [h2,w2] = size(pyr{i+1}); 
     
    %find maxima 
    mx = find_extrema(pyr{i},thresh,radius);                        %find maxima at current scale level 
    mx2 = round((mx-1)/scl) + 1;                                    %find coords in level above         
    mx_above = neighbor_max(pyr{i},pyr{i+1},mx,mx2,radius2);        %do neighbor comparison in scale space above 
 
    if i>1 
        mx2 = round((mx-1)*scl) + 1;                                %find coords in level below      
        mx_below = neighbor_max(pyr{i},pyr{i-1},mx,mx2,radius2);    %do comparison in scale below 
     
        maxima{i} = plist(mx, mx_below & mx_above);                 %get coord list for retained maxima and minima 
    else 
        maxima{i} = plist(mx, mx_above); 
    end 
        
    %find minima 
    %if i==11, 
    %    keyboard; 
    %end; 
     
    mx = find_extrema(-pyr{i},thresh,radius);                       %find minima at current scale level 
    mx2 = round((mx-1)/scl) + 1;                                    %find coords in level above         
    mx_above = neighbor_max(-pyr{i},-pyr{i+1},mx,mx2,radius2);      %do neighbor comparison in scale space above 
    if i>1 
        mx2 = round((mx-1)*scl) + 1;                                %find coords in level below      
        mx_below = neighbor_max(-pyr{i},-pyr{i-1},mx,mx2,radius2);  %do comparison in scale below 
     
        mxtemp = plist(mx, mx_below & mx_above);                    %get coord list for retained maxima and minima 
    else 
        mxtemp = plist(mx, mx_above); 
    end 
   
 
    maxima{i} = [maxima{i}; mxtemp];                                %combine maxima and minima into list for return 
     
     
    %display results if desired 
    if disp_flag > 0                                                 
        figure 
        if img_flag == 0 
            tmp=resample_bilinear(img,himg/h); 
            imagesc(tmp); 
            colormap gray; 
            show_plist(maxima{i},mcolor(mod(i-1,7)+1,:),'+'); 
        else 
            imagesc(pyr{i}); 
            colormap gray; 
            show_plist(maxima{i},mcolor(mod(i-1,7)+1,:),'+'); 
        end 
    end     
end 
     
%////////////////////////////////////////////////////////////////////////////////////////////// 
% 
%   Compare a vector of pixels with its neighbors in another scale  
% 
%////////////////////////////////////////////////////////////////////////////////////////////// 
 
function v = neighbor_max(img1,img2,i,i2,radius)                    % i and i2 are column vectors of r,c coords 
     
    if (size(i2,1))==0 | size(img2,1)<11 | size(img2,2)<11 
        v=zeros(length(i),1); 
    else 
         
        [h,w] = size(img1); 
        [h2,w2] = size(img2); 
     
        [y,x]=meshgrid(-20:20,-20:20);                              %create set of offsets within radius  
        z = (x.^2+y.^2)<=radius^2; 
        [y,x]=find(z); 
        x=x-21; y=y-21; 
     
        radius=ceil(radius); 
     
        bound = ones(size(i2,1),2)*[h2-radius 0;0 w2-radius];        %create boundary listing 
        i2 = i2 - ((i2 > bound).*(i2-bound+1));                      %test bounds to make all points within image 
        i2 = i2 + ((i2 < radius+1).*(radius-i2+1)); 
         
         
        i2 = vec(i2,h2);                                             %create indices from x,y coords 
        i = vec(i,h); 
     
        p = img1(i); 
        res = ones(length(i),1); 
     
        for j=1:length(x)                                            %check against all points within radius 
            itest = i2 + x(j) + h2*y(j); 
            p2 = img2(itest); 
            res = res & (p>=p2); 
        end 
 
         
        v = res;                                                     %store results in binary vector 
    end 
     
%////////////////////////////////////////////////////////////////////////////////////////////// 
 
function v = vec(points,h) 
 
    y = points(:,1); 
    x = points(:,2); 
     
    v = y + (x-1)*h;  %create index vectors 
     
%////////////////////////////////////////////////////////////////////////////////////////////// 
 
function p = plist(points, flags) 
 
    p = points(find(flags),:);