EM+GMM基于EM和GMM算法的目标轨迹跟踪和异常行为识别matlab仿真

Posted fpga和matlab

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了EM+GMM基于EM和GMM算法的目标轨迹跟踪和异常行为识别matlab仿真相关的知识,希望对你有一定的参考价值。

1.软件版本

matlab2013b

2.系统概述

3.部分源码

clc;
clear;
close all;
warning off;
addpath 'func\\'

%训练数据
Is = imread('Track_file\\tower1.bmp');
load Track_file\\track1.mat
load Track_file\\track2.mat

figure(1);
subplot(121);
imshow(Is);
hold on;
X1=[];
Y1=[];
T1=[];
W1=[];
H1=[];
for i = 1:length(track1)
    i
    if isempty(track1(1,i).x) == 1
       X = 0;
       Y = 0;
    else
       X = cell2mat(track1(1,i).x(1,1));
       Y = cell2mat(track1(1,i).y(1,1));
       X1=[X1;cell2mat(track1(1,i).x(1,1))];
       Y1=[Y1;cell2mat(track1(1,i).y(1,1))];
       T1=[T1;cell2mat(track1(1,i).time(1,1))];
       W1=[W1;cell2mat(track1(1,i).w(1,1))];
       H1=[H1;cell2mat(track1(1,i).h(1,1))];        
    end
    plot(X,Y,'c.');
    hold on;
end
title('Track1');

subplot(122);
imshow(Is);
hold on;
X2=[];
Y2=[];
T2=[];
W2=[];
H2=[];
for i = 1:length(track2)
    i
    if isempty(track2(1,i).x) == 1
       X = 0;
       Y = 0;
    else
       X = cell2mat(track2(1,i).x(1,1));
       Y = cell2mat(track2(1,i).y(1,1));
       X2=[X2;cell2mat(track2(1,i).x(1,1))];
       Y2=[Y2;cell2mat(track2(1,i).y(1,1))];
       T2=[T2;cell2mat(track2(1,i).time(1,1))];
       W2=[W2;cell2mat(track2(1,i).w(1,1))];
       H2=[H2;cell2mat(track2(1,i).h(1,1))];  
    end
    plot(X,Y,'c.');
    hold on;
end
title('Track2');
 

tao=20;
X1s=[];
Y1s=[];
T1s=[];
W1s=[];
H1s=[];
X_samp  = [];
Y_samp  = [];
T_samp  = [];
W_samp  = [];
H_samp  = [];
X_samps1 = [];
Y_samps1 = [];
T_samps1 = [];
W_samps1 = [];
H_samps1 = [];
for i = 1:length(track1)
    i
    if isempty(track1(1,i).x) == 1
    X1si=0;
    Y1si=0;
    T1si=0;
    W1si=0;
    H1si=0;
    else
    X1si=cell2mat(track1(1,i).x(1,1));
    Y1si=cell2mat(track1(1,i).y(1,1));
    T1si=cell2mat(track1(1,i).time(1,1));
    W1si=cell2mat(track1(1,i).w(1,1));
    H1si=cell2mat(track1(1,i).h(1,1));
    end
    if i <= tao;
       X_samp = (X1si);
       Y_samp = (Y1si);
       T_samp = 0.00001*ones(length(X_samp),1);
       W_samp = (W1si);
       H_samp = (H1si);
    else
       X_samp = (X1si);
       Y_samp = (Y1si);
       T_samp = 0.00001*ones(length(X_samp),1);
       Ls           = min(length(X1si),length(X1si-tao));
       for j = 1:Ls
           T_samp(j,1) = sqrt((X1si(j) - X1si-tao(j))^2 + (Y1si(j) - Y1si-tao(j))^2)/tao;
       end
       W_samp = (W1si);
       H_samp = (H1si);
    end
    X_samps1 = [X_samps1;X_samp];
    Y_samps1 = [Y_samps1;Y_samp];
    T_samps1 = [T_samps1;T_samp];
    W_samps1 = [W_samps1;W_samp];
    H_samps1 = [H_samps1;H_samp];
end


