matlab实现sift提取特征点

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);