基于颜色直方图的粒子滤波目标跟踪MATLAB实现

来源:互联网 发布:mysql join on关键词 编辑:程序博客网 时间:2024/05/19 13:22

       视频运动目标估计是一个非线性、非高斯的过程,而粒子滤波目标跟踪由于其对于非线性和非高斯的独特特点,已被广泛应用于视频目标跟踪领域。下面主要是基于颜色直方图的粒子滤波目标跟踪结果和程序,颜色通道为hsv通道,量化bin为16 16 16。先上结果:



下面是代码:

      

function main()mov=VideoReader('00002.avi');    % 更换视频目录即可Startframe=300;% -----------------------------------粒子滤波变量初始化-------------------------------------N=99;num=2;                      % num为要跟踪的目标个数  , 粒子个数delta=0.7;                      % 状态转移矩阵A = [1 delta 0 0 0 0  ; 0 1 0 0 0 0  ; 0 0 1 delta 0 0 ; 0 0 0 1 0 0  ; 0 0 0 0 1 0  ; 0 0 0 0 0 1];B = [delta^2/2 0 0 0 0 0;0 delta 0 0 0 0;0 0 delta^2/2 0 0 0;0 0 0 delta 0 0;0 0 0 0 1 0;0 0 0 0 0 1];% img1=read(mov,Startframe);% 手动选取目标,并输出目标框的中心坐标以及包含目标框宽高的变量   图1变量初始化[temp1,rect1,center1,loc1]=Get_Target1(img1,Startframe,num);            % loc1存放的是图1中的目标框的中心点坐标以及框的宽高[X01,XX1,XX01,XX0_temp1,x_temp1,hist1]=state_initial1(loc1,temp1,N,num);% 粒子状态初始化%% 读取序列图像% 初始化绘制结果的矩阵Position1=cell(1,num);Estmation_x1=cell(1,num);Estmation_y1=cell(1,num);w1=cell(1,num);for i=1:num    % 初始化权值矩阵    w1{1,i}=ones(1,N);endsave_hist1=cell(1,3);        % 这个数组为后面处理遮挡后更新hist用for frame=(Startframe+1):(Startframe+500)    Im1 =read(mov, frame);    for k=1:num        % 跟踪处理        [XX1{1,k},value1(k),w1{1,k},Estmation_x1{1,k},Estmation_y1{1,k},loc1{1,k}]=TrackingProcessing1(Im1,hist1{1,k},A,B,XX1{1,k},loc1{1,k},N,w1{1,k});        %状态更新        [hist1{1,k},XX01{1,k},x_temp1{1,k},Position1{1,k},save_hist1{1,k}]=UPdate1(Im1,value1(k),Estmation_x1{1,k},Estmation_y1{1,k},loc1{1,k},rect1{1,k},frame,Startframe,XX01{1,k},x_temp1{1,k},hist1{1,k},save_hist1{1,k});        %Resample        [XX1{1,k},w1{1,k}]=resample_particles1(XX1{1,k},w1{1,k},N);        % 重采样后粒子权重归一化        %   w1{1,k}=ones(1,N);    end    % 显示结果    figure(1);    imshow(Im1);set(gca, 'position', [0 0 1 1 ]);axis normal;    DrawResult1(num,Position1,Estmation_x1,Estmation_y1,XX1,frame+1);    title('Camera1');    AFrame=getframe(gcf);    %存储调整过大小的图片 imwrite(AFrame.cdata,strcat('F:\粒子滤波相关\HSV求直方图的\视频的\HSV求直方图的\2017.5.18备份\1\',num2str(frame),'.jpg'))end    function [temp,rect,center,location]=Get_Target1(rgb1,Startframe,num)        %------这个函数文件是针对于读取图片的-------        % 手动选取目标,并输出目标框的中心坐标以及包含目标框        % 宽高的变量                temp=cell(1,num);        rect=cell(1,num);        center=cell(1,num);        location=cell(1,num);        % is=num2str(Startframe);        % number = '0000';        % number(end-length(is)+1:end)=is;        % filename11=[ file_path1 'frame_' number '.jpg'];        % rgb1 = imread(filename11);                for kk=1:num            figure(1),            imshow(rgb1);            [temp{1,kk},rect{1,kk}]=imcrop(rgb1);            % 得到目标框中西的坐标            center{1,kk}=[rect{1,kk}(1)+rect{1,kk}(3)/2,rect{1,kk}(2)+rect{1,kk}(4)/2];            location{1,kk}=[center{1,kk}(1),center{1,kk}(2),rect{1,kk}(3),rect{1,kk}(4)];        end        close(gcf);    end    function [X0,XX,XX0,XX0_temp,x_temp,hist1]=state_initial1(location,temp,N,num)        % 粒子状态初始化        X0=cell(1,num);        XX=cell(1,num);        XX0=cell(1,num);        XX0_temp=cell(1,num);        x_temp=cell(1,num);        hist1=cell(1,num);        for ii=1:num            X0{1,ii}=[location{1,ii}(1),0,location{1,ii}(2),0,size(temp{1,ii},1),size(temp{1,ii},2)];            XX{1,ii}=zeros(6,N);            for k=1:N                XX{1,ii}(:,k)=X0{1,ii}';            end            XX0{1,ii}=X0{1,ii};            XX0_temp{1,ii}=zeros(1,6);            x_temp{1,ii}=X0{1,ii};            %%计算权值矩阵            [hist1{1,ii},m_wei1{1,ii}]=calc_hist1(temp{1,ii});        end    end    function [XX,value,w,Estmation_x,Estmation_y,location]=TrackingProcessing1(Im1,hist1,A,B,XX,location,N,w)        % Im1            要跟踪的图像        % hist1          图像中目标的直方图        % A               状态转移矩阵        % B               转移矩阵        % XX             所有粒子初始状态        % location    目标框中心位置坐标及目标框的宽高        % N              粒子个数        % w              粒子权重                [HistDist,XX,w]=Cale_HistDist1(Im1,hist1,A,B,XX,location,N,w);   % 计算每个粒子与目标框的相似度以及权重                % 得到相似度最大的粒子的位置及最大相似度值        w=w/sum(w);                                      % 权值归一化        [Estmation_x,Estmation_y]=estimation1(XX,w);     % 加权估计目标中心位置        BounderProcess1(Estmation_x,Estmation_y,location,Im1);         % 防止出界        Position=[Estmation_x-location(3)/2,Estmation_y-location(4)/2,location(3),location(4)];        temp1=imcrop(Im1,round(Position));        [hist_temp1,~]=calc_hist1(temp1);        % 计算两个直方图的相似度        value=Cal_D_Between_Hist(hist_temp1,hist1);    end    function [ Estmation_x,Estmation_y]=estimation1(XX,w)        x_new=XX(1,:);        y_new=XX(3,:);                %利用加权平均权值估计目标新位置        Estmation_x=x_new*w';        Estmation_y=y_new*w';    end    function BounderProcess1(Estmation_x,Estmation_y,location,Im)        % 边界处理,防止粒子状态更新得到的新位置在图像外面        if Estmation_x-location(3)/2<1            Estmation_x=location(3)/2+1;        end        if Estmation_x-location(3)/2>size(Im,2)            Estmation_x=location(3)/2+size(Im,2);        end        if Estmation_y-location(4)/2<1            Estmation_y=location(4)/2+1;        end        if Estmation_y-location(4)/2>size(Im,1)            Estmation_y=location(4)/2+size(Im,1);        end    end    function value=Cal_D_Between_Hist(hist1,hist2)        Sum1=sum(hist1);        Sum2=sum(hist2);        Sumup = sqrt(hist1.*hist2);        SumDown = sqrt(Sum1*Sum2);        Sumup = sum(Sumup);        value=1-sqrt(1-Sumup/SumDown);    end    function [hist,k_m_wei] = calc_hist1(temp)        % 转为hsv通道        temp=round(255*(rgb2hsv_mex(double(temp))));                %hist为绘制直方图        [a,b,c]=size(temp);        y(1)=a/2;        y(2)=b/2;        k_m_wei=zeros(a,b);%权值矩阵        h=y(1)^2+y(2)^2 ;%带宽        %%计算权值矩阵        for i=1:a            for j=1:b                dist=(i-y(1))^2+(j-y(2))^2;                m_wei=dist/h;                if m_wei<=1                    k_m_wei(i,j)=1-m_wei*m_wei;             % 此种颜色直方图的计算是参考  基于分块粒子滤波的目标跟踪算法                else                                        % 的研究  刘晓龙的硕士论文中的                    k_m_wei(i,j)=0;                end            end        end        C=1/sum(sum(k_m_wei));    % 归一化系数        %%计算目标权值直方图        hist=zeros(1,4096);                for i=1:a            for j=1:b  %rgb颜色空间量化为16*16*16bins                q_h=floor(double(temp(i,j,1))/16);%fix是趋于0取整函数                q_s=floor(double(temp(i,j,2))/16);                q_v=floor(double(temp(i,j,3))/16);                q_temp=q_h*16+q_s*256+1*q_v;%设置每个像素点红色、绿色、蓝色分量所占的比重                hist(q_temp+1)= hist(q_temp+1)+k_m_wei(i,j);%计算直方图中每个像素点所占的权重                %(q_temp+1)=hist(q_temp+1)*C;            end        end        hist=hist*C;    end    function [hist1,XX0,x_temp,Position1,save_hist]=UPdate1(Im,value,Estmation_x,Estmation_y,location,rect,frame,Startframe,XX0,x_temp,hist,save_hist)        % 该函数实现对粒子状态的更新        % 输入:    value                主程序里面求出的相似度最大的值        %               Estmation_x    表示估计当前帧的粒子的目标框中心位置横坐标        %               Estmation_y    表示估计当前帧的粒子的目标框中心位置纵坐标        %               XX                   主程序里面求出的每个粒子的状态矩阵        %               XX0                 主程序里面初始的最初粒子状态        %               XX0_temp         6*1的全零矩阵        %               x_temp            主程序里面初始的最初粒子状态        %               frame              当前帧的图像        %               startframe       当前帧所在帧数        %               location          存放抠取的图像的宽高即坐标的数组        %               rect                 主函数里面存放抠取的图像的宽高即坐标的数组        % 输出:    vv                    这个结构体保存的是画矩形框所用的左上角的坐标及矩形框的宽高        %               X0                   初始状态的更新        %               XX0                 初始状态的更新                XX0_temp=zeros(1,6);        if value>0.8            % 如果相似度大于0.7,直接用估计的目标位置作为下一帧的初始目标            rect(1)=Estmation_x-location(3)/2;            rect(2)=Estmation_y-location(4)/2;            temp1=imcrop(Im,round(rect));            [hist1,~]=calc_hist1(temp1);                        XX0(1)=Estmation_x;            XX0(3)=Estmation_y;            XX0_temp(2)=( XX0_temp(2)+XX0(1)-x_temp(1) );            XX0_temp(4)=( XX0_temp(4)+XX0(3)-x_temp(3) );            XX0(2)=XX0_temp(2)/(frame-Startframe);            XX0(4)=XX0_temp(4)/(frame-Startframe);            x_temp=XX0;            Position1=[rect(1),rect(2),location(3),location(4)];        else            % 否则结合上一帧目标的位置以及速度估计下一帧的初始位置                        XX0(1)= x_temp(1)+x_temp(2);                      % 上一帧粒子的位置加上粒子的速度(即一帧的位移长度)            XX0(3)= x_temp(3)+x_temp(4);            XX0_temp(2)=(XX0_temp(2)+XX0(1)-x_temp(1));       % 重新估计粒子的速度            XX0_temp(4)=(XX0_temp(4)+XX0(3)-x_temp(3));            XX0(2)=XX0_temp(2)/(frame-Startframe);            XX0(4)=XX0_temp(4)/(frame-Startframe);            x_temp=XX0;            % 若value小于阈值 判断目标已经被遮挡  输出下面这个            % 输出目标位置            Position1(1,1)=Estmation_x-location(3)/2;            Position1(1,2)=Estmation_y-location(4)/2;            Position1(1,3)=location(3);            Position1(1,4)=location(4);                        % 得到目标位置对应的区域的直方图            temp1=imcrop(Im,round(Position1));            [hist_temp1,~]=calc_hist1(temp1);            save_hist=[save_hist;hist_temp1];                        % 对hist1直接赋值遮挡前的直方图            hist1=hist;            % 如果满足下面条件 更改hist1赋值            % 本程序里面更新hist的条件是:            % 如果判断了目标已经被遮挡,保存被遮挡前的直方图,然后依次保存连续三帧的            % 估计的目标区域的直方图,若连续两帧的直方图的相似度大于0.7并且第三个估计            % 的直方图与遮挡时保存的直方图的相似度大于0.5,即判断目标已经跟踪回来,用            % 第三个估计的直方图代替之前保存的遮挡时保存的直方图                        if size((save_hist),1)==3                tmp1=Cal_D_Between_Hist(save_hist(1,:),save_hist(2,:));                tmp2=Cal_D_Between_Hist(save_hist(2,:),save_hist(3,:));                tmp3=Cal_D_Between_Hist(save_hist(3,:),hist);                if tmp1>0.75&tmp2>0.75&tmp3>0.6                    hist1=save_hist(3,:);                end                save_hist=[];            end        end    end    function [XX,w_new]=resample_particles1(XX,w,N)        % 该函数实现对粒子的重采样        % 输入:    w                粒子权重        %               XX              主程序里面求出的每个粒子的状态矩阵        %               N                粒子个数        % 输出:    w                粒子权重        %               XX              主程序里面求出的每个粒子的状态矩阵        XX_resample=XX;        w_new=w;        % 判断是否需要重采样        % Neff表示有效的粒子个数(也即是否需要重采样的阈值)        Neff=1/sum(w.*w);        if Neff<50            for i = 1 : N                u = rand; % uniform random number between 0 and 1                qtempsum = 0;                for j = 1 : N                    qtempsum = qtempsum + w(j);                    if qtempsum >= u                        XX_resample(:,i)=XX(:,j);                        w_new(i)=w(i);                        break;                    end                end            end            % 重采样结果输出            XX=XX_resample;            w_new=ones(1,N)/N;        end    end    function DrawResult1(num,Position1,Estmation_x,Estmation_y,XX,frame)        % function DrawResult1(num,Position1,Estmation_x,Estmation_y,XX,frame,EpipolarLine)        % 定义不同颜色        Color=[1 0 0;            0 1 0;            0 0 1            1 1 0;            1 0 1;            0 1 1;            1 1 1];        for k=1:num            hold on;            rectangle('position',Position1{1,k},'EdgeColor',Color(k,:));            hold on;            plot(XX{1,k}(1,:),XX{1,k}(3,:),'r.')            hold on;            plot(Estmation_x{1,k},Estmation_y{1,k},'g+','MarkerSize',10);        end        text(10,20,int2str(frame),'FontSize',15,'Color','r');    endend

    代码里面的rgb转hsv可用MATLAB自带函数或者在点击打开链接这个网址里面下载包含混编的rgb转hsv的全部打包程序。
 
原创粉丝点击