【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颜色转换为gray
img=im2double(img); %变为双精度  im2double Convert image to double precision.
origin=img;
toc

%% Scale-Space Extrema Detection
tic
% 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.

>> D=cell(1,3)

D =

    []    []    []

for i=1:octave
D(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');
end
toc

Keypoint Localistaion

%% Keypoint Localistaion
% search each pixel in the DoG map to find the extreme point
tic
interval=level-1;
number=0;
for i=2:octave+1
    number=number+(2^(i-octave)*colum)*(2*row)*interval;
end
extrema=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
    end
end
idx= 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 on
plot(ry,rx,'r+');