%--------------------------------------------------------------------------
% ENSEEIHT - 2SN MM - Traitement des donnˇes audio-visuelles
% TP6 - Restauration d'images
% exercice_0 : dˇbruitage avec mod¸le de Tikhonov (niveaux de gris)
%--------------------------------------------------------------------------

clear
close all
clc

% Mise en place de la figure pour affichage :
taille_ecran = get(0,'ScreenSize');
L = taille_ecran(3);
H = taille_ecran(4);

% Lecture de l'image :
u0 = double(imread('Images/lena_avec_bruit.bmp'));
image = zeros(size(u0));
[nb_lignes,nb_colonnes,nb_canaux] = size(u0);
u_max = max(u0(:));

% Affichage de l'image :
% subplot(1,2,1)
figure;
	imagesc(max(0,min(1,u0/u_max)),[0 1])
	colormap gray
	axis image off
	title('Image bruitee')
    imwrite(max(0,min(1,u0/u_max)), "saves/exo0bis_bruit.png")
    
% Vectorisation de u0 :
nb_pixels = nb_lignes * nb_colonnes;

for c=1:nb_canaux

    ut = u0(:,:,c);
    ut = reshape(ut, nb_pixels, 1);
    
    % Opˇrateur gradient :
    e = ones(nb_pixels,1);
    Dx = spdiags([-e e],[0 nb_lignes],nb_pixels,nb_pixels);
    Dx(end-nb_lignes+1:end,:) = 0;
    Dy = spdiags([-e e],[0 1],nb_pixels,nb_pixels);
    Dy(nb_lignes:nb_lignes:end,:) = 0;
    
    % Second membre b du syst¸me :
    b = ut;
    
    % Matrice A du syst¸me :
    Lap = -Dx'*Dx -Dy'*Dy;
    lambda = 2;				% Poids de la rˇgularisation
    A = speye(nb_pixels) - lambda*Lap;
    
    % Rˇsolution du syst¸me A*u = b :
    u = A\b;
    image(:,:,c) = max(0,min(1,reshape(u,[nb_lignes nb_colonnes])/u_max));
end


% Affichage de l'image restaurˇe :
% subplot(1,2,2)
figure;
	imagesc(image, [0 1])
	colormap gray
	axis image off
	title('Image restauree')
    imwrite(image, "saves/exo0bis_clean.png")