%ee368b project fall2000
%Laplacian pyramid with many filters
%use switch_reduce.m and expand2.m
%switch_reduce.m uses filter4x4.m which uses comp_diff.m

clear; 
close all;

%g=imread('lena.tiff','tif');
g=imread('zone.tif','tif');
imagesize=size(g);
stepsize=[32 16 8 4];

g=double(g);
figure;
imshow(uint8(g));
title('g');

%filter and subsample image, and find which filter is used in different area in the image
[g1,met1]=switch_reduce(g);
figure;
imshow(uint8(g1));
title('g1');
figure;
imshow(uint8(met1));
title('method');

[g2,met2]=switch_reduce(g1);
figure;
imshow(uint8(g2));
title('g2');
figure;
imshow(uint8(met2));
title('method');

[g3,met3]=switch_reduce(g2);
figure;
imshow(uint8(g3));
title('g3');
figure;
imshow(uint8(met3));
title('method');

%interpolate using the same interpolator as the nonlinear_pyramid
g11=expand2(g1);
figure;
imshow(uint8(g11));
title('g11');

g21=expand2(g2);

g31=expand2(g3);

%compute difference images
L0=g-g11+128;
figure;
imshow(uint8(L0));
title('L0');

L1=g1-g21+128;
figure;
imshow(uint8(L1));
title('L1');

L2=g2-g31+128;
figure;
imshow(uint8(L2));
title('L2');

%quantize difference images
C3=stepsize(4)*round(g3/stepsize(4));
C2=stepsize(3)*round(L2/stepsize(3));
C1=stepsize(2)*round(L1/stepsize(2));
C0=stepsize(1)*round(L0/stepsize(1));

%find frequency of pixel values
p3=hist(C3(:),min(C3(:)):stepsize(4):max(C3(:)));
p2=hist(C2(:),min(C2(:)):stepsize(3):max(C2(:)));
p1=hist(C1(:),min(C1(:)):stepsize(2):max(C1(:)));
p0=hist(C0(:),min(C0(:)):stepsize(1):max(C0(:)));

%find probability of pixel values
p3=p3/length(C3(:));
p2=p2/length(C2(:));
p1=p1/length(C1(:));
p0=p0/length(C0(:));

%compute entropy
H3=-sum(nonzeros(p3).*log2(nonzeros(p3)));
H2=-sum(nonzeros(p2).*log2(nonzeros(p2)));
H1=-sum(nonzeros(p1).*log2(nonzeros(p1)));
H0=-sum(nonzeros(p0).*log2(nonzeros(p0)));

%compute bit-rate
bitrate=(H3*length(C3(:))+H2*length(C2(:))+H1*length(C1(:))+H0*length(C0(:)))/length(C0(:));

%reconstruct image
r2=C2+expand2(C3)-128;
r1=C1+expand2(r2)-128;
r0=C0+expand2(r1)-128;
figure;
imshow(uint8(r0));
title('reconstructed image');

%compute distortion and PSNR
D=sum(sum((g-r0).^2))/length(r0(:));
PSNR=10*log10(255^2/D);

imwrite(uint8(g1),'switch_g1_lena.tif','tif','compression','none');
imwrite(uint8(met1),'switch_met1_lena.tif','tif','compression','none');
imwrite(uint8(r0),'switch_r0_lena.tif','tif','compression','none');
