图像sift特征提取实验

%%% Test sift Dog clear;close all; Img=imread('F:\乱七八糟\photo\im24.jpg'); img=double(rgb2gray(Img)); [M,N]=size(img); S=3; sigma=0.5;thresh=0.03;gama=10; % param Ht=cell(S+3,1); Gauss=cell(S+3,1); Dog=cell(S+2,1); for i=1:S+3 Ht{i}=fspecial('gaussian', 6, sigma*2^((i-1)/S)); Gauss{i}=imfilter(img,Ht{i},'replicate'); end for i=1:S+2 Dog{i}=Gauss{i+1}-Gauss{i}; end %% 遍历搜索极值 % max&min Improc=cell(S+2,1); % 进行处理的中间变量 for i=1:S+2 Improc{i}=zeros(M+2,N+2); Improc{i}(2:M+1,2:N+1)=Dog{i}; % 扩充 Improc{i}(1:M+2,1)=Improc{i}(1:M+2,2); Improc{i}(1:M+2,end)=Improc{i}(1:M+2,end-1); Improc{i}(1,2:N+1)=Improc{i}(2,2:N+1); Improc{i}(end,2:N+1)=Improc{i}(end-1,2:N+1); end indic=zeros(M,N,S); for i=2:S+1 for j=2:M+1 for k=2:N+1 if Improc{i}(j,k)>Improc{i}(j+1,k) && Improc{i}(j,k)>Improc{i}(j-1,k) && Improc{i}(j,k)>Improc{i}(j,k+1) && Improc{i}(j,k)>Improc{i}(j,k-1) && Improc{i}(j,k)>Improc{i}(j+1,k+1) && Improc{i}(j,k)>Improc{i}(j-1,k-1) && Improc{i}(j,k)>Improc{i}(j+1,k-1) && Improc{i}(j,k)>Improc{i}(j-1,k+1)... && Improc{i}(j,k)>Improc{i+1}(j,k) && Improc{i}(j,k)>Improc{i+1}(j+1,k) && Improc{i}(j,k)>Improc{i+1}(j-1,k) && Improc{i}(j,k)>Improc{i+1}(j,k+1) && Improc{i}(j,k)>Improc{i+1}(j,k-1) && Improc{i}(j,k)>Improc{i+1}(j+1,k+1) && Improc{i}(j,k)>Improc{i+1}(j-1,k-1) && Improc{i}(j,k)>Improc{i+1}(j+1,k-1) && Improc{i}(j,k)>Improc{i+1}(j-1,k+1)... && Improc{i}(j,k)>Improc{i-1}(j,k) && Improc{i}(j,k)>Improc{i-1}(j+1,k) && Improc{i}(j,k)>Improc{i-1}(j-1,k) && Improc{i}(j,k)>Improc{i-1}(j,k+1) && Improc{i}(j,k)>Improc{i-1}(j,k-1) && Improc{i}(j,k)>Improc{i-1}(j+1,k+1) && Improc{i}(j,k)>Improc{i-1}(j-1,k-1) && Improc{i}(j,k)>Improc{i-1}(j+1,k-1) && Improc{i}(j,k)>Improc{i-1}(j-1,k+1) indic(j-1,k-1,i-1)=1; end if Improc{i}(j,k)<Improc{i}(j+1,k) && Improc{i}(j,k)<Improc{i}(j-1,k) && Improc{i}(j,k)<Improc{i}(j,k+1) && Improc{i}(j,k)<Improc{i}(j,k-1) && Improc{i}(j,k)<Improc{i}(j+1,k+1) && Improc{i}(j,k)<Improc{i}(j-1,k-1) && Improc{i}(j,k)<Improc{i}(j+1,k-1) && Improc{i}(j,k)<Improc{i}(j-1,k+1)... && Improc{i}(j,k)<Improc{i+1}(j,k) && Improc{i}(j,k)<Improc{i+1}(j+1,k) && Improc{i}(j,k)<Improc{i+1}(j-1,k) && Improc{i}(j,k)<Improc{i+1}(j,k+1) && Improc{i}(j,k)<Improc{i+1}(j,k-1) && Improc{i}(j,k)<Improc{i+1}(j+1,k+1) && Improc{i}(j,k)<Improc{i+1}(j-1,k-1) && Improc{i}(j,k)<Improc{i+1}(j+1,k-1) && Improc{i}(j,k)<Improc{i+1}(j-1,k+1)... && Improc{i}(j,k)<Improc{i-1}(j,k) && Improc{i}(j,k)<Improc{i-1}(j+1,k) && Improc{i}(j,k)<Improc{i-1}(j-1,k) && Improc{i}(j,k)<Improc{i-1}(j,k+1) && Improc{i}(j,k)<Improc{i-1}(j,k-1) && Improc{i}(j,k)<Improc{i-1}(j+1,k+1) && Improc{i}(j,k)<Improc{i-1}(j-1,k-1) && Improc{i}(j,k)<Improc{i-1}(j+1,k-1) && Improc{i}(j,k)<Improc{i-1}(j-1,k+1) indic(j-1,k-1,i-1)=1; end end end end %% 插值求亚像素极值点——根据像素点筛选 % 由于要计算x,y,sigma(z) 3个变量的偏微分，分别计算差分矩阵。 hy=[-0.5;0;0.5]; hx=[-0.5,0,0.5]; fx=zeros(M+2,N+2,S); fy=zeros(M+2,N+2,S); fxx=zeros(M+2,N+2,S); fyy=zeros(M+2,N+2,S); fxy=zeros(M+2,N+2,S); fyx=zeros(M+2,N+2,S); fz=zeros(M+2,N+2,S); fzy=zeros(M+2,N+2,S); fzx=zeros(M+2,N+2,S); fzz=zeros(M+2,N+2,S); fxz=zeros(M+2,N+2,S); fyz=zeros(M+2,N+2,S); for i=1:S fx(:,:,i)=imfilter(Improc{i+1},hx,'replicate'); fy(:,:,i)=imfilter(Improc{i+1},hy,'replicate'); fxx(:,:,i)=imfilter(fx(:,:,i),hx,'replicate'); fxy(:,:,i)=imfilter(fx(:,:,i),hy,'replicate'); fyx(:,:,i)=imfilter(fy(:,:,i),hx,'replicate'); fyy(:,:,i)=imfilter(fy(:,:,i),hy,'replicate'); fz(:,:,i)=(Improc{i+2}-Improc{i})/2; fzy(:,:,i)=imfilter(fz(:,:,i),hy,'replicate'); fzx(:,:,i)=imfilter(fz(:,:,i),hx,'replicate'); if i==1 fzz(:,:,i)=(fz(:,:,i+1)-Improc{i+1}/2)/2; fxz(:,:,i)=(fx(:,:,i+1)-imfilter(Improc{i},hx,'replicate'))/2; fyz(:,:,i)=(fy(:,:,i+1)-imfilter(Improc{i},hy,'replicate'))/2; elseif i==S fzz(:,:,i)=(-Improc{i+1}/2-fz(:,:,i-1))/2; fxz(:,:,i)=(imfilter(Improc{i+2},hx,'replicate')-fx(:,:,i-1))/2; fyz(:,:,i)=(imfilter(Improc{i+2},hy,'replicate')-fy(:,:,i-1))/2; else fzz(:,:,i)=(fz(:,:,i+1)-fz(:,:,i-1))/2; fxz(:,:,i)=(fx(:,:,i+1)-fx(:,:,i-1))/2; fyz(:,:,i)=(fy(:,:,i+1)-fy(:,:,i-1))/2; end end % 对每个极值点像素，插值计算极值点偏差x^和响应|D(x^)|。 Indic=zeros(size(indic,1),size(indic,2)); for i=1:S Ind=find(indic(:,:,i)==1); [I,J]=ind2sub(size(indic(:,:,i)),Ind); % 还需各+1，因为少了4条边。 D_x=zeros(length(Ind),1); for j=1:length(Ind) F=[fxx(I(j)+1,J(j)+1,i),fxy(I(j)+1,J(j)+1,i),fxz(I(j)+1,J(j)+1,i); fyx(I(j)+1,J(j)+1,i),fyy(I(j)+1,J(j)+1,i),fyz(I(j)+1,J(j)+1,i); fzx(I(j)+1,J(j)+1,i),fzy(I(j)+1,J(j)+1,i),fzz(I(j)+1,J(j)+1,i)]; Finv=inv(F); x_x=-Finv*[fx(I(j)+1,J(j)+1,i);fy(I(j)+1,J(j)+1,i);fz(I(j)+1,J(j)+1,i)]; D_x(j)=Improc{i+1}(I(j)+1,J(j)+1)+0.5*[fx(I(j)+1,J(j)+1,i),fy(I(j)+1,J(j)+1,i),fz(I(j)+1,J(j)+1,i)]*x_x; if max(x_x)>0.5 || D_x(j)/max(max(img))<thresh indic(I(j),J(j),i)=0; end end Indic=Indic | indic(:,:,i); end disp(sum(sum(Indic))); % Hessian矩阵删除边缘效应，fxx、fxy、fyy等已求出。 for i=1:S Ind=find(indic(:,:,i)==1); [I,J]=ind2sub(size(indic(:,:,i)),Ind); % 还需各+1，因为少了4条边。 for j=1:length(Ind) H=[fxx(I(j)+1,J(j)+1,i),fxy(I(j)+1,J(j)+1,i);fxy(I(j)+1,J(j)+1,i),fyy(I(j)+1,J(j)+1,i)]; Tr=trace(H); Det=det(H); if Tr^2/Det >= (gama+1)^2/gama indic(I(j),J(j),i)=0; end end Indic=Indic | indic(:,:,i); end disp(sum(sum(Indic))); %% 绘出关键点的位置 imgshow=drawkeypoints(img,J+1,I+1,15); imgshow=uint8(imgshow); imshow(imgshow);