SIFT检测特征点之生成128维描述子

接着上一篇 现在要对每个特征点生成128维描述符 这一步的理论还是来自于论文以及这位大神http://blog.csdn.net/v_JULY_v/article/details/6245939 总结下论文上和这位大神所说：先对每个特征点取8x8邻域 然后把这8x8区域又分成四个4x4区域 每个4x4区域生成一个有8个方向向量信息的种子 这样对一个特征点就有4x8=32维向量信息了 但是发明SIFT算法的这个外国爷爷建议我们不这样 建议我们把8x8区域分成16个2x2区域 每个2x2区域生成一个8个方向向量的种子
那么一个特征点就有了16x8=128维向量信息了 就是我题目所说生成了128维描述子

但是我还是按照简单的来吧 我们实际中没必要那么遵守规则 还是按照那位大神的每4x4区域生成一个种子 对每个特征点生成32维向量就可以了

接着上一篇的程序继续写：

[Y,X]=find(E11==1| E11==-1);
[m,n]=size(X);
allvector1=zeros(m,32); %表示m个特征点 每个特征点有32维信息
binnew1=zeros(1,8); %第一个4x4区域生成的种子 也就是8个方向向量信息
binnew2=zeros(1,8); %第二个4x4区域生成的8维向量信息
binnew3=zeros(1,8); %第三个4x4区域生成的8维向量信息
binnew4=zeros(1,8); %第四个4x4区域生成的8维向量信息
Xnew=zeros(1,m); %将坐标轴旋转到特征点主方向后的坐标
Ynew=zeros(1,m);
%将坐标旋转到主方向
for t=1:m
arfa=2*pi*allbin(t)*10/360;
Xnew(t)=X(t)*cos(arfa)-Y(t)*sin(arfa);
Ynew(t)=X(t)*sin(arfa)+Y(t)*cos(arfa);
if(Xnew(t)<5)
Xnew(t)=6;
end
if(Ynew(t)<5)
Ynew(t)=6;
end
if(Xnew(t)>(m1-5))
Xnew(t)=m1-6;
end
if(Ynew(t)>n1-5)
Ynew(t)=n1-6;
end
end
%下面对m个特征点生成32维描述子
for t=1:m
for i=-4:-1
for j=1:4
mag=sqrt((D12(ceil(Xnew(t)+i+1),ceil(Ynew(t)+j))-D12(ceil(Xnew(t)+i-1),ceil(Ynew(t)+j)))^2+(D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j+1))-D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j-1)))^2);
ori=atan2(D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j+1))-D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j-1)),D12(ceil(Xnew(t)+i+1),ceil(Ynew(t)+j))-D12(ceil(Xnew(t)+i-1),ceil(Ynew(t)+j)));
w=exp(-(i^2+j^2)/(2*0.4^2));
n=ceil(8*(ori+pi)/(2*pi));
if(n>8);
n=n-8;
end
binnew1(n)=binnew1(n)+w*mag;
end
end
for i=1:4
for j=1:4
mag=sqrt((D12(ceil(Xnew(t)+i+1),ceil(Ynew(t)+j))-D12(ceil(Xnew(t)+i-1),ceil(Ynew(t)+j)))^2+(D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j+1))-D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j-1)))^2);
ori=atan2(D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j+1))-D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j-1)),D12(ceil(Xnew(t)+i+1),ceil(Ynew(t)+j))-D12(ceil(Xnew(t)+i-1),ceil(Ynew(t)+j)));
w=exp(-(i^2+j^2)/(2*0.4^2));
n=ceil(8*(ori+pi)/(2*pi));
if(n>8);
n=n-8;
end
binnew2(n)=binnew2(n)+w*mag;
end
end
for i=-4:-1
for j=-4:-1
mag=sqrt((D12(ceil(Xnew(t)+i+1),ceil(Ynew(t)+j))-D12(ceil(Xnew(t)+i-1),ceil(Ynew(t)+j)))^2+(D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j+1))-D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j-1)))^2);
ori=atan2(D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j+1))-D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j-1)),D12(ceil(Xnew(t)+i+1),ceil(Ynew(t)+j))-D12(ceil(Xnew(t)+i-1),ceil(Ynew(t)+j)));
w=exp(-(i^2+j^2)/(2*0.4^2));
n=ceil(8*(ori+pi)/(2*pi));
if(n>8);
n=n-8;
end
binnew3(n)=binnew3(n)+w*mag;
end
end
for i=1:4
for j=-4:-1
mag=sqrt((D12(ceil(Xnew(t)+i+1),ceil(Ynew(t)+j))-D12(ceil(Xnew(t)+i-1),ceil(Ynew(t)+j)))^2+(D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j+1))-D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j-1)))^2);
ori=atan2(D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j+1))-D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j-1)),D12(ceil(Xnew(t)+i+1),ceil(Ynew(t)+j))-D12(ceil(Xnew(t)+i-1),ceil(Ynew(t)+j)));
w=exp(-(i^2+j^2)/(2*0.4^2));
n=ceil(8*(ori+pi)/(2*pi));
if(n>8);
n=n-8;
end
binnew4(n)=binnew4(n)+w*mag;
end
end
allvector1(t,:)=[binnew1 binnew2 binnew3binnew4];
binnew1=zeros(1,8);
binnew2=zeros(1,8);
binnew3=zeros(1,8);
binnew4=zeros(1,8);
end
到此 对第一幅图像处理已经结束 因为我是用SIFT进行鱼眼图像的拼接 所以我先对第一幅鱼眼图像进行提取有效区域 然后校正 再进行SIFT的高斯卷积 然后就是SIFT的DOG差分空间 然后就是得到极值 再就是去除两种特征点 然后再生成梯度直方图 最后就是生成128维（理论上）描述子 这样对第一幅鱼眼图像的处理已经结束

