% Usage: [x2,P] = jpegprogressive(x,r,p)
% Simple JPEG compression/progressive transmission demo.  Input 
% argument x is an M by N array of numbers (can be uint8 or double,
% it's converted to a double below) assumed to encode grayscale image 
% on 0 to 255 range.  Argument "r" is for scaling the quantization 
% matrix as in book. Uses quantization matrix from book.
% Argument "p" is a sort of frequency limit that uses ONLY DCT 2D
% frequencies which satisfy k + l <= p to reconstruct. Choose 
% 0<=p<=14. With p=14 this is equivalent to "jpegdemo". With p=0
% we have reconstruction from DC coefficients only.

% Output is M by N array "x2", simulated compressed/decompressed 
% image as an array of doubles, and also "P", fraction of non-zero 
% block DCT coefficients (as a crude measure of compression).

function [x2,P] = jpegprogressive(x,r,p)

T = r*[16 11 10 16 24 40 51 61; 12 12 14 19 26 58 60 55; 
    14 13 16 24 40 57 69 56; 14 17 22 29 51 87 80 62;
    18 22 37 56 68 109 103 77; 24 35 55 64 81 194 113 92;
    49 64 78 87 103 121 120 101; 72 92 95 98 121 100 103 99];

 m=8; n=8; %Really, the blocksize m, n can be anything.
 
 %Extend image x so dimensions divisible by m, n.
 [M0,N0] = size(x); M=M0; N=N0;
 blocksM = ceil(M/m); blocksN = ceil(N/n);
 M = m*blocksM; N = n*blocksN;
 for k=(M0+1):M %Extend each column
     x(k,:) = x(M0,:);
 end
 for k=(N0+1):N %Extend each row
     x(:,k) = x(:,N0);
 end
 
 Y = zeros(M,N); %Storage for block DCT
 
 %Arrays for 8x8 DCT; straight matrix multiply is faster than dct.
 C1 = dct(eye(m));
 C2 = dct(eye(n))';
 
 %Construct mask B for zeroing out all but desired DCT coefficients, to
 %simulate progressive transmission
 B = ones(m,n);
 for i=1:m
     q=max(p+3-i,1);
     B(i,q:n) = 0.0;
 end;
 
 %Block DCT code, quantization
 for i=1:blocksM
     for j=1:blocksN
         
         %Pick off each 8 x 8 block, put in "Z"
         e = (i-1)*m+1;
         f = (j-1)*n+1;
         Z = double(x(e:(e+m-1),(f:f+n-1)));
         
         %Compute DCT of block
         Q = C1*Z*C2;
         
         %Quantize (quant. matrix already scaled), install in Y.
         Q2 = T.*round(Q./T);
         Q2 = Q2.*B; %Use only coefficients that satisfy resolution limits.
         Y(e:(e+m-1),(f:f+n-1)) = Q2;
         
         end;
     end;
     
 P = sum(sum(sign(abs(Y))))/M/N;

 %Inverse block DCT
 Y = double(blkidct2(Y,m,n));
 
 %Pick off original image dimensions
 x2 = double(Y(1:M0,1:N0));