clc
clear
close all
%%
%讀取視頻
mov = VideoReader('test.avi');
N=mov.NumFrames;         %讀取視頻長(zhǎng)度
%%
% 幀設(shè)置
frame = read(mov,1);           % 讀取第一幀作為背景
frame_bw = rgb2gray(frame);          % 將背景轉(zhuǎn)換為灰度圖像
frame_size = size(frame);             %取幀大小
width = frame_size(2);
height = frame_size(1);
frame_g = zeros(height, width);              %前景
backdrop_bw = zeros(height, width);           %背景
%%
%視頻變量
C = 3;                                  % 組成混合高斯的單高斯數(shù)目 (一般3-5)
M = 3;                                  % 
D = 2.5;                                % 閾值（一般2.5個(gè)標(biāo)準(zhǔn)差）
alpha = 0.01;                           % learning rate 學(xué)習(xí)率決定更新速度(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 初始化標(biāo)準(zhǔn)差(for new components) var = 36 in paper
w = zeros(height,width,C);              % initialize weights array 初始化權(quán)值數(shù)組
mean = zeros(height,width,C);           % pixel means 像素均值
sd = zeros(height,width,C);             % pixel standard deviations 像素標(biāo)準(zhǔn)差
u_diff = zeros(height,width,C);         % difference of each pixel from mean 與均值的差
p = alpha/(1/C);                        % initial p variable 參數(shù)學(xué)習(xí)率(used to update mean and sd)
rank = zeros(1,C);                      % rank of components (w/sd)
%%
%初始化均值和權(quán)值
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之間的隨機(jī)數(shù))
            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);       %轉(zhuǎn)換為灰度圖像
    frame_bw=medfilt2(frame_bw);

    % 計(jì)算像素差值
    for m=1:C
        u_diff(:,:,m) = abs(double(frame_bw) - double(mean(:,:,m)));
    end
    %更新每個(gè)像素的背景模型
    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;  % 設(shè)置匹配記號(hào)
                    %更新權(quán)值，均值帖汞，標(biāo)準(zhǔn)差和參數(shù)學(xué)習(xí)率
                    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 權(quán)值減小
                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
            % 如果沒(méi)有匹配的模型則創(chuàng)建新模型
            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,:);             % 計(jì)算優(yōu)先級(jí)
            rank_ind = [1:1:C];
            % 計(jì)算前景
            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