将多尺度SIFT改为单尺度上应用

最新推荐文章于 2022-05-28 11:44:09 发布

moonlightran

最新推荐文章于 2022-05-28 11:44:09 发布

阅读量2k

点赞数

文章标签： Harris-SIFT

本文链接：https://blog.csdn.net/moonlightran/article/details/25035931

版权

在网上看到好多Harris-SIFT特征提取的文章，可惜找不到程序，因此本人痛下决心，决定自己改写一份。

我们都知道，SIFT是在多尺度上提取特征点，生成特征描述子，然后再匹配。

在多尺度上提取的一个特征点有三个信息：分别是位置信息、尺度信息、方向信息；

将SIFT改在单尺度上应用，提取的特征点就只能包含位置信息和方向信息了，然后通过这两个信息来生成特征描述子。

一般情况下，我们可以根据自己需要的算法提取出特征点，这时候后提取的特征点仅仅拥有位置信息，如果利用区域灰度匹配的话，我们通常会以特征点为中心选择一块区域去做匹配。现在我们在单尺度上改为应用SIFT中生成特征描述子的方法来匹配。

首先，根据我们自己的算法提取出想要的特征点以后，特征点具有位置信息，下来我们该对特征点赋予方向信息：本人依据多尺度SIFT自己改写的在单尺度下的特征点赋予主方向的程序为：

frames=MajorOrientation(mag,ang,x,y)

其中x,y 为我们提取特征的位置坐标，mag,ang为以位置坐标为中心的某领域的方向和角度矩阵，在调用MajorOrientation前首先的求出来；

求mag和ang的程序如下：其中img为包含特征点的图片

function [mag,ang]=MagAndAng(img)
%input:
% img:input image;
%output:
% mag: the magnitudes;
% ang: the angles;

%the fucniton is to calculate the magintudes and angles and save it in mag
%and ang;
%step 1:compute x and y derivates using pixel difference;
% %compute x and y derivates using pixel difference;
% diff_x = 0.5*(img(2:(end-1),3:(end))-img(2:(end-1),1:(end-2)));
% diff_y = 0.5*(img(3:(end),2:(end-1))-img(1:(end-2),2:(end-1)));
dx_filter = [-0.5 0 0.5];
dy_filter = dx_filter';
gradient_x = imfilter(img, dx_filter);
gradient_y = imfilter(img, dy_filter);
gradient_x=double(gradient_x);
gradient_y=double(gradient_y);
% %compute the magnitude of the gradient;
magnitudes=sqrt(gradient_x.^2+gradient_y.^2);

% %compute the orientation of the gradient;
% angles=atan2( gradient_x, gradient_y );
angles=mod(atan(gradient_y ./ (eps + gradient_x)) + 2*pi, 2*pi);
% angles(find(angles == pi)) = -pi;
% Store the magnitude of the gradient and the orientation of the gradient
% mag=magnitudes(2:end-1,2:end-1);
% ang=angles(2:end-1,2:end-1);
mag=magnitudes;
ang=angles;

end

根据mag和ang计算主方向的程序如下：

function frames=MajorOrientation(mag,ang,x,y)
% The next step of the algorithm is to assign orientations to the keypoints
% that have been located. This is done by looking for peaks in histograms of
% gradient orientations in regions surrounding each keypoint. A keypoint may be
% assigned more than one orientation. If it is, then two identical descriptors
% are added to the database with different orientations.
%input:
%mag:幅度
%ang:角度
%x;特征点x坐标
%y:特征点y坐标
%output:
%frames:位置信息及主方向和若干个从方向信息;
win_factor = 1.5 ;
NBINS = 36; %直方图36柱
histo = zeros(1, NBINS); %初始化直方图数组
sigmaw = win_factor *6.4 ; %高斯加权因子
W = floor(3.0* 1.5*6.4); %高斯加权窗口直径
[M, N] = size(mag); % M 是图像的高度, N 是图像的宽度; num_level 是该组尺度空间的层数

xp=x;
yp=y;
magnitudes=mag;
angles=ang;

for xs = xp - max(W, xp-1): min((N - 2), xp + W)
for ys = yp - max(W, yp-1) : min((M-2), yp + W)
dx = (xs - x);
dy = (ys - y);
if dx^2 + dy^2 <= W^2 % 点在高斯加权圆内
wincoef = exp(-(dx^2 + dy^2)/(2*sigmaw^2));
bin = round( NBINS * angles(ys, xs)/(2*pi) + 0.5); %最近取整，确定所在方向直方图的柱

histo(bin) = histo(bin) + wincoef * magnitudes(ys, xs); %用高斯加权后的向量幅度做直方图值累加
end

end
end
% figure(2);
% bar(histo);
frames=[];
theta_max = max(histo); %找到直方图峰值
theta_indx = find(histo> 0.8 * theta_max); %关键点方向索引，大于80%峰值的角度
for i = 1: size(theta_indx, 2)
theta = 2*pi * theta_indx(i) / NBINS;
frames = [frames, [xp yp theta]'];
end

end

通过函数MajorOrientation传出的参数frames 我们可以看到，frames包含每个特征点的位置信息xp,yp,还有包含主方向信息theta;

接下来我们就依据位置和主方向信息在单尺度上生成特征描述子，其程序如下：

function desc=DESC(frames,mag,angles)
nbp=4;
desc=zeros(1,nbp*nbp*8);

o=1;
s=1;
[M,N]=size(mag);
sigm0=3.2;
sbp=3*sigm0;
W=floor(sqrt(2)*sbp*(nbp+1)/2+0.5);
xp=fix(frames(1,1));
yp=fix(frames(2,1));
theta0=frames(3,1);
sin0=sin(theta0);
cos0=cos(theta0);
index=0;
histo=zeros(1,nbp*nbp*8);
for xs = xp - min(W, xp-1): min((M - 1), xp + W)
for ys = yp - min(W, yp-1) : min((N-1), yp + W)
dx=xp-xs;
dy=yp-ys;
if dx^2+dy^2<W^2
theta=angles(xs,ys);
theta=mod((theta-theta0),2*pi);
nx=(dy*cos0+dx*sin0)/sbp;
ny=(-dy*sin0+dx*cos0)/sbp;
nt=theta/pi*4;
wsigma = nbp/2;
mags=mag(xs,ys);
wc=exp(-(nx*nx + ny*ny)/(2.0 * wsigma * wsigma));
for cx=-1.5:1:1.5
for cy=-1.5:1:1.5

if abs(nx-cx)<1&&abs(ny-cy)<1
ki=cx+2.5;
kj=cy+2.5;
ni1=ceil(nt);
ni2=floor(nt);
d=(1-abs(nx-cx))*(1-abs(ny-cy));
index=(((ki-1)*nbp+kj)-1)*8;
if ni1~=ni2

dt1=1-(ni1-nt);
dt2=1-(-ni2+nt);
index1=index+mod(ni1,8)+1;
index2=index+ni2+1;
histo(index1)=histo(index1)+wc*d*dt1*mags;
histo(index2)=histo(index2)+wc*d*dt2*mags;
else
histo(index+ni1+1)=histo(index+ni1+1)+wc*d*mags;
end
end
end
end
end
end
end
desc=histo;

desc = desc ./ norm(desc);
indx = desc > 0.2;
desc(indx) = 0.2;
desc = desc ./ norm(desc);