📄 sift.m
字号:
function [image,count,Keypoint,featureAll] = SIFT(imagefile)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SIFT 算法测试程序
% identify the feature descriptors
% program by super xiaoweige
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 12.1=(10+1)^2/10
% 10 = (8+1)^2/8;
% 7.2= (5+1)^2/5;
% 4.5= (2+1)^2/2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%image = imread('E:\MatlabWork\sift\scene.pgm');
image=imread(imagefile);
if isrgb(image)
image = rgb2gray(image);
end
image = im2double(image);
[ix,iy] = size(image);
%检测可能的感兴趣点。即检测尺度空间中的关于尺度的极值点。即找出相邻高斯滤波后的极值点坐标和尺度（即sigma）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% double the input image using linear interpolation prior to building the first level of the pyramid
%imageT = double(zeros(2*ix-1,2*iy-1));
%for i = 1:ix
% for j = 1:iy
% imageT(2*i-1,2*j-1) = image(i,j);
% end
%end
%for i = 1:ix-1
% for j = 1:iy
% imageT(2*i,2*j-1) = (imageT(2*i-1,2*j-1)+imageT(2*i+1,2*j-1))/2;
% end
%end
%for i = 1:2*ix-1
% for j = 1:iy-1
% imageT(i,2*j) = (imageT(i,2*j-1)+imageT(i,2*j+1))/2;
% end
%end
imageT = image;
[sizeX,sizeY] = size(imageT);
imageT = imageT-min(imageT(:)); % transform the image to [0 1]
imageT = imageT/max(imageT(:));
FirstSigma = 0.8;
Gau = fspecial('gaussian',[9 9],FirstSigma); %prior smoothing using sigma = 1.414
imageT=imfilter(imageT,Gau); %与高斯模板做卷积
count = 0;
Keypoint = zeros(1,6);
featureAll= zeros(1,128);
RunNum = 0;%RunNum = -1;
ScaDoG = zeros(6,1);
while (sizeX>20 & sizeY>20)
ix=sizeX
iy=sizeY
IMG = zeros(sizeX,sizeY,6);
DoG = zeros(sizeX,sizeY,5); % DoG
IMG(:,:,1) = imageT;
%gwid = max(5,fix(Sig));
ScaDoG(1) = FirstSigma*2^RunNum;
for i=1:1:5
%Sigma = Sig+SigDiff*(i-3);
%GauWid = [2*gwid+1 2*gwid+1]; % 可变宽度的高斯滤波器 和不可变宽度的高斯滤波器. [2*i+1 2*i+1]
Sigma = 2^(i/5);
Gau = fspecial('gaussian',[9 9],Sigma);
IG=imfilter(imageT,Gau); %与高斯模板做卷积
IMG(:,:,i+1) = IG;
ScaDoG(i+1) = Sigma*ScaDoG(1);
end
clear Gau;
clear IG;
%DoG(:,:,1) = IMG(:,:,1) - imageT;
for i = 1:1:5;
DoG(:,:,i) = IMG(:,:,i+1)-IMG(:,:,i); % DoG
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% identify the keypoint
%L = DoG;
countT = 1;
countTT = 1;
for k = 3:4
for i = 11:1:ix-10
for j = 11:1:iy-10
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D = DoG(i,j,k);
DoGTemp = DoG(i-1:i+1,j-1:j+1,k);
if D == max(DoGTemp(:))
DoGTemp1 = DoG(i-1:i+1,j-1:j+1,k-1);
if D > max(DoGTemp1(:))
DoGTemp2 = DoG(i-1:i+1,j-1:j+1,k+1);
if D > max(DoGTemp2(:))
countT = countT+1;
Dxx = DoG(i+1,j,k)+DoG(i-1,j,k)-2*DoG(i,j,k);
Dyy = DoG(i,j+1,k)+DoG(i,j-1,k)-2*DoG(i,j,k);
Dxy = DoG(i+1,j+1,k)+DoG(i,j,k)-DoG(i+1,j,k)-DoG(i,j+1,k);
r = (Dxx+Dyy)^2/(Dxx*Dyy-Dxy^2);
if r<12.1
countTT=countTT+1;
Dx = DoG(i+1,j,k)-DoG(i,j,k);
Dy = DoG(i,j+1,k)-DoG(i,j,k);
Dz = DoG(i,j,k+1)-DoG(i,j,k);
Dxz = DoG(i+1,j,k+1)+DoG(i,j,k)-DoG(i+1,j,k)-DoG(i,j,k+1);
Dyz = DoG(i,j+1,k+1)+DoG(i,j,k)-DoG(i,j+1,k)-DoG(i,j,k+1);
Dzz = DoG(i,j,k+1)+DoG(i,j,k-1)-2*DoG(i,j,k);
xhTemp = [Dxx Dxy Dxz;Dxy Dyy Dyz;Dxz Dyz Dzz];
xHead = inv(xhTemp)*[Dx Dy Dz]';
if abs(xHead(1))<1 & abs(xHead(2))<1 & abs(xHead(3))<0.5
DxHead = DoG(i,j,k) + 1/2*[Dx Dy Dz]*xHead;
%if abs(DxHead) >0.03
count = count+1;
key = [i+xHead(1),j+xHead(2),ScaDoG(k)+xHead(3),0,0,0];
%key = [scaleX scaleY scale 0 0 0];
[featureA,key] = descriptor(key,IMG(:,:,k),count,RunNum);
Keypoint(count,:) = key;
featureAll(count,:) = featureA;
%end
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif D==min(DoGTemp(:))
DoGTemp1 = DoG(i-1:i+1,j-1:j+1,k-1);
if D < min(DoGTemp1(:))
DoGTemp2 = DoG(i-1:i+1,j-1:j+1,k+1);
if D < min(DoGTemp2(:))
countT = countT+1;
Dxx = DoG(i+1,j,k)+DoG(i-1,j,k)-2*DoG(i,j,k);
Dyy = DoG(i,j+1,k)+DoG(i,j-1,k)-2*DoG(i,j,k);
Dxy = DoG(i+1,j+1,k)+DoG(i,j,k)-DoG(i+1,j,k)-DoG(i,j+1,k);
r = (Dxx+Dyy)^2/(Dxx*Dyy-Dxy^2);
if r<12.1
countTT=countTT+1;
Dx = DoG(i+1,j,k)-DoG(i,j,k);
Dy = DoG(i,j+1,k)-DoG(i,j,k);
Dz = DoG(i,j,k+1)-DoG(i,j,k);
Dxz = DoG(i+1,j,k+1)+DoG(i,j,k)-DoG(i+1,j,k)-DoG(i,j,k+1);
Dyz = DoG(i,j+1,k+1)+DoG(i,j,k)-DoG(i,j+1,k)-DoG(i,j,k+1);
Dzz = DoG(i,j,k+1)+DoG(i,j,k-1)-2*DoG(i,j,k);
xhTemp = [Dxx Dxy Dxz;Dxy Dyy Dyz;Dxz Dyz Dzz];
xHead = inv(xhTemp)*[Dx Dy Dz]';
if abs(xHead(1))<1 & abs(xHead(2))<1 & abs(xHead(3))<0.5
DxHead = DoG(i,j,k) + 1/2*[Dx Dy Dz]*xHead;
%if abs(DxHead) >0.03
count = count+1;
key = [i+xHead(1),j+xHead(2),ScaDoG(k)+xHead(3),0,0,0];
%key = [scaleX scaleY scale 0 0 0];
[featureA,key] = descriptor(key,IMG(:,:,k),count,RunNum);
Keypoint(count,:) = key;
featureAll(count,:) = featureA;
%end
end
end
end
end
end
end %for j
end %for i
end % for k = 2:4
countT
countTT
imageT = imageT(1:2:end,1:2:end);
[sizeX,sizeY] = size(imageT);
clear IMG;
clear DoG;
RunNum = RunNum +1;
end % while(sizeX+sizeY)/2 >
count
function [featureA,key] = descriptor(key,IM,icount,RunNum)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%for i = countPrevious+1:count
x = fix(key(1));
y = fix(key(2));
%k = key(3);
key(4) = sqrt((IM(x+1,y)-IM(x-1,y))^2+(IM(x,y+1)-IM(x,y-1))^2);
key(6) = atan2(IM(x,y+1)-IM(x,y-1),IM(x+1,y)-IM(x-1,y));
Orient = zeros(1,4); % 先设置关键点只有四个方向，即只包含与坐标轴平行的方向。上下左右，先看看效果，再考虑8个以上方向
for m = x-1:1:x+1
for n = y-1:1:y+1
theta = atan2(IM(m,n+1)-IM(m,n-1),IM(m+1,n)-IM(m-1,n));
thetaT = fix((theta + pi)/(pi/4)); % 将-pi到 pi的theta 转化成代表其区间整数
switch thetaT
case {0,7}
Orient(1) = Orient(1)+1;
case {1,2}
Orient(2) = Orient(2)+1;
case {3,4}
Orient(3) = Orient(3)+1;
case {5,6}
Orient(4) = Orient(4)+1;
otherwise
% 4
% 1 3
% 2
end
end
end % for m = x-1:1:x+1
[OriMax,OriI] = max(Orient);
key(5) = OriI; %
%T(:,i) = Orient';
key(1) = key(1)*(2^RunNum); % must be checked later.
key(2) = key(2)*(2^RunNum);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% feature descriptor & descriptor representation
ori = key(5);
igra = double(zeros(18,18)); % %因为需要16*16的特征矩阵来计算特征的128维表示.又因为需要计算每个点的梯度方向,所以,向各个方向放大一格
% igra 只是将IM矩阵中特征点附近对应的数值按照特征点方向的要求复制出来. 矩阵的X 方向与特征点的梯度方向一致
[sizeX,sizeY] = size(IM);
% 运算中先对图象的大小做了限定,不会有图片的宽度或长度少于20.
if x>=10 & x<=sizeX-8 & y>=10 & y<=sizeY-8
igra = IM(x-9:x+8,y-9:y+8);
end
%elseif x>=10 & x<=sizeX-8 & y>=10 & y>sizeY-8
%igraT = IM(x-9:x+8,y-9:end);
%igra(:,1:size(igraT,2)) = igraT;
%elseif x>=10 & x<=sizeX-8 & y<10 & y<=sizeY-8
switch ori
case {1}
igra = rot90(igra,2);
case {2}
rot90(igra);
case {3}
case {4}
igra = rot90(igra,3);
otherwise
end % switch ori
gra = double(zeros(16,16));
graMag = double(zeros(16,16));
DescriptorSig = 2;
for m = 1:16
for n = 1:16
graMag(m,n) = sqrt((igra(m+2,n+1)-igra(m,n+1))^2+(igra(m+1,n+2)-igra(m+1,n))^2);%*(1/(2*pi*DescriptorSig^2)*exp(-((m-8.5)^2 + (n-8.5)^2)/(2*DescriptorSig^2)));
gra(m,n) = atan2(igra(m+1,n+2)-igra(m+1,n),igra(m+2,n+1)-igra(m,n+1));
end
end
%featureA = double(zeros(1,128));
for m = 1:4
for n = 1:4
feaOri8 = double(zeros(1,8));
for p = (m-1)*4+1:(m-1)*4+4
for q = (n-1)*4+1:(n-1)*4+4
feaOriT = fix((gra(p,q)+pi)/(pi/8));
switch feaOriT
case {0,15}
feaOri8(1) = feaOri8(1) + graMag(p,q);
case {1,2}
feaOri8(2) = feaOri8(2) + graMag(p,q);
case {3,4}
feaOri8(3) = feaOri8(3) + graMag(p,q);
case {5,6}
feaOri8(4) = feaOri8(4) + graMag(p,q);
case {7,8}
feaOri8(5) = feaOri8(5) + graMag(p,q);
case {9,10}
feaOri8(6) = feaOri8(6) + graMag(p,q);
case {11,12}
feaOri8(7) = feaOri8(7) + graMag(p,q);
case {13,14}
feaOri8(8) = feaOri8(8) + graMag(p,q);
otherwise
end
end
end % for p = (m-1)*4+1:(m-1)*4+4
feaOriAll(:,(m-1)*4+n) = feaOri8';
end
end %for m = 1:4
featureA = feaOriAll;
featureA = reshape(featureA,1,128);
if sum(featureA.^2) > 0
featureA = featureA / sqrt(sum(featureA.^2)); % 特征归一化
end
%featureAll(icount,:) = featureA;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -