clear all, close all

I = im2double(imread('cameraman.tif'));
figure, imshow(I); title('original image')
H = fspecial('motion',20,40);
% H = fspecial('gaussian',[19,19],3);

fftH=psf2otf(H,size(I));
I_blur= imfilter(I,H,'conv','circular');


% deconvolution directe

sigma=[0,0.000001,0.0001,0.001,0.01,0.1];

I_blurnoise0= imnoise(I_blur,'gaussian',0,0);
figure, imshow(I_blurnoise0); title('blur, no noise');
I_deconvdirect0 = ifft2(fft2(I_blurnoise0)./fftH);
figure, imshow(I_deconvdirect0); title('direct deconvolution 0');
MSEa(1)=mean((I_deconvdirect0(:)-I(:)).^2);
disp(['MSE: ', num2str(MSEa(1))])

I_blurnoise1= imnoise(I_blur,'gaussian',0,0.000001^2);
figure, imshow(I_blurnoise1); title('blur + noise 1');
I_deconvdirect1 = ifft2(fft2(I_blurnoise1)./fftH);
figure, imshow(I_deconvdirect1); title('direct deconvolution 1');
MSEa(2)=mean((I_deconvdirect1(:)-I(:)).^2);
disp(['MSE: ', num2str(MSEa(2))])

I_blurnoise2= imnoise(I_blur,'gaussian',0,0.0001^2);
figure, imshow(I_blurnoise2); title('blur + noise 2');
I_deconvdirect2 = ifft2(fft2(I_blurnoise2)./fftH);
figure, imshow(I_deconvdirect2); title('direct deconvolution 2');
MSEa(3)=mean((I_deconvdirect2(:)-I(:)).^2);
disp(['MSE: ', num2str(MSEa(3))])

I_blurnoise3= imnoise(I_blur,'gaussian',0,0.001^2);
figure, imshow(I_blurnoise3); title('blur + noise 3');
I_deconvdirect3 = ifft2(fft2(I_blurnoise3)./fftH);
figure, imshow(I_deconvdirect3); title('direct deconvolution 3');
MSEa(4)=mean((I_deconvdirect3(:)-I(:)).^2);
disp(['MSE: ', num2str(MSEa(4))])

I_blurnoise4= imnoise(I_blur,'gaussian',0,0.01^2);
figure, imshow(I_blurnoise4); title('blur + noise 4');
I_deconvdirect4 = ifft2(fft2(I_blurnoise4)./fftH);
figure, imshow(I_deconvdirect4); title('direct deconvolution 4');
MSEa(5)=mean((I_deconvdirect4(:)-I(:)).^2);
disp(['MSE: ', num2str(MSEa(5))])

I_blurnoise5= imnoise(I_blur,'gaussian',0,0.1^2);
figure, imshow(I_blurnoise5); title('blur + noise 5');
I_deconvdirect5 = ifft2(fft2(I_blurnoise5)./fftH);
figure, imshow(I_deconvdirect5); title('direct deconvolution 5');
MSEa(6)=mean((I_deconvdirect5(:)-I(:)).^2);
disp(['MSE: ', num2str(MSEa(6))])

figure, semilogx(sigma,MSEa,'-x'); xlabel('sigma'), ylabel('MSE')
axis([0,1e-1,0,.1])
legend('direct deconv');


% deconvolution Wiener

powspec = abs(fft2(I)).^2;
% essayer aussi (on ne connait pas I):
% powspec = abs(fft2(I_blurnoise0)).^2;
% selon la thÃ©orie, la MSE doit Ãªtre moins bonne dans ce cas

I_deconvwiener0 = deconvwnr(I_blurnoise0,H,numel(I)*0^2./powspec);
figure, imshow(I_deconvwiener0); title('Wiener deconvolution 0');
MSEb(1)=mean((I_deconvwiener0(:)-I(:)).^2);
disp(['MSE wnr 0: ', num2str(MSEb(1))])

I_deconvwiener1 = deconvwnr(I_blurnoise1,H,numel(I)*0.000001^2./powspec);
figure, imshow(I_deconvwiener1); title('Wiener deconvolution 1');
MSEb(2)=mean((I_deconvwiener1(:)-I(:)).^2);
disp(['MSE wnr 1: ', num2str(MSEb(1))])

I_deconvwiener2 = deconvwnr(I_blurnoise2,H,numel(I)*0.0001^2./powspec);
figure, imshow(I_deconvwiener2); title('Wiener deconvolution 2');
MSEb(3)=mean((I_deconvwiener2(:)-I(:)).^2);
disp(['MSE wnr 2: ', num2str(MSEb(3))])

I_deconvwiener3 = deconvwnr(I_blurnoise3,H,numel(I)*0.001^2./powspec);
figure, imshow(I_deconvwiener3); title('Wiener deconvolution 3');
MSEb(4)=mean((I_deconvwiener3(:)-I(:)).^2);
disp(['MSE wnr 3: ', num2str(MSEb(4))])

I_deconvwiener4 = deconvwnr(I_blurnoise4,H,numel(I)*0.01^2./powspec);
figure, imshow(I_deconvwiener4); title('Wiener deconvolution 4');
MSEb(5)=mean((I_deconvwiener4(:)-I(:)).^2);
disp(['MSE wnr 4: ', num2str(MSEb(5))])

