純MATLAB版本 SIFT代碼

先貼幾個連接：

http://blog.csdn.net/abcjennifer/article/details/7639681  Rachel-Zhang的

http://blog.csdn.net/manji_lee/article/details/8922474   

David G. Lowe的兩篇標誌性的文章分別叫 Object recognition from local scale-invariant features 和

Distinctive Image Features from Scale-Invariant Keypoints

其實我就是想在博客裏保存下網上下載的純MATLAB版本的sift代碼

 

%% 功能：提取灰度圖像的尺度不變特徵（SIFT特徵）
% 輸入:
% im - 灰度圖像，該圖像的灰度值在0到1之間（注意：應首先對輸入圖像的灰度值進行歸一化處理）
% octaves - 金字塔的組數：octaves (默認值爲4). 
% intervals - 該輸入參數決定每組金字塔的層數(默認值爲 2).
% object_mask - 肯定圖像中尺度不變特徵點的搜索區域，若是沒有特別指出，則算法將搜索整個圖像
% contrast_threshold - 對比度閾值(默認值爲0.03).
% curvature_threshold - 曲率閾值（默認值爲10.0）.
% interactive - 函數運行顯示標誌，將其設定爲1,則顯示算法運行時間和過程的相關信息；
% 若是將其設定爲2，則僅顯示最終運行結果(default = 1).
% 
% 輸出:
% pos - Nx2 矩陣，每一行包括尺度不變特徵點的座標(x,y) 
% scale - Nx3 矩陣，每一行包括尺度不變特徵點的尺度信息
% (第一列是尺度不變特徵點所在的組，
% 第二列是其所在的層, 
% 第三列是尺度不變特徵點的sigma).
% orient - Nx1 向量，每一個元素是特徵點的主方向，其範圍在 [-pi,pi)之間.
% desc - N x 128 矩陣，每一行包含特徵點的特徵向量.

%% 步驟0：載入圖像後歸一化，並擴充成 2xN-1 大小
% 用於完全搞懂SIFT的MATLAB版本代碼的一個測試文件
clear all;clc;close all;
% im = rgb2gray(imread('nest.png'));
im = imread('nest.png');
im = scjMatMap(im,0,1);

antialias_sigma = 0.5;
signal = im;

[X Y] = meshgrid( 1:0.5:size(signal,2), 1:0.5:size(signal,1) );
signal = interp2( signal, X, Y, '*linear' ); % im被擴充成  2 x N - 1  了

subsample = [0.5]; %  降採樣率
%% 步驟1：生成高斯和差分高斯(DOG)金字塔
%  這兩個金字塔的數據分別存儲在名爲
%  gauss_pyr{orient,interval}和DOG_pyr{orient,interval}
%  的元胞數組中。高斯金字塔含有s+3層，差分高斯金字塔含有s+2層。

preblur_sigma = sqrt(sqrt(2)^2 - (2*antialias_sigma)^2);

g = gaussian_filter( preblur_sigma );

gauss_pyr{1,1} = conv2( g, g, signal, 'same' );
clear signal;  %  完成了初始化工做

% 第一組第一層的sigma
initial_sigma = sqrt( (2*antialias_sigma)^2 + preblur_sigma^2 );

octaves = 4;                      % 有4層金字塔   
intervals = 2;                    % 每一層有2張圖片    s+3
% 之因此+3，是由於：生成DOG時會減小一個圖片，在計算極值的時候
% 第一個和最後一個圖片不能使用，因此要+3.
% 爲了保證差分高斯金字塔DOG在每一層至少有一張圖片可用來尋找極值點
% 必須讓高斯金字塔有+3張圖片
% 最後會有octaves層，每層intervals+3張圖片

% 記錄每一層和每個尺度的sigma
absolute_sigma = zeros(octaves,intervals+3);
absolute_sigma(1,1) = initial_sigma * subsample(1);

% 記錄產生金字塔的濾波器的尺寸和標準差
filter_size = zeros(octaves,intervals+3);
filter_sigma = zeros(octaves,intervals+3);

