SIFT关键代码.docx

有关SIFT函数的关键代码(2010-06-08 13:35:18) 转载标签： 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 )if exist(octaves) octaves = 4;endif exist(intervals) intervals = 2;endif exist(object_mask) object_mask = ones(size(im);endif size(object_mask) = size(im) object_mask = ones(size(im);endif exist(contrast_threshold) contrast_threshold = 0.02;endif exist(curvature_threshold) curvature_threshold = 10.0;endif exist(interactive) interactive = 1;end% Check that the image is normalized to 0,1if( (min(im(:) 1) ) fprintf( 2, Warning: image not normalized to 0,1.n );end% 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 );endtic;antialias_sigma = 0.5;if antialias_sigma = 0 signal = im;else g = gaussian_filter( antialias_sigma );%至少比0.5大 if exist(corrsep) = 3 signal = corrsep( g, g, im ); else signal = conv2( g, g, im, same ); endendsignal = im;X Y = meshgrid( 1:0.5:size(signal,2), 1:0.5:size(signal,1) );signal = interp2( signal, X, Y, *linear );subsample = 0.5; % subsampling rate for doubled image is 1/2if interactive = 1 fprintf( 2, Prebluring image.n );endpreblur_sigma = sqrt(sqrt(2)2 - (2*antialias_sigma)2);if preblur_sigma = 0 gauss_pyr1,1 = signal;else g = gaussian_filter( preblur_sigma ); if exist(corrsep) = 3 gauss_pyr1,1 = corrsep( g, g, signal ); else gauss_pyr1,1 = conv2( g, g, signal, same ); endendclear signalpre_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_sigma2 );% Keep track of the absolute sigma for the octave and scaleabsolute_sigma = zeros(octaves,intervals+3);absolute_sigma(1,1) = initial_sigma * subsample(1);% Keep track of the filter sizes and standard deviations used to generate the pyramidfilter_size = zeros(octaves,intervals+3);filter_sigma = zeros(octaves,intervals+3);% Generate the remaining levels of the geometrically sampled gaussian and DOG pyramidsif interactive = 1 fprintf( 2, Expanding the Gaussian and DOG pyramids.n );endtic;for octave = 1:octaves if interactive = 1 fprintf( 2, tProcessing octave %d: image size %d x %d subsample %.1fn, octave, size(gauss_pyroctave,1,2), size(gauss_pyroctave,1,1), subsample(octave) ); fprintf( 2, ttInterval 1 sigma %fn, absolute_sigma(octave,1) ); end sigma = initial_sigma; g = gaussian_filter( sigma ); filter_size( octave, 1 ) = length(g); filter_sigma( octave, 1 ) = sigma; DOG_pyroctave = zeros(size(gauss_pyroctave,1,1),size(gauss_pyroctave,1,2),intervals+2);%M*N*Z for interval = 2:(intervals+3)%从第二层开始方便差分DOG计算 % Compute the standard deviation of the gaussian filter needed to % produce the % next level of the geometrically sampled pyramid. Here, sigma_i+1 = k*sigma. % By definition of successive convolution, the required blurring sigma_f to % produce sigma_i+1 from sigma_i is: % % sigma_i+12 = sigma_f,i2 + sigma_i2 % (k*sigma_i)2 = sigma_f,i2 + sigma_i2 % % therefore: % % sigma_f,i = sqrt(k2 - 1)sigma_i % % where k = 2(1/intervals) to span the octave, so: % % sigma_f,i = sqrt(2(2/intervals) - 1)sigma_i sigma_f = sqrt(2(2/intervals) - 1)*sigma; g = gaussian_filter( sigma_f ); sigma = (2(1/intervals)*sigma;%得到下一个sigma % Keep track of the absolute sigma absolute_sigma(octave,interval) = sigma * subsample(octave); % Store the size and standard deviation of the filter for later use filter_size(octave,interval) = length(g); filter_sigma(octave,interval) = sigma; if exist(corrsep) = 3 gauss_pyroctave,interval = corrsep( g, g, gauss_pyroctave,interval-1 ); else%卷积后获得当前层的高斯金子塔 gauss_pyroctave,interval = conv2( g, g, gauss_pyroctave,interval-1, same ); end%获得DOG金字塔 DOG_pyroctave(:,:,interval-1) = gauss_pyroctave,interval - gauss_pyroctave,interval-1; if interactive = 1 fprintf( 2, ttInterval %d sigma %fn, interval, absolute_sigma(octave,interval) ); end end if octave = 1 fprintf( 2, Pryamid processing time %.2f seconds.n, pyr_time );end% Display the gaussian pyramid when in interactive modeif interactive = 2 sz = zeros(1,2); sz(2) = (intervals+3)*size(gauss_pyr1,1,2); for octave = 1:octaves sz(1) = sz(1) + size(gauss_pyroctave,1,1); end pic = zeros(sz); y = 1; for octave = 1:octaves%显示所有组所有层的图像 x = 1; sz = size(gauss_pyroctave,1); for interval = 1:(intervals + 3)pic(y:(y+sz(1)-1),x:(x+sz(2)-1) = gauss_pyroctave,interval; x = x + sz(2); end y = y + sz(1); end fig = figure; clf; imshow(pic); resizeImageFig( fig, size(pic), 0.25 ); fprintf( 2, The gaussian pyramid (0.25 scale).nPress any key to continue.n ); pause; close(fig)end% Display the DOG pyramid when in interactive modeif interactive = 2 sz = zeros(1,2); sz(2) = (intervals+2)*size(DOG_pyr1(:,:,1),2); for octave = 1:octaves sz(1) = sz(1) + size(DOG_pyroctave(:,:,1),1); end pic = zeros(sz); y = 1; for octave = 1:octaves x = 1; sz = size(DOG_pyroctave(:,:,1); for interval = 1:(intervals + 2)pic(y:(y+sz(1)-1),x:(x+sz(2)-1) = DOG_pyroctave(:,:,interval); x = x + sz(2); end y = y + sz(1); end fig = figure; clf; imshow(pic); resizeImageFig( fig, size(pic), 0.25 ); fprintf( 2, The DOG pyramid (0.25 scale).nPress any key to continue.n ); pause; close(fig)end%已经可能到这边xx = 1 -2 1 ;yy = xx;xy = 1 0 -1; 0 0 0; -1 0 1 /4;%| Dxx Dxy |%| |%| Dxy Dyy |%| |% Coordinates of keypoints after each stage of processing for display% in interactive mode.raw_keypoints = ;contrast_keypoints = ;curve_keypoints = ;% Detect local maxima in the DOG pyramidif interactive = 1 fprintf( 2, Locating keypoints.n );endtic;loc = cell(size(DOG_pyr); % boolean maps of keypointsfor octave = 1:octaves if interactive = 1 fprintf( 2, tProcessing octave %dn, octave ); end for interval = 2:(intervals+1) keypoint_count = 0; contrast_mask = abs(DOG_pyroctave(:,:,interval) = contrast_threshold; lococtave,interval = zeros(size(DOG_pyroctave(:,:,interval); if exist(corrsep) = 3 edge = 1; else edge = ceil(filter_size(octave,interval)/2); end for y=(1+edge):(size(DOG_pyroctave(:,:,interval),1)-edge) for x=(1+edge):(size(DOG_pyroctave(:,:,interval),2)-edge) % Only check for extrema where the object mask is 1 if object_mask(round(y*subsample(octave),round(x*subsample(octave) = 1 if( (interactive = 2) | (contrast_mask(y,x) = 1) ) tmp = DOG_pyroctave(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(:) ) % The point is a local extrema of the DOG image. Store its coordinates for % displaying keypoint location in interactive mode. raw_keypoints = raw_keypoints; x*subsample(octave) y*subsample(octave); if abs(DOG_pyroctave(y,x,interval) = contrast_threshold % The DOG image at the extrema is above the contrast threshold. Store % its coordinates for displaying keypoint locations in interactive mode. contrast_keypoints = contrast_keypoints; raw_keypoints(end,:);%最后一行加入对比度 % Compute the entries of the Hessian matrix at the extrema location. Dxx = sum(DOG_pyroctave(y,x-1:x+1,interval) .* xx); Dyy = sum(DOG_pyroctave(y-1:y+1,x,interval) .* yy); Dxy = sum(sum(DOG_pyroctave(y-1:y+1,x-1:x+1,interval) .* xy); % Compute the trace and the determinant of the Hessian. Tr_H = Dxx + Dyy; Det_H = Dxx*Dyy - Dxy2; % Compute the ratio of the principal curvatures. curvature_ratio = (Tr_H2)/Det_H; if (Det_H = 0) & (curvature_ratio = 1 fprintf( 2, tt%d keypoints found on interval %dn, keypoint_count, interval ); end endendkeypoint_time = toc;if interactive = 1 fprintf( 2, Keypoint location time %.2f seconds.n, keypoint_time );end% Display results of extrema detection and keypoint filtering in interactive mode.if interactive = 2 fig = figure; clf; imshow(im); hold on; plot(raw_keypoints(:,1),raw_keypoints(:,2),y+); resizeImageFig( fig, size(im), 2 ); fprintf( 2, DOG extrema (2x scale).nPress any key to continue.n ); pause; close(fig); fig = figure; clf; imshow(im); hold on; plot(contrast_keypoints(:,1),contrast_keypoints(:,2),y+); resizeImageFig( fig, size(im), 2 ); fprintf( 2, Keypoints after removing low contrast extrema (2x scale).nPress any key to continue.n ); pause; close(fig); fig = figure; clf; imshow(im); hold on; plot(curve_keypoints(:,1),curve_keypoints(:,2),y+); resizeImageFig( fig, size(im), 2 ); fprintf( 2, Keypoints after removing edge points using principal curvature filtering (2x scale).nPress any key to continue.n ); pause; close(fig);endclear raw_keypoints contrast_keypoints curve_keypointsg = gaussian_filter( 1.5 * absolute_sigma(1,intervals+3) / subsample(1) );zero_pad = ceil( length(g) / 2 );if interactive = 1 fprintf( 2, Computing gradient magnitude and orientation.n );endtic;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) diff_x = 0.5*(gauss_pyroctave,interval(2:(end-1),3:(end)-gauss_pyroctave,interval(2:(end-1),1:(end-2); diff_y = 0.5*(gauss_pyroctave,interval(3:(end),2:(end-1)-gauss_pyroctave,interval(1:(end-2),2:(end-1); mag = zeros(size(gauss_pyroctave,interval); mag(2:(end-1),2:(end-1) = sqrt( diff_x . 2 + diff_y . 2 ); mag_pyroctave,interval = zeros(size(mag)+2*zero_pad); mag_pyroctave,interval(zero_pad+1):(end-zero_pad),(zero_pad+1):(end-zero_pad) = mag; grad = zeros(size(gauss_pyroctave,interval); grad(2:(end-1),2:(end-1) = atan2( diff_y, diff_x ); grad(find(grad = pi) = -pi; grad_pyroctave,interval = zeros(size(grad)+2*zero_pad); grad_pyroctave,interval(zero_pad+1):(end-zero_pad),(zero_pad+1):(end-zero_pad) = grad; endendclear mag gradgrad_time = toc;if interactive = 1 fprintf( 2, Gradient calculation time %.2f seconds.n, grad_time );endnum_bins = 36;hist_step = 2*pi/num_bins;hist_orient = -pi:hist_step:(pi-hist_step);% Initialize the positions, orientations, and scale information% of the keypoints to emtpy matrices.pos = ;orient = ;scale = ;if interactive = 1 fprintf( 2, Assigining keypoint orientations.n );endtic;for octave = 1:octaves if interactive = 1 fprintf( 2, tProcessing octave %dn, octave ); end for interval = 2:(intervals + 1) if interactive = 1 fprintf( 2, ttProcessing interval %d , interval ); end 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(lococtave,interval)+2*zero_pad); loc_pad(zero_pad+1):(end-zero_pad),(zero_pad+1):(end-zero_pad) = lococtave,interval; iy ix=find(loc_pad=1); for k = 1:length(iy) x = ix(k); y = iy(k); wght = g.*mag_pyroctave,interval(y-hf_sz):(y+hf_sz),(x-hf_sz):(x+hf_sz); grad_window = grad_pyroctave,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) % Compute the diference of the orientations mod pi diff = mod( grad_window - hist_orient(bin) + pi, 2*pi ) - pi; % Accumulate the histogram bins orient_hist(bin)=orient_hist(bin)+sum(sum(wght.*max(1 - abs(diff)/hist_step,0); %orient_hist(bin)=orient_hist(bin)+sum(sum(wght.*(abs(diff) = hist_step); end % Find peaks in the orientation histogram using nonmax suppression. 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 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 % Store the keypoint position, orientation, and scale information 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; % Get the next peak peaks(ipeak) = 0; peak_val ipeak = max(peaks); end end if interactive = 1 fprintf( 2, (%d keypoints)n, keypoint_count ); end endendclear loc loc_padorient_time = toc;if interactive = 1 fprintf( 2, Orientation assignment time %.2f seconds.n, orient_time );end% Display the keypoints with scale and orientation in interactive mode.if interactive = 2 fig = figure; clf; imshow(im); hold on; display_keypoints( pos, scale(:,3), orient, y ); resizeImageFig( fig, size(im), 2 ); fprintf( 2, Final keypoints with scale and orientation (2x scale).nPress any key to continue.n ); pause; close(fig);endorient_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 = ;if interactive = 1 fprintf( 2, Computing keypoint feature descriptors for %