% Napredne metode digitalne obradbe signala - Projekt
% Compressed sensing - Potpuna rekonstrukcija uz smanjeni broj uzoraka
% Davor Bonaci, davor.bonaci@fer.hr

%% 1. Initialization
clear all;
close all hidden;
pack;

N = 1024;
M = 600;

%% 2. Generation of the matrixes
%psi_inv = dctmtx(N)';
%psi = dctmtx(N);
psi_inv = WaveMtx(N, 3, 'sym4');
psi = pinv(psi_inv);

N2 = size(psi);
N2 = N2(2);

randn('state', 0);
fi = 1/N .* randn(M, N);

theta = fi * psi;

%% 3. Read input image
im_rgb = imread('lenna_oko.tiff');
im_gray = rgb2gray(im_rgb);
image = double(reshape(im_gray, 1, prod(size(im_gray))));
figure, imshow(im_gray);
len = length(image);

%% 4. Compress and decompress
image_est = [];
for i = 1 : floor(len / N)
    signal = image((i - 1) * N + 1 : i*N);
    
    % Performing measurements
    y = fi * signal.';

    % Reconstructing initial signal
    cvx_begin
        variable s_est(N2, 1);
        minimize( norm(s_est, 1) );
        subject to
            theta * s_est == y;
    cvx_end

    signal_est = (psi * s_est).';
    
    image_est = [ image_est signal_est ];
end

%% 6. Verification
im_gray_est = uint8(round(reshape(image_est, sqrt(len), sqrt(len))));

max_error_signal = max(abs(signal - signal_est))

figure,
subplot(2, 1, 1), subimage(im_gray),
subplot(2, 1, 2), subimage(im_gray_est);
