目标检测算法

function MovingTargetDetectionByMMI() %Moving Target Detection

%Based on Maximun Mutual Information %

%EDIT BY PSL@CSU

%QQ：547423688

%Email：anyuezhiji@https://www.360docs.net/doc/364941471.html,

%2012-06-01

%读文件

Im1=imread('001.jpg');

Im2=imread('002.jpg');

Im3=imread('003.jpg');

Im1=rgb2gray(Im1);

Im2=rgb2gray(Im2);

Im3=rgb2gray(Im3);

tic;

d12=GetDifferenceImg(Im2,Im1);

d23=GetDifferenceImg(Im2,Im3);

d=d12.*d23;

se =;

for i=1:4

d = imfilter(d,se);

end

for i=1:2

d = medfilt2(d,);

end

%%d=abs((d12-d23).^0.7);

d=uint8(d/max(max(d))*255);

level = graythresh(d);

BW = im2bw(d,level);

s=regionprops(BW,'BoundingBox'); figure(1)

subplot(2,2,1);

imshow(uint8(d12/max(max(d12))*255));

title('参考帧与前一帧的差值')

subplot(2,2,2);

imshow(uint8(d23/max(max(d23))*255));

title('参考帧与后一帧的差值')

subplot(2,2,3);

imshow(BW);

title('由前后帧得出的差值')

subplot(2,2,4);

imshow(Im2);

%imshow(d);

rectangle('Position',s(1).BoundingBox,'Curvature',,'LineWidth',2,'LineStyle','--','EdgeColor', 'r')

title('参考帧与检测结果')

%求相邻两帧重合部分差值主函数

function outImg=GetDifferenceImg(R,F)

=dwt2(R,'db1');

=dwt2(F,'db1');

CA1=uint8(CA1);

CA2=uint8(CA2);

fprintf('\n------PSO start\n');

=PSO(CA1,CA2);

while mi<1.2

=PSO(CA1,CA2);

end

fprintf('tx:%f ty:%f ang:%f mi:%f\n',pa(1),pa(2),pa(3),mi);

fprintf('------PSO end\n\n');

%pa=;

fprintf('------Powell start\n');

mi_old=0;

while abs(mi-mi_old)>0.01

mi_old=mi;

=powell(R,F,pa);

end

fprintf('------Powell end\n\n');

time=toc;

fprintf('tx:%.4f ty:%.4f ang:%.2f mi:%f\n',pa(1),pa(2),pa(3),mi);

fprintf('time:%f\n',time);

outImg=GetDifference(pa(1),pa(2),pa(3),R,F);

%figure(6);imshow(outImg);

%求相邻两帧重合部分差值

function outImg=GetDifference(tx,ty,ang,R,F)

=size(R);

%

R=im2double(R);

F=im2double(F);

theta=ang*pi/180; %旋转角度转弧度

cx=floor(n/2); %旋转的中心点

cy=floor(m/2);

outImg=zeros(m,n);

for j=1:m

for i=1:n

%参考图像在浮动图像平移后的对应点

% x=i-tx; %列

% y=j-ty; %

x=(i-cx)*cos(theta)-(j-cy)*sin(theta)+cx-tx;

y=(i-cx)*sin(theta)+(j-cy)*cos(theta)+cy-ty;

x1=floor(x);

y1=floor(y);

rval=R(j,i);

%图像重合部分求差

