%%%% 2D Fourier Transform

y=im2double(imread('image_lena512.png'));
figure(1), imshow(y),title('original image');

Y=fft2(y);

%Y=fftshift(fft2(y));

% Y(1,1) = 0;

figure(4), imshow(real(Y),[]),title('real part of Fourier Transform'),impixelinfo;
figure(5), imshow(imag(Y),[]),title('imag part of Fourier Transform'),impixelinfo;

y=ifft2(Y);
figure(6), imshow(y,[]);


%% the basis element
Y = zeros(512);
u = 5
v = 100
Y(u , v) = 1;
% 
% u=1
% v=1
% Y(u,v)=1;


figure(1),imshow(Y,[]),title(['the u=',num2str(u), 'v=',num2str(v), 'space domain basis element'] )

y=ifft2(Y);

figure(2),imshow(real(y),[]),title(['the u=',num2str(u), 'v=',num2str(v), 'Fourier domain basis element'] )
saveas(2,['fourier_',num2str(u),'_',num2str(v),'.png'],'png');
% y=ifft2(fftshift(Y));


% figure(6), imshow(real(y),[]);



%% low pass
sigma_noise = 0;

y = im2double(imread('image_lena512.png'));
y = y + sigma_noise * randn(size(y)); 
figure(1), imshow(y , []), title('original image');

MASK = ones(size(y));

% MASK=zeros(size(y));
bounds=220;

MASK(1 : bounds , :) = 0;
MASK(end - bounds : end , :) = 0;
MASK(:  ,1 : bounds) = 0;
MASK(: , end - bounds : end) = 0;

% MASK=MASK';

%
%MASK=fspecial('gaussian',size(y,1),10); %gaussian
figure(2), imshow(MASK,[]),title('mask')

Y = fft2(y);

Y = fftshift(fft2(y));

% Y_low= Y.*MASK;
Y_low= fftshift(Y .* MASK);

% Y_low= (Y.*MASK);
y_low=real(ifft2(Y_low));

figure(3), imshow(y_low,[]), title(['lowpassed lena']);


%what about the filter?
mask=fftshift(ifft2(MASK));
figure(4), imshow(real(mask),[]),title('low pass filter'),impixelinfo


%% Deconvolution

%% Out of Focus Blur
h_size=7;
sigma_gauss=2;
h = fspecial('gaussian',h_size,sigma_gauss)
figure(3),bar3(h),title('kernel')


%% Point Spread Function

h_size=64;
sigma_gauss=10;
h = fspecial('gaussian',h_size,sigma_gauss);

% display
figure(3),imshow(fftshift(real(fft2(h)))),title('Fourier Transform of the Kernel')
figure(4),imshow(h,[]),title('kernel')

%% Parametric Blur PSF is given by two parameters 
y=im2double(imread('image_lena512.png'));

h=ones(1,21);
h=h/sum(h(:));
z = conv2(y,h,'same');


% display
figure(1), imshow(y),title('original')
figure(2), imshow(z),title('observation')
figure(3), bar3(h),title('kernel')




% 
len_h=15;
theta_h=30;
h = fspecial('motion',len_h,theta_h);

z = conv2(y,h,'same');


% display
figure(1), imshow(y),title('original')
figure(2), imshow(z),title('observation')
figure(3), bar3(h),title('kernel')



%% deblurring
% 
y = im2double(imread('image_lena512.png'));
sigma_noise = 0.00001;
h = fspecial('gaussian' , 30 , 100);

% h = fspecial('motion',30,10);

[yN , xN] = size(y);
size_z_1 = yN;
size_z_2 = xN;
[ghy , ghx] = size(h);

% BLURRING has to be performed by periodic padding of the image out of the image boundaries.
% therefore we compute blurred image directly in fourier domain, by pixel-wise multiplication 
% of the image coefficients against a PSF (which is padded to have the same size of the original image)
big_v = zeros(yN,xN); % pad PSF with zeros to whole image domain, and centers it.
big_v(1 : ghy , 1 : ghx) = h;
big_v = circshift(big_v , -round([(ghy - 1) / 2 (ghx - 1) / 2]));
V = fft2(big_v); % Frequency response of the PSF
% performs blurring (convolution is obtained by product in frequency domain)
y_blur = real(ifft2(V .* fft2(y))) + sigma_noise * randn(size(y)); 

Y_blur = fft2(y_blur);
% RI= conj(V)./( (abs(V).^2) +epsStar^2);
% RI(1,1)=conj(V(1,1))./(abs(V(1,1)).^2);
epsStar = 0.000001
z_ri=real(ifft2(Y_blur .* conj(V) ./ (abs(V).^2 +epsStar)));


figure(1), imshow(y,[]),title('original')
figure(2), imshow(y_blur,[]),title('observation')
figure(7), imshow(z_ri,[]),title('restored image');
figure(8), imshow(abs(z_ri-y),[]),title([' Absolute error, max ',num2str(max(abs(y(:)-z_ri(:))))]);


%%
% 
% [yN,xN]=size(y);
% [ghy,ghx]=size(h);
% 
% big_h=zeros(yN,xN);
% big_h(1:ghy,1:ghx)=h;
% big_h=circshift(big_h,-round([(ghy-1)/2 (ghx-1)/2]));
% 
% % Frequency response of the PSF
% H=fft2(big_h);
% % performs blurring (convolution is obtained by product in frequency domain)
% y_blur=real(ifft2(H.*fft2(y))) + sigma_noise*randn(size(y));
% 
% % inverse
% % RI=1./H;
% 
% epsStar =0;
% RI= conj(H)./( (abs(H).^2) +epsStar );
% 
% Z_RI=RI.*fft2(y_blur);
% 
% y_hat=real(ifft2(Z_RI));


%% wavelet Transform
y=im2double(imread('image_lena512.png'));
sigma_noise=0.1;
z=y+sigma_noise*randn(size(y));

[cA1,cH1,cV1,cD1] = dwt2(z , 'db4');

T = sigma_noise*sqrt(2*log(size(y,1)*size(y,2)));

WT=[cA1,cH1;cV1,cD1];

figure(1), imshow(z,[])
figure(2), imshow(WT,[])

 cA1_hat=cA1;
% cA1_hat(find(abs(cA1)<T))=0;

cH1_hat=cH1;
cH1_hat(find(abs(cH1)<T))=0;

cV1_hat=cV1;
cV1_hat(find(abs(cV1)<T))=0;

cD1_hat=cD1;
cD1_hat(find(abs(cD1)<T))=0;

WT_hat=[cA1_hat,cH1_hat;cV1_hat,cD1_hat];

y_hat = idwt2(cA1_hat,cH1_hat,cV1_hat,cD1_hat,'db4');


figure(2), imshow(y_hat,[])
figure(3), imshow(y,[])