function [ pos, scale, orient, desc ] = SIFT( im, octaves, intervals, object_mask, contrast_threshold, curvature_threshold, interactive )
% [ pos, scale, orient, desc ] = SIFT( im, octaves, intervals, object_mask, contrast_threshold, curvature_threshold, interactive )
%
% Apply David Lowe's Scale Invariant Feature Transform (SIFT) algorithm
% to a grayscale image. This algorithm takes a grayscale image as input
% and returns a set of scale- and rotationally-invariant keypoints
% allong with their corresponding feature descriptors.
%
% This implementation is based on:
% [1] David G. Lowe, "Distinctive Image Features from Sacle-Invariant Keypoints",
% accepted for publicatoin in the International Journal of Computer
% Vision, 2004.
% [2] David G. Lowe, "Object Recognition from Local Scale-Invariant Features",
% Proc. of the International Conference on Computer Vision, Corfu,
% September 1999.
%
% Input:
% im - the input image, with pixel values normalize to lie betwen [0,1].%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%attention
%输入图像
% octaves - the number of octaves to search for keypoints (default = 4).
% 八度的个数
% intervals - the number of geometrically sampled intervals to divide
% each octave into when searching for keypoints (default = 2).
% 尺度采样数(组),比如在同一个ocave下三个scale中产生极值点,本尺度+上一尺度+下一尺度,产生一个采样数(组)
% object_mask - a binary mask specifying the location of the object in
% the image to search for keypoints on. If not specified, the whole
% image is searched.
%检测区域
% contrast_threshold - the threshold on the contrast of the DOG extrema
% before classifying them as keypoints (default = 0.03).
%对比度阈值--用于去除不明显的特征点
% curvature_threshold - the upper bound on the ratio between the principal
% curvatures of the DOG extrema before classifying it as a keypoint
% (default = 10.0).
%曲率阈值---用于去除边缘上的特征点
% interactive - >= 1 displays progress and timing information,
% >= 2 displays intermediate results of the algorihtm (default = 1).
%调试参数 1:显示进程和耗时。2:显示中间结果
% Output:
% pos - an Nx2 matrix containing the (x,y) coordinates of the keypoints
% stored in rows.
% 特征点的位置(x,y)序列
% scale - an Nx3 matrix with rows describing the scale of each keypoint (i.e.,
% first column specifies the octave, second column specifies the interval, and
% third column specifies sigma).
%特征点的尺度(octave,interval,sigma)序列
% orient - a Nx1 vector containing the orientations of the keypoints [-pi,pi).
%特征点的方向序列
% desc - an Nx128 matrix with rows containing the feature descriptors
% corresponding to the keypoints.
%SIFT描述子序列
% Thomas F. El-Maraghi
% May 2004
% assign default values to the input variables
if ~exist('octaves','var')
octaves = 4;
end
if ~exist('intervals','var')
intervals = 2;
end
if ~exist('object_mask','var')
object_mask = ones(size(im));
end
if size(object_mask) ~= size(im)
object_mask = ones(size(im));
end
if ~exist('contrast_threshold','var')
contrast_threshold = 0.02;%注意,这里没有用0.03,而是0.02
end
if ~exist('curvature_threshold','var')
curvature_threshold = 5.0;
end
if ~exist('interactive','var')
interactive = 1;
end
% Check that the image is normalized to [0,1]
if( (min(im(:)) < 0) | (max(im(:)) > 1) )
fprintf( 2, 'Warning: image not normalized to [0,1].\n' );%1 for standard output (the screen) or 2 for standard error
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%先进行标准差0.5的高斯滤波(实际上是相当于没有进行的),然后线性插值产生subsample=0.5的图像,就是把图像放大一倍
% Blur the image with a standard deviation of 0.5 to prevent aliasing
% and then upsample the image by a factor of 2 using linear interpolation.
% Lowe claims that this increases the number of stable keypoints by
% a factor of 4.
if interactive >= 1
fprintf( 2, 'Doubling image size for first octave...\n' );
end
tic;%开启一个秒表计数器
antialias_sigma = 0.5;
if antialias_sigma == 0
signal = im;
else
g = gaussian_filter( antialias_sigma );%产生一个1D高斯滤波器
if exist('corrsep') == 3 %3 if corrsep is a MEX-file on MATLAB's search path
signal = corrsep( g, g, im );
else
signal = conv2( g, g, im, 'same' );%MATLAB标准函数,convolves im first with the vector g along the rows and then with the vector g along the columns.
% 'same' - returns the central part of the convolution that is the same size as A.
end
end
if(signal==im)
fprintf(2,'signal==im\n');
else
fprintf(2,'signal!=im\n');
end
signal = im;%!!!!这里有问题:进行这次赋值后,相当于没有进行上面的高斯滤波
[X Y] = meshgrid( 1:0.5:size(signal,2), 1:0.5:size(signal,1) ); %SIZE(X,1) returns the number of rows.对于m*n维矩阵,插值产生的矩阵是(2m-1)*(2n-1)维的
signal = interp2( signal, X, Y, '*linear' );%使用时,必须在方法的名字前加上一个星号'*',如interp1(x, Y, xx ,'*cubic')
%函数原型是ZI = INTERP2(X,Y,Z,XI,YI)。ZI = INTERP2(Z,XI,YI) assumes X=1:N and Y=1:M where [M,N]=SIZE(Z).
subsample = [0.5]; % subsampling rate for doubled image is 1/2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%下一步产生高斯金字塔和DOG金字塔。如果采样数ineval=s,则需要产生s+3幅gaussion图像,s+2幅DOG图像
% The next step of the algorithm is to generate the gaussian and difference-of-
% gaussian (DOG) pyramids. These pyramids will be stored as two cell arrays(元胞数组),
% gauss_pyr{orient(octave),interval} and DOG_pyr{orient(octave),interval}, respectively. In order
% to detect keypoints on s intervals per octave, we must generate s+3 blurred
% images in the gaussian pyramid. This is becuase s+3 blurred images generates
% s+2 DOG images, and two images are needed (one at the highest and one lowest scales
% of the octave) for extrema detection.
% Generate the first image of the first octave of the gaussian pyramid
% by preblurring the doubled image with a gaussian with a standard deviation
% of 1.6. This choice for sigma is a trade off between repeatability and
% efficiency.
if interactive >= 1
fprintf( 2, 'Prebluring image...\n' );
end
preblur_sigma = sqrt(sqrt(2)^2 - (2*antialias_sigma)^2);% 1.0
if preblur_sigma == 0
gauss_pyr{1,1} = signal;
else
g = gaussian_filter( preblur_sigma );%产生一个1D高斯滤波器,sigma^2 = 1.0
if exist('corrsep') == 3
gauss_pyr{1,1} = corrsep( g, g, signal );
else
gauss_pyr{1,1} = conv2( g, g, signal, 'same' );
end
end
clear signal
pre_time = toc;
if interactive >= 1
fprintf( 2, 'Preprocessing time %.2f seconds.\n', pre_time );
end
% The initial blurring for the first image of the first octave of the pyramid.
initial_sigma = sqrt( (2*antialias_sigma)^2 + preblur_sigma^2 );%sqrt(2)
% Keep track of the absolute sigma for the octave and scale
absolute_sigma = zeros(octaves,intervals+3);
absolute_sigma(1,1) = preblur_sigma * subsample(1);
% Keep track of the filter sizes and standard deviations used to generate the pyramid
filter_size = zeros(octaves,intervals+3);%滤波器个数
filter_sigma = zeros(octaves,intervals+3);
% Generate the remaining levels of the geometrically sampled gaussian and DOG pyramids
if interactive >= 1
fprintf( 2, 'Expanding the Gaussian and DOG pyramids...\n' );
end
tic;
for octave = 1:octaves
if interactive >= 1
fprintf( 2, '\tProcessing octave %d: image size %d x %d subsample %.1f\n', octave, size(gauss_pyr{octave,1},2), size(gauss_pyr{octave,1},1), subsample(octave) );
fprintf( 2, '\t\tInterval 1 sigma %f\n', absolute_sigma(octave,1) );
end
sigma = initial_sigma;
g = gaussian_filter( preblur_sigma );
filter_size( octave, 1 ) = length(g);
filter_sigma( octave, 1 ) = preblur_sigma;% 注意,对于每一组中的第一帧图像都没有用高斯函数平滑,平滑从第二帧开始
DOG_pyr{octave} = zeros(size(gauss_pyr{octave,1},1),size(gauss_pyr{octave,1},2),intervals+2);