% ex01            Example 1: compression of example image, different methods

%----------------------------------------------------------------------
% Copyright (c) 2010.  Karl Skretting.  All rights reserved.
% Stavanger University, Signal Processing Group
% Mail:  karl.skretting@uis.no   Homepage:  http://www.ux.uis.no/~karlsk/
% 
% HISTORY:
% Ver. 1.0  15.03.2010  KS: m-file made
% Ver. 1.1  01.02.2011  KS: added JP2 option
%----------------------------------------------------------------------

if ~exist('imageFile', 'var')
    imageFile = 'lena.bmp';
end
if ~exist('makeFigure', 'var')
    makeFigure = true;
end
if ~exist('makeTable', 'var')
    makeTable = true;
end
if ~exist('transforms', 'var')
    transforms = {'JPEG', 'DCT', 'LOT', 'M97', 'SPIHT', 'JP2'};
end

if ~exist(imageFile, 'file')
    fprintf('Image %s does not exist.',imageFile);
    return
end

A = double(imread(imageFile))-128;
[N,M] = size(A);
peak = 255;
res = cell(numel(transforms),1);

for tno = 1:numel(transforms)
    
    if strcmpi(transforms{tno}, 'JP2')
        if verLessThan('matlab', '7.11.0')
            disp(' ');
            disp('Matlab imwrite with JP2 option not available in this version.');
            if strcmpi(imageFile, 'lena.bmp')
                disp('Results for lena.bmp from another computer is used.');
                res{tno} = struct('bits',[26208, 33464, 50968, 64032, 80240, 100064, ...
                    126728, 157728, 192664, 217072, 257504, 313824, 388160, 452240, ...
                    519904, 595472, 695344], ...
                    'PSNR',[31.0695, 32.1500, 34.0647, 35.1917, 36.2943, 37.2754, ...
                    38.4545, 39.5063, 40.3786, 41.0092, 41.8209, 42.8968, 44.0517, ...
                    45.0905, 46.3231, 47.7135, 49.3444]  );
            elseif strcmpi(imageFile, 'barbara.bmp')
                disp('Results for barbara.bmp from another computer is used.');
                res{tno} = struct('bits',[26880, 32936, 50968, 64040, 79896, ...
                    101336, 127568, 157512, 192248, 215728, 258384, 314192, ...
                    389488, 452112, 520808, 592416, 692440], ...
                    'PSNR',[24.7801, 25.4151, 27.1745, 28.2991, 29.3373, 30.6678, ...
                    32.0822, 33.3990, 34.7408, 35.5528, 37.0757, 38.5875, 40.4192, ...
                    41.7770, 43.0677, 44.2891, 46.2447] );
            elseif strcmpi(imageFile, 'boat.bmp')
                disp('Results for boat.bmp from another computer is used.');
                res{tno} = struct('bits',[26872, 33168, 50392, 64008, 80288, ...
                    101384, 127536, 156256, 191208, 216824, 257328, 313552, ...
                    387936, 451736, 520304, 592488, 694032], ...
                    'PSNR',[26.6568, 27.4076, 29.0232, 30.0227, 30.9728, ...
                    32.1109, 33.2005, 34.1659, 35.1261, 35.7873, 36.6386, ...
                    37.7957, 39.0926, 40.4254, 41.9366, 43.3093, 45.1097]  );
            else
                res{tno} = struct('bits',8*numel(A), 'PSNR', inf);
            end
        else
            disp(' ');
            disp('Use Matlab imwrite with JP2 option.');
            CR = [80 60 40 32 25 20 16 13 10.67 9.5 8 6.6 5.33 4.6 4 3.5 3];
            res{tno} = struct('bits',zeros(1,numel(CR)), 'PSNR',zeros(1,numel(CR)) );
            i = 0;
            for cr = CR
                i = i+1;
                imwrite(uint8(A+128),'tempfile.jp2','jp2','CompressionRatio',cr);
                h = dir('tempfile.jp2');
                res{tno}.bits(i) = (h.bytes)*8;  % bits
                Ar = double(imread('tempfile.jp2'));
                R = double(A+128)-Ar;
                res{tno}.PSNR(i) = 10*log10(((numel(A)*peak^2)/sum(sum(R.*R))));
            end
        end
        
    elseif strcmpi(transforms{tno}, 'JPEG')
        disp(' ');
        disp('Use Matlab imwrite with JPEG option.');
        quality = [10:5:90,93];
        res{tno} = struct('bits',zeros(1,numel(quality)), 'PSNR',zeros(1,numel(quality)) );
        i = 0;
        for q = quality
            i = i+1;
            imwrite(uint8(A+128),'tempfile.jpg','JPEG','Quality',q);
            h = dir('tempfile.jpg');
            res{tno}.bits(i) = (h.bytes)*8;  % bits
            Ar = double(imread('tempfile.jpg'));
            R = double(A+128)-Ar;
            res{tno}.PSNR(i) = 10*log10(((numel(A)*peak^2)/sum(sum(R.*R))));
        end
    
    elseif strcmpi(transforms{tno}, 'SPIHT')
       if strcmpi(imageFile, 'lena.bmp')
             disp('Results for lena.bmp and SPIHT found earlier.');
            % The SPIHT results were found using the Matlab files from
            % http://www.cipr.rpi.edu/research/SPIHT/spiht3.html
            res{tno} = struct('bits',[26222, 39329, 52436, 65544, 78651, 91758, 104865, 117972, 131080, ...
                144187, 157294, 170401, 183508, 196616, 209723, 222830, 235937, 249044, 262152, ...
                288366, 314580, 340795, 367009, 393224, 419438, 445652, 471867, 498081, 524296], ...
                'PSNR',[30.9116, 32.5325, 33.9122, 34.7792, 35.6850, 36.5298, 37.1707, 37.6565, 38.1067, 38.5569, 38.9265, ...
                39.2193, 39.7848, 40.1233, 40.4369, 40.7175, 40.9311, 41.1792, 41.4044, 41.8447, 42.2400, 42.5826, ...
                42.8745, 43.2606, 43.9558, 44.3019, 44.6520, 44.9937, 45.3059]  );
       else
           res{tno} = compressSPIHT(A, imageFile, 0.2:0.1:2.1);
           if ~isfield(res{tno}, 'bits'); res{tno}.bits = 8*size(A); end;
           if ~isfield(res{tno}, 'PSNR'); res{tno}.PSNR = inf; end;
       end
       
    else    % results for more transforms
        del = fliplr([2.5,3,4,5,6,8,10,13,16,20,25,30,35,40,50,60]);
        X = myim2col(A, 't',transforms{tno}, 's',8);
        res{tno} = struct('bits',zeros(1,numel(del)), ...
            'PSNR',zeros(1,numel(del)), ...
            'del', del);
        for i=1:numel(del)
            disp(' ');
            disp(['Use ',transforms{tno},' and delta = ',num2str(del(i))]);
            Z = uniquant(X, del(i), del(i), 2000);
            Xr = uniquant(Z, del(i), del(i));  % inverse quantizing
            Ar = mycol2im(Xr, 't',transforms{tno}, 's',8, 'imsize',[M,N]);
            R = A-Ar;
            res{tno}.PSNR(i) = 10*log10((((M*N)*peak^2)/sum(sum(R.*R))));
            xC = myreshape(Z, 'm',1, 'sDC',[M/8,N/8]);
            res{tno}.bits(i) = estimateBits(xC, 'Huff06', 2);
        end
    end
