clear;clf;

% Input and display the binary image
N = 256;
I = readbin('lenna.256',N,N);
figure(1)
imshow(I,gray(256))

% Normalize the image: zero mean and unit energy
mI = mean2(I);
NI = I-mI;
ENI = norm(NI,'fro');
x = NI/ENI;

% Blue the image, corrupt the image using WGN and ddisplay it
% h is the blurring filter, and sigma is the noise std
h = ones(4,4)/16;
sigma = sqrt(4*10^(-7));

Xf = fft2(x);
Hf = fft2(h,N,N);
y = real(ifft2(Hf.*Xf))+sigma*randn(N,N);

figure(2)
Ih = y*ENI+mI;
imshow(Ih,gray(256))

% Estimate the power spectrum of the original image
Yf = fft2(y);
% Pxf = Yf*Yf'/N^2;
% Pxf = Xf*Xf'/N^2;
Pxf = Xf.^2/N^2;

% Regularized inverse filter, when alpha = 1, it is the Wiener filter
%alpha = 1;
alpha = 0.27;
Gf = ((Hf.')').*abs(Pxf)./(abs(Hf.^2).*abs(Pxf)+alpha*sigma^2);

% Restorated image without denoising
EXf = Gf.*Yf;
EI = real(ifft2(EXf))*ENI+mI;
figure(3)
imshow(EI,gray(256))

% Wavelet denoising using Daubechies tap 8 wavelet transform
level = 6;
[LO_D,HI_D,LO_R,HI_R] = wfilters('db4');
h = LO_R;

[Idh] = round(denoise(EI,h,0,[0 3.0 0 1 0 0]));
Idh = Idh.*(0<Idh & Idh<256)+(Idh>255)*255;
figure(4)
imshow(Idh,gray(256))

return