I_deconvwiener5 = deconvwnr(I_blurnoise5,H,numel(I)*0.1^2./powspec);
figure, imshow(I_deconvwiener5); title('Wiener deconvolution 5');
MSEb(6)=mean((I_deconvwiener5(:)-I(:)).^2);
disp(['MSE wnr 5: ', num2str(MSEb(6))])

figure, semilogx(sigma,MSEa,'-x',sigma,MSEb,'-x'); xlabel('sigma'), ylabel('MSE')
axis([0,1e-1,0,.1])
legend('direct deconv','wiener deconv');


% regularisation

I_deconvreg0 = deconvreg(I_blurnoise0,H,numel(I)*0^2);
figure, imshow(I_deconvreg0); title('Regularization deconvolution 0');
MSEc(1)=mean((I_deconvreg0(:)-I(:)).^2);
disp(['MSE reg 0: ', num2str(MSEc(1))])

I_deconvreg1 = deconvreg(I_blurnoise1,H,numel(I)*0.000001^2);
figure, imshow(I_deconvreg1); title('Regularization deconvolution 1');
MSEc(2)=mean((I_deconvreg1(:)-I(:)).^2);
disp(['MSE reg 1: ', num2str(MSEc(3))])

I_deconvreg2 = deconvreg(I_blurnoise2,H,numel(I)*0.0001^2);
figure, imshow(I_deconvreg2); title('Regularization deconvolution 2');
MSEc(3)=mean((I_deconvreg2(:)-I(:)).^2);
disp(['MSE reg 2: ', num2str(MSEc(3))])

I_deconvreg3 = deconvreg(I_blurnoise3,H,numel(I)*0.001^2);
figure, imshow(I_deconvreg3); title('Regularization deconvolution 1');
MSEc(4)=mean((I_deconvreg3(:)-I(:)).^2);
disp(['MSE reg 3: ', num2str(MSEc(4))])

I_deconvreg4 = deconvreg(I_blurnoise4,H,numel(I)*0.01^2);
figure, imshow(I_deconvreg4); title('Regularization deconvolution 4');
MSEc(5)=mean((I_deconvreg4(:)-I(:)).^2);
disp(['MSE reg 4: ', num2str(MSEc(5))])

I_deconvreg5 = deconvreg(I_blurnoise5,H,numel(I)*0.1^2);
figure, imshow(I_deconvreg5); title('Regularization deconvolution 5');
MSEc(6)=mean((I_deconvreg5(:)-I(:)).^2);
disp(['MSE reg 5: ', num2str(MSEc(6))])

figure, semilogx(sigma,MSEa,'-x',sigma,MSEb,'-x',sigma,MSEc,'-x'); xlabel('sigma'), ylabel('MSE')
axis([0,1e-1,0,.05])
legend('direct deconv','wiener deconv','regul deconv');


% TV 

p=2000;  % essayez différentes valeurs de p

out = deconvtv(I_blurnoise0,H,p);
I_deconvtv0=out.f;
figure, imshow(I_deconvtv0); title('TV deconvolution 0');
MSEd(1)=mean((I_deconvtv0(:)-I(:)).^2);
disp(['MSE tv 0: ', num2str(MSEd(1))])

out = deconvtv(I_blurnoise1,H,p);
I_deconvtv1=out.f;
figure, imshow(I_deconvtv1); title('TV deconvolution 1');
MSEd(2)=mean((I_deconvtv1(:)-I(:)).^2);
disp(['MSE tv 1: ', num2str(MSEd(2))])

out = deconvtv(I_blurnoise2,H,p);
I_deconvtv2=out.f;
figure, imshow(I_deconvtv2); title('TV deconvolution 2');
MSEd(3)=mean((I_deconvtv2(:)-I(:)).^2);
disp(['MSE tv 2: ', num2str(MSEd(3))])

out = deconvtv(I_blurnoise3,H,p);
I_deconvtv3=out.f;
figure, imshow(I_deconvtv3); title('TV deconvolution 3');
MSEd(4)=mean((I_deconvtv3(:)-I(:)).^2);
disp(['MSE tv 3: ', num2str(MSEd(4))])

out = deconvtv(I_blurnoise4,H,p);
I_deconvtv4=out.f;
figure, imshow(I_deconvtv4); title('TV deconvolution 4');
MSEd(5)=mean((I_deconvtv4(:)-I(:)).^2);
disp(['MSE tv 4: ', num2str(MSEd(5))])

out = deconvtv(I_blurnoise5,H,p/15);
I_deconvtv5=out.f;
figure, imshow(I_deconvtv5); title('TV deconvolution 5');
MSEd(6)=mean((I_deconvtv5(:)-I(:)).^2);
disp(['MSE tv 5: ', num2str(MSEd(6))])

figure, semilogx(sigma,MSEa,'-x',sigma,MSEb,'-x',sigma,MSEc,'-x',sigma,MSEd,'-x'); 
xlabel('sigma'), ylabel('MSE')
axis([0,1e-1,0,.05])
legend('direct deconv','wiener deconv','regul deconv','TV deconv');