因为我要拼接的鱼眼图像有三幅：s1，s2和s3






就是这三幅鱼眼图像
顺序我好像发反了

所以我才说CSDN发图好麻烦 不说了 说多了都是泪 反正就是第三幅应该拼接以后位于最坐标 然后第一幅图在中间 然后第二幅图拼接后位于最右边 如果拼接程序编出来 应该就是这样的效果

这三幅是原图 然后进行校正 校正后的三幅图如下 原谅我没有用3D的校正方法 也没有进行插值 因为目前球面透视投影校正和经纬度坐标校正我编出来显示出来简直亮瞎了我的狗眼 不忍直视 所以等我弄完拼接再回去慢慢编3D的校正方法 之前的插值方法并不适合所有的图像 所以我没有插值







接下来我把这三幅图SIFT检测后的效果图也贴出来：我实在受不了之前的红色的星号表示的特征点了 我觉得我之前脑袋被雷劈了 我还是把它们换成绿色的菱形吧 好看些







下面每个特征点主方向画在图上：







为了下一步的拼接方便 我把第一幅鱼眼图生成的32维描述子放在allvector1里面 把第二幅的描述子放在allvector2里 把第三幅的描述子放在allvector3里 因为代码都一样的 只是图换了而已 所以就不贴代码了 参照前面就可以了。

为了方便 我把整个SIFT检测特征点编程M文件 这样不用每次都写长长的代码了 直接调用就好了

第一幅：

A1=imread('F:\orl_zhifangtu\s1.jpg');

T=40;

[A1,R1]=kuaisusaomiao(A1,T); %快速扫描算法提取鱼眼有效区域

C1=jingduzuobiao(A1,T,R1); %鱼眼经度坐标校正

[E11,H1,D12,D11,D13]=primerpoints(C1); %得到候选特征点

XY1=last_points(E11,D12,D11,D13); %去除两种不要的特征点 即精确特征点

allbin1=tiduzhifangtu(D12,XY1,H1); %对每个特征点生成梯度直方图 找到主方向

allvector1=miaoshuzi(allbin1,XY1,D12,H1); %生成128维描述子 并归一化

[jj,kk]=size(XY1); %如果想显示特征点出来 可以加上这一小段 如果想显示每个特征点的方向幅值也可以 之前已写过 这里就不重写了

>> X=[];

>> Y=[];

>> for t=1:jj

X=[X XY1(t,1)];

Y=[Y XY1(t,2)];

end

>> imshow(H1)

>> hold on

>> plot(X,Y,'g.')

第二幅：

>> A2=imread('F:\orl_zhifangtu\s2.jpg');

T=40;

[A2,R2]=kuaisusaomiao(A2,T);

C2=jingduzuobiao(A2,T,R2);

[E22,H2,D22,D21,D23]=primerpoints(C2);

XY2=last_points(E22,D22,D21,D23);

allbin2=tiduzhifangtu(D22,XY2,H2);

allvector2=miaoshuzi(allbin2,XY2,D22,H2);

[jj,kk]=size(XY2);

>> X=[];

>> Y=[];

>> for t=1:jj

X=[X XY2(t,1)];

Y=[Y XY2(t,2)];

end

>> figure

>> imshow(H2)

>> hold on

>> plot(X,Y,'g.')

第三幅：

>> A3=imread('F:\orl_zhifangtu\s3.jpg');

