SIFT检测特征点之找到主方向
来源：互联网 发布：sas编程 编辑：程序博客网 时间：2024/05/01 11:22
接着前面的，前面分别得到了 ：高斯卷积，DOG差分尺度空间，检测极值点，去除两种不要的特征点（精确特征点），接下来就是第五步计算每个特征点的梯度mag和方向ori生成梯度直方图 这部分的理论可以参考
http://blog.csdn.net/v_JULY_v/article/details/6245939
http://blog.csdn.net/abcjennifer/article/details/7639681 只不过 在编的过程中 发现 这两个理论有的细节部分具体怎么实现还是讲得不够清楚 因为只有把细节都弄明白了 才能真得编好程序 可能这两位大神都觉得我们疑惑的细节问题太简单了 所以没写出来 好吧 只能自己去慢慢理解 摸索 调试
在此感谢陌生朋友ccworm帮我解惑 让我更加理解了生成梯度直方图这部分 才能正确编程实现 谢谢啦 如果他是男生 祝他越长越帅 如果她是女生 祝她越来越漂亮
言归正传 因为之前写的程序里已经得到装着特征点位置的E11矩阵 以及特征点D12图像
接着之前的程序写：
在以半径为rad的圆内遍历像素计算梯度的幅值mag和方向ori并进行高斯权重的累加
[X,Y]=find(E11==1 | E11==-1);
rad=1.5*0.4; %0.4是之前就已经写的D12的sigma 有个opencv程序这里的半径还乘以了一个3倍 我不知道为什么 我依旧按照论文上说的 只乘以一个1.5倍
[m,n]=size(X);
bin=zeros(1,36); %36个方向 即 36个柱子bin 即每个特征点领域的直方图 第一行全部装直方图的幅值 第二行用来装幅值对应的下标
max_mag=zeros(1,m); %这个用来放每个特征点主方向时对应的最大梯度值
allbin=zeros(1,m); %这个是用来装每个特征点的主方向的 我直接放下标 因为下标可以转换成角度 乘以一个10度就好了
for t=1:m
for i=-rad:rad
for j=-rad:rad
if(ceil(X(t)+i+1)>m1 || ceil(Y(t)+j+1)>n1)
continue;
end
mag=sqrt((D12(ceil(X(t)+i+1),ceil(Y(t)+j))-D12(ceil(X(t)+i-1),ceil(Y(t)+j)))^2+(D12(ceil(X(t)+i),ceil(Y(t)+j+1))-D12(ceil(X(t)+i),ceil(Y(t)+j-1)))^2);
ori=atan2(D12(ceil(X(t)+i),ceil(Y(t)+j+1))-D12(ceil(X(t)+i),ceil(Y(t)+j-1)),D12(ceil(X(t)+i+1),ceil(Y(t)+j))-D12(ceil(X(t)+i-1),ceil(Y(t)+j)));
w=exp(-(i^2+j^2)/(2*0.4^2)); %高斯权重
n=ceil(36*(ori+pi)/(2*pi)); %最终出来下标是1,2,3.....，35,36
if(n>36);
n=n-36;
end
bin(n)=bin(n)+w*mag;
end
end
[max_mag(t),allbin(t)]=max(bin);（将每个特征点的主方向（即最大幅值对应的直方图的下标）放在allbin里保存起来）
bin=zeros(1,36);
end
如果要对特征点的梯度直方图进行高斯平滑 可以自己在每个特征点直方图后面加上一次或者几次平滑（不过这个次数我不知道怎么确定 没查到讲这个的）
其实也可以不平滑 没太大影响 上面的程序就可以找到主方向
为了看效果 我对第一个特征点的直方图进行平滑并显示出来 如下所示
没进行高斯平滑之前 显示直方图如下所示：
i=1:36;
>> bar(i,bin,'g');
allbin(1)
ans =
21
即第一个特征点的主方向是21x10=210度 和上图显示出来的结果吻合 自己可以把下标换成弧度的角度
我先在这个结果的基础上设置7次平滑看看：
for i=1:7
bin(1)=0.25*bin(36)+0.5*bin(1)+0.25*bin(2);
for i=2:35
bin(i)=0.25*bin(i-1)+0.5*bin(i)+0.25*bin(i+1);
end
bin(36)=0.25*bin(35)+0.5*bin(36)+0.25*bin(1);
end
i=1:36;
>> bar(i,bin,'g');
可以看到还是第21个方向上幅值最大 即主方向仍是210度
接下来 在原图中将主方向和梯度值用向量的形式画出来：
[Y,X]=find(E11==1 | E11==-1); %这里千万别弄反 find函数第一个找到的是Y坐标 我之前弄反了 得到的是一个乍一看鬼一样的东西（就是倒过来的样子 四不像）
[m,n]=size(X);
imshow(H)
hold on
for t=1:m
len=round(max_mag(t)*5);
hlen=round(max_mag(t)*0.75);
blen=len-hlen;
theta=2*pi*allbin(t)*10/360;
ytermi=round(len*(-sin(theta)))+Y(t);
xtermi=round(len*(cos(theta)))+X(t);
xA=round(blen*(cos(theta-pi/18)))+X(t);
yA=round(blen*(-sin(theta-pi/18)))+Y(t);
xB=round(blen*(cos(theta+pi/18)))+X(t);
yB=round(blen*(-sin(theta+pi/18)))+Y(t);
line([X(t),xtermi],[Y(t),ytermi]);
line([xA,xtermi],[yA,ytermi]);
line([xB,xtermi],[yB,ytermi]);
end
hold off
怎么样 是不是很漂亮呢 虽然这些方向乱七八糟的 但在我心里 这依旧美得天昏地暗日月无光 美得我都不忍直视了
好了 不开玩笑了 接下来的工作就是生成描述子然后几幅图匹配了 加油！宇宙无敌超级狂拽酷炫屌炸天帅女侠！
我用图像处理女神图再试一下：
把全部代码重新贴一遍：
C=imread('F:\orl_zhifangtu\lena.png');
>> C=rgb2gray(C);
H=C;
C=double(C);
[m1,n1,h1]=size(C);
k=2^(1/3);
threshold=3;
h11=fspecial('gaussian',[5 5],0.3);
C11=imfilter(C,h11,'conv');
h12=fspecial('gaussian',[5 5],0.4);
C12=imfilter(C,h12,'conv');
h13=fspecial('gaussian',[5 5],0.5);
C13=imfilter(C,h13,'conv');
h14=fspecial('gaussian',[5 5],0.6);
C14=imfilter(C,h14,'conv');
h15=fspecial('gaussian',[5 5],0.7);
C15=imfilter(C,h15,'conv');
h16=fspecial('gaussian',[5 5],0.8);
C16=imfilter(C,h16,'conv');
D11=C11-C12;
D12=C12-C13;
D13=C13-C14;
D14=C14-C15;
D15=C15-C16;
E11=zeros(m1,n1);
for i=3:m1-2
for j=3:n1-2
if(D12(i,j)>D12(i-1,j-1) && D12(i,j)>D12(i,j-1) && D12(i,j)>D12(i+1,j-1) && D12(i,j)>D12(i-1,j) && D12(i,j)>D12(i+1,j) && D12(i,j)>D12(i-1,j+1) && D12(i,j)>D12(i,j+1) && D12(i,j)>D12(i+1,j+1))
if(D12(i,j)>D11(i,j) && D12(i,j)>D11(i-1,j-1) && D12(i,j)>D11(i,j-1) && D12(i,j)>D11(i+1,j-1) && D12(i,j)>D11(i-1,j) && D12(i,j)>D11(i+1,j) && D12(i,j)>D11(i-1,j+1) && D12(i,j)>D11(i,j+1) && D12(i,j)>D11(i+1,j+1))
if(D12(i,j)>D13(i,j) && D12(i,j)>D13(i-1,j-1) && D12(i,j)>D13(i,j-1) && D12(i,j)>D13(i+1,j-1) && D12(i,j)>D13(i-1,j) && D12(i,j)>D13(i+1,j) && D12(i,j)>D13(i-1,j+1) && D12(i,j)>D13(i,j+1) && D12(i,j)>D13(i+1,j+1))
if(D12(i,j)>threshold)
E11(i,j)=1;
end
end
end
elseif(D12(i,j)<D12(i-1,j-1) && D12(i,j)<D12(i,j-1) && D12(i,j)<D12(i+1,j-1) && D12(i,j)<D12(i-1,j) && D12(i,j)<D12(i+1,j) && D12(i,j)<D12(i-1,j+1) && D12(i,j)<D12(i,j+1) && D12(i,j)<D12(i+1,j+1))
if(D12(i,j)<D11(i,j) && D12(i,j)<D11(i-1,j-1) && D12(i,j)<D11(i,j-1) && D12(i,j)<D11(i+1,j-1) && D12(i,j)<D11(i-1,j) && D12(i,j)<D11(i+1,j) && D12(i,j)<D11(i-1,j+1) && D12(i,j)<D11(i,j+1) && D12(i,j)<D11(i+1,j+1))
if(D12(i,j)<D13(i,j) && D12(i,j)<D13(i-1,j-1) && D12(i,j)<D13(i,j-1) && D12(i,j)<D13(i+1,j-1) && D12(i,j)<D13(i-1,j) && D12(i,j)<D13(i+1,j) && D12(i,j)<D13(i-1,j+1) && D12(i,j)<D13(i,j+1) && D12(i,j)<D13(i+1,j+1))
if(D12(i,j)<-threshold)
E11(i,j)=-1;
end
end
end
end
end
end
[X,Y]=find(E11==1 | E11==-1);
[M,N]=size(X);
r=10;
for i=1:M
dxx=D12(X(i),Y(i)+1)+D12(X(i),Y(i)-1)-2*D12(X(i),Y(i));
dyy=D12(X(i)+1,Y(i))+D12(X(i)-1,Y(i))-2*D12(X(i),Y(i));
dxy=(D12(X(i)+1,Y(i)+1)-D12(X(i)+1,Y(i)-1)-D12(X(i)-1,Y(i)+1)+D12(X(i)-1,Y(i)-1))/4;
tr=dxx+dyy;
det=dxx*dyy-dxy*dxy;
if(det<=0)
E11(X(i),Y(i))=0;
end
if(tr*tr/det>(r+1/r+1))
E11(X(i),Y(i))=0;
end
end
[X,Y]=find(E11==1 | E11==-1);
[p,q]=size(X);
for i=1:p
dx=(D12(X(i),Y(i)+1)-D12(X(i),Y(i)-1))/2;
dy=(D12(X(i)+1,Y(i))-D12(X(i)-1,Y(i)))/2;
ds=(D13(X(i),Y(i))-D11(X(i),Y(i)))/2;
dD=[dx,dy,ds]';
dxx=D12(X(i),Y(i)+1)+D12(X(i),Y(i)-1)-2*D12(X(i),Y(i));
dyy=D12(X(i)+1,Y(i))+D12(X(i)-1,Y(i))-2*D12(X(i),Y(i));
dxy=(D12(X(i)+1,Y(i)+1)-D12(X(i)+1,Y(i)-1)-D12(X(i)-1,Y(i)+1)+D12(X(i)-1,Y(i)-1))/4;
dss=D13(X(i),Y(i))+D11(X(i),Y(i))-2*D12(X(i),Y(i));
dxs=(D13(X(i),Y(i)+1)-D13(X(i),Y(i)-1)-D11(X(i),Y(i)+1)+D11(X(i),Y(i)-1))/4;
dys=(D13(X(i)+1,Y(i))-D13(X(i)-1,Y(i))-D11(X(i)+1,Y(i))+D11(X(i)-1,Y(i)))/4;
Hiss=[dxx,dxy,dxs;dxy,dyy,dys;dxs,dys,dss];
XX=-inv(Hiss)*dD;
DD=D13(X(i),Y(i))+0.5*dD'*XX;
if(DD<0.03)
E11(X(i),Y(i))=0;
end
end
[X,Y]=find(E11==1 | E11==-1);
rad=1.5*0.4;
[m,n]=size(X);
bin=zeros(1,36);
max_mag=zeros(1,m);
allbin=zeros(1,m);
for t=1:m
for i=-rad:rad
for j=-rad:rad
if(ceil(X(t)+i+1)>m1 || ceil(Y(t)+j+1)>n1)
continue;
end
mag=sqrt((D12(ceil(X(t)+i+1),ceil(Y(t)+j))-D12(ceil(X(t)+i-1),ceil(Y(t)+j)))^2+(D12(ceil(X(t)+i),ceil(Y(t)+j+1))-D12(ceil(X(t)+i),ceil(Y(t)+j-1)))^2);
ori=atan2(D12(ceil(X(t)+i),ceil(Y(t)+j+1))-D12(ceil(X(t)+i),ceil(Y(t)+j-1)),D12(ceil(X(t)+i+1),ceil(Y(t)+j))-D12(ceil(X(t)+i-1),ceil(Y(t)+j)));
w=exp(-(i^2+j^2)/(2*0.4^2));
n=ceil(36*(ori+pi)/(2*pi));
if(n>36);
n=n-36;
end
bin(n)=bin(n)+w*mag;
end
end
[max_mag(t),allbin(t)]=max(bin);
bin=zeros(1,36);
end
[Y,X]=find(E11==1 | E11==-1);
[m,n]=size(X);
imshow(H)
hold on
for t=1:m
len=round(max_mag(t)*5);
hlen=round(max_mag(t)*0.75);
blen=len-hlen;
theta=2*pi*allbin(t)*10/360;
ytermi=round(len*(-sin(theta)))+Y(t);
xtermi=round(len*(cos(theta)))+X(t);
xA=round(blen*(cos(theta-pi/18)))+X(t);
yA=round(blen*(-sin(theta-pi/18)))+Y(t);
xB=round(blen*(cos(theta+pi/18)))+X(t);
yB=round(blen*(-sin(theta+pi/18)))+Y(t);
line([X(t),xtermi],[Y(t),ytermi]);
line([xA,xtermi],[yA,ytermi]);
line([xB,xtermi],[yB,ytermi]);
end
hold off
把特征点也显示出来看看 我换成绿色 绿色好看些
[X22,Y22]=find(E11==1 | E11==-1);
figure
imshow(H)
hold on
plot(Y22,X22,'g.')
我觉得还是红色的'x'星号表示特征点 好看点 这个绿色的菱形。。。有密集恐惧症的不要看 不然会像我一样不敢看第二眼 头皮发麻的赶脚。。。
这是我在大神群里找到的opencv的SIFT代码 关于检测特征点以及生成128维描述子讲解得非常详细 虽然
http://blog.csdn.net/v_JULY_v/article/details/6245939
http://blog.csdn.net/abcjennifer/article/details/7639681
虽然这两位大神都给出了opencv的代码 但是对代码的解释太少了 大部分是讲整个框架的理论以及步骤
而下面这个代码是对代码的每一行 每个变量 每个函数 基本上都有解释
http://download.csdn.net/detail/wd1603926823/8804967
- SIFT检测特征点之找到主方向
- SIFT特征点检测
- SIFT特征点方向确定
- SIFT检测特征点之去除边缘不稳定噪声点
- SIFT检测特征点之生成128维描述子
- SIFT学习--确定特征点的方向
- SIFT算法：确定特征点方向
- SIFT第三步：找到特征点
- 鱼眼拼接之SIFT特征点检测之得到候选特征点
- 目标特征检测之SIFT特征
- SIFT特征中关键点检测技术研究
- 特征点检测学习_1(sift算法)
- 特征点检测学习_1(sift算法)
- 特征点检测学习_1(sift算法)
- 【opencv】 特征点检测 sift和surf
- 特征点检测学习_1(sift算法)
- 特征点检测学习_1(sift算法)
- 特征点检测学习_1(sift算法)
- apache2: Could not reliably determine the server's fully qualified domain name, using 127.0.1.1 for
- 属性参数
- hibernate的“ No suitable driver”错误
- Gson解析复杂的json数据
- 键盘VK键值列表
- SIFT检测特征点之找到主方向
- 在centos6上搭建git服务器
- 【牛腩新闻发布系统 一】如何高效写代码
- 使用PHP制作 简易员工管理系统之七（MVC实现用户信息增、删、改、查）
- KVC
- 开始刷leetcode day40：Set Matrix Zeroes
- 技术 blog20 之 Java正则表达式 by 八戒
- Android之各种Adapter
- 我的Java开发学习之旅------>Java使用Fork/Join框架来并行执行任务
