clc
clear
close all
%%
%读取视频
mov = VideoReader('test.avi');
N=mov.NumFrames;         %读取视频长度
%%
% 帧设置
frame = read(mov,1);           % 读取第一帧作为背景
frame_bw = rgb2gray(frame);          % 将背景转换为灰度图像
frame_size = size(frame);             %取帧大小
width = frame_size(2);
height = frame_size(1);
frame_g = zeros(height, width);              %前景
backdrop_bw = zeros(height, width);           %背景
%%
%视频变量
C = 3;                                  % 组成混合高斯的单高斯数目 (一般3-5)
M = 3;                                  % 
D = 2.5;                                % 阈值（一般2.5个标准差）
alpha = 0.01;                           % learning rate 学习率决定更新速度(between 0 and 1) (from paper 0.01)
thresh = 0.25;                          % foreground threshold 前景阈值(0.25 or 0.75 in paper)
sd_init = 6;                            % initial standard deviation 初始化标准差(for new components) var = 36 in paper
w = zeros(height,width,C);              % initialize weights array 初始化权值数组
mean = zeros(height,width,C);           % pixel means 像素均值
sd = zeros(height,width,C);             % pixel standard deviations 像素标准差
u_diff = zeros(height,width,C);         % difference of each pixel from mean 与均值的差
p = alpha/(1/C);                        % initial p variable 参数学习率(used to update mean and sd)
rank = zeros(1,C);                      % rank of components (w/sd)
%%
%初始化均值和权值
pixel_depth = 8;                        % 8-bit resolution 像素深度为8位
pixel_range = 2^pixel_depth -1;         % pixel range 像素范围2的7次方0—255（# of possible values)
for i=1:height
    for j=1:width
        for k=1:C
            mean(i,j,k) = rand*pixel_range;     % means random (0-255之间的随机数)
            w(i,j,k) = 1/C;                     % weights uniformly dist
            sd(i,j,k) = sd_init;                % initialize to sd_init
            
        end
    end
end

%%
%处理帧
for n = 1:N
    frame = read(mov,n);       %  读%第n帧
    frame_bw = rgb2gray(frame);       %转换为灰度图像
    frame_bw=medfilt2(frame_bw);

    % 计算像素差值
    for m=1:C
        u_diff(:,:,m) = abs(double(frame_bw) - double(mean(:,:,m)));
    end
    %更新每个像素的背景模型
    for i=1:height
        for j=1:width
            match = 0;
            for k=1:C                       
                if (abs(u_diff(i,j,k)) <= D*sd(i,j,k))       %像素匹配了模型
                    match = 1;  % 设置匹配记号
                    %更新权值，均值，标准差和参数学习率
                    w(i,j,k) = (1-alpha)*w(i,j,k) + alpha;
                    p = alpha/w(i,j,k);                  
                    mean(i,j,k) = (1-p)*mean(i,j,k) + p*double(frame_bw(i,j));
                    sd(i,j,k) =   sqrt((1-p)*(sd(i,j,k)^2) + p*((double(frame_bw(i,j)) - mean(i,j,k)))^2);
                else                                   
                    w(i,j,k) = (1-alpha)*w(i,j,k);      % weight slighly decreases 权值减小
                end
            end
            backdrop_bw(i,j)=0;
            for k=1:C
                backdrop_bw(i,j) = backdrop_bw(i,j)+ mean(i,j,k)*w(i,j,k);  %更新背景
            end
            % 如果没有匹配的模型则创建新模型
            if (match == 0)
                [min_w, min_w_index] = min(w(i,j,:));  
                mean(i,j,min_w_index) = double(frame_bw(i,j));
                sd(i,j,min_w_index) = sd_init;
            end
            rank = w(i,j,:)./sd(i,j,:);             % 计算优先级
            rank_ind = [1:1:C];
            % 计算前景
            frame_g(i,j) = 0;
            while ((match == 0)&&(k<=M))
                    if (abs(u_diff(i,j,rank_ind(k))) <= D*sd(i,j,rank_ind(k)))
                        frame_g(i,j) = 0;     %black = 0
                    else
                        frame_g(i,j) = frame_bw(i,j);     
                    end
                k = k+1;
            end
        end
        
    end
    
    figure(1),subplot(1,2,1),imshow(frame)    %显示输入图像
%     subplot(1,3,2),imshow(uint8(backdrop_bw))    %显示背景图像
    subplot(1,2,2),imshow(uint8(frame_g))     %显示前景图像
end