if(x1>=1&&x1=1&&y1

dy=y1-y;dx=x1-x;

%双线性插值

fval=(F(y1+1,x1)-F(y1,x1))*dy+(F(y1,x1+1)-F(y1,x1))*dx+(F(y1+1,x1+1)+F (y1,x1)-F(y1,x1+1)-F(y1+1,x1))*dy*dx+F(y1,x1);

outImg(j,i)=abs((rval-fval).^0.7*exp(-min()/20));

%outImg(j,i)=abs((rval-fval).^2.5/(min()).^0.2);

end

end

end

%outImg=uint8(outImg/max(max(outImg))*255);

%双线性插值求互信息

function out=BI_mi(tx,ty,ang,R,F)

=size(R);

hist=zeros(256,256);

ha = zeros(1,256);

hb = zeros(1,256);

%归一化到256级灰度

% if max(max(r))~=min(min(r)) %max(max(a))结果是A中最大的元素，max(A)结果一个行向量，元素分别是A的每个列向量的最大的元素

% r = (r-min(min(r)))/(max(max(r))-min(min(r)));

% else

% r = zeros(M,N);

% end

%

% if max(max(f))-min(min(f))

% f = (f-min(min(f)))/(max(max(f))-min(min(f)));

% else

% f = zeros(M,N);

% end

%

% r = double(int16(r*255))+1;

% f = double(int16(f*255))+1;

R=R+1;

F=F+1;

theta=ang*pi/180; %旋转角度转弧度

cx=floor(n/2); %旋转的中心点

cy=floor(m/2);

%求联合概率密度

for j=1:m

for i=1:n

%参考图像在浮动图像平移后的对应点

% x=i-tx; %列

% y=j-ty; %

x=(i-cx)*cos(theta)-(j-cy)*sin(theta)+cx-tx;

y=(i-cx)*sin(theta)+(j-cy)*cos(theta)+cy-ty;

x1=floor(x);

y1=floor(y);

rval=R(j,i);

%图像重合部分求差

if(x1>=1&&x1=1&&y1

dy=y1-y;dx=x1-x;

%双线性插值

fval=(F(y1+1,x1)-F(y1,x1))*dy+(F(y1,x1+1)-F(y1,x1))*dx+(F(y1+1,x1+1)+F (y1,x1)-F(y1,x1+1)-F(y1+1,x1))*dy*dx+F(y1,x1);

hist(fval,rval)=hist(fval,rval)+1;

end

end

end

%下面求边缘概率密度

for i=1:256

ha(i)=sum(hist(i,: ));

hb(i)=sum(hist(:,i));

end

%调用互信息函数

out=MI(hist,ha,hb);

%下面是求互信息的函数

function mi=MI(hist,ha,hb)

%

hsum = sum(sum(hist));

index = find(hist~=0);

p=hist/hsum;

Hab=sum(sum(-p(index).*log(p(index)))); %联合熵hsum = sum(sum(ha));

index = find(ha~=0);

p = ha/hsum;

Ha = sum(sum(-p(index).*log(p(index)))); %边缘熵hsum = sum(sum(hb));

index = find(hb~=0);

p = hb/hsum;

Hb = sum(sum(-p(index).*log(p(index)))); %边缘熵mi = Ha+Hb-Hab;

%粒子群算法

function =PSO(R,F)

%粒子群算法

%初始化

D=3; %维数

ps=30; %种群规模

VRmin=ones(D,1)*-20; %最小速度

VRmax=ones(D,1)*20; %最大速度

VR=;

%minmax = 1;

pos=40*rand(ps,D)-20; %随机产生初始位置

vel=8*rand(ps,D)-4; %产生随机速度

% ps=15; %种群规模

% VRmin=ones(D,1)*-10; %最小速度

% VRmax=ones(D,1)*10; %最大速度

% VR=;

% %minmax = 1;

% pos=20*rand(ps,D)-10; %随机产生初始位置

% vel=4*rand(ps,D)-2; %产生随机速度

%一些参数

maxinterations=20; %最大迭代次数

iw=1; %固定权重

iw1 = 0.9; % 最大惯性权重

iw2 = 0.4;

iwe = 15;

ac1=2;

ac2=2;

flagg=1;

% ergrd=1e-5;

% ergrdep=5; %

% mv=4;%

ergrd=1e-4;

ergrdep=5; %

mv=4;%

%初始个体极值

pbest=pos;

%求初始全局极值

for i=1:ps

p=pos(i,: ); %第i个粒子位置

out(i)=BI_mi(p(1),p(2),p(3),R,F); %求函数值

end

pbestval=out; %每个粒子当前函数值

=max(pbestval); %全局最优函数值

gbest=pbest(idx,: ); %全局极值

tr(1)=gbestval; %保存当前全局最优函数值% start PSO iterative procedures

cnt=0; % counter used for updating display according to df in the options

cnt2=0; % counter used for the stopping subroutine based on error convergence %开始迭代

for i=1:maxinterations

for j=1:ps

if flagg==1 % randomization control, one random set for each particle at ea ch epoch

rannum1=rand(1);

rannum2=rand(1);

end

p=pos(j,: ); %第i个粒子位置

out(j)=BI_mi(p(1),p(2),p(3),R,F); %求函数值

e(j)=out(j);

%更新pbest

if pbestval(j)<=e(j);%%%====

pbestval(j)=e(j);

pbest(j,: )=pos(j,: );

end

%更新gbest

=max(pbestval);

if gbestval<=iterbestval %%%===

gbestval=iterbestval;

gbest=pbest(idx1,: );

end

tr(i+1)=gbestval;

te=i;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%

%%%%%%更新速度，位置

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%

if i<=iwe %%%===

iwt(i)=((iw2-iw1)/(iwe-1))*(i-1)+iw1; %惯性权重

else

iwt(i)=iw2;

end

% iwt(i)=1;

%%%%%%%%%%速度%%%%%%%%%%%%%%%%%%%%

if flagg==2 % 粒子的每个参数的随机数不同

for dimcnt=1: D

rannum1=rand(1);

rannum2=rand(1);

vel(j,dimcnt)=iwt(i)*vel(j,dimcnt)...

+ac1*rannum1*(pbest(j,dimcnt)-pos(j,dimcnt))...

+ac2*rannum2*(gbest(1,dimcnt)-pos(j,dimcnt));

end

else % 粒子的每个参数的随机数相同

vel(j,: )=iwt(i)*vel(j,: )...

+ac1*rannum1*(pbest(j,: )-pos(j,: ))...

+ac2*rannum2*(gbest(1,: )-pos(j,: ));

end

% 固定权重

% vel(j,: )=iw*vel(j,: )...

% +ac1*rannum1*(pbest(j,: )-pos(j,: ))...

% +ac2*rannum2*(gbest(1,: )-pos(j,: ));

%%%%%%%%%%位置%%%%%%%%%%%%%%%%%%%%%%% pos(j,: )=pos(j,: )+vel(j,: );

%%%%%%%%%速度和位置范围％％％％％％％％％

for dimcnt=1: D

if vel(j,dimcnt)>mv

vel(j,dimcnt)=mv;

end

if vel(j,dimcnt)<-mv

vel(j,dimcnt)=-mv;

end

if pos(j,dimcnt)>=VR(dimcnt,2)

pos(j,dimcnt)=VR(dimcnt,2);

end

if pos(j,dimcnt)<=VR(dimcnt,1) %%%===

pos(j,dimcnt)=VR(dimcnt,1);

end

end

end %结束一次粒子循环

temp=gbest';

fprintf('%f,%f,%f,%f\n',temp(1),temp(2),temp(3),gbestval);

Y(i)=gbestval;

X(i)=i;

% 收敛条件

tmp1=abs(tr(i)-gbestval);

if tmp1>ergrd

cnt2=0;

elseif tmp1<=ergrd %%%===

cnt2=cnt2+1;

if cnt2>=ergrdep

break

end

end

end %迭代结束

fprintf('total interations:%d\n',i);

%OUT=;

pa=gbest;

mi=gbestval;

%POWELL优化算法

function =powell(R,F,x)

len=5; %搜索区间

itmax=30; %最大循环次数

e=1e-3; %允许误差

%方向矢量矩阵存放d1,d2,d3三个方向矢量

D=;

%起始点

%x0=;

x0=x;

fx0=BI_mi(x0(1),x0(2),x0(3),R,F);

%fx0=pv_mi(x0(1),x0(2),-x0(3),R,F);

%循环搜索

for k=0:itmax

%从起始点出发，沿d1方向搜索，得到最大值fx1,对应点x1 d1=D(1,: );

=oneDimSearch(R,F,x0,d1,len);

fprintf('%f, %f, %f, %f\n',x1(1),x1(2),x1(3),fx1);

%从x1出发，沿d2方向搜索

d2=D(2,: );

=oneDimSearch(R,F,x1,d2,len);

fprintf('%f, %f, %f, %f\n',x2(1),x2(2),x2(3),fx2);

%从x2出发，沿d3方向搜索

d3=D(3,: );

=oneDimSearch(R,F,x2,d3,len);

%加速方向

dn=x3-x0;

y=sum(dn.*dn);

fprintf('%f\n',y);

fprintf('%f, %f, %f, %f\n',x3(1),x3(2),x3(3),fx3);

if y<=e %满足结束条件，退出循环%%%===

pa=x3;

mi=fx3;

return;

end

%调整搜索方向

%计算目标值下降最多的那个方向

cha=;

=max(cha);

%

xe=2*x3-x0;

fe=BI_mi(xe(1),xe(2),xe(3),R,F);

%fe=pv_mi(xe(1),xe(2),-xe(3),R,F);

%这里求极大值

if (fe<=fx0) || (2*(fx0-2*fx3+fe)*(fx0-fx3-maxcha)*(fx0-fx3-maxcha) >= (fx0-fe)* (fx0-fe)*maxcha)

%不引入新的方向%%%===

x0=x3; %下次搜索的新起点

fx0=fx3;

else %引进新的方向

%以xn为起点沿dn进行搜索，得到下次搜索的新起点

=oneDimSearch(R,F,x3,dn,len);

%替换方向

D(4,: )=dn;

for i=j0:3

D(i,: )=D(i+1,: );

end

end

end %end for

pa=x3;

mi=fx3;

%一维搜索。从-len到len，找出最大值的点

function =oneDimSearch(R,F,X,direction,len)

%一维brent搜索求函数最小值

CGOLD =0.3819660;

a=-len;

b=len;

%a=a0;

%b=b0;

tol=0.01;

e=0.0;

ITMAX=100;

v=a+CGOLD*(b-a);

Xv=X+direction*v;

%fv=pv_mi(Xv(1),Xv(2),-Xv(3),R,F);

fv=BI_mi(Xv(1),Xv(2),Xv(3),R,F);

w=v;

x=v;

fw=fv;

fx=fv;

for k=0:ITMAX

xm=0.5*(a+b); %中点

tol1=tol*abs(x);

tol2=2*tol1;

%if abs(x-xm)=tol;

if abs(x-xm)<=(tol2-0.5*(b-a)) %满足精度，退出%%%=== step=x;

Y=X+direction*step;

fY=fx;

break;

end

if (abs(e) > tol1)

r=(x-w)*(fx-fv);

q=(x-v)*(fx-fw);

p=(x-v)*q-(x-w)*r; %分子

q=2*(q-r); %分母

if (q>0)

p=-p;

end

q=abs(q);

e=d;

if (abs(p) >= abs(0.5*q*etemp) || p<=q*(a-x) || p >= q*(b-x)) %使用黄金分割%%%===

if x >= xm

e=a-x;

else

e=b-x;

end

d=CGOLD*e;

else %抛物线内插

d=p/q;

u=x+d; %u

if(u-a if x > xm;

d=-abs(tol1);

else

d=abs(tol1);

end

end

end

else %使用黄金分割

if x >= xm

e=a-x;

else

e=b-x;

end

d=CGOLD*e;

end %if (abs(e) > tol1)

if abs(d)>=tol1

u=x+d;

else

u=x+abs(tol1);

else

u=x-abs(tol1);

end

end

Xu=X+direction*u;

% fu=pv_mi(Xu(1),Xu(2),-Xu(3),R,F); %函数值 fu=BI_mi(Xu(1),Xu(2),Xu(3),R,F); %函数值

%更新参数

if fu>=fx

if u>=x

a=x;

else

b=x;

end

v=w;w=x;x=u;

fv=fw;fw=fx;fx=fu;

else

if u

a=u;

else

b=u;

end

if (fu>=fw || w==x)

v=w;w=u;

fv=fw;fw=fu;

elseif (fu>=fv || v==x || v==w)

v=u;

fv=fu;

end

end

end %end for

用于检测的原始序列如下：

前一帧 001.jpg

参考帧002.jpg

后一帧003.jpg

运行结果如下：

------PSO start

6.219116,3.410710,-4.782166,1.132898 3.181979,1.134846,-3.718395,1.227941 2.302928,-0.910530,-2.792701,1.258413 4.141745,1.269771,-0.407340,1.409782 4.141745,1.269771,-0.407340,1.409782 2.782469,0.226964,-0.933429,1.415868

2.782469,0.226964,-0.933429,1.415868

2.782469,0.226964,-0.933429,1.415868

2.782469,0.226964,-0.933429,1.415868

2.782469,0.226964,-0.933429,1.415868

3.591808,0.754384,-0.656160,1.481511

3.591808,0.754384,-0.656160,1.481511

3.473467,0.846358,-0.247672,1.483317

3.810920,0.631524,-0.892638,1.489102

3.857352,0.603985,-0.596614,1.490289

3.857352,0.603985,-0.596614,1.490289

3.862678,0.614558,-0.632967,1.493891

3.862678,0.614558,-0.632967,1.493891

3.862678,0.614558,-0.632967,1.493891

3.862678,0.614558,-0.632967,1.493891

total interations:20

tx:3.862678 ty:0.614558 ang:-0.632967 mi:1.493891 ------PSO end

------Powell start

7.719056, 0.614558, -0.632967, 1.298058

7.719056, 1.904456, -0.632967, 1.466426

16.537888

7.719056, 1.904456, -0.681950, 1.463247

8.229286, 2.035531, -0.686927, 1.481208

8.229286, 1.980443, -0.686927, 1.482608

0.015995

8.223876, 1.978633, -0.686858, 1.482462

8.057189, 1.978633, -0.686858, 1.482155

8.057189, 1.967980, -0.686858, 1.481889

0.027765

8.057581, 1.968111, -0.686863, 1.481945

8.060448, 1.968111, -0.686863, 1.481762

8.060448, 1.967566, -0.686863, 1.481860

0.000009

8.060564, 1.967605, -0.686865, 1.481848

8.056211, 1.967605, -0.686865, 1.481832

8.056211, 1.980282, -0.686865, 1.482319

0.000180

8.056211, 1.980282, -0.687108, 1.482356

------Powell end

tx:8.0562 ty:1.9803 ang:-0.69 mi:1.482356

time:5.523760

------PSO start

-4.218613,-1.831505,-2.592978,1.198497

-3.096033,0.220999,-1.548051,1.355306

-3.096033,0.220999,-1.548051,1.355306

-3.788383,0.085873,-1.170346,1.356069

-2.132791,0.115566,-1.059226,1.400370

-2.536139,-0.302278,-0.703993,1.402590

-3.122589,0.042784,-0.841495,1.411028

-2.294945,0.111933,-0.833426,1.440319

-2.601845,0.432572,-0.042504,1.558405

-2.601845,0.432572,-0.042504,1.558405

-2.601845,0.432572,-0.042504,1.558405

-2.601845,0.432572,-0.042504,1.558405

-2.623005,0.577600,0.167668,1.558405

-2.623005,0.577600,0.167668,1.558405

total interations:14

tx:-2.623005 ty:0.577600 ang:0.167668 mi:1.558405 ------PSO end

------Powell start

-5.479949, 0.577600, 0.167668, 1.500895

-5.479949, 1.084894, 0.167668, 1.552012

8.420517

-5.479949, 1.084894, 0.135345, 1.551926

-4.800852, 1.084894, 0.135345, 1.625440 -4.800852, 1.146578, 0.135345, 1.639162 0.464979

-4.800852, 1.146578, 0.134208, 1.638934 -4.657180, 1.145903, 0.134221, 1.638934 -4.657180, 1.144253, 0.134221, 1.638942 0.018618

-4.671874, 1.142918, 0.134245, 1.639685 -4.679493, 1.142918, 0.134245, 1.639685 -4.679493, 1.144981, 0.134245, 1.638942 0.000028

-4.676647, 1.145240, 0.134241, 1.639855 -4.670146, 1.145240, 0.134241, 1.639855 -4.670146, 1.145272, 0.134241, 1.639788 0.000042

-4.670146, 1.145272, 0.134258, 1.639799 ------Powell end

tx:-4.6701 ty:1.1453 ang:0.13 mi:1.639799 time:10.290607