test_dwt.m

SPIHT coding implementation. using wavelet transforms as a key tool. compression and decoding is don

function test_DWT(soubor, wavelet, level, bpp, zobr)
% DWT testing function
% parameters:     soubor - input image file
%                 wavelet - wavelet name or 'cdf97'
%                 level - transform depth
%                 bpp - bits per pixel quantifier
%                 zobr - result 
%                            0 - only input + output images
%                            1 - more details

% residuals coefficient
mult = 10;

% init
image = floor(double(imread(soubor)));
size_x = size(image,1);

% image -> DWT domain
if strcmp(wavelet, 'cdf97')
  tile = waveletcdf97(image, level);
else
  mode = 'per';
  dwtmode(mode);
  type = wavelet;
  [Lo_D,Hi_D,Lo_R,Hi_R] = wfilters(type);
  [tile] = makeDwtTile(image, level, Lo_D, Hi_D);
end

% SPIHT encoder
% bpp -> bits count
bpp_full = 8;
bytes = ceil((bpp/bpp_full)*(size_x^2));

eff = bytes;
bits = bytes*8;

% SPIHT coder 
[max_pass, passess, bitstream, timeel] = encodeSPIHT(tile, bits, level);
disp(['ENCODER: Encoding stopped at ' num2str(max_pass - passess) ' after ' num2str(toc) 's, ' num2str(ceil(size(bitstream,2)/8)+1) 'B transmitted']);
timeel = timeel + toc;
disp(['ENCODER: Total time elapsed: ' num2str(timeel) 's.']);  

% SPIHT decoder
[max_pass, passess, tile2, timeel2] = decodeSPIHT(bitstream);
disp(['DECODER: Decoding stopped at ' num2str(max_pass - passess) ' after ' num2str(toc) 's']);
timeel2 = timeel2 + toc;
disp(['DECODER: Total time elapsed: ' num2str(timeel2) 's.']);

% DWT domain -> image
if strcmp(wavelet, 'cdf97')
  recon = waveletcdf97(tile2, -level);
else
  recon = reconstructDwtTile(tile2, Lo_R, Hi_R, level);
end

% residuals computation
residuals = abs(recon-image);
mult = 10;
residuals = residuals .* mult;
% PSNR + rounding
sum_1 = 0;

for i=1:size(image,1)
    for j=1:size(image,1)
        sum_1 = sum_1 + (image(i,j) - recon(i,j))^2;
    end
end

MSE = sum_1/(size_x^2);
RMSE = sqrt(MSE);
PSNR = 20*log10(255/RMSE)*100;
PSNR = round(PSNR);
PSNR = PSNR/100

% show results
str_orig = ['Original, dimensions=' num2str(size_x) 'x' num2str(size_x) 'px, size=' num2str(size_x^2) 'B'];
str_tile = ['DWT decomposition, wavelet=' wavelet ', level=' num2str(level)];
str_dec = ['DWT-SPIHT result, size=' num2str(ceil(size(bitstream,2)/8)+1) 'B (1:' num2str(round((size_x^2)/(size(bitstream,2)/8))) ' / ' num2str(bpp) 'bpp), PSNR=' num2str(PSNR) 'dB'];
str_res = ['Absolute difference (values magnified by ' num2str(mult) 'x)'];

if zobr == 1
    figure(5);
    
    subplot(1,2,1);
    imagesc(tile);
    colormap('gray');
    title(str_tile);
    
    subplot(1,2,2);
    imagesc(residuals, [0 255]);
    colormap('gray');
    title(str_res);
    
    figure(6);
    mesh(tile);
    title('3D graph of the decomposition image, level=3');
    
end
figure(1);
subplot(1,2,1);
imagesc(image, [0 255]);
colormap('gray');
title(str_orig);
    
subplot(1,2,2);
imagesc(recon, [0 255]);
colormap('gray');
title(str_dec);