% demo_recon_time_pcs.m
%   video encoder and deocder using parallel compressed sensing with
%   zigzag-scan-based permutation
%
% Written: Hao Fang
% Created: Apr. 2013

%=======================================
clear all;
close all;
clc;

path(path, '../Videos');
path(path, '../Tools');
path(path, '../Tools/YUV');
path(path, '../Tools/ZIGZAG_PERM');
path(path, '../Tools/ZIGZAG_SCAN');
path(path, '../ParallelCS');

% akiyo_qcif
% carphone_qcif
% claire_qcif
% coastguard_qcif
% foreman_qcif
% salesman_qcif
filename = 'akiyo_qcif';
yuvformat = 'qcif';

%======================================

switch yuvformat
    case 'cif'
        dims = [352 288];
        block_size = 16;
    case 'qcif'
        dims = [176 144];
        block_size = 8;
end
numfrm = 2;
startfrm = 0;
[Y, U, V] = yuv_import(sprintf('%s.yuv', filename), dims, numfrm, startfrm, yuvformat);

X_Ref = dct2(double(Y{1})-128); % range of Y-128: [-128, 127]
max_X_Ref = max(abs(X_Ref(:)));
X_Ref = X_Ref ./ max_X_Ref; % normalize since l1magic won't work for too large max(abs(x))
X_Non_Ref = double(Y{1}-Y{2});
max_X_Non_Ref = max(abs(X_Non_Ref(:)));
X_Non_Ref = X_Non_Ref ./ max_X_Non_Ref;

% parameters
msr_ratio = 0.8;
N = dims(2);
ave_CR = 0.1:0.1:0.5;
iter = 200;
ave_rec_time = zeros(1, length(ave_CR));
ave_PSNR_Ref = zeros(1, length(ave_CR));
ave_PSNR_Non_Ref = zeros(1, length(ave_CR));
for jj = 1:length(ave_CR)
    CR_Ref = ave_CR(jj)*2*msr_ratio;
    CR_Non_Ref = ave_CR(jj)*2*(1-msr_ratio);
    K_Ref = ceil(N*CR_Ref);
    K_Non_Ref = ceil(N*CR_Non_Ref);
    
    PSNR_Ref = zeros(1, iter);
    PSNR_Non_Ref = zeros(1, iter);
    rec_time = zeros(1, iter);
    for ii = 1:iter
        A_Ref = sqrt(1/K_Ref) * randn(K_Ref, N);
        A_Ref = orth(A_Ref')';
        A_Non_Ref = sqrt(1/K_Non_Ref) * randn(K_Non_Ref, N);
        A_Non_Ref = orth(A_Non_Ref')';
        
        % parallel cs encoding
        Z_Ref = parallel_cs_encoder(X_Ref, A_Ref, @perm_zigzag);
        Z_Non_Ref = parallel_cs_encoder(X_Non_Ref, A_Non_Ref);
        
        % parallel cs decoding
        t = cputime;
        rec_X_Ref = parallel_cs_decoder_l1magic(Z_Ref, A_Ref, @invperm_zigzag);
        rec_X_Non_Ref = parallel_cs_decoder_l1magic(Z_Non_Ref, A_Non_Ref);
        rec_Y_Ref = mod(idct2(rec_X_Ref.*max_X_Ref) + 128, 255);
        rec_Y_Non_Ref = rec_Y_Ref + rec_X_Non_Ref.*max_X_Non_Ref;
        rec_time(ii) = cputime - t;
        
        PSNR_Ref(ii) = psnr(double(Y{1}), double(rec_Y_Ref));
        PSNR_Non_Ref(ii) = psnr(double(Y{2}), double(rec_Y_Non_Ref));
    end
    ave_rec_time(jj) = mean(rec_time);
    ave_PSNR_Ref(jj) = mean(PSNR_Ref);
    ave_PSNR_Non_Ref(jj) = mean(PSNR_Non_Ref);
end