% Course in Wavelets, TU Chemnitz, 2012

C=imread('Lenna.png');

A=double(C(:,:,1));
[m,n]=size(A);
s=0.05;
r=4;
sigma=s*(max(A(:))-min(A(:)));
An=A+sigma*randn(m,n);

figure(1);
image(A);
colormap(gray(256));
title('Original');

figure(2);
image(An);
colormap(gray(256));
title('Noisy');

% Single level DWT
figure(3)
dwtmode('per');
[C,H,V,D] = dwt2(A,'db4');
imagesc([C,H;V,D]);
colormap(gray(256));
title('Single level DWT');

% Combined plots
figure(4);
l=4;
[c,S]=wavedec2(An,l,'db4');

C=reshape(c(1:S(1,1)*S(1,2)),S(1,1),S(1,2));
offset=1+S(1,1)*S(1,2);

for j=2:l+1
  s1=S(j,1);
  s2=S(j,2);
  C=[C,...
     reshape(c(offset:offset+s1*s2-1),s1,s2);...
     reshape(c(offset+s1*s2:offset+2*s1*s2-1),s1,s2),...
     reshape(c(offset+2*s1*s2:offset+3*s1*s2-1),s1,s2)];
  offset=offset+3*s1*s2;
end;

imagesc(C);
colormap(gray(256));
title('DWT');

% Filtering
figure(5)
r=1;
cc=zeros(m*n,1);
cc(1:m*n/4^r)=c(1:m*n/4^r);

Af=waverec2(cc,S,'db4');

image(Af);
colormap(gray(256));
title('Filtered');

% Thresholding
figure(6)
s=0.01;
c(abs(c)<s*max(abs(c)))=0;

At=waverec2(c,S,'db4');

image(At);
colormap(gray(256));
title('Thresholded');
