load SurvTestFrames.mat

% M = [];
% 
% for i=1:size(ag,3)
%     imshow(ag(:,:,i));
%     M = [M; getframe()];
% end

BG1 = mean(ag,3);
BGstd = std(ag,0,3);
FG1 = abs((ag - repmat(BG1,[1 1 size(ag,3)]))) > 2*repmat(BGstd,[1 1 size(ag,3)]);

% Mfg = [];
% for i=1:size(ag,3)
%     imshow(FG1(:,:,i));
%     Mfg = [Mfg; getframe()];
% end

v = sum(sum(FG1,1),2); v=v(:)';

bg2ind = find(v<mean(v));
BG2 = mean(ag(:,:,bg2ind),3);

chmask = mod(repmat([1:320],240,1)+repmat([1:240]',1,320), 2);

uiwait(msgbox('Following movie shows a) original, b) foreground mask, c) original with foreground mask, and d) vertical bar indicator of percentage of foreground pixels','','non-modal'));

maxverbar = size(ag,1)*3;
percfg = round(maxverbar*v/max(v));

for i=1:size(ag,3)
    imshow([ [ag(:,:,i) ; FG1(:,:,i) ; ag(:,:,i).*(1-FG1(:,:,i)) + chmask.*FG1(:,:,i)] [zeros(maxverbar - percfg(i),20); ones(percfg(i),20)]]);
    
    pause(0.03);
end

uiwait(msgbox('Following images are BG1, BG2, and the difference (not equalized, but emphasized: (BG2-BG1)*4+0.5)\nWe can clearly see the "ghost"','','non-modal'));

close all;
subplot(3,1,1); imshow(BG1);
subplot(3,1,2); imshow(BG2);
subplot(3,1,3); imshow((BG2-BG1)*4+0.5);

FG2 = abs((ag - repmat(BG2,[1 1 size(ag,3)]))) > 2*repmat(BGstd,[1 1 size(ag,3)]);
v2 = sum(sum(FG2,1),2); v2=v2(:)';

uiwait(msgbox('Vectors v (blue) and v2 (green); v2 indicates less movement when there is no movement, and more movement when there is actual movement','','non-modal'));

plot([v; v2]','o:');

uiwait(msgbox('Following movie shows a) foreground mask #1, b) original with foreground mask #1, c) foreground mask #2, d) original with foreground mask #2','','non-modal'));


for i=1:size(ag,3)
    imshow([ [FG1(:,:,i)  ag(:,:,i).*(1-FG1(:,:,i)) + chmask.*FG1(:,:,i)] ; [FG2(:,:,i) ag(:,:,i).*(1-FG2(:,:,i)) + chmask.*FG2(:,:,i)] ]);
    
    pause(0.03);
end

BG3 = sum((1-FG1).*ag,3)./sum(1-FG1,3);

uiwait(msgbox('Following image is BG3 (individual foreground pixels are discarded, not the whole images)','','non-modal'));

imshow(BG3);

bg4ind = find(v<100);
BG4 = mean(ag(:,:,bg4ind),3);

uiwait(msgbox('Following image is BG4 (similar to BG2, but much smaller threshold)','','non-modal'));

imshow(BG4);

uiwait(msgbox('Finally: [BG1 BG2; BG3 BG4]','','non-modal'));

subplot(2,2,1); imshow(BG1)
subplot(2,2,2); imshow(BG2)
subplot(2,2,3); imshow(BG3)
subplot(2,2,4); imshow(BG4);


% BG5 = zeros(size(ag,1),size(ag,2));
% binc = 0.005:0.01:0.995;
% for i=1:size(ag,1)
%     for j=1:size(ag,2)
%         h = hist(ag(i,j,:),binc);
%         h1 = h;
%         w1 = 3;
%         for k=1:length(binc)
%             i1=max(1,k-w1);
%             i2=min(length(binc),k+w1);
%             h1(k) = sum(h(i1:i2))/(i2-i1+1);
%         end
%         [hmax, imax] = max(h1);
%         BG5(i,j) = binc(imax); 
%     end
%     i
% end
