【运动检测】基于最大互信息运动目标检测含Matlab源码

最新推荐文章于 2022-10-11 22:11:46 发布

原创最新推荐文章于 2022-10-11 22:11:46 发布 · 292 阅读

3 ·

CC 4.0 BY-SA版权

文章标签：

#matlab #图像处理 #计算机视觉

图像处理专栏收录该内容

871 篇文章

订阅专栏

本文介绍了一种基于MATLAB的运动目标检测方法，通过改进的帧差算法、自适应阈值分割、滤波去噪，并结合逻辑运算和形态学处理，实现对运动物体轮廓的高效提取。关键步骤包括PSO优化算法的应用，展示了完整代码和仿真结果，以及相关参考文献。

1 简介

基于MATLAB工具对运动目标的检测进行了研究,利用改进的帧差算法对获取的图像做帧差并将帧差图像做自适应阈值分割,然后自适应滤波去燥,最后利用逻辑与运算将两幅二值化图像合并成一幅图像,利用形态学腐蚀和膨胀的检测方法,提取运动物体的轮廓达到检测运动目标的目的.

2 完整代码

function  MovingTargetDetectionByMMI() %Moving Target Detection%Based on Maximun Mutual Information%%读文件  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 =[0     0     0     1     0     0     0     0     0     1     2     1     0     0     0     1     2     4     2     1     0     1     2     4     8     4     2     1     0     1     2     4     2     1     0     0     0     1     2     1     0     0     0     0     0     1     0     0     0];for i=1:4d = imfilter(d,se);endfor i=1:2d = medfilt2(d,[4 4]);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',[0,0],'LineWidth',2,'LineStyle','--','EdgeColor','r')  title('参考帧与检测结果') %求相邻两帧重合部分差值主函数function outImg=GetDifferenceImg(R,F)[CA1,~,~,CD1]=dwt2(R,'db1'); [CA2,~,~,CD2]=dwt2(F,'db1'); CA1=uint8(CA1); CA2=uint8(CA2); fprintf('\n------PSO start\n'); [pa,mi]=PSO(CA1,CA2); while mi<1.2    [pa,mi]=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=[0,0,0]; fprintf('------Powell start\n'); mi_old=0;while abs(mi-mi_old)>0.01    mi_old=mi;    [pa,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) [m,n]=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<n&&y1>=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([rval fval])/20));            %outImg(j,i)=abs((rval-fval).^2.5/(min([rval fval])).^0.2);        end    end end %outImg=uint8(outImg/max(max(outImg))*255);%双线性插值求互信息 function out=BI_mi(tx,ty,ang,R,F) [m,n]=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<n&&y1>=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;        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 [pa,mi]=PSO(R,F) %粒子群算法  %初始化 D=3;                  %维数 ps=30;                %种群规模 VRmin=ones(D,1)*-20;  %最小速度 VRmax=ones(D,1)*20;   %最大速度 VR=[VRmin,VRmax]; %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=[VRmin,VRmax]; % %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;                        %每个粒子当前函数值 [gbestval,idx]=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 each 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         [iterbestval,idx1]=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=[gbest';gbestval]; pa=gbest; mi=gbestval; %POWELL优化算法 function [pa,mi]=powell(R,F,x) len=5;   %搜索区间 itmax=30; %最大循环次数 e=1e-3;   %允许误差  %方向矢量矩阵存放d1,d2,d3三个方向矢量 D=[1 0 0;     0 1 0;     0 0 1]; %起始点 %x0=[-10,-12,0]; 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,:);     [x1,fx1,step]=oneDimSearch(R,F,x0,d1,len);     fprintf('%f,  %f,  %f,  %f\n',x1(1),x1(2),x1(3),fx1);     %从x1出发，沿d2方向搜索     d2=D(2,:);       [x2,fx2,step]=oneDimSearch(R,F,x1,d2,len);     fprintf('%f,  %f,  %f,  %f\n',x2(1),x2(2),x2(3),fx2);     %从x2出发，沿d3方向搜索     d3=D(3,:);     [x3,fx3,step]=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=[fx0-fx1,fx1-fx2,fx2-fx3];     [maxcha,j0]=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进行搜索，得到下次搜索的新起点     [x0,fx0,step]=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 [Y,fY,step]=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);         etemp=e;         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<tol2 || b-u<tol2)  %u靠近边界 %%%===                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         if d>=0             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<x         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