X2s=[];
Y2s=[];
T2s=[];
W2s=[];
H2s=[];
X_samp  = [];
Y_samp  = [];
T_samp  = [];
W_samp  = [];
H_samp  = [];
X_samps2 = [];
Y_samps2 = [];
T_samps2 = [];
W_samps2 = [];
H_samps2 = [];
for i = 1:length(track2)
    i
    if isempty(track2(1,i).x) == 1
    X2si=0;
    Y2si=0;
    T2si=0;
    W2si=0;
    H2si=0;
    else
    X2si=cell2mat(track2(1,i).x(1,1));
    Y2si=cell2mat(track2(1,i).y(1,1));
    T2si=cell2mat(track2(1,i).time(1,1));
    W2si=cell2mat(track2(1,i).w(1,1));
    H2si=cell2mat(track2(1,i).h(1,1));
    end
    if i <= tao;
       X_samp = (X2si);
       Y_samp = (Y2si);
       T_samp = 0.00001*ones(length(X_samp),1);
       W_samp = (W2si);
       H_samp = (H2si);
    else
       X_samp = (X2si);
       Y_samp = (Y2si);
       T_samp = 0.00001*ones(length(X_samp),1);
       Ls           = min(length(X2si),length(X2si-tao));
       for j = 1:Ls
           T_samp(j,1) = sqrt((X2si(j) - X2si-tao(j))^2 + (Y2si(j) - Y2si-tao(j))^2)/tao;
       end
       W_samp = (W2si);
       H_samp = (H2si);
    end
    X_samps2 = [X_samps2;X_samp];
    Y_samps2 = [Y_samps2;Y_samp];
    T_samps2 = [T_samps2;T_samp];
    W_samps2 = [W_samps2;W_samp];
    H_samps2 = [H_samps2;H_samp];
end
 

%原始数据的训练,用set2的数据来训练
feature1    = [X_samps1,Y_samps1,T_samps1,W_samps1,H_samps1];
feature2    = [X_samps2,Y_samps2,T_samps2,W_samps2,H_samps2];
feature     = [feature2(5000:7000,:)];%
tmps        = feature;
[Px, model] = gmm(tmps,2); 
% save EmGmm_model.mat model
 
