% read the image and from input file
f = im2double(rgb2gray(imread('Vegetation-028.jpg')));
[m, n] = size(f);


%% Ex 1,2 
% Compute sparse matrices Dx_tilde, Dy_tilde such that
% f * Dx_tilde = f_x
%    and
% Dy_tilde * f = f_y:

% create a row-vector containing n many ones.
e = ones(n, 1);

% Put the values of e along the lower diagonal (identified in MATLAB as -1) and
% Put the values of -e along the actual diagonal (0)
% Dx_tilde is of dimension nxn
Dx_tilde = spdiags([e -e], -1:0, n, n); 
% Take care of the Boundary conditions.
Dx_tilde(:, end) = 0;
% compute the x-derivative f_x (Output f_x is a matrix)
f_x = f*Dx_tilde;

% Put the values of -e along the actual diagonal (0) and
% Put the values of e along the upper diagonal (1)
% Dx_tilde is of dimension nxn
e = ones(m, 1);
Dy_tilde = spdiags([-e e], 0:1, m, m);
% boundary conditions
Dy_tilde(end, end) = 0;
% compute the y-derivative f_y (Output f_y is a matrix)
f_y = Dy_tilde*f;


%%
% Visualize f_x and f_y
figure;
imshow(f_x, []);
figure;
imshow(f_y, []);

%% Ex 3
% Apply the formula from the exercise sheet.
Dx = kron(Dx_tilde', speye(m));
Dy = kron(speye(n), Dy_tilde);


%% Ex 4
% vectorize the image.
vec_f = f(:);
% vectorize the x-derivative f_x
vec_f_x = f_x(:);
% According to the formula from the exercise sheet Dx*vec(f) should be the 
% same as vec(f_x).
diff = norm(Dx*vec_f - vec_f_x(:));
['Difference should be zero: ' num2str(diff)]

% Same holds for the y-derivative f_y
diff = norm(Dy*f(:) - f_y(:));
['Difference should be zero: ' num2str(diff)]

%% Ex 5
% Assemble a color gradient using the kronecker product kron and cat and the matrices Dx and Dy.
% (See Wikipedia for a definition of the kronecker product.)
Nabla_c = cat(1, kron(speye(3), Dx), kron(speye(3), Dy));

%% Ex 6
% Read the image from the file
f1_c = double(imread('Vegetation-028.jpg'));

% Compute the color gradient by applying Nabla_c to the vectorized color
% image. The result is a vectorized n*m x 3 x 2 Tensor containing the
% 3 x 2 Color gradients for each pixel.
Nabla_c_f1_vec = Nabla_c * f1_c(:);

% to compute color TV we at each pixel need to compute the Frobenius
% norm of the 3 x 2 Color gradients and sum the values.
% Since the Frobenius norm of a matrix is equal to the 2-norm of the
% vectorized matrix we reshape Nabla_c_f1 to a m*n x 6 matrix...
Nabla_c_f1_mat = reshape(Nabla_c_f1_vec, m*n, 6);

% ...and for each row (corresponding to a pixel) compute the 2-norm of the
% row and...
sq1 = sqrt(sum(Nabla_c_f1_mat.^2, 2));
% ...finally sum over all rows.
TV_c_f1 = sum(sq1);
['TV of Vegetation-028.jpg is ' num2str(TV_c_f1/(m*n))]

% Same for the other image.
f2_c = double(imread('Vegetation-043.jpg'));
Nabla_c_f2 = Nabla_c * f2_c(:);
sq2 = sqrt(sum(reshape(Nabla_c_f2, m*n, 6).^2, 2));
TV_c_f2 = sum(sq2);
['TV of Vegetation-043.jpg is ' num2str(TV_c_f2/(m*n))]

% Same for a constant image.
f3_c = 255 * ones(m, n, 3);
Nabla_c_f3 = Nabla_c * f3_c(:);
sq3 = sqrt(sum(reshape(Nabla_c_f3, m*n, 6).^2, 2));
TV_c_f3 = sum(sq3);
['TV of a constant image should be zero ' num2str(TV_c_f3/(m*n))]