% 生成高斯和差分高斯金字塔
for octave = 1:octaves
   sigma = initial_sigma;
   g = gaussian_filter(sigma);
   filter_size(octave,1) = length(g);
   filter_sigma(octave,1) = sigma;
   DOG_pyr{octave} = ...
       zeros(size(gauss_pyr{octave,1},1),size(gauss_pyr{octave,1},2),intervals+2);
   for interval = 2:(intervals+3)     
      sigma_f = sqrt(2^(2/intervals) - 1)*sigma;
      g = gaussian_filter( sigma_f );
      sigma = (2^(1/intervals))*sigma;
      
      % 記錄sigma的值
      absolute_sigma(octave,interval) = sigma * subsample(octave);
      
      % 記錄濾波器的尺寸和標準差
      filter_size(octave,interval) = length(g);
      filter_sigma(octave,interval) = sigma;
      
      gauss_pyr{octave,interval} = conv2(g,g,gauss_pyr{octave,interval-1}, 'same' );    
      DOG_pyr{octave}(:,:,interval-1) = ...
          gauss_pyr{octave,interval} - gauss_pyr{octave,interval-1};                 
   end      
   if octave < octaves  
      sz = size(gauss_pyr{octave,intervals+1});
      [X Y] = meshgrid( 1:2:sz(2), 1:2:sz(1) );
      gauss_pyr{octave+1,1} = interp2(gauss_pyr{octave,intervals+1},X,Y,'*nearest'); 
      absolute_sigma(octave+1,1) = absolute_sigma(octave,intervals+1);
      subsample = [subsample subsample(end)*2];
   end      
end
%% 步驟2：查找差分高斯金字塔中的局部極值，並經過曲率和照度進行檢驗
contrast_threshold = 0.02;
curvature_threshold = 10.0;

curvature_threshold = ((curvature_threshold + 1)^2)/curvature_threshold;

% 二階微分核
xx = [ 1 -2  1 ];
yy = xx';
xy = [1 0 -1; 
      0 0  0; 
      -1 0 1]/4;

raw_keypoints = [];
contrast_keypoints = [];
curve_keypoints = [];

object_mask = ones(size(im));

% 在高斯金字塔中查找局部極值
loc = cell(size(DOG_pyr)); 
for octave = 1:octaves
   for interval = 2:(intervals+1)
      keypoint_count = 0;
      contrast_mask = abs(DOG_pyr{octave}(:,:,interval)) >= contrast_threshold;
      loc{octave,interval} = zeros(size(DOG_pyr{octave}(:,:,interval)));
        
      edge = ceil(filter_size(octave,interval)/2);
     
      for y=(1+edge):(size(DOG_pyr{octave}(:,:,interval),1)-edge)         
         for x=(1+edge):(size(DOG_pyr{octave}(:,:,interval),2)-edge)
            if object_mask(round(y*subsample(octave)),round(x*subsample(octave))) == 1    
                interactive = 1;
                if( (interactive >= 2) | (contrast_mask(y,x) == 1) )  % 注意！            
                  % 經過空間核尺度檢測最大值和最小值
                  tmp = DOG_pyr{octave}((y-1):(y+1),(x-1):(x+1),(interval-1):(interval+1));  
                  pt_val = tmp(2,2,2);
                  if( (pt_val == min(tmp(:))) | (pt_val == max(tmp(:))) )
                     % 存儲對灰度大於對比度閾值的點的座標
                     raw_keypoints = [raw_keypoints; x*subsample(octave) y*subsample(octave)];
                     if abs(DOG_pyr{octave}(y,x,interval)) >= contrast_threshold
                        contrast_keypoints = [contrast_keypoints; raw_keypoints(end,:)];     
                        % 計算局部極值的Hessian矩陣
                        Dxx = sum(DOG_pyr{octave}(y,x-1:x+1,interval) .* xx);
                        Dyy = sum(DOG_pyr{octave}(y-1:y+1,x,interval) .* yy);
                        Dxy = sum(sum(DOG_pyr{octave}(y-1:y+1,x-1:x+1,interval) .* xy));
                        % 計算Hessian矩陣的直跡和行列式.
                        Tr_H = Dxx + Dyy;
                        Det_H = Dxx*Dyy - Dxy^2;
                        % 計算主曲率.
                        curvature_ratio = (Tr_H^2)/Det_H;
                        if ((Det_H >= 0) & (curvature_ratio < curvature_threshold))
                           % 存儲主曲率小於閾值的的極值點的座標（非邊緣點）
                           curve_keypoints = [curve_keypoints; raw_keypoints(end,:)];
                           % 將該點的位置的座標設爲1，並計算點的數量.
                           loc{octave,interval}(y,x) = 1;
                           keypoint_count = keypoint_count + 1;
                        end
                     end                  
                  end
               end               
            end
         end         
      end

         fprintf( 2, '%d keypoints found on interval %d\n', keypoint_count, interval );

   end
end
clear raw_keypoints contrast_keypoints curve_keypoints;
%% 步驟3：計算特徵點的主方向. 
% 在特徵點的一個區域內計算其梯度直方圖
g = gaussian_filter( 1.5 * absolute_sigma(1,intervals+3) / subsample(1) );
zero_pad = ceil( length(g) / 2 );

% 計算高斯金字塔圖像的梯度方向和幅值
mag_thresh = zeros(size(gauss_pyr));
mag_pyr = cell(size(gauss_pyr));
grad_pyr = cell(size(gauss_pyr));

