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


f0 = reshape(A,n2,1);
h=[1 1 1; 1 2 1; 1 1 1]./10;
K=costruisci(h, n,n);
%K=K'*K;
gnf = K*f0;
subplot(2,3,2)
imagesc(reshape(gnf,n,n));
title('Immagine sfocata')
pause

noise=randn(n2,1)*0.01;
g=gnf+noise;
subplot(2,3,3)
imagesc(reshape(g,n,n));
title('Immagine sfocata+rumore')
pause

finv=K\gnf;
subplot(2,3,5)
imagesc(reshape(finv, n,n));
title('Imm. ric. senza noise');
pause

fpinv=K\g;
subplot(2,3,6)
imagesc(reshape(fpinv, n, n));
title('Imm ric. con noise');


figure
mu=linspace(10^-6,0.1);
nf0=norm(f0);
ng=norm(g);
I=eye(n2);
b=K'*g;

D=norm(K*f0-g)
for i=1:100
    R=K'*K+mu(i)*I;
    fmu=R\b;
    Err(i)=norm(fmu-f0);
    Dis(i)=norm(K*fmu-g);
    En(i)=norm(fmu);
end
subplot(1,3,1)
plot(mu, Err);
title('err. ricostr.')
pause

subplot(1,3,2)
plot(mu, Dis, mu, D, 'b');
title('discrep');
pause

subplot(1,3,3)
plot(mu, En, mu, nf0, 'r');
title('Energia')
pause

[Erropt, iop]=min(Err);

muopt=mu(iop)
Eopt=Erropt
R=K'*K+muopt*I;
fmuopt=R\b;
figure;
colormap gray
subplot(2,2,1)
imagesc(reshape(fmuopt,n,n));
title('con mu_opt')
pause


id=find(Dis<=D);
md=mu(id(end))
R=K'*K+md*I;
fd=R\b;
subplot(2,2, 2)
imagesc(reshape(fd, n,n));
title('con mu_discr');
Ed=norm(fd-f0)
pause

E=nf0;

ie=find(En<=E);
me=mu(ie(1));
R=K'*K+me*I;
fe=R\b;
subplot(2,2, 3)
imagesc(reshape(fe, n,n));
title('con mu_energia');
Ee=norm(fe-f0)
pause


M=D^2/E^2
R=K'*K+M*I;
fM=R\b;
subplot(2,2, 4)
imagesc(reshape(fM, n,n));
title('con mu_Miller');
Em=norm(fM-f0)
pause