T=40;

[A3,R3]=kuaisusaomiao(A3,T);

C3=jingduzuobiao(A3,T,R3);

[E33,H3,D32,D31,D33]=primerpoints(C3);

XY3=last_points(E33,D32,D31,D33);

allbin3=tiduzhifangtu(D32,XY3,H3);

allvector3=miaoshuzi(allbin3,XY3,D32,H3);

[jj,kk]=size(XY3);

>> X=[];

>> Y=[];

>> for t=1:jj

X=[X XY3(t,1)];

Y=[Y XY3(t,2)];

end

>> figure

>> imshow(H3)

hold on

>> plot(X,Y,'g.')

下面是各个函数：

第一个函数：

function [B,x,y,R]=kuaisusaomiao(A,T)

% The function is meant to get yuyan picture from a primer picture.It seems

%like everyline_everyrow_sweep but have a little difference.A is the primer

%picture and T is the yuzhi .B is the output picture we get.

[m,n,k]=size(A);

for i=1:m-1

for j=1:n

I(i,j)=0.59*A(i,j,1)+0.11*A(i,j,2)+0.3*A(i,j,3);

if(I(i,j)>=T)

I(i+1,j)=0.59*A(i+1,j,1)+0.11*A(i+1,j,2)+0.3*A(i+1,j,3);

if(I(i+1,j)>=T)

bottom=j;

break;

end

end

end

end

for i=m:-1:2

for j=1:n

I(i,j)=0.59*A(i,j,1)+0.11*A(i,j,2)+0.3*A(i,j,3);

if(I(i,j)>=T)

I(i-1,j)=0.59*A(i-1,j,1)+0.11*A(i-1,j,2)+0.3*A(i-1,j,3);

if(I(i-1,j)>=T)

top=i;

break;

end

end

end

end

for j=1:n-1

for i=1:m

I(i,j)=0.59*A(i,j,1)+0.11*A(i,j,2)+0.3*A(i,j,3);

if(I(i,j)>=T)

I(i,j+1)=0.59*A(i,j+1,1)+0.11*A(i,j+1,2)+0.3*A(i,j+1,3);

if(I(i,j+1)>=T)

right=j;

break;

end

end

end

end

for j=n:-1:2

for i=1:m

I(i,j)=0.59*A(i,j,1)+0.11*A(i,j,2)+0.3*A(i,j,3);

if(I(i,j)>=T)

I(i,j-1)=0.59*A(i,j-1,1)+0.11*A(i,j-1,2)+0.3*A(i,j-1,3);

if(I(i,j-1)>=T)

left=j;

break;

end

end

end

end

x=(left+right)/2;

y=(top+bottom)/2;

R=max(abs((right-left)/2),abs((bottom-top)/2));

B=imcrop(A,[left,top,2*R,2*R]);

第二个：

function C=jingduzuobiao(A,T,R)

%The function is meant to correct the wrong yuyan picture.A is the primer

%wrong yuyan picture,[x y] is the center of the A. B is the correct yuyan

%picture we get.The function's name is jingduzuobiaojiaozheng suanfa.

[m,n,k]=size(A);

C=[];

x=n/2;

y=m/2;

for u=1:m

for v=1:n

i=u;

j=round(sqrt(R^2-(y-u)^2)*(v-x)/R+x);

if(R^2-(y-u)^2<0)

continue;

end

C(u,v,1)=A(i,j,1);

C(u,v,2)=A(i,j,2);

C(u,v,3)=A(i,j,3);

end

end

C=uint8(C);

第三个：

function [E11,H1,D12,D11,D13]=primerpoints(C)

%The function is meant to get the primer character points of SIFT.

C=rgb2gray(C);

H1=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

第四个：

function XY1=last_points(E11,D12,D11,D13)

%The function is get the newest and the most correct character points of

%SIFT.E11 (in the right) is the primer points and E11(in the right) is the

%last points.D12,D11,D13 are the DOG space_scales.

