close all; clear
n = 30;
n2 = n*n;

A = zeros(n,n);
A(4:12,10:20) = .5;
A(:,5:6) = .6;
A(20:30,7:15) = .4;
A(10:17,20:27) = rand(8,8);
A(25,25) = 1;
subplot(2,3,1)
imagesc(A)
colormap gray
title('Immagine originale')
pause

h=ones(3);
h(2, 2)=3;


K=costruisci(h, n,n);
K=K*K';
g=K*A(:); %immagine sfuocata
subplot(2, 3, 2);
imagesc(reshape(g, n,n));
title('immagine sfocata');
pause

f=K\g;
subplot(2,3,3)
imagesc(reshape(f, n,n));
title('immagine ricostruita');
pause

r=randn(n^2, 1)*0.1; %rumore gausssiano bianco con varianza 0.1
gr=g+r;

subplot(2, 3, 4);
imagesc(reshape(gr, n, n));
title('sfuocata+rumore');
pause

subplot(2,3, 5)
fr=K\gr;
imagesc(reshape(fr, n, n));
title('inversione con rumore');
pause

[V,D]=eig(full(K));
c=condest(K);
disp(['numero di condizionamento: ', num2str(c)]);

figure(2);
d=diag(D);
[dd, i]=sort(abs(d));
DD=diag(d(i));
VV=V(:, i);
subplot(1, 4, 1)
semilogy(dd);
title('modulo degli autovalori');
pause
%K1=VV*inv(DD)*VV';
r1=VV'*r;
subplot(1,4, 2)
semilogy(abs(r1));
title('rumore');
pause
gg=VV'*g;
subplot(1,4,3);
semilogy(abs(gg));
title('imm. sf. nella base degli autovet.');
pause

ff=VV'*A(:);
subplot(1,4,4);
semilogy(abs(ff));
title('imm. or. nella base degli autovet.');
pause


p=400;
Dinv=inv(DD(p:n^2, p:n^2)); %prendo solo gli autovalori piu' grandi in modulo
Kinv=VV(:, p:n^2)*Dinv*VV(:, p:n^2)';
fr=Kinv*gr;
figure(1);
subplot(2,3, 6)
imagesc(reshape(fr, 30, 30));
title('imagine ric con meno autoval.');