end

if makeFigure % make the figure
    col = 'bgrcmkmbgrckbgrcmk';
    figure(1); clf; hold on; grid on;
    if strcmpi(imageFile, 'lena.bmp')
        axis([0,2,30,46]); x2=0.9; y2=31; 
    elseif strcmpi(imageFile, 'barbara.bmp')
        axis([0,2,26,42]); x2=1.05; y2=27; 
    elseif strcmpi(imageFile, 'boat.bmp')
        axis([0,2,28,42]); x2=1.05; y2=28.5; 
    end
    x1=0.6; dy=1.1; dx=0.03;
    for tno = 1:numel(res);   % or [1,3,4,2,5]
        plot(res{tno}.bits/(M*N), res{tno}.PSNR, [col(tno),'-']);
        y1 = interp1(res{tno}.bits/(M*N), res{tno}.PSNR, x1, 'linear');
        plot([x1,x2], [y1, y2], [col(tno),'.-']);
        h = text(x2+dx, y2, [transforms{tno},', ',num2str(y1,'%5.2f')]);
        set(h, 'BackgroundColor', [1,1,1]);
        y2 = y2+dy;
    end
    title(['Image compression of ',imageFile]);
    xlabel(['Bit rate.    (',datestr(now),')']);
    ylabel('PSNR');
end

if makeTable
    xtab = [0.25, 0.5, 0.75, 1.0, 1.25, 1.50, 1.75, 2.0];
    linje = ' %10s & %5.2f & %5.2f & %5.2f & %5.2f & %5.2f & %5.2f & %5.2f & %5.2f \n';
    fprintf(linje, 'bit rate', xtab); 
    for tno = 1:numel(res)
        ytab = zeros( size(xtab) );
        for i = 1:numel(xtab) 
            if (res{tno}.bits(1)/(M*N) <= xtab(i)) && (res{tno}.bits(end)/(M*N) >= xtab(i))    
                ytab(i) = interp1(res{tno}.bits/(M*N), res{tno}.PSNR, xtab(i), 'linear');
            end
        end
        fprintf(linje, transforms{tno}, ytab);
    end 
end


