%%  Separace statickeho a dynamickeho obsahu ve videu
%   Pristup 2: Principal Component Analysis pomoci SVD
%   Autor: Vadym Serediuk

%% Inicializace
clc
clear

% Zadani hodnosti L slozky
r = 1;

% Vyber vstupniho videa:
inputVideoCase = 7;
switch inputVideoCase
    case 1, inputVideo = 'test';
    case 2, inputVideo = 'escalator';
    case 3, inputVideo = 'rowing';
    case 4, inputVideo = 'shop';
    case 5, inputVideo = 'traffic';
    case 6, inputVideo = 'pedestrian';
    case 7, inputVideo = 'pedestrian_noisy';
end

inputVideoFolder = 'input_video';
inputVideoFull = [inputVideo, '.avi'];

videoOutputTrigger = false;         % <- ulozeni jednotlivych videi ve formatu AVI
videoCompilationTrigger = true;     % <- ulozeni kompilaci 4 videi ve formatu AVI
videoStructTrigger = true;          % <- ulozeni matic videi do MAT souboru (pro dalsi zpracovani)
logOutputTrigger = true;            % <- vytvoreni log souboru s detaily vypoctu

outputFolder =  'output_PCA';
outputA =       'PCA_A.avi';
outputL =       'PCA_L.avi';
outputS =       'PCA_S.avi';
outputM =       'PCA_M.avi';
outputComp =    'PCA.avi';
outputStruct =  'PCA.mat';
outputLog =     'PCA_log.txt';

fprintf('Start.\n')
tStart = tic;

%% Import videa
if inputVideoCase <= 5 % nahravani AVI souboru
    [ A, vidWidth, vidHeight ] = importVideoToColumn( ...
        [inputVideoFolder, '\', inputVideoFull] );
else % nahravani simulovaneho videa z MAT souboru
    load([inputVideoFolder, '\pedestrian.mat']);
    switch inputVideoCase
        case 6, A = video.arrays.A * 255;
        case 7, A = video.arrays.A_noise * 255;
    end
    Lorig = video.arrays.L * 255; % L originalni
    Sorig = video.arrays.S * 255; % S originalni
    
    vidWidth = video.parameters.width;
    vidHeight = video.parameters.height;
end

[pixelCount, frameCount] = size(A);

%% Vypocet
% SVD rozklad a nasledny vyber r prvnich singularnich cisel
[U, Sigma, V] = svd(A, 'econ');
L = U(:,1:r)*Sigma(1:r,1:r)*V(:,1:r)';
S = A - L;

% Maska
M = createMask(S);

% Normalizace a vyrovnani jasu matice S (pro graficky vystup)
Str = shiftBrightness (S);

%% Hodnost L a ridkost S, M
rank_L = rank(L);
sparsity_S = length(find(S == 0)) / length(S(:));
sparsity_M = length(find(M == 0)) / length(M(:));

fprintf('rank(L) = %2.0f (rank(A) = %2.0f)\n', rank_L, min(pixelCount, frameCount))
fprintf('sparsity(S) = %2.4f\n', sparsity_S)
fprintf('sparsity(M) = %2.4f\n', sparsity_M)

%% Signal-to-noise Ratio (pouze pro simulovana data)
if inputVideoCase == 6 || inputVideoCase == 7
    SNR_L = 20 * log10( norm(Lorig, 'fro') / norm(Lorig - L, 'fro'));
    
    fprintf('SNR(L) = %2.4f dB\n', SNR_L)
end

%% Ulozeni videi
if videoOutputTrigger == true
    makeVideoFromColumn( A, [outputFolder, '\', inputVideo, '_', outputA] , vidWidth, vidHeight );
    makeVideoFromColumn( L, [outputFolder, '\', inputVideo, '_', outputL] , vidWidth, vidHeight );
    makeVideoFromColumn( Str, [outputFolder, '\', inputVideo, '_', outputS] , vidWidth, vidHeight );
    makeVideoFromColumn( M, [outputFolder, '\', inputVideo, '_',outputM] , vidWidth, vidHeight );
end

if videoCompilationTrigger == true
    makeCompilation( A, L, Str, M, [outputFolder, '\', inputVideo, '_', outputComp] , vidWidth, vidHeight, 1 );
end

%% Ulozeni promennych do MAT souboru
if videoStructTrigger == true
    video = struct;
    video.arrays = struct;
    video.parameters = struct;
    
    video.arrays.A = A;
    video.arrays.L = L;
    video.arrays.S = S;
    video.arrays.Str = Str;
    video.arrays.M = M;
    
    video.parameters.width = vidWidth;
    video.parameters.height = vidHeight;
    video.parameters.frames = frameCount;
    
    save([outputFolder, '\', inputVideo, '_', outputStruct], 'video');
end

%% Log soubor
totalElapsedTime = toc(tStart);

if logOutputTrigger == true
    % inicializace
    fid = fopen( [outputFolder, '\', inputVideo, '_', outputLog], 'w' );
    fprintf( fid, 'Start.\r\n' );
    
    % import videa
    fprintf( fid, 'Imported video:\r\n');
    fprintf( fid, '\tName: \t\t%s \r\n', inputVideo );
    fprintf( fid, '\tVideo size: \t%1.0fx%1.0f px \r\n', vidWidth, vidHeight );
    fprintf( fid, '\tFrame count:\t%1.0f \r\n', frameCount );
    fprintf( fid, '\tMatrix size:\t%1.0fx%1.0f \r\n', pixelCount, frameCount );
    
    % parametry
    fprintf(fid, 'Parameters, constants:\r\n' );
    fprintf( fid, '\tr = rank(L) \t= %1.0f\r\n', r );
    
    % vypocet
    fprintf( fid, 'Computation ... \r\n' );
    fprintf( fid, '\tObservation matrix singular values: \r\n');
    for i = 1:size(Sigma,1)
        fprintf(fid, '\t%3.0f:%12.4f\r\n', i, Sigma(i,i));
    end
    
    fprintf(fid, '\trank(L) = %2.0f (rank(A) = %2.0f)\r\n', rank_L, min(pixelCount, frameCount));
    fprintf(fid, '\tsparsity(S) = %2.4f\r\n', sparsity_S);
    fprintf(fid, '\tsparsity(M) = %2.4f\r\n', sparsity_M);
    if inputVideoCase == 6 || inputVideoCase == 7  
        fprintf(fid, '\tSNR(L) = %2.4f dB\r\n', SNR_L);
    end
    
    % finished
    fprintf(fid, 'Finished! Total elapsed time is %2.2f seconds.', ...
        totalElapsedTime);
    fclose all;
end

%% Konec
fprintf( 'Konec! Celkovy cas %1.2f s.\n', totalElapsedTime )

%% Ukazka ortonormalni baze
% k = 50;
% % sing = U(:,k)*Sigma(k,k)*V(:,k)';
% sing = U(:,k)*V(:,k)';
% sing = (sing - min(sing(:)))/ (max(sing(:)) - min(sing(:)));
% imwrite(  imresize( sing, [200, 200])  , 'base50.jpg')

%% Ukazka singularnich cisel
% figure
% for i = 2:size(Sigma,1)
%     plot([i-1 i], [Sigma(i-1,i-1) Sigma(i,i)], '-k', 'Linewidth', 1.5)
%     hold on
% end
% grid on
% title('Graf velikosti singularnich cisel ve videu c. 4')
% xlabel('i')
% ylabel('sigma_i')


