sumimg = zeros(240, 320);
%Create a matrix to store the sum of the images.
for a=1:size(ag, 3)
sumimg = plus(sumimg, ag(:, :, a));
%Populate the matrix of the sum of images.
end
BG1 = sumimg ./ size(ag, 3);
% Calculate the background image.
BGstd = std(ag, 0, 3);
%Calculate the standard deviation.
vector_v = zeros(1, size(ag, 3));
%Create a vector to find the number of foreground pixels.
for b=1:size(ag,3)
    diffimg = minus(abs(ag(:, :, b) - BG1), (BGstd .* 2));
    %Calculate the difference image.
    for c=1:240
        for d=1:320
            if diffimg(c, d)<0
                diffimg(c, d)=0;
            else
                diffimg(c, d)=1;
            end
        end
        frgimg(:, :, b) = diffimg;
        %Fill the difference image matrix for each image.
        ones_vec = sum(sum(diffimg));
        %Find the number of ones in each difference image.
    end
            vector_v(b) = ones_vec;
            %Find the vector that contains the number of foreground pixels
            %in each frame.
end
avgv = sum(vector_v) ./ size(ag, 3);
%Calculate the average number of foreground pixels.
newvec_v = vector_v - avgv;
%Find the vector that shows the deviation from the average foreground pixel
%value.
for e=1:size(ag, 3)
    if newvec_v(b)<0
        newsample(:, :, e)=ag(:, :, e);
%The new sample with reduced foreground activity.
    else
        newsample(:, :, e)=zeros(240, 320);
    end
end
sumimg1 = zeros(240, 320);
%Create a matrix to store the sum of the images.
for f=1:size(ag, 3)
    sumimg1 = plus(sumimg1, newsample(:, :, f));
end
%Populate the matrix of the sum of images.
newvec = zeros(1, size(ag, 3));
for g=1:size(ag, 3)
    if newvec_v(g)<0;
        newvec(g) = 0;
    else
        newvec(g) = 1;
    end
end
lngthsmpl = size(ag, 3) - sum(newvec);
%Finding the number of images with lower foreground activity than average.
BG2 = sumimg1 ./ lngthsmpl;
%Calculate the background image with frames of reduced foreground activity.
%So far all we have done is taken two different average images to get the 
%foreground. Taking the next step we move to the spatial filtering methods
%for the images.
%We start with convolution, and use the Laplacian to get the gradient.
%I used two Laplacians, one with positively weighed center, and 
%the other with negative center. The difference comes on whether we add
%the convoluted image or subtract it from the original image. But due to 
%memory problems I have limited the code to include only one of them.
L1 = [0 1 0; 1 -4 1; 0 1 0];
%L2 = [-1 -1 -1; -1 8 -1; -1 -1 -1];
%The convolution can be done in two ways. We can either apply the Laplacian
%to the average image we obtained above, or we can apply the Laplacian to 
%all the frames and then take the average. I will show both the results, 
%for this case, but after repeated executions, I think that applying the 
%operators to the average background image gives better results. I don't 
%know if this is true in general, or it is happening with these images only
%or there is some flaw in my code!
convimg1 = zeros(240, 320);
convimg2 = zeros(240, 320);
convimg3 = zeros(240, 320);
for h=1:size(ag, 3)
   agnew(:, :, h) = conv2(ag(:, :, h), L1, 'same');
   %Calculate the convoluted frames.
   agnewnormal(:, :, h) = (agnew(:, :, h) - min(agnew(:))) / (max(agnew(:))-min(agnew(:)));
   %Normalize.
   convimg12(:, :, h) = plus(ag(:, :, h), agnew(:, :, h));
   %Find the new image as the sum of original image and convolved image.
   convimg13(:, :, h) = plus(ag(:, :, h), agnewnormal(:, :, h));
   %Do the same with normalized matrix.
end
for i=1:size(ag, 3)
    convimg1 = plus(convimg1, agnew(:, :, i));
    convimg2 = plus(convimg2, convimg12(:, :, i));
    convimg3 = plus(convimg3, convimg13(:, :, i));
