【Matlab】SIFT_feature代码学习

%this code is the Matlab implimentation of David G. Lowe,%"Distinctive image features from scale-invariant keypoints,"%International Journal of Computer Vision, 60, 2 (2004), pp. 91-110.%this code should be used only for academic research.

tic  % tic and toc为time clear;clc;row=256;colum=256;img=imread('timg1.jpg');  %读文件img=imresize(img,[row,colum]); %读尺寸img=rgb2gray(img); %rgb颜色转换为grayimg=im2double(img); %变为双精度  im2double Convert image to double PRecision.origin=img; toc

%% Scale-Space Extrema Detectiontic% original sigma and the number of actave can be modified. the larger sigma0, the more quickly-smooth images越大的sigma0，越快的平滑图片sigma0=sqrt(2);octave=3;              %6*sigma*k^(octave*level)<=min(m,n)/(2^(octave-2))level=3;D=cell(1,octave);     %    cell(M,N) or cell([M,N]) is an M-by-N cell array of empty matrices.

for i=1:octaveD(i)=mat2cell(zeros(row*2^(2-i)+2,colum*2^(2-i)+2,level),row*2^(2-i)+2,colum*2^(2-i)+2,level);% mat2cell Break matrix up into a cell array of matrices.end

% first image in first octave is created by interpolating the original one.   %上采样插值获得新图形temp_img=kron(img,ones(2));    % ones(N) is an N-by-N matrix of ones.

% kron(X,Y) is the Kronecker tensor product of X and Y.    The result is a large matrix formed by taking all possible    products between the elements of X and those of Y. For    example, if X is 2 by 3, then kron(X,Y) is        [ X(1,1)*Y  X(1,2)*Y  X(1,3)*Y         X(2,1)*Y  X(2,2)*Y  X(2,3)*Y ]

temp_img=padarray(temp_img,[1,1],'replicate');

%    B = padarray(A,PADSIZE,METHOD,DIRECTION) pads array A using the    specified METHOD.  METHOD can be one of these strings:

    'replicate'   Repeats border elements of A.

figure(2)subplot(1,2,1);imshow(origin)   

%进行上采样，重复各个像素的值并且重复边缘像素

%create the DoG pyramid.for i=1:octave %层    temp_D=D{i};    for j=1:level          scale=sigma0*sqrt(2)^(1/level)^((i-1)*level+j); % 和level有关 1.5874 1.7818 2.0000 2.2449 2.5198 2.8284 3.1748 3.5636 4.0000        p=(level)*(i-1);        figure(1);        subplot(octave,level,p+j);        f=fspecial('gaussian',[1,floor(6*scale)],scale);        L1=temp_img;        if(i==1&&j==1)        L2=conv2(temp_img,f,'same');        L2=conv2(L2,f','same');        temp_D(:,:,j)=L2-L1;        imshow(uint8(255 * mat2gray(temp_D(:,:,j))));  % mat2gray Convert matrix to intensity image.        L1=L2;        else        L2=conv2(temp_img,f,'same');        L2=conv2(L2,f','same');        temp_D(:,:,j)=L2-L1;        L1=L2;        if(j==level)            temp_img=L1(2:end-1,2:end-1);        end        imshow(uint8(255 * mat2gray(temp_D(:,:,j))));        end    end    D{i}=temp_D;    temp_img=temp_img(1:2:end,1:2:end);    temp_img=padarray(temp_img,[1,1],'both','replicate');endtoc

Keypoint Localistaion

%% Keypoint Localistaion% search each pixel in the DoG map to find the extreme pointticinterval=level-1;number=0;for i=2:octave+1    number=number+(2^(i-octave)*colum)*(2*row)*interval;endextrema=zeros(1,4*number);flag=1;for i=1:octave    [m,n,~]=size(D{i});    m=m-2;    n=n-2;    volume=m*n/(4^(i-1));    for k=2:interval              for j=1:volume            % starter=D{i}(x+1,y+1,k);            x=ceil(j/n);            y=mod(j-1,m)+1;            sub=D{i}(x:x+2,y:y+2,k-1:k+1);            large=max(max(max(sub)));            little=min(min(min(sub)));            if(large==D{i}(x+1,y+1,k))                temp=[i,k,j,1];                extrema(flag:(flag+3))=temp;                flag=flag+4;            end            if(little==D{i}(x+1,y+1,k))                temp=[i,k,j,-1];                extrema(flag:(flag+3))=temp;                flag=flag+4;            end        end    endendidx= extrema==0;extrema(idx)=[];toc[m,n]=size(img);x=floor((extrema(3:4:end)-1)./(n./(2.^(extrema(1:4:end)-2))))+1;y=mod((extrema(3:4:end)-1),m./(2.^(extrema(1:4:end)-2)))+1;ry=y./2.^(octave-1-extrema(1:4:end));rx=x./2.^(octave-1-extrema(1:4:end));figure(2)subplot(1,2,2);imshow(origin)hold onplot(ry,rx,'r+');