% figure;
% imshow(Is);
% hold on;
% [~, cls_ind] = max(Px',[],1); 
% for tt = 1:1
%     tt
%     for jj = 1:length(cls_ind);
%         if cls_ind(jj)==1
%            plot(feature(jj,1),feature(jj,2),'b.');
%         end
%         if cls_ind(jj)==2
%            plot(feature(jj,1),feature(jj,2),'r.'); 
%         end
%         hold on;
%     end
%     title('异常行为识别');
% end
clc;
clear;
close all;
warning off;
addpath 'func\\'




%初始背景的提取
%注意,原始的视频非常大,这里需要截取其中需要处理的帧数
Bg_Length   = 20;
Start_Frame = Bg_Length+1;
End_Frame   = 999;
Simulation_frame = min(800,End_Frame);
avi1        = aviread('vedio\\4.AVI',[Start_Frame:End_Frame]);
vedio1      = avi1.cdata;
avi2        = aviread('vedio\\4.AVI',[1:Bg_Length]);
vedio2      = avi2.cdata;
vedio       = [vedio2,vedio1];
%步骤一:结合高斯混合模型的帧差法进行目标的提取
%这个步骤主要是目标提取,并获得特征序列两个部分
RR                  = 240;%处理视频大小
CC                  = 320;
K                   = 4;                   %组件
Alpha               = 0.002;               %适应权重速度
Rho                 = 0.01;                %适应权重速度协方差
Deviation_sq        = 49;                  %阈值用于查找匹配 
Variance            = 2;                   %初始方差为新放置组件
Props               = 0.00001;             %最初为新放置 
Back_Thresh         = 0.6;                 %体重的比例必须占背景模型
Comp_Thresh         = 10;                  %滤掉连接组件的较小的尺寸
SHADOWS             = [0.7,0.25,0.85,0.95]; 
frameNum_Original   = length(vedio);


CRGB = 3;
D    = RR * CC;
Temp = zeros(RR,CC,CRGB,frameNum_Original,'uint8');
 
for tt = 1:frameNum_Original
  im                       = vediott;
  Temp(:,:,:,tt)           = im;
end
Temp = reshape(Temp,size(Temp,1)*size(Temp,2),size(Temp,3),size(Temp,4));      

Mus                 = zeros(D,K,CRGB);
Mus(:,1,:)          = double(Temp(:,:,1));
Mus(:,2:K,:)        = 255*rand([D,K-1,CRGB]);
Sigmas              = Variance*ones(D,K,CRGB); 
Weights             = [ones(D,1),zeros(D,K-1)];
Squared             = zeros(D,K);    
Gaussian            = zeros(D,K);     
Weight              = zeros(D,K);
background          = zeros(RR,CC);
Shadows             = zeros(RR,CC);
Images0             = zeros(RR,CC);        
Images1             = zeros(RR,CC);  
Images2             = zeros(RR,CC);   
background_Update   = zeros(RR,CC,CRGB); 
 

Xft        = [];
Yft        = [];
Timeft     = [];
Wft        = [];
Hft        = [];
Smin       = 50;
Smax       = 150;
PAlpha_min = 0.001;
PAlpha_max = 0.004;
Pdf1       = 0;
Pdf2       = 0;
tao        = 20;
Xx         = [];
flag       = 0;
Xf_save    = cell(End_Frame,1);
Yf_save    = cell(End_Frame,1);
Timef_save = cell(End_Frame,1);
Wf_save    = cell(End_Frame,1);
Hf_save    = cell(End_Frame,1);
S          = 0;
load EmGmm_model.mat;

for tt = 1:Simulation_frame
    disp('当前帧数');
    tt
    %由分析得到的反馈信号
    %由分析得到的反馈信号
    %Minimum Object Size 反馈,目标大小检测
    if tt > 1
       Ss = Smin*Pdf1 + Smax*(1-Pdf1);
    else
       Ss = Smin; 
    end
    %Background Learning Rate 反馈,背景更新学习率
    if tt > 1
       Alpha = PAlpha_min*Pdf2 + PAlpha_max*(1-Pdf2);
    else
       Alpha = PAlpha_min; 
    end
 
    %基于高斯混合模型的帧差分背景提取算法
    %基于高斯混合模型的帧差分背景提取算法
    image = Temp(:,:,tt);
    for kk = 1:K   
        Datac         = double(Temp(:,:,tt))-reshape(Mus(:,kk,:),D,CRGB);
        Squared(:,kk) = sum((Datac.^ 2)./reshape(Sigmas(:,kk,:),D,CRGB),2); 
    end
    [junk,index] = min(Squared,[],2); 
    Gaussian                                                = zeros(size(Squared));
    Gaussian(sub2ind(size(Squared),1:length(index),index')) = ones(D,1);
    Gaussian                                                = Gaussian&(Squared<Deviation_sq);
    %参数更新
    Weights = (1-Alpha).*Weights+Alpha.*Gaussian;
    for kk = 1:K
        pixel_matched   = repmat(Gaussian(:,kk),1,CRGB);
        pixel_unmatched = abs(pixel_matched-1);
        Mu_kk           = reshape(Mus(:,kk,:),D,CRGB);
        Sigma_kk        = reshape(Sigmas(:,kk,:),D,CRGB);
        Mus(:,kk,:)     = pixel_unmatched.*Mu_kk+pixel_matched.*(((1-Rho).*Mu_kk)+(Rho.*double(image)));
        Mu_kk           = reshape(Mus(:,kk,:),D,CRGB); 
        Sigmas(:,kk,:)  = pixel_unmatched.*Sigma_kk+pixel_matched.*(((1-Rho).*Sigma_kk)+repmat((Rho.* sum((double(image)-Mu_kk).^2,2)),1,CRGB));       
    end
    replaced_gaussian   = zeros(D,K); 
    mismatched          = find(sum(Gaussian,2)==0);       
    for ii = 1:length(mismatched)
        [junk,index]                            = min(Weights(mismatched(ii),:)./sqrt(Sigmas(mismatched(ii),:,1)));
        replaced_gaussian(mismatched(ii),index) = 1;
        Mus(mismatched(ii),index,:)             = image(mismatched(ii),:);
        Sigmas(mismatched(ii),index,:)          = ones(1,CRGB)*Variance;
        Weights(mismatched(ii),index)           = Props;  
    end
    Weights         = Weights./repmat(sum(Weights,2),1,K);
    active_gaussian = Gaussian+replaced_gaussian;
    %背景分割 
    [junk,index]    = sort(Weights./sqrt(Sigmas(:,:,1)),2,'descend');
    bg_gauss_good   = index(:,1);
    linear_index    = (index-1)*D+repmat([1:D]',1,K);
    weights_ordered = Weights(linear_index);
    for kk = 1:K
        Weight(:,kk)= sum(weights_ordered(:,1:kk),2);
    end
    bg_gauss(:,2:K) = Weight(:,1:(K-1)) < Back_Thresh;
    bg_gauss(:,1)   = 1;           
    bg_gauss(linear_index)     = bg_gauss;
    active_background_gaussian = active_gaussian & bg_gauss;
    foreground_pixels          = abs(sum(active_background_gaussian,2)-1);
    foreground_map             = reshape(sum(foreground_pixels,2),RR,CC);
    Images1(:,:,tt)            = foreground_map;   
    objects_map                = zeros(size(foreground_map),'int32');
    object_sizes               = [];
    Obj_pos                    = [];
    new_label                  = 1;
    %计算连通区域
    [label_map,num_labels]     = bwlabel(foreground_map,8);
    
    for label = 1:num_labels 
       object      = (label_map == label);
       object_size = sum(sum(object));
       if(object_size >= Comp_Thresh)
          objects_map             = objects_map + int32(object * new_label);
          object_sizes(new_label) = object_size;
          [X,Y]                   = meshgrid(1:CC,1:RR);    
          object_x                = X.*object;
          object_y                = Y.*object;
          Obj_pos(:,new_label)    = [sum(sum(object_x)) / object_size;
				                     sum(sum(object_y)) / object_size];
          new_label               = new_label + 1;
       end
    end
    num_objects = new_label - 1;
    %去除阴影
    index                       = sub2ind(size(Mus),reshape(repmat([1:D],CRGB,1),D*CRGB,1),reshape(repmat(bg_gauss_good',CRGB,1),D*CRGB,1),repmat([1:CRGB]',D,1));
    background                  = reshape(Mus(index),CRGB,D);
    background                  = reshape(background',RR,CC,CRGB); 
    background                  = uint8(background);
    background_Update           = background;
    background_hsv              = rgb2hsv(background);
    image_hsv                   = rgb2hsv(vediott);
    for i = 1:RR
        for j = 1:CC      
            if (objects_map(i,j))&&...
               (abs(image_hsv(i,j,1)-background_hsv(i,j,1))<SHADOWS(1))&&...
               (image_hsv(i,j,2)-background_hsv(i,j,2)<SHADOWS(2))&&...
               (SHADOWS(3)<=image_hsv(i,j,3)/background_hsv(i,j,3)<=SHADOWS(4))
               Shadows(i,j) = 1;  
            else
               Shadows(i,j) = 0;  
            end               
        end    
    end
    %运动目标检测,二值图
    Images0                  = objects_map;
    objecs_adjust_map        = Shadows;
    Images3                  = objecs_adjust_map;   
    %运动目标检测-最后处理结果
    [res3,res4,Xf,Yf,Timef,Wf,Hf] = func_fangk(Images0,vediott,tt,Bg_Length,Ss);
 
    %根据上面的结果,提取特征向量
    %根据上面的结果,提取特征向量
    Lens  = 800;
    if length(Xft) < Lens
        Xft   = [Xft;Xf];
        Yft   = [Yft;Yf];
        Timeft= [Timeft;Timef];
        Wft   = [Wft;Wf];
        Hft   = [Hft;Hf];
    else
        Xft   = [Xft;Xf];
        Yft   = [Yft;Yf];
        Timeft= [Timeft;Timef];
        Wft   = [Wft;Wf];
        Hft   = [Hft;Hf];
        Xft(1:end-800)   = [];
        Yft(1:end-800)   = [];
        Timeft(1:end-800)= [];
        Wft(1:end-800)   = [];
        Hft(1:end-800)   = [];
    end
    
    
    
    figure(2);
    subplot(121)
    imshow(vediott);
    title('原始图像');
    subplot(122)
    imshow(uint8(res3),[]);
    hold on;
    if  tt > Bg_Length
        if length(Xft) > Lens
           VX = Xft(length(Xft)-Lens:length(Xft));
           VY = Yft(length(Xft)-Lens:length(Xft));
           plot(Xft(length(Xft)-Lens:length(Xft)),Yft(length(Xft)-Lens:length(Xft)),'c.');
        else
           VX = Xft;
           VY = Yft;
           plot(Xft,Yft,'c.'); 
        end
    end
    title('运动目标检测结果');
 
 
    %计算Pdf1
    tmps    = [Wf,Hf];
    [Ns,Ds] = size(tmps);
    Pdf1s   = min(calc_prob(tmps,Ns,2,model.Miu(:,4:5),model.Sigma(4:5,4:5,:),Ds));
    Pdf1    = mean(Pdf1s);
    %计算Pdf2
    if mod(tt,Lens) <= 20
       tmps = 0;
    else
       tmps = sqrt((Xft(mod(tt,Lens))-Xft(mod(tt,Lens)-20))^2+(Yft(mod(tt,Lens))-Yft(mod(tt,Lens)-20))^2)/20;
    end
    [Ns,Ds] = size(tmps);
    Pdf2s = min(calc_prob(tmps,Ns,2,model.Miu(:,5),model.Sigma(5,5,:),Ds)); 
    Pdf2  = mean(Pdf2s);
    
    pause(0.01);
    %将数据保存,用于后期处理
 
    
    Xf_savett=    Xf;
    Yf_savett=    Yf;
    Timef_savett= Timef;
    Wf_savett=    Wf;
    Hf_savett=    Hf;
    Pimagestt=    Images0;
 
end
 
save TEMP\\R4.mat  Xf_save Yf_save Timef_save Wf_save Hf_save Pimages

 







clc;
clear;
close all;
warning off;
addpath 'func\\'

%通过前面的目标跟踪,这里我们主要对最后的结果进行显示
%本质上C这个程序是B里面的,但是由于涉及到的变量较多,如果再增加一些显示的变量,仿真及其缓慢,而且会出现OUT OF Memory
%所以这里单独对结论进行显示


%初始背景的提取
%注意,原始的视频非常大,这里需要截取其中需要处理的帧数
Bg_Length   = 20;
Start_Frame = Bg_Length+1;
End_Frame   = 999;
Simulation_frame = min(800,End_Frame);


load EmGmm_model.mat;
 
load TEMP\\R4.mat
avi1        = aviread('vedio\\4.AVI',[Start_Frame:End_Frame]);
vedio1      = avi1.cdata;
avi2        = aviread('vedio\\4.AVI',[1:Bg_Length]);
vedio2      = avi2.cdata;
vedio       = [vedio2,vedio1];
tao     = 20;
X_samp  = [];
Y_samp  = [];
T_samp  = [];
W_samp  = [];
H_samp  = [];
X_samps = [];
Y_samps = [];
T_samps = [];
W_samps = [];
H_samps = [];
X_samps2 = [];
Y_samps2 = [];
T_samps2 = [];
W_samps2 = [];
H_samps2 = [];
for tt = 1:Simulation_frame;
    tt
    Images0 = Pimagestt;
    [res3,res4,Xf,Yf,Timef,Wf,Hf] = func_fangk2(Images0,vediott,tt,Bg_Length);
    
    %下面开始进行异常类型识别
    %重新采样
    X_sampstt = Xf;
    Y_sampstt = Yf;
    T_sampstt = Timef;
    W_sampstt = Wf;
    H_sampstt = Hf;
    
     if tt <= tao;
       X_samp = Xf;
       Y_samp = Yf;
       T_samp = Timef;
       W_samp = Wf;
       H_samp = Hf;
    else
       X_samp = Xf;
       Y_samp = Yf;
       T_samp = 0.001*ones(length(X_samp),1);
       Ls           = min(length(X_sampstt),length(X_sampstt-tao));
       for j = 1:Ls
           T_samp(j,1) = sqrt((X_sampstt(j) - X_sampstt-tao(j))^2 + (Y_sampstt(j) - Y_sampstt-tao(j))^2)/tao;
       end
       W_samp = Wf;
       H_samp = Hf;
    end
    
    X_samps2 = [X_samps2;X_samp];
    Y_samps2 = [Y_samps2;Y_samp];
    T_samps2 = [T_samps2;T_samp];
    W_samps2 = [W_samps2;W_samp];
    H_samps2 = [H_samps2;H_samp];
end


feature = [X_samps2,Y_samps2,T_samps2,W_samps2,H_samps2];
[Ns,Ds] = size(feature);  
Px      = calc_prob(feature,Ns,2,model.Miu,model.Sigma,Ds);
[~, cls_ind] = max(Px',[],1); 


%最后完异常检测效果
XV = [];
YV = [];
for tt = 1:1
    tt
    Images0 = Pimagestt;
    [res3,res4,Xf,Yf,Timef,Wf,Hf] = func_fangk2(Images0,vediott,tt,Bg_Length);

    figure(1);
    subplot(121)
    imshow(vediott);
    title('原始图像');
    subplot(122)
    imshow(uint8(res3),[]);
    hold on;
    for jj = 1:length(cls_ind);
        if length(find(cls_ind==1)) > 1.2*length(find(cls_ind==2)) 
           if cls_ind(jj)==1
              plot(X_samps2(jj),Y_samps2(jj),'b.');
              XV(jj,1) = [X_samps2(jj)];
              YV(jj,1) = [Y_samps2(jj)];
              XV(jj,2) = [1];
              YV(jj,2) = [1];
           else
              plot(X_samps2(jj),Y_samps2(jj),'r.'); 
              XV(jj,1) = [X_samps2(jj)];
              YV(jj,1) = [Y_samps2(jj)];
              XV(jj,2) = [2];
              YV(jj,2) = [2];
           end
        end
        if length(find(cls_ind==1)) <= 1.2*length(find(cls_ind==2)) 
           if cls_ind(jj)==1
              plot(X_samps2(jj),Y_samps2(jj),'r.');
              XV(jj,1) = [X_samps2(jj)];
              YV(jj,1) = [Y_samps2(jj)];
              XV(jj,2) = [2];
              YV(jj,2) = [2];
           else
              plot(X_samps2(jj),Y_samps2(jj),'b.');  
              XV(jj,1) = [X_samps2(jj)];
              YV(jj,1) = [Y_samps2(jj)];
              XV(jj,2) = [1];
              YV(jj,2) = [1];
           end
        end        
 
        hold on;
    end
    title('异常行为识别');
    pause(0.0025);
end

%异常检测动画效果,将系统最后处理的结果进行显示
S1 = 0;
E1 = 0;
for tt = 1:Simulation_frame;
    tt
    Images0 = Pimagestt;
    [res3,res4,Xf,Yf,Timef,Wf,Hf] = func_fangk2(Images0,vediott,tt,Bg_Length);
    figure(2);
    subplot(121)
    imshow(vediott);
    title('原始图像');
    subplot(122)
    imshow(uint8(res3),[]);
    hold on;
    S1 = E1 + 1;
    E1 = E1 + length(Xf);
    Xtmp = XV(1:E1,1);
    Ytmp = YV(1:E1,1);
    Flag = XV(1:E1,2);
    for j = 1:length(Flag)
        if Flag(j) == 1
           plot(Xtmp(j),Ytmp(j),'b.');  
           hold on;
        end
        if Flag(j) == 2
           plot(Xtmp(j),Ytmp(j),'r.');  
           hold on;
        end
        hold on;
    end
    hold on;
    title('异常行为检测');
    pause(0.01);
end






3.仿真结果

A10-19

以上是关于EM+GMM基于EM和GMM算法的目标轨迹跟踪和异常行为识别matlab仿真的主要内容,如果未能解决你的问题,请参考以下文章

05 EM算法 - 高斯混合模型 - GMM

用EM思想估计GMM(高斯混合聚类)

高斯混合模型(GMM)及EM算法的初步理解

EM算法及其应用GMM/pLSA/LDA

EM算法:GMM训练算法

GMM与EM共舞