for octave = 1:octaves
   for interval = 2:(intervals+1)      
      
      % 計算x,y的差分
      diff_x = 0.5*(gauss_pyr{octave,interval}(2:(end-1),3:(end))-gauss_pyr{octave,interval}(2:(end-1),1:(end-2)));
      diff_y = 0.5*(gauss_pyr{octave,interval}(3:(end),2:(end-1))-gauss_pyr{octave,interval}(1:(end-2),2:(end-1)));
      
      % 計算梯度幅值
      mag = zeros(size(gauss_pyr{octave,interval}));      
      mag(2:(end-1),2:(end-1)) = sqrt( diff_x .^ 2 + diff_y .^ 2 );
      
      % 存儲高斯金字塔梯度幅值
      mag_pyr{octave,interval} = zeros(size(mag)+2*zero_pad);
      mag_pyr{octave,interval}((zero_pad+1):(end-zero_pad),(zero_pad+1):(end-zero_pad)) = mag;      
      
      % 計算梯度主方向
      grad = zeros(size(gauss_pyr{octave,interval}));
      grad(2:(end-1),2:(end-1)) = atan2( diff_y, diff_x );
      grad(find(grad == pi)) = -pi;
      
      % 存儲高斯金字塔梯度主方向
      grad_pyr{octave,interval} = zeros(size(grad)+2*zero_pad);
      grad_pyr{octave,interval}((zero_pad+1):(end-zero_pad),(zero_pad+1):(end-zero_pad)) = grad;
   end
end

clear mag grad
%% 步驟4：肯定特徵點的主方向
% 方法：經過尋找每一個關鍵點的子區域內梯度直方圖的峯值（注：每一個關鍵點的主方向能夠有不止一個）

% 將灰度直方圖分爲36等分，每隔10度一份
num_bins = 36;
hist_step = 2*pi/num_bins;  % 步進值
hist_orient = [-pi:hist_step:(pi-hist_step)];  % 步進方向

% 初始化關鍵點的位置、方向和尺度信息
pos = [];
orient = [];
scale = [];

% 給關鍵點肯定主方向
for octave = 1:octaves
   for interval = 2:(intervals + 1)           
      keypoint_count = 0;
      
      % 構造高斯加權掩模
      g = gaussian_filter( 1.5 * absolute_sigma(octave,interval)/subsample(octave) );
      hf_sz = floor(length(g)/2);
      g = g'*g;      
      
      loc_pad = zeros(size(loc{octave,interval})+2*zero_pad);
      loc_pad((zero_pad+1):(end-zero_pad),(zero_pad+1):(end-zero_pad)) = loc{octave,interval};
      
      [iy ix]=find(loc_pad==1);
      for k = 1:length(iy)

         x = ix(k);
         y = iy(k);
         wght = g.*mag_pyr{octave,interval}((y-hf_sz):(y+hf_sz),(x-hf_sz):(x+hf_sz));
         grad_window = grad_pyr{octave,interval}((y-hf_sz):(y+hf_sz),(x-hf_sz):(x+hf_sz));
         orient_hist=zeros(length(hist_orient),1);
         for bin=1:length(hist_orient)     
            diff = mod( grad_window - hist_orient(bin) + pi, 2*pi ) - pi;
            orient_hist(bin) = ...
                orient_hist(bin)+sum(sum(wght.*max(1 - abs(diff)/hist_step,0)));
         end
         
         % 運用非極大抑制法查找主方向直方圖的峯值
         peaks = orient_hist;        
         rot_right = [ peaks(end); peaks(1:end-1) ];
         rot_left = [ peaks(2:end); peaks(1) ];         
         peaks( find(peaks < rot_right) ) = 0;
         peaks( find(peaks < rot_left) ) = 0;
         
         % 提取最大峯值的值和其索引位置 
         [max_peak_val ipeak] = max(peaks);
         
         % 將大於等於最大峯值80% 的直方圖的也肯定爲特徵點的主方向
         peak_val = max_peak_val;
         while( peak_val > 0.8*max_peak_val )
	
            % 最高峯值最近的三個柱值經過拋物線插值精確獲得
            A = [];
            b = [];
            for j = -1:1
               A = [A; (hist_orient(ipeak)+hist_step*j).^2 (hist_orient(ipeak)+hist_step*j) 1];
	            bin = mod( ipeak + j + num_bins - 1, num_bins ) + 1;
               b = [b; orient_hist(bin)];
            end
            c = pinv(A)*b;
            max_orient = -c(2)/(2*c(1));
            while( max_orient < -pi )
               max_orient = max_orient + 2*pi;
            end
            while( max_orient >= pi )
               max_orient = max_orient - 2*pi;
            end            
            
            % 存儲關鍵點的位置、主方向和尺度信息
            pos = [pos; [(x-zero_pad) (y-zero_pad)]*subsample(octave) ];
            orient = [orient; max_orient];
            scale = [scale; octave interval absolute_sigma(octave,interval)];
            keypoint_count = keypoint_count + 1;   

            peaks(ipeak) = 0;
            [peak_val ipeak] = max(peaks);
         end         
      end
         fprintf( 2, '(%d keypoints)\n', keypoint_count );          
   end
