基于高斯混合模型的目标检测和背景提取算法matlab仿真,含仿真操作录像

clc; clear; close all; warning off; addpath 'func\' path = [cd '/Save_Temp/RGB/']; files = dir([path '*.' 'jpg']); files = sort({files.name}); T = size(files, 2); % Number of frames [HEIGHT, WIDTH, C] = size(imresize(imread([path files{1}]),[128,100])); % Get the dimension of one picture D = HEIGHT * WIDTH; image_sequence = zeros(HEIGHT, WIDTH, C, T, 'uint8'); % Used for image dispaly data = zeros(HEIGHT, WIDTH, C, T, 'uint8'); % load data from the image sequence for tt = 1:T im = imresize(imread([path files{tt}]),[128,100]); data(:, :, :, tt) = im; image_sequence(:, :, :, tt) = im; end % Turn each picture into a column vector of color triplets data = reshape(data, size(data, 1) * size(data, 2), size(data, 3), ... size(data, 4)); %-------------------------------------------------------------------------- %-------------------Initialise mix_gaussian parameters--------------------- K = 3; % how many components does the mixture have ALPHA = 0.01; % How fast to adapt the component weights RHO = 0.01; % How fast to adapt the component means and covariances DEVIATION_SQ_THRESH = 7^2; % Threshold used to find matching versus unmatching components. INIT_VARIANCE = 3; % Initial variance for newly placed components INIT_MIXPROP = 0.00001; % Initial prior for newly placed components BACKGROUND_THRESH = 0.8; % Percentage of weight that must be accounted for by background models. COMPONENT_THRESH = 10; % Filter out connected components of smaller size %-------------------------------------------------------------------------- %----------------------------model initialise------------------------------ % Place the mean of first Gaussian on the value of first observed pixel. % Place all other Gaussians mean randomly. Mus = zeros(D, K, C); Mus(:, 1, :) = double(data(:, :, 1)); Mus(:, 2:K, :) = rand([D, K-1, C]) * 255; % Genenrally, R��G��B component are independent and have the same variance % We just initialise the variance of the Gaussian with a relatively high value. Sigmas = ones(D, K, C) * INIT_VARIANCE; % Give each Gaussian except the first one zero weight. This ensures that uninitialised % Gaussians (i.e. everything but the first) will be replaced first time we have the chance. Weights = [ones(D, 1), zeros(D, K-1)]; %-------------------------------------------------------------------------- %---------------------------------- iteration ----------------------------- deviations_squared = zeros(D, K); % the deviation from the mean of the Gaussian matched_gaussian = zeros(D, K); % All gaussians that matched a pixel accumulated_weights = zeros(D, K); background = zeros(HEIGHT, WIDTH, T); shadow_mark = zeros(HEIGHT, WIDTH); % Used for image dispaly foreground_map_sequence = zeros(HEIGHT, WIDTH, T); foreground_with_map_sequence = zeros(HEIGHT, WIDTH, T); foreground_map_adjust_sequence = zeros(HEIGHT, WIDTH, T); background_sequence = zeros(HEIGHT, WIDTH, C, T); for tt = 1:T tt image = data(:, :, tt); %----------------------Light and shadow correction--------------------- % Gama correction % ratio = sum(sum(background_hsv(:,:,3)))/(sum(sum(background_hsv(:,:,3)))... % + sum(sum(background_hsv(:,:,3)))); % image_hsv(:,:,3) = image_hsv(:,:,3).^ratio; % image_adjust = hsv2rgb(image_hsv) * 255; %---------------------------------------------------------------------- %---------------------------match definition--------------------------- % Compute the squared Mahalanobis distance. This only works if % Sigmas represents a diagonal covariance matrix in the form % of a vector of diagonal entries. for kk = 1:K % Center the image around the means of the kk-th Gaussian % data_centered = double(reshape(image_adjust, D, C)- reshape(Mus(:, kk, :), D, C)); data_centered = double(data(:, :, tt)) - reshape(Mus(:, kk, :), D, C); deviations_squared(:, kk) = sum((data_centered .^ 2) ./ reshape(Sigmas(:, kk, :), D, C), 2); end % Now see if any of the Gaussians matches to within DEVIATION_SQ_THRESH deviations. % The best match for a pixel has min deviation. [junk, index] = min(deviations_squared, [], 2); % From index, whose elements are indeces of the minimum element of each % row, compute an indicator matrix that has a 1 in the corresponding % position of the min element and a 0 otherwise. matched_gaussian = zeros(size(deviations_squared)); matched_gaussian(sub2ind(size(deviations_squared), 1:length(index), index')) = ones(D, 1); matched_gaussian = matched_gaussian & (deviations_squared < DEVIATION_SQ_THRESH); %---------------------------------------------------------------------- %---------------------------Parameters Update-------------------------- % Now update the weights. Increment weight for the selected Gaussian (if any), % and decrement weights for all other Gaussians. Weights = (1 - ALPHA) .* Weights + ALPHA .* matched_gaussian; % Adjust Mus and Sigmas for matching distributions. for kk = 1:K pixel_matched = repmat(matched_gaussian(:, kk), 1, C); pixel_unmatched = abs(pixel_matched - 1); % Inverted and mutually exclusive Mu_kk = reshape(Mus(:, kk, :), D, C); Sigma_kk = reshape(Sigmas(:, kk, :), D, C); Mus(:, kk, :) = pixel_unmatched .* Mu_kk + ... pixel_matched .* (((1 - RHO) .* Mu_kk) + ... (RHO .* double(image))); % Get updated Mus; Sigmas is still unchanged Mu_kk = reshape(Mus(:, kk, :), D, C); Sigmas(:, kk, :) = pixel_unmatched .* Sigma_kk + ... pixel_matched .* (((1 - RHO) .* Sigma_kk) + ... repmat((RHO .* sum((double(image) - Mu_kk) .^ 2, 2)), 1, C)); end % Maintain an indicator matrix of those components that were replaced because no component matched. replaced_gaussian = zeros(D, K); % Find those pixels which have no Gaussian that matches mismatched = find(sum(matched_gaussian, 2) == 0); % A method that works well: Replace the component we % are least confident in. This includes weight in the choice of % component. for ii = 1:length(mismatched) [junk, index] = min(Weights(mismatched(ii), :) ./ sqrt(Sigmas(mismatched(ii), :, 1))); % Mark that this Gaussian will be replaced replaced_gaussian(mismatched(ii), index) = 1; % With a distribution that has the current pixel as mean Mus(mismatched(ii), index, :) = image(mismatched(ii), :); % And a relatively wide variance Sigmas(mismatched(ii), index, :) = ones(1, C) * INIT_VARIANCE; % Also set the weight to be relatively small Weights(mismatched(ii), index) = INIT_MIXPROP; end % Now renormalise the weights so they still sum to 1 Weights = Weights ./ repmat(sum(Weights, 2), 1, K); active_gaussian = matched_gaussian + replaced_gaussian; %---------------------------------------------------------------------- %--------------------------background Segment-------------------------- % Find maximum weight/sigma per row. [junk, index] = sort(Weights ./ sqrt(Sigmas(:, :, 1)), 2, 'descend'); % Record indeces of those each pixel's component we are most confident % in, so that we can display a single background estimate later. While % our model allows for a multi-modal background, this is a useful % visualisa