end
convimg111 = convimg1 ./ size(ag, 3);
%Average convolved image.
convimg112 = convimg2 ./ size(ag, 3);
%Average transformed image with no normalization.
convimg113 = convimg3 ./ size(ag, 3);
% Average transformed image with normalization.
% Now applying the convolution to the average background image.
convimg14 = conv2(BG2, L1, 'same');
%Applying the Laplacian to the background image.
convimg15 = (convimg14 - min(convimg14(:)))/(max(convimg14(:)) - min(convimg14(:)));
%Normalize.
convimg114 = plus(convimg14, BG2);
%Output image.
convimg115 = plus(convimg15, BG2);
%Normalized output image.
convimg16 = edge(BG2, 'prewitt');
%Finding the Prewitt edge.
convimg116 = plus(convimg16, BG2);
%Adding the edge to the background.
A = [1 1 1 1; 1 1 1 1; 1 1 1 1; 1 1 1 1];
%Taking an averaging filter.
A1 = A .* (1/16);
%Normalizing.
convimg117 = conv2(convimg116, A1, 'same');
%Applying the averaging mask.
convimg118 = BG2 .* convimg117;
%Taking the product image.
convimg119 = plus(BG2, convimg118);
%Adding the result to the background image.
convimg120 = minus(BG2, convimg118);
%Subtracting the result from the background image.
%The following process gave the best result for the movie. I'm not sure if
%it was because of the structure of the movie or the way I applied the
%transforms!
%High Pass Filtering followed by the application of the median filter.
for j=1:size(ag, 3)
    aghp(:, :, j) = minus(ag(:, :, j), BG1);
    %Background image subtracted from each frame.
end
imgtoshow = zeros(240, 320);
for k=1:size(ag, 3)
    imgtoshow = plus(imgtoshow, aghp(:, :, k));
    %Sum of the frames from where the background image has been subtracted.
end
hpimg1 = imgtoshow ./ size(ag, 3);
%Average of those frames.
hpimg2 = plus(imgtoshow, BG1);
%High pass image.
%Median Filter
%I will be using the hybrid median filter (Nieminen et al, 1987) as well as
%the normal median filter supplied in MATLAB
%The hybrid median filter works on the subwindows(5*5)
% |E - D - C|
% |- E D C -|
% |F F A B B|
% |- G H I -|
% |G - H - I|
% and the center A is replaced the value A'
% A' = median(A,median(A,B,D,F,H),med(A,C,E,G,I))
PIM = padarray(hpimg2,[2,2],'replicate');
%Padded image, with hpimg2 being padded to ensure the application of a 5 by 5
%filter.
dim = size(PIM);
%Assign the size of the padded image for the counter variable.
for m = (1+2):(dim(1,1)-2)
    for n = (1+2):(dim(1,2)-2)
        % parameters for filtering
        a = PIM(m,n);
        b1 = PIM(m,n+1);
        b2 = PIM(m,n+2);
        c1 = PIM(m-1,n+1);
        c2 = PIM(m-2,n+2);
        d1 = PIM(m-1,n);
        d2 = PIM(m-2,n);
        e1 = PIM(m-1,n-1);
        e2 = PIM(m-2,n-2);
        f1 = PIM(m,n-1);
        f2 = PIM(m,n-2);
        g1 = PIM(m+1,n-1);
        g2 = PIM(m+2,n-2);
        h1 = PIM(m+1,n);
        h2 = PIM(m+2,n);
        i1 = PIM(m+1,n+1);
        i2 = PIM(m+2,n+2);
        para1 = [a,b1,b2,d1,d2,f1,f2,h1,h2];
        %Making a subsequence.
        para2 = [a,c1,c2,e1,e2,g1,g2,i1,i2];
        %Making the second subsequence.
        para3 = [a,median(para1),median(para2)];
        %Applying the modified median form.
        MIM(m-2,n-2) = median(para3);
        %Output filtered image
    end
end
M = medfilt2(hpimg2,[5,5]);
%Filtered image by MATLAB's median filter for comparison
figure
subplot(4,4,1, 'align')
imshow(BG1);
title('Average Background Image')
subplot(4,4,2, 'align')
imshow(BG2);
title('With Reduced Foreground Activity')
subplot(4,4,3, 'align')
imshow(convimg111);
title('Laplacian Applied to All Frames')
subplot(4,4,4, 'align')
imshow(convimg112);
title('Laplacian Applied to All Frames:Final')
subplot(4,4,5, 'align')
imshow(convimg113);
title('Laplacian Applied to All Frames:Normalized')
subplot(4,4,6, 'align')
imshow(convimg114);
title('Laplacian Applied to Background')
subplot(4,4,7, 'align')
imshow(convimg115);
title('Laplacian Applied to Background:Normalized')
subplot(4,4,8, 'align')
imshow(convimg116);
title('Image with Prewitt Edge')
subplot(4,4,9, 'align')
imshow(convimg117);
title('4 by 4 Averaging Filter')
subplot(4,4,10, 'align')
imshow(convimg118);
title('Sum Image')
subplot(4,4,11, 'align')
imshow(convimg119);
title('Difference Image')
subplot(4,4,12, 'align')
imshow(hpimg1);
title('Result of High Pass Filtering')
subplot(4,4,13, 'align')
imshow(hpimg2);
title('Final High Pass Filter Image')
subplot(4,4,14, 'align')
imshow(M);
title('Normal Median Filter')
subplot(4, 4, 15, 'align')
imshow(MIM);
title('Hybrid Median Filter')