[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

[Y,X]=find(E11==1 | E11==-1);

[m11,n11]=size(Y);

XY1=[m11,2];

for tt=1:m11

XY1= [XY1;X(tt),Y(tt)];

end

第五个：

function allbin=tiduzhifangtu(D12,XY1,H1)

%The function is meant to get zhifangtu of every character point of XY1.

X=XY1(:,1);

Y=XY1(:,2);

rad=1.5*0.4;

[m,n]=size(X);

bin=zeros(1,36);

max_mag=zeros(1,m);

allbin=zeros(1,m);

[m1,n1]=size(H1);

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

第六个：

function allvector1=miaoshuzi(allbin,XY1,D12,H1)

%The function is meant to get miaoshuzi of every character point ,then get

%the whole miaoshuzi of a picture.

[m1,n1]=size(H1);

X=XY1(:,1);

Y=XY1(:,2);

[m,n]=size(X);

allvector1=zeros(m,32);

binnew1=zeros(1,8);

binnew2=zeros(1,8);

binnew3=zeros(1,8);

binnew4=zeros(1,8);

Xnew=zeros(1,m);

Ynew=zeros(1,m);

for t=1:m

arfa=2*pi*allbin(t)*10/360;

Xnew(t)=X(t)*cos(arfa)-Y(t)*sin(arfa);

Ynew(t)=X(t)*sin(arfa)+Y(t)*cos(arfa);

if(Xnew(t)<5)

Xnew(t)=6;

end

if(Ynew(t)<5)

Ynew(t)=6;

end

if(Xnew(t)>(m1-5))

Xnew(t)=m1-6;

end

if(Ynew(t)>n1-5)

Ynew(t)=n1-6;

end

end

for t=1:m

for i=-4:-1

for j=1:4

mag=sqrt((D12(ceil(Xnew(t)+i+1),ceil(Ynew(t)+j))-D12(ceil(Xnew(t)+i-1),ceil(Ynew(t)+j)))^2+(D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j+1))-D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j-1)))^2);

ori=atan2(D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j+1))-D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j-1)),D12(ceil(Xnew(t)+i+1),ceil(Ynew(t)+j))-D12(ceil(Xnew(t)+i-1),ceil(Ynew(t)+j)));

w=exp(-(i^2+j^2)/(2*0.4^2));

n=ceil(8*(ori+pi)/(2*pi));

if(n>8);

n=n-8;

end

binnew1(n)=binnew1(n)+w*mag;

end

end

for i=1:4

for j=1:4

mag=sqrt((D12(ceil(Xnew(t)+i+1),ceil(Ynew(t)+j))-D12(ceil(Xnew(t)+i-1),ceil(Ynew(t)+j)))^2+(D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j+1))-D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j-1)))^2);

ori=atan2(D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j+1))-D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j-1)),D12(ceil(Xnew(t)+i+1),ceil(Ynew(t)+j))-D12(ceil(Xnew(t)+i-1),ceil(Ynew(t)+j)));

w=exp(-(i^2+j^2)/(2*0.4^2));

n=ceil(8*(ori+pi)/(2*pi));

if(n>8);

n=n-8;

end

binnew2(n)=binnew2(n)+w*mag;

end

end

for i=-4:-1

for j=-4:-1

mag=sqrt((D12(ceil(Xnew(t)+i+1),ceil(Ynew(t)+j))-D12(ceil(Xnew(t)+i-1),ceil(Ynew(t)+j)))^2+(D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j+1))-D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j-1)))^2);

ori=atan2(D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j+1))-D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j-1)),D12(ceil(Xnew(t)+i+1),ceil(Ynew(t)+j))-D12(ceil(Xnew(t)+i-1),ceil(Ynew(t)+j)));

w=exp(-(i^2+j^2)/(2*0.4^2));

n=ceil(8*(ori+pi)/(2*pi));

if(n>8);

n=n-8;

end

binnew3(n)=binnew3(n)+w*mag;

end

end

for i=1:4

for j=-4:-1

mag=sqrt((D12(ceil(Xnew(t)+i+1),ceil(Ynew(t)+j))-D12(ceil(Xnew(t)+i-1),ceil(Ynew(t)+j)))^2+(D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j+1))-D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j-1)))^2);

ori=atan2(D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j+1))-D12(ceil(Xnew(t)+i),ceil(Ynew(t)+j-1)),D12(ceil(Xnew(t)+i+1),ceil(Ynew(t)+j))-D12(ceil(Xnew(t)+i-1),ceil(Ynew(t)+j)));

w=exp(-(i^2+j^2)/(2*0.4^2));

n=ceil(8*(ori+pi)/(2*pi));

if(n>8);

n=n-8;

end

binnew4(n)=binnew4(n)+w*mag;

end

end

allvector1(t,:)=[binnew1 binnew2 binnew3 binnew4];

binnew1=zeros(1,8);

binnew2=zeros(1,8);

binnew3=zeros(1,8);

binnew4=zeros(1,8);

end

[m1,n1]=size(allvector1);

sum1=0;

for t=1:m1

for k=1:32

sum1=sum1+allvector1(t,k);

end

for k=1:32

allvector1(t,k)=allvector1(t,k)/sum1;

end

end