%% MCHI, cv9, file 1 (cv9_1) - Photography digitalization
%
% Data compression using SVD - compress a BW picture
%
% Author:       Martin Isoz
% Organisation: ICT Prague
% Date:         13. 11. 2014 (update 04.11.2016)
%
% License: This code is published under MIT License, please do not abuse
% it.
%

%% Clear the environment
close all
clear all
% clc

% rng default


%% Set up wanter image relative accuracy
relE = 1e-2;

%% Load new image
% Picture = imread('sergeant.png','png');                                     %reading image
% Picture = imread('crown.jpg','jpg');                                        %reading image
% Picture = imread('small-black-white-triangles.jpg','jpg');                  %reading image
% Picture = imread('facebookLogo.png','png');                                 %reading image
Picture = imread('hradcany.jpg','jpg');                                     %reading image
% imshow(Picture);

%% Convert the picture to black and white
PictureBW = im2bw(Picture);
% PictureBW = rgb2gray(Picture);                                              %or to grayscale

%% Compare original and the new picture
% figure
% subplot(121)
% imshow(PictureBW)
% title('Black and white')
% subplot(122)
% imshow(Picture)
% title('Original')

disp('The original picture is represented by matrix of a size:')
disp(size(Picture))

disp('The black and white picture is represented by matrix of a size:')
disp(size(PictureBW))

%% Get singular values
% and of images, start of math
A   = double(PictureBW);                                                    %we will work with the black and white picture
% Note: matlab has problems with black and white picture matrix
% represantation, by default, it considers it as a type logical, which
% would cause problems in further calculations. Thus, we need to remind it,
% that it is in fact a double precision matrix


% preparatory calculations
% Note: it is usefull to identify the smaller of the A dimensions
sizeA = size(A);                                                            %get the dimensions of the input matrix
[~,aMin] = min(sizeA);                                                      %identify the minimal dimension of A

% get the smallest possible matrix (small matrix -> SM)
if aMin == 1;
    SM = A'*A;                                                              %m<n => A'*A is of type (n,n)
else
    SM = A*A';                                                              %m>n => A*A' is of type (m,m)
end

ESM= round(1e4*eig(SM))./1e4;                                               %get eigenvalues of the resulting matrix
% Note: from numerical calculation of eigenvalues, I get stuck with some
% wrongly deteremined (small) imaginary parts of them, so I use rounding to
% 4 decimal places to convert all the eigenvalues to real (and I in fact
% ommit the realy small ones...)

sVal = sqrt(ESM);                                                            %singular values are root square of eigenvalues
% Note: eig returns eigenvalues sorted in ascending order, so to keep it
% consistent with our notes, we need to reverse order of sVal
sVal = sort(sVal,1,'descend');

sVal = sVal(sVal~=0);                                                       %reduce sVal to its non-zero entries

disp('In the studied picture, we found')
disp(numel(sVal)')
disp('non-zero singular values.')
disp('Those values are:')
disp(sVal)

%% Find how many members of singular decomposition of matrix A we need

% define relative error of new image compared to the old one
eK = @(k) 1 - sqrt(sum(sVal(1:k).^2)/sum(sVal.^2));

% find how many members I need to keep
k = 1;
while eK(k) > relE && k < numel(sVal)-1;                                    %second condition is there to prevent 'index out of bounds' error
    k = k+1;
end

% truncate the data (vector sVal)
sVal = sVal(1:k);

% output to the screen
disp('To achieve the wanted relative error of')
disp(relE)
disp('we need to keep')
disp(k)
disp('members of singular decomposition of A.')
disp('Thus we need to calculate first')
disp(k)
disp('eigenvectors of ATA and AAT corresponding to the')
disp(k)
disp('biggest singular values of ATA.')

%% Calculate the needed members
% Note: for the clarity, we will base our calculations on the full SVD

% prepare options for the calculation
OPTS.issym = 'true';                                                        %ATA and AAT matrices are symmetric
OPTS.isreal= 'true';                                                        %those matrices are real
OPTS.v0 = ones(size(SM,1),1);                                              %initial vector, supply ones
% OPTS.v0 = round(rand(size(ATA,1),1));
OPTS.v0 = OPTS.v0./norm(OPTS.v0);
OPTS.disp = 0;
OPTS.p  = 2*k;
% for the matrix ATA
[VSM,~,FLAG1] = eigs(SM,[],k,'LM',OPTS);

% get the other relevant matrix
invS = diag(ones(numel(sVal),1)./sVal);                                     %get the inverse of S matrix
if aMin == 1;
    VATA = VSM;
    UAAT = A*VATA*invS;
else
    UAAT = VSM;
    VATA = A'*UAAT*invS;
end


%% Get the approximation matrix
B   = zeros(size(A));                                                       %this time, I actually need this
for i = 1:k
    B   = B + sVal(i).*UAAT(:,i)*VATA(:,i)';
end

%% Get a new image out of it and compare it
PictureCP = B;
PictureCPBW = im2bw(B);

% graphical output
figure('Units','Pixels','Position',[100 50 1000 700],'Visible','on');
subplot(221)
imshow(Picture)
title('Original')
subplot(222)
imshow(PictureBW)
title('Black and White')
subplot(223)
imshow(PictureCP)
title('Compressed')
subplot(224)
imshow(PictureCPBW)
title('Compressed Black and White')
