demoems.m

混合高斯概率密度模型,其参数估计可以通过期望最大化( EM) 迭代算法获得。

%EM-shift demo
%a simple color-histogram tracking demo script
%
%Z.Zivkovic, B.Krose
%An EM-like algorithm for color-histogram-based object tracking"
%IEEE Conference CVPR, June, 2004
%
%University of Amsterdam, The Netherlands
%Author:Zoran Zivkovic, www.zoranz.net
%Date: 17-5-2005
clear
sAvi='e:\EM算法\em\mix2.avi';%set your own avi file name
nNImageStart=120;nNImageEnd=270;%select frames

%choose the tracker
%sTracker='MS'
dScale=0.1;
sTracker='EMS'
beta=1.2;
nIterationMax=6;

%%%%%%%%%%%%%%%%%%%
%Initialization
%%%%%%%%%%%%%%%%%%%
%manually select object
data=aviread('mix2.avi',nNImageStart);
im0=data.cdata;
figure(1),imMask=roipoly(im0);%roipoly is from the Image processing toolbox
%shape=>ellipse
[ry,rx]=find(imMask);m0=mean([rx,ry])';V0=cov([rx,ry]);
[K]=GaussKernel(V0);%make kernel of appropriate size
%make object histogram
data=GetAffineRegion(im0,m0,V0,K.rX,K.rY,K.sigmas);
nBinsD=4;histO=ColorHist(data,nBinsD,K.rK);
%state vector
s0=MeanVarToS(m0,V0);nStates=length(s0);s=s0;
%dynamic model s=A*s+W
A=eye(nStates);
W=[3^2  0   0   0   0;
    0   3^2 0   0   0;
    0   0   1^2 0   0;
    0   0   0   1^2 0;
    0   0   0   0   1^2];
lambda=0.2;%observation model p~exp(-0.5*(1-rho)/lambda^2)
Pk=100^2*eye(nStates);%some high initial value for the variance

%%%%%%%%%%%%%%%%%%%
%Tracking
%%%%%%%%%%%%%%%%%%%
for iFrame=nNImageStart:nNImageEnd
    data=aviread(sAvi,iFrame);im0=data.cdata;%read an image
    %different trackers
    switch sTracker
        case 'EMS'
            %Z.Zivkovic, B.Krose
            %An EM-like algorithm for color-histogram-based object tracking"
            %IEEE Conference CVPR, June, 2004
            %predict
            s=A*s;
            Pk=A*Pk*A'+W;
            %find max
            [xt,V]=SToMeanVar(s);%prediction as starting point
            for iIteration=1:nIterationMax
                [xt,V]=EMShift(im0,xt,V,histO,K,beta);
            end
            z=MeanVarToS(xt,V);
            R=FitGauss(im0,z,histO,K)*lambda^2;
            %Kalman filter
            KFgain=Pk*inv(Pk+R);
            s=s+KFgain*(z-s);
            Pk=(eye(nStates)-KFgain)*Pk;
        case 'MS'
            %D. Comaniciu, V. Ramesh, P. Meer
            %Kernel-Based Object Tracking,
            %IEEE Trans. Pattern Analysis Machine Intell., Vol. 25, No. 5, 564-575, 2003
            %predict
            s=A*s;
            Pk=A*Pk*A'+W;
            %Find the local max using the prediction as the start point
            [xt,V]=SToMeanVar(s);%prediction as starting point
            for iIteration=1:nIterationMax
                xt=EMShift(im0,xt,V,histO,K,beta);%mean shift does not use the second moment
            end
            rho0=Similarity(im0,xt,V,histO,K);
            z0=MeanVarToS(xt,V);
            %measure V by trying a number of different scales
            rD=dScale*[0 0 0 0 0;0 0 s(3) 0 s(5); 0 0 -s(3) 0 -s(5)]';
            %try bigger
            [xt,V]=SToMeanVar(s+rD(:,2));
            for iIteration=1:nIterationMax
                xt=EMShift(im0,xt,V,histO,K,beta);%mean shift does not use the second moment
            end
            rho1=Similarity(im0,xt,V,histO,K);
            z1=MeanVarToS(xt,V);
            %try smaller
            [xt,V]=SToMeanVar(s+rD(:,3));
            for iIteration=1:nIterationMax
                xt=EMShift(im0,xt,V,histO,K,beta);%mean shift does not use the second moment
            end
            rho2=Similarity(im0,xt,V,histO,K);
            z2=MeanVarToS(xt,V);
            %choose the best
            [dummy,ri]=max([rho0 rho1 rho2]);
            if (ri(1)==1) z=z0; end;
            if (ri(1)==2) z=z1; end;
            if (ri(1)==3) z=z2; end;
            
            R=FitGauss(im0,z,histO,K)*lambda^2;
            %since we can be far from a local minima for the V
            %the local curvature is not good estimate for the V uncertainty !!! 
            %-use some fixed value - mean value from the EMS for example
            R(3,3)=1;R(4,4)=1;R(5,5)=1;
            %Kalman filter equations
            KFgain=Pk*inv(Pk+R);
            s=s+KFgain*(z-s);
            Pk=(eye(nStates)-KFgain)*Pk;    
    end

    %show
    figure(1),imshow(im0);hold on;%imshow is from the Image processing toolbox
    h=PlotEllipse(s);set(h,'LineWidth',4,'Color',[1 0.7 1],'LineStyle','--');
    drawnow;hold off
end