function MovingTargetDetectionByMMI() %Moving Target Detection%Based on Maximun Mutual Information %%EDIT BY PSL@CSU%QQ：547423688%Email：anyuezhiji@%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:4d = imfilter(d,se);endfor i=1:2d = 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);endfprintf('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.01mi_old=mi;=powell(R,F,pa);endfprintf('------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:mfor 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<n&&y1 style="margin: 0px; word-break: break-all; word-wrap: break-word;">=1&&y1<m)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);endendend%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:mfor 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<n&&y1 style="margin: 0px; word-break: break-all; word-wrap: break-word;">=1&&y1<m)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;endendend%下面求边缘概率密度for i=1:256ha(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:psp=pos(i,: ); %第i个粒子位置out(i)=BI_mi(p(1),p(2),p(3),R,F); %求函数值endpbestval=out; %每个粒子当前函数值=max(pbestval); %全局最优函数值gbest=pbest(idx,: ); %全局极值tr(1)=gbestval; %保存当前全局最优函数值% start PSO iterative procedurescnt=0; % counter used for updating display according to df in the optionscnt2=0; % counter used for the stopping subroutine based on error convergence %开始迭代for i=1:maxinterationsfor j=1:psif flagg==1 % randomization control, one random set for each particle at ea ch epochrannum1=rand(1);rannum2=rand(1);endp=pos(j,: ); %第i个粒子位置out(j)=BI_mi(p(1),p(2),p(3),R,F); %求函数值e(j)=out(j);%更新pbestif pbestval(j)<=e(j);%%%====pbestval(j)=e(j);pbest(j,: )=pos(j,: );end%更新gbest=max(pbestval);if gbestval<=iterbestval %%%===gbestval=iterbestval;gbest=pbest(idx1,: );endtr(i+1)=gbestval;te=i;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%更新速度，位置%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%if i<=iwe %%%===iwt(i)=((iw2-iw1)/(iwe-1))*(i-1)+iw1; %惯性权重elseiwt(i)=iw2;end% iwt(i)=1;%%%%%%%%%%速度%%%%%%%%%%%%%%%%%%%%if flagg==2 % 粒子的每个参数的随机数不同for dimcnt=1: Drannum1=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));endelse % 粒子的每个参数的随机数相同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: Dif vel(j,dimcnt)>mvvel(j,dimcnt)=mv;endif vel(j,dimcnt)<-mvvel(j,dimcnt)=-mv;endif pos(j,dimcnt)>=VR(dimcnt,2)pos(j,dimcnt)=VR(dimcnt,2);endif pos(j,dimcnt)<=VR(dimcnt,1) %%%===pos(j,dimcnt)=VR(dimcnt,1);endendend %结束一次粒子循环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>ergrdcnt2=0;elseif tmp1<=ergrd %%%===cnt2=cnt2+1;if cnt2>=ergrdepbreakendendend %迭代结束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:3D(i,: )=D(i+1,: );endendend %end forpa=x3;mi=fx3;%一维搜索。