%%  Separace statickeho a dynamickeho obsahu ve videu
%   Pristup 3: Principal Component Pursuit pomoci metody Augmented Lagrangian Multiplier (ALM)
%   Autor: Vadym Serediuk

%% Inicializace
clc
clear

% Vyber vstupniho videa:
inputVideoCase = 1;
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_PCP';
outputA =       'PCP_A.avi';
outputL =       'PCP_L.avi';
outputS =       'PCP_S.avi';
outputM =       'PCP_M.avi';
outputComp =    'PCP.avi';
outputStruct =  'PCP.mat';
outputLog =     'PCP_log.txt';

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

%% 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
% Nastaveni parametru:
% mistar -  korekce mi        ... mistar = 0.01 ~ 0.1 (mensi mi -> vetsi skok na zacatku, ale rychlejsi konvergence, mensi a kratsi 'kopec' - kratka etapa diverdence na zacatku vypoctu)
    mistar = 0.01; 
% miIn -    inkrement mi      ... miIn = 1.00 ~ 1.02   (mirne urychleni konvergence, prilis velka hodnota znehodnocuje kvalitu separace)
    miIn = 1.01;
% lambdastar - korekce lambda ... lambdastar = 0.5 ~ 1.5 (vetsi lambda -> vice agresivni minimalizace ridkosti S, vyssi hodnost L)    
    lambdastar = 0.7;        
% delta -   krit ukonceni     ... delta >= 1e-6, jinak bude konvergovat roky
    delta = 1e-6;
% Maximalni pocet iteraci     ... maxIter < 500
    maxIter = 500;

% Vypocet 
tComp = tic;

mi = mistar * pixelCount * frameCount / (4 * norm(A, 1));
lambda = lambdastar / sqrt( max(pixelCount, frameCount) );

normALS = 0;                            % norm_F(A-L-S)
normA_delta = delta * norm( A, 'fro' ); % delta * norm_F(A)

logMat = zeros(maxIter, 6); % pro ukladani hodnot z iteraci
S = zeros( size(A) );
Y = zeros( size(A) );

fprintf('Iterovani ... \n');
fprintf('Ukonceni pri normF(A-L-S) < delta*normF(A) nebo pri dosazeni %1.0f iteraci.\n', maxIter);
fprintf('   #    normALS    d.normA  desc per it     mi     time\n')
for i = 1 : maxIter
    iterTime = tic;
    
    L = svt( A - S + mi^(-1)*Y , 1/mi );
    S = soft( A - L + mi^(-1)*Y, lambda/mi );
    Y = Y + mi*(A - L - S);
    mi = miIn * mi;
    
    normALSold = normALS; % ukladani predchozi hodnoty norm_F(A-L-S) pro porovnani
    normALS = norm( A-L-S, 'fro' );
    
    iterDur = toc(iterTime);
    
    % Ulozeni hodnot iterace do log matice
    logMat(i,:) = [i, normALS, normA_delta, normALSold-normALS, mi, iterDur];
    
    % Vypis hodnot na obrazovku
    if      (i < 10) || (mod(i,5) == 0)
        fprintf('%4.0f %10.4f %10.4f %10.4f %10.6f %6.2f\n', i, normALS, normA_delta, normALSold - normALS, mi, iterDur)
    end
    
    % Podminka ukonceni cyklu
    if normALS <= normA_delta, break, end
end

durCalc = toc(tComp);
durIterationAverage = durCalc/i;

fprintf( 'Konec vypoctu, cas %1.2fs.\n', durCalc )
fprintf( 'Prumerny cas jedne iterace %1.2fs.\n', durIterationAverage )

%% Postprocessing
% 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

%% Statistika pameti
mem = memory;
memoryMaxPossibleArray = floor( mem.MaxPossibleArrayBytes/(2^20) );
memoryAvailableAllArrays = floor( mem.MemAvailableAllArrays/(2^20) );
memoryUsedMATLAB = floor( mem.MemUsedMATLAB/(2^20) );

%% Statistika casu vypoctu
totalElapsedTime = toc(tStart);
totalHour = floor( totalElapsedTime / 3600 );
totalMin = floor( (totalElapsedTime - 3600*totalHour) / 60 );
totalSec = totalElapsedTime - 3600*totalHour - 60*totalMin;

%% Log soubor
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, '\tmi^star \t= %1.6f\r\n', mistar );
    fprintf( fid, '\tmiIn \t= %1.6f\r\n', miIn );
    fprintf( fid, '\tlambda^star \t= %1.6f\r\n', lambdastar );
    fprintf( fid, '\tlambda \t= %1.6f\r\n', lambda );
    fprintf( fid, '\tdelta \t= %1.0e\r\n', delta );
    fprintf( fid, '\tmaxIter \t= %1.0f\r\n', maxIter );
    
    % vypocet
    fprintf( fid, 'Computation ... ' );
    fprintf( fid, 'Terminating condition: normF(A-L-S) < delta*normF(A)\r\n' );
    fprintf( fid, '\t# it. \tnormF(A-L-S)\tdelta*normF(A)\tdescent per it. \t\tmi \t\t time\r\n' );
    for i = 1 : maxIter
        if logMat(i,1) == 0
            fprintf(fid, '\tConverged!!!\r\n');
            break
        end
        fprintf( fid, '\t%4.0f \t%12.4f \t%12.4f \t%12.4f \t%12.6f %8.2f\r\n', logMat(i,:) );
    end
    if (i == maxIter), fprintf(fid, '\tIteration limit reached!!!\r\n'); end
    fprintf(fid, '\tCompleted. Elapsed time is %2.2f seconds.\r\n', durCalc);
    fprintf(fid, '\tAverage iteration time is %1.2f seconds.\r\n', durIterationAverage);

    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
    
    % memory
    fprintf(fid, 'Memory stats:\r\n');
    fprintf(fid, '\tMaximum possible array: \t\t%1.0f MB\r\n', memoryMaxPossibleArray);
    fprintf(fid, '\tMemory available for all arrays: \t%1.0f MB\r\n', memoryAvailableAllArrays);
    fprintf(fid, '\tMemory used by MATLAB: \t\t%1.0f MB\r\n', memoryUsedMATLAB);
    
    % finished
    fprintf(fid, 'Finished! Total elapsed time is %2.2f seconds (%1.0fh %1.0fm %1.2fs).', ...
        totalElapsedTime, totalHour, totalMin, totalSec);
    fclose all;   
end

%% KONEC PROGRAMU
fprintf( 'Konec! celkovy cas %1.0fh %1.0fm %1.2fs.\n', totalHour, totalMin, totalSec )


