%function_3

function C = recover_reg(B, psf)
% function C = recover_reg(B, psf) % In input:
% B, matrice dell’immagine sfocata
% psf, matrice della PSF
% In output:
% C, matrice con l’immagine restaurata
%RECOVER_GRADIENT_REG
[m,n,~]=size(B);
% Costante per la regolarizzazione
alpha=0.01;
% PSF "flippata", corrispondente alla moltiplicazione per la matrice trasposta.
psfr = psf(end:-1:1,end:-1:1);
% Right hand side
b = B;
for j = 1 : size(b, 3)
b(:,:,j) = conv2(B(:,:,j), psfr, 'same');
end
% Operatore lineare PSF * PSF come funzione inline
function y = Op(x)
V=reshape(x,m,n,3);
for jj=1:3
V(:,:,jj) = conv2(conv2(V(:,:,jj), psf, 'same'), psfr, 'same');
end
y=reshape(V,m*n*3,1)+alpha*x;
end
% Risolvo il sistema
v=reshape(b,m*n*3,1);
v=pcg(@Op, v, 1.e-10, 500);
C=reshape(v,m,n,3);
end