end
clear loc loc_pad
%% 步驟5：特徵向量生成
orient_bin_spacing = pi/4;
orient_angles = [-pi:orient_bin_spacing:(pi-orient_bin_spacing)];

grid_spacing = 4;
[x_coords y_coords] = meshgrid( [-6:grid_spacing:6] );
feat_grid = [x_coords(:) y_coords(:)]';
[x_coords y_coords] = meshgrid( [-(2*grid_spacing-0.5):(2*grid_spacing-0.5)] );
feat_samples = [x_coords(:) y_coords(:)]';
feat_window = 2*grid_spacing;

desc = [];

for k = 1:size(pos,1)
   x = pos(k,1)/subsample(scale(k,1));
   y = pos(k,2)/subsample(scale(k,1));   
   
   % 將座標軸旋轉爲關鍵點的方向，以確保旋轉不變性
   M = [cos(orient(k)) -sin(orient(k)); sin(orient(k)) cos(orient(k))];
   feat_rot_grid = M*feat_grid + repmat([x; y],1,size(feat_grid,2));
   feat_rot_samples = M*feat_samples + repmat([x; y],1,size(feat_samples,2));
   
   % 初始化特徵向量.
   feat_desc = zeros(1,128);
   
   for s = 1:size(feat_rot_samples,2)
      x_sample = feat_rot_samples(1,s);
      y_sample = feat_rot_samples(2,s);
      
      % 在採樣位置進行梯度插值
      [X Y] = meshgrid( (x_sample-1):(x_sample+1), (y_sample-1):(y_sample+1) );
      G = interp2( gauss_pyr{scale(k,1),scale(k,2)}, X, Y, '*linear' );  % 耗時太長
      G(find(isnan(G))) = 0;
      diff_x = 0.5*(G(2,3) - G(2,1));
      diff_y = 0.5*(G(3,2) - G(1,2));
      mag_sample = sqrt( diff_x^2 + diff_y^2 );
      grad_sample = atan2( diff_y, diff_x );
      if grad_sample == pi
         grad_sample = -pi;
      end      
      
      % 計算x、y方向上的權重
      x_wght = max(1 - (abs(feat_rot_grid(1,:) - x_sample)/grid_spacing), 0);
      y_wght = max(1 - (abs(feat_rot_grid(2,:) - y_sample)/grid_spacing), 0); 
      pos_wght = reshape(repmat(x_wght.*y_wght,8,1),1,128);
      
      diff = mod( grad_sample - orient(k) - orient_angles + pi, 2*pi ) - pi;
      orient_wght = max(1 - abs(diff)/orient_bin_spacing,0);
      orient_wght = repmat(orient_wght,1,16);         
      
      % 計算高斯權重
      g = exp(-((x_sample-x)^2+(y_sample-y)^2)/(2*feat_window^2))/(2*pi*feat_window^2);
      
      feat_desc = feat_desc + pos_wght.*orient_wght*g*mag_sample;
   end
   
   % 將特徵向量的長度歸一化，則能夠進一步去除光照變化的影響.
   feat_desc = feat_desc / norm(feat_desc);
   
   feat_desc( find(feat_desc > 0.2) ) = 0.2;
   feat_desc = feat_desc / norm(feat_desc);
   
   % 存儲特徵向量.
   desc = [desc; feat_desc];
   tmp = mod(k,25);
   if ( tmp == 0 )
      fprintf( 2, '.' );
   end
end

% 調整採樣誤差
sample_offset = -(subsample - 1);
for k = 1:size(pos,1)
   pos(k,:) = pos(k,:) + sample_offset(scale(k,1));
end

if size(pos,1) > 0
	scale = scale(:,3);
end

hh = display_keypoints( pos, scale, orient);

調用了

function g = gaussian_filter(sigma)
 
% 功能：生成一維高斯濾波器
% 輸入：
% sigma – 高斯濾波器的標準差
% 輸出：
% g – 高斯濾波器

sample = 7.0/2.0;
n = 2*round(sample * sigma)+1; 
x = 1:n;
x = x - ceil(n/2);
 
g = exp(-(x.^2)/(2*sigma^2))/(sigma*sqrt(2*pi));

 

原做者是Xiaochuan ZHAO 郵箱：zhaoxch@bit.edu.cn

 

我去除了他的代碼中一些可有可無的部分。