function []=SGimages(file_id,output_id,num_comps,alpha,T,init_w,init_std,min_p_std,max_p_std) 
%SG   Video background subtraction using Stauffer and Grimson [1] method
%   SG(FILE_ID,OUTPUT_ID,NUM_COMPS,ALPHA,T,INIT_W,INIT_STD,MIN_P_STD,MAX_P_STD)
%   INPUT:
%   - FILE_ID: directory\input;
%   - OUTPUT_ID: directory in which the FG output frames will be located;
%   - NUM_COMPS: number of Gaussian components;
%   - ALPHA: learning rate (range: [0,1]);
%   - T: FG threashold (range: [0,1]);
%   - INIT_W: initial weight for the Gaussian components;
%   - INIT_STD: initial variance for the Gaussian components;
%   - MIN_P_STD: minimal std allowed for the Gaussians;
%   - MAX_P_STD: maximal  std allowed for the Gaussians;
%
% OUTPUT:
%   SG saves in output_id directory two files for each frame processed: a
%   (r x c) .bmp image which depicts the FG extracted, and a .mat file which
%   contains: 
%   - Track: (r x c) the FG mask extracted;
%   - w: (r*c x number of Gauss components) the weights associated to each Gaussian component, for each pixel;
%   - mu:(r*c x number of Gauss components) the mean values associated to each Gaussian component;
%   - p_std:(r*c x number of Gauss components) the std values associated to each Gaussian component;
%   - c_matched:(r*c x 1) the index of the Gaussian component matched; 
%
% [1] C. Stauffer and W. Grimson. Adaptive background mixture
% models for real-time tracking. In Int. Conf. Computer Vision
% and Pattern Recognition, volume 2, 1999.

%%% 13/02/2008 - Marco Cristani 

%----------------------
% DATA ALLOCATION
%----------------------


files               =   dir(strcat(file_id,'*.bmp'));
num_images          =   size(files);
[r,c,channels]      =   size(imread(strcat(file_id,files(1).name)));
num_obs             =   r*c;

w   =   single(zeros(num_obs,num_comps)); % Weights,
mu  =   single(ones(num_obs,num_comps)).*-1000; % means,
p_std =   single(ones(num_obs,num_comps)).*100; % variances of the Gaussian components;


c_matched   =   single(zeros(num_obs,1)); % Index of the component which matches with the current observation;
tmp_novel_c =   single(zeros(num_obs,1)); % Indexes of novel components put in the mixtures;
 
Track       =   single(zeros(r,c));       % FG mask;
Im_v        =   single(zeros(r*c,1));

%----------------------
% OUTPUT SETUP
%----------------------

screen      =   get(0,'screensize'); 
size_fig    =   floor(screen(3)/3); % Size of the visual output;
h_original = figure; 
set(h_original,'Position',[10,screen(4)-size_fig-80,size_fig,size_fig],'name','Input frame for BG subtraction');
colormap gray

h_activity = figure; 
set(h_activity,'Position',[size_fig+10,screen(4)-size_fig-80,size_fig,size_fig],'name','FG mask');
colormap gray;

h_contour = figure; 
set(h_contour,'Position',[2*size_fig+10,screen(4)-size_fig-80,size_fig,size_fig],'name','FG contour');


%----------------------
% MAIN CYCLE
%----------------------

for t = 1:num_images
    fprintf('%d|%d)---------------------------------------------\n ',t,num_images);
    Track       =   Track*0;
   
    figure(h_original);
    if channels ==3
    Tmp         =  rgb2gray(imread(strcat(file_id,files(t).name)));
    else
        Tmp         =  imread(strcat(file_id,files(t).name));
    end
    imagesc(Tmp,[0,255]);
    title(strcat('Input frame n:',num2str(t))); axis image
    Im_v        =   single(Tmp(:)); 
    tmp_novel_c =   tmp_novel_c*0;
    c_matched   =   eval_Mix_Gauss(Im_v(:),mu,p_std,num_comps);
            %  Check for selected pixel (locations) the component that matches with the current related
            %  observation;
    if  ~isempty(c_matched)
        pos                                     =   find(c_matched==0);
        num_unmatch                             =   length(pos);
        tmp_novel_c                             =   sum((w(pos,:)>0),2)+1; % Number of used components;
        tmp_novel_c(tmp_novel_c==num_comps+1)   =   num_comps; % Position of the novel components in the mixture;
        pos1D                                   =   couple2single(pos,tmp_novel_c,num_obs,num_comps);
        w(pos1D)                                =   init_w;
        p_std(pos1D)                            =   init_std;
        mu(pos1D)                               =   Im_v(pos); 
        c_matched(pos)                          =   tmp_novel_c; % Update the index vector indicating the matching components;
    end
    [w,mu,p_std]=update_parameters(Im_v(:),w,mu,p_std,c_matched,alpha,num_comps,min_p_std,max_p_std);
    for g=1:num_obs
        Sum(g,1)=sum(w(g,1:(c_matched(g)-1)),2); % Sum the weights of the components up the matched one;
    end
    Track(:)=Sum>T; % Mark in the FG mask the FG pixels;
    [w,mu,p_std]=order_mixture(w,mu,p_std,num_obs,num_comps); % Order the mixtures for w/p_std;

    Track = reshape(Track,r,c); % Final FG;
    
    % -- Lezione 9 - 5 marzo 2008 --

    Track = bwareaopen(Track,30,8);
    Labels = bwlabel(Track,8);
    desc = regionprops(Labels,'Area');
    massimo = -inf;
    if size(desc,1)>=1
        for obj=1:size(desc,1)
            if desc(obj).Area >massimo
                massimo = desc(obj).Area;
                ind_massimo = obj;
            end
        end
        MYarea(t)=desc(ind_massimo).Area;
        Contorno(t).mat = contour_following(Labels==ind_massimo);
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        Qx= Contorno(t).mat(1:4:end,1)'; % nodi in x;
        Qy= Contorno(t).mat(1:4:end,2)'; % nodi in y;
        nB = length(Qx);
        L = nB;
        iter=1;
        xs = zeros(nB*100+1,1);
        ys = zeros(nB*100+1,1);
        %xs(iter)=0;ys(iter)=0;
        s=[0:0.01:nB];
        for n=1:nB
        xs= xs+Qx(n).*qb_spline(mod(s-n,nB));
        ys= ys+Qy(n).*qb_spline(mod(s-n,nB));
        end
        figure(h_contour); set(gcf,'name','Curva 2D chiusa')
        scatter(fliplr(ys),fliplr(xs)); hold on
        scatter(fliplr(Qy),fliplr(Qx),10,[1,0,0],'marker','O','markerfacecolor',[1,0,0])
        axis([0,120,0,150])
        axis ij
        hold off
        
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        figure(h_activity); colormap jet
        imagesc(Labels,[0,10]); 
    end
    
    
    

    drawnow;
    pause(0.001);
    if t<10
        fill_n = '000';
    elseif t<100
        fill_n = '00';
    elseif t<1000
        fill_n = '0';
    else
        fill_n = '';
    end
    output_name_no_ext = strcat(output_id,'\',fill_n,num2str(t));
    imwrite(Track,strcat(output_name_no_ext,'.bmp'),'bmp');
    % Save in output_id a .bmp copy of the current FG frame computed;
    save(strcat(output_name_no_ext,'.mat'),'Track','w','mu','p_std','c_matched');
    % Save in output_id some information on the current FG computation;
end
figure; plot(MYarea);
   % -- FINE Lezione 9 - 5 marzo 2008 --

fprintf('END!!!\n');