%%%%%%%%%%%%%%%%% Echo image
im=imread('lena.tif');
im=double(im);
m=zeros(size(im));
m(1,1)=0.5;
m(8,13)=0.5;
M=fft2(m);

bim=abs((ifft2(fft2(im).*(M))));
showImage(bim);
%imwrite(uint8(badim1),'BadImages/echo.jpg','jpg');

bim=readImage('BadImages/echo.tif');
M=fft2(m);
h=find(abs(M)<0.01);
M(h)=1;
rim=abs((ifft2(fft2(bim)./(M))));
showImage(rim);

% pim=rim;
% pim=scaleIm(pim,0,255);
%imwrite(uint8(pim),'BadImages/junk.jpg','jpg');

%%%%%%% prior  term1
noise=randn(size(bim));
nbim=bim+5.*noise;
showImage(nbim);
M=fft2(m);
rim=abs((ifft2(     fft2(nbim).*conj(M)./(conj(M).*M+0.0001))));
showImage(rim);

% pim=rim;
% pim=scaleIm(pim,0,255);
% imwrite(uint8(pim),'BadImages/junk.jpg','jpg');

%%%%%%% prior  term2
% Use same image as above
% noise=randn(size(bim));
% nbim=bim+5.*noise;
% imagesc(nbim);
M=fft2(m);
[U,V]=meshgrid(-128:127,-128:127);
U = fftshift(U);
V = fftshift(V);
R=U.^2+V.^2;
rim=abs((ifft2(     fft2(nbim).*conj(M)./(conj(M).*M+0.0001*sqrt(R)))));
showImage(rim);





%%%%%%%%%%%% moving image
k=zeros(size(im));
k(1,1:12)=1/12;
K=fft2(k);
bim=abs((ifft2(fft2(im).*(K))));
showImage(bim);
%writeImage(round(bim),'BadImages/moving.tif');

bim=readimage('BadImages/moving.tif');
K=fft2(k);
h=find(abs(K)<0.01);
K(h)=1;
rim=abs((ifft2(fft2(bim)./(K))));
showImage(rim);


%%%%%%%%%%%%%%%%%%%%% cyclic noise 
[x,y]=meshgrid(1:256);
s=sin(2*pi*20.*x/256+2*pi*5.*y/256);
imagesc(s);
badim3=im./2+50+60.*s;
imagesc(badim3);
badim3=scaleIm(badim3,0,255);
imwrite(uint8(badim3),'Images/blena.jpg','jpg');

%writeImage(round(badim3),'BadImages/noise1.tif');

bim=imread('BadImages/noise1.tif');
B=fft2(bim);
imagesc(log(abs(fftshift(B)))+1);

B(6,21)=0;
B(252,237)=0;
imagesc(abs(ifft2(B)));
imwrite(uint8(badim3),'Images/blena.jpg','jpg');

% different technique for denoising
h=find(abs(B)>1.e6);
dc=B(1,1);
B(h)=0;
B(1,1)=dc;
imagesc(abs(ifft2(B)));

%%%%%%%%%%%%%%%%%% more cyclic noise
[x,y]=meshgrid(1:256);
s1=0.5*sin(2*pi*12.*x/256+2*pi*7.*y/256)+0.5;
s2=0.5*sin(2*pi*8.*x/256+2*pi*15.*y/256)+0.5;
badim4=im./2+125.*(s1+s2)/2;
imagesc(badim4);
%writeImage(round(badim4),'BadImages/noise2.tif');

bim=readimage('BadImages/noise2.tif');
B=fft2(bim);
h=find(abs(B)>1.e6);
dc=B(1,1);
B(h)=0;
B(1,1)=dc;
r=abs(ifft2(B));
imagesc(r);

%%%%%%%% multiplicative noise
[x,y]=meshgrid(1:256);
s=exp(0.5*sin(2*pi*12.*x/256+2*pi*7.*y/256)+0.5);
imagesc(s);
badim5=im.*s./2.8;
imagesc(badim5);
%writeImage(round(badim5),'BadImages/mulNoise.tif');

bim=double(imread('BadImages/mulNoise.tif'));
B=fft2(log(bim+1));
h=find(abs(B)>1.6e4);
dc=B(1,1);
B(h)=0;
B(1,1)=dc;
r=abs(ifft2(B));
res=exp(r);
imagesc(res);


%%%%%%%%%%%%%%%%% blurring
[x,y]=meshgrid(1:256);
k=zeros(256,256);
d=[1 2 3 5 3 2 1]/17;
dd=conv2(d',d);
k(1:7,1:7)=dd;
K=fft2(k);
B=fft2(im).*K;
badim6=abs(ifft2(B));
imagesc(badim6);
%writeImage(round(badim6),'BadImages/blur.tif');


bim=double(imread('BadImages/blur.tif'));
imagesc(log(abs(fftshift(K))+1))

h=find(abs(K)<0.03);
K(h)=1;
imagesc(log(abs(fftshift(K))+1))

B=fft2(bim);
R=B./K;
r=abs(ifft2(R));
imagesc(r);

