load clown.mat;
[Nx Ny] = size(X);
figure(1); clf, colormap('gray'); image(X); pause,

% definiamo la risposta al campione unitario del sistema che descrive l'alterazione da sfuocamento:
%
h = zeros(size(X));
R = 4;  % raggio dello sfuocamento
for x=1:Nx
  for y=1:Ny
    if sqrt((x-Nx/2)^2+(y-Ny/2)^2)<=R
      h(x,y) = 1 / (pi*R^2);
    end
  end
end

% ora otteniamo la figura sfuocata (rumorosa):
%
%Xr = conv2(X,h,'same');    %save Xr.mat Xr;
load Xr.mat;
figure(2); clf, colormap('gray'); image(Xr); pause,


% ora applichiamo il filtraggio inverso
%
H = fft2(h);
% creo un limite inferiore "liminf" ai valori di H :
liminf = 0.1;
for x=1:Nx
  for y=1:Ny
    ratio = 1;
    if abs(H(x,y)) < 1e-12
      H(x,y) = liminf;
    elseif abs(H(x,y)) < liminf
      ratio = liminf / abs(H(x,y));
    end
    if ratio > 1
      H(x,y) = H(x,y) * ratio ;
    end
  end
end
% applico il filtro inverso mediante FFT :
Xe = fftshift(real(ifft2( fft2(Xr) ./ H )));
figure(3); clf, colormap('gray'); image(Xe);

