matlab实现混合高斯模型运动物体检测算法

clear all;

Source_Video = VideoReader('viptraffic.avi'); % read the video for processing

Out_Back = VideoWriter('GMM_Background'); Out_Back.FrameRate = 30; open(Out_Back);

Out_Dect = VideoWriter('GMM_Foreground'); Out_Dect.FrameRate = 30; open(Out_Dect);

% ----------------------- 帧尺寸变量 -----------------------

height = Source_Video.Height; width = Source_Video.Width;

% --------------------- 模型变量 -----------------------------------

C = 3; % 混合高斯模型个数(typically 3-5)

D = 2.5; % 正偏差阈值

alpha = 0.01; % 学习率 (between 0 and 1) (from paper 0.01)

thres = 0.75; % 前景阈值 (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)

S = strel('square', 3); % sturcture unit for mophology

% --------------------- initialize component means and weights ---------

pixel_depth = 8; % 8-bit resolution pixel_range = 2^pixel_depth -1; % pixel range

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

% --------------------- process frames -----------------------------------

while(hasFrame(Source_Video))

fr = readFrame(Source_Video);

fr_bw = rgb2gray(fr); % convert frame to grayscale

fg = zeros(height, width); % temporary foreground for each frame

bg_bw = zeros(height, width); % temporary background for each frame

% calculate difference of pixel values from mean for m = 1 : C

u_diff(:,:,m) = abs(double(fr_bw) - double(mean(:,:,m))); end

% update gaussian components for each pixel 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)) % pixel matches component |fr_bw(i,j) - u(i,j,m)| <= 2.5\\sigma(i,j,m)

match = 1; % variable to signal component match

% update weights, mean, sd, p

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(fr_bw(i,j));

sd(i,j,k) = sqrt((1 - p) * (sd(i,j,k)^2) + p * ((double(fr_bw(i,j)) - mean(i,j,k)))^2);

else % pixel doesn't match component

w(i,j,k) = (1 - alpha) * w(i,j,k); % weight slighly decreases

end

end

% normalize and update the weights in each step w_vect = squeeze(w(i,j,:));

w(i,j,:) = w(i,j,:) / sum(w_vect);

% calculate component rank

rank = w(i,j,:) ./ sd(i,j,:); % if no components match, create new component if (match == 0)

% replace the the distribution with lowest confidence [min_confid, min_confid_index] = min(rank); mean(i,j,min_confid_index) = double(fr_bw(i,j)); sd(i,j,min_confid_index) = sd_init; w(i,j,min_confid_index) = min(w(i,j,:));

% calculate moving object for foreground fg(i,j) = fr_bw(i,j);

end

% calculate mixtures of Gaussian for the background confid = w(i,j,:) ./ sd(i,j,:);

[con_sorted, w_index] = sort(confid, 'descend'); w_sum = 0; for k = 1 : C

w_sum = w_sum + w(i,j,k);

bg_bw(i,j) = bg_bw(i,j) + mean(i,j,w_index(k)) * w(i,j,k); if w_sum >= thres

break; end end

end end

% use morphology to file object fg = 255 * (fg > 0); fg = imdilate(fg, S);

% display frame, background and moving object figure(1);

subplot(3,1,1), imshow(fr);

subplot(3,1,2), imshow(uint8(bg_bw)); subplot(3,1,3), imshow(uint8(fg));

% write the background and moving object to videos writeVideo(Out_Back, uint8(bg_bw)); writeVideo(Out_Dect, uint8(fg)); end

% close video structure close(Out_Back); close(Out_Dect);

本文来源：https://www.bwwdw.com/article/m4ax.html