%ee368b project 11/14/2000
%Burt/Adelson Laplacian pyramid coder
%use function reduce.m and expand.m

clear; 
close all;

g=imread('lena.tiff','tif');
%g=imread('zone.tiff','tif');
imagesize=size(g);
stepsize=[32 16 8 4]; %quantizer stepsize
a=.4; %parameter for filter
what=[1/4-a/2 1/4 a 1/4 1/4-a/2];   %filter
w=what'*what;

g=[g(1,1) g(1,1:imagesize(2));
   g(1:imagesize(1),1) g]; %resize the image from 512x512 to 513x513 for appropriate filtering
g=double(g);
figure;
imshow(uint8(g));
title('g');

%generate the next 3 levels of Gaussian pyramid
g1=reduce(g,w);
figure;
imshow(uint8(g1));
title('g1');

g2=reduce(g1,w);
figure;
imshow(uint8(g2));
title('g2');

g3=reduce(g2,w);
figure;
imshow(uint8(g3));
title('g3');

%interpolate images in the Gaussian pyramid
g11=expand(g1,w);
figure;
imshow(uint8(g11));
title('g11');

g21=expand(g2,w);
g22=expand(g21,w);
figure;
imshow(uint8(g22));
title('g22');

g31=expand(g3,w);
g32=expand(g31,w);
g33=expand(g32,w);
figure;
imshow(uint8(g33));
title('g33');

%generate difference images in the Laplacian pyramid
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');

%pass the difference images to the quantizers
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 the distribution after quantization
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(:)));

%histogram of each difference image
figure;
subplot(2,2,1); hist(C3(:),min(C3(:)):stepsize(4):max(C3(:)));
subplot(2,2,2); hist(C2(:),min(C2(:)):stepsize(3):max(C2(:)));
subplot(2,2,3); hist(C1(:),min(C1(:)):stepsize(2):max(C1(:)));
subplot(2,2,4); hist(C0(:),min(C0(:)):stepsize(1):max(C0(:)));

%normalize pixel frequency to probability
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)));

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

%reconstruct image
r2=C2+expand(C3,w)-128;
r1=C1+expand(r2,w)-128;
r0=C0+expand(r1,w)-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(g),'ba_original_lena.tif','tif','compression','none');
%imwrite(uint8(g1),'ba_g1_lena.tif','tif','compression','none');
%imwrite(uint8(g2),'ba_g2_lena.tif','tif','compression','none');
%imwrite(uint8(g3),'ba_g3_lena.tif','tif','compression','none');
%imwrite(uint8(L0),'ba_l0_lena.tif','tif','compression','none');
%imwrite(uint8(L1),'ba_l1_lena.tif','tif','compression','none');
%imwrite(uint8(L2),'ba_l2_lena.tif','tif','compression','none');
%imwrite(uint8(C3),'ba_r3_lena.tif','tif','compression','none');
%imwrite(uint8(r2),'ba_r2_lena.tif','tif','compression','none');
%imwrite(uint8(r1),'ba_r1_lena.tif','tif','compression','none');
%imwrite(uint8(r0),'ba_r0_lena.tif','tif','compression','none');
