jpeg.m

simple jpeg compression and decompression

function jpeg
% This is a JPEG encoding & decoding program of still image.
% it does not use level shifting.
% Discrete Cosine transform (DCT) is performed both by classical & Chen's
% Flowgraph methods. Predefined JPEG quantization array & Zigzag order are
% used here. 'RUN', 'LEVEL' coding is used instead of  Huffman coding.
% Compression ratio is compared for each DCT method. Effect of coarse and fine quantization is
% also examined. The execution time of 2 DCT methods is also checked.
% In addition, most energatic DCT coefficients are also applied to examine
% the effect in MatLab 7.4.0 R2007a. Input is 9 gray scale pictures &
% output is 9*9=81 pictures to compare. Blocking effect is obvious.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %-------------------------- Initialization -----------------------------
    % JPEG default quantization array
    Q_8x8 =uint8([
            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 104 113  92
            49 64 78 87 103 121 120 101
            72 92 95 98 112 100 103 99]);

    % I am interested in the energy of the dct coefficients, I will use 
    % this matrix (found from long observation) to select dct coefficients.
    % lowest number -> highest priority
    dct_coefficient_priority_8x8 =[
           1   2   6   7  15  16  28  29;
           3   5   8  14  17  27  30  43;
           4   9  13  18  26  31  42  44;
          10  12  19  25  32  41  45  54;
          11  20  24  33  40  46  53  55;
          21  23  34  39  47  52  56  61;
          22  35  38  48  51  57  60  62;
          36  37  49  50  58  59  63  64];
    % if we decide to take 10 coefficients with the most energy, we will assign 
    % 99 to ignore the other coefficients and remain with a matrix of 8x8

    %This suitable Zigzag order is formed from the  JPEG standard
    ZigZag_Order = uint8([
            1  9  2  3  10 17 25 18
            11 4  5  12 19 26 33 41
            34 27 20 13 6  7  14 21 
            28 35 42 49 57 50 43 36 
            29 22 15 8  16 23 30 37
            44 51 58 59 52 45 38 31 
            24 32 39 46 53 60 61 54 
            47 40 48 55 62 63 56 64]);

    % Finding the reverse zigzag order (8x8 matrix)
    reverse_zigzag_order_8x8 = zeros(8,8);
    for k = 1:(size(ZigZag_Order,1) *size(ZigZag_Order,2)) 
        reverse_zigzag_order_8x8(k) = find(ZigZag_Order== k); 
    end;
                    
    Compressed_image_size=0;
    %---------------------------------------------------------------------

    
close all;

for Image_Index = 8:8    % the whole program will be tested for 9 images (0->8)
   figure;          %keep the input-output for each image seperately

   
    %--------------------------load a picture ----------------------------
    switch Image_Index
    case {0,1}, input_image_128x128 = im2double( imread( sprintf( '%d.tif',Image_Index ),'tiff' ) );
    otherwise, input_image_128x128 = im2double( imread( sprintf( '%d.tif',Image_Index),'jpeg' ) );
    end   
    %-------------------------- ------------------------------------------

    
    %---------------- show the input image -------------------------------
    subplot(3,3,1);
    imshow(input_image_128x128);
    title( sprintf('original image #%d',Image_Index) );
    %---------------------------------------------------------------------

    for Quantization_Quality = 0:1 % 0 -> coarse quantization, 1 -> fine quantization

    for DCT_type = 0:1    % 0 -> classic DCT, 1 -> Flowgraph fast DCT

    for chosen_number_of_dct_coefficient = 1:63:64    % select 1 or 64 dct coefficients

        
    %---------------- choose energetic DCT coefficients ------------------
    % I will use this matrix to choose only the wanted number of dct coefficients
    % the matrix is initialized to zeros -> zero coefficient is chosen at the beginning
    coef_selection_matrix_8x8 = zeros(8,8);

    % this loop will choose 1 dct coefficients each time
    for l=1:chosen_number_of_dct_coefficient 
        % find the most energetic coefficient from the mean_matrix
        [y,x] = find(dct_coefficient_priority_8x8==min(min(dct_coefficient_priority_8x8)));
    
        % select specific coefficients by location index y,x for the image to be compressed
        coef_selection_matrix_8x8(y,x) = 1;
    
        % set it as 99 for the chosen dct coefficient, so that in the next loop, we will choose the "next-most-energetic" coefficient
        dct_coefficient_priority_8x8(y,x) = 99;
    end
    
    % replicate the selection matrix for all the parts of the dct transform
    selection_matrix_128x128 = repmat( coef_selection_matrix_8x8,16,16 );
    %---------------------------------------------------------------------
    
    
    tic ;    % start mark for elapsed time for encoding & decoding
    
    
    %------------------------- Forward DCT -------------------------------        
    % for each picture perform a 2 dimensional dct on 8x8 blocks.    
    if DCT_type==0
        dct_transformed_image = Classic_DCT(input_image_128x128) .* selection_matrix_128x128;
    else
        dct_transformed_image = image_8x8_block_flowgraph_forward_dct(input_image_128x128) .* selection_matrix_128x128;
    end
    %---------------------------------------------------------------------

    
    %---------------- show the DCT of image ------------------------------
% one can use this portion to show DCT coefficients of the image
%    subplot(2,2,2);
%    imshow(dct_transformed_image);
%    title( sprintf('8x8 DCT of image #%d',Image_Index) );
    %---------------------------------------------------------------------
    
    
    %normalize dct_transformed_image by the maximum coefficient value in dct_transformed_image 
    Maximum_Value_of_dct_coeffieient = max(max(dct_transformed_image));
    dct_transformed_image = dct_transformed_image./Maximum_Value_of_dct_coeffieient;
    
    %integer conversion of dct_transformed_image 
    dct_transformed_image_int = im2uint8( dct_transformed_image ); 
    
    
    %-------------------- Quantization -----------------------------------
    % replicate the 'Q_8x8' for at a time whole (128x128) image quantization
    if Quantization_Quality==0
        quantization_matrix_128x128 = repmat(Q_8x8,16,16 ); %for coarse quantization
    else
        quantization_matrix_128x128 = repmat(uint8(ceil(double(Q_8x8)./40)),16,16 );  %for fine quantization
    end
    
    %at a time whole image (128x128) quantization
    quantized_image_128x128 =  round(dct_transformed_image_int ./quantization_matrix_128x128) ; %round operation should be done here for lossy quantization
    %---------------------------------------------------------------------
    

    % Break 8x8 block into columns
    Single_column_quantized_image=im2col(quantized_image_128x128, [8 8],'distinct');

    
    %--------------------------- zigzag ----------------------------------
    % using the MatLab Matrix indexing power (specially the ':' operator) rather than any function
    ZigZaged_Single_Column_Image=Single_column_quantized_image(ZigZag_Order,:);    
    %---------------------------------------------------------------------


    %---------------------- Run Level Coding -----------------------------
    % construct Run Level Pair from ZigZaged_Single_Column_Image
    run_level_pairs=uint8([]);
    for block_index=1:256    %block by block - total 256 blocks (8x8) in the 128x128 image
        single_block_image_vector_64(1:64)=0;
        for Temp_Vector_Index=1:64
            single_block_image_vector_64(Temp_Vector_Index) = ZigZaged_Single_Column_Image(Temp_Vector_Index, block_index);  %select 1 block sequentially from the ZigZaged_Single_Column_Image
        end
        non_zero_value_index_array = find(single_block_image_vector_64~=0); % index array of next non-zero entry in a block
        number_of_non_zero_entries = length(non_zero_value_index_array);  % # of non-zero entries in a block

    % Case 1: if first ac coefficient has no leading zeros then encode first coefficient
        if non_zero_value_index_array(1)==1,  
           run=0;   % no leading zero
            run_level_pairs=cat(1,run_level_pairs, run, single_block_image_vector_64(non_zero_value_index_array(1)));
        end

    % Case 2: loop through each non-zero entry    
        for n=2:number_of_non_zero_entries, 
            % check # of leading zeros (run)
            run=non_zero_value_index_array(n)-non_zero_value_index_array(n-1)-1;
            run_level_pairs=cat(1, run_level_pairs, run, single_block_image_vector_64(non_zero_value_index_array(n)));
        end
        
    % Case 3: "End of Block" mark insertion
        run_level_pairs=cat(1, run_level_pairs, 255, 255);
    end
    %---------------------------------------------------------------------
    
    
    Compressed_image_size=size(run_level_pairs);        % file size after compression
    Compression_Ratio = 20480/Compressed_image_size(1,1);



% % %  -------------------------------------------------------------------
% % %  -------------------------------------------------------------------
% % %                DECODING
% % %  -------------------------------------------------------------------
% % %  -------------------------------------------------------------------

    

    %---------------------- Run Level Decoding ---------------------------
    % construct  ZigZaged_Single_Column_Image from Run Level Pair 
    c=[];
    for n=1:2:size(run_level_pairs), % loop through run_level_pairs
        % Case 1 & Cae 2 
        % concatenate zeros according to 'run' value
        if run_level_pairs(n)<255 % only end of block should have 255 value
            zero_count=0;
            zero_count=run_level_pairs(n);
            for l=1:zero_count    % concatenation of zeros accouring to zero_count
                c=cat(1,c,0);   % single zero concatenation
            end
            c=cat(1,c,run_level_pairs(n+1)); % concatenate single'level' i.e., a non zero value
       
        % Case 3: End of Block decoding
        else
            number_of_trailing_zeros= 64-mod(size(c),64);
            for l= 1:number_of_trailing_zeros    % concatenate as much zeros as needed to fill a block
                c=cat(1,c,0); 
            end
        end
    end
    %---------------------------------------------------------------------
    

    %---------------------------------------------------------------------
    %    prepare the ZigZaged_Single_Column_Image vector (each column represents 1 block) from the
    %    intermediate concatenated vector "c"
    for i=1:256
        for j=1:64
            ZigZaged_Single_Column_Image(j,i)=c(64*(i-1)+j);
        end
    end
    %---------------------------------------------------------------------
     
    
    %--------------------------- reverse zigzag --------------------------
    %reverse zigzag procedure using the matrix indexing capability of MatLab (specially the ':' operator)
    Single_column_quantized_image = ZigZaged_Single_Column_Image(reverse_zigzag_order_8x8,:);
    %---------------------------------------------------------------------
    

   %image matrix construction from image column
    quantized_image_128x128 = col2im(Single_column_quantized_image,   [8 8],   [128 128],   'distinct');

    
    %-------------------- deQuantization ---------------------------------
    dct_transformed_image =  quantized_image_128x128.*quantization_matrix_128x128 ;
    %---------------------------------------------------------------------

    
    %-------------------------- Inverse DCT ------------------------------
    % restore the compressed image from the given set of coeficients
    if DCT_type==0
        restored_image = image_8x8_block_inv_dct( im2double(dct_transformed_image).*Maximum_Value_of_dct_coeffieient  ); %Maximum_Value_of_dct_coeffieient is used for reverse nornalization
    else
        restored_image = image_8x8_block_flowgraph_inverse_dct( im2double(dct_transformed_image).*Maximum_Value_of_dct_coeffieient  ); %Maximum_Value_of_dct_coeffieient is used for reverse nornalization
    end
    %---------------------------------------------------------------------

    
    elapsed_time = toc;    % time required for both enconing & decoding

    
    %-------------------------- Show restored image ----------------------
    subplot(3,3, Quantization_Quality*2^2+ DCT_type*2+ floor(chosen_number_of_dct_coefficient/64)+2);
    imshow( restored_image );
    
    if DCT_type == 0 
        if Quantization_Quality == 0
            title( sprintf('coarse quantize\nClassic DCT\nRestored image with %d coeffs\nCompression ratio %.2f\nTime %f',chosen_number_of_dct_coefficient,Compression_Ratio,elapsed_time) );
        else
            title( sprintf('fine quantize\nclassic DCT\nRestored image with %d coeffs\nCompression ratio %.2f\nTime %f',chosen_number_of_dct_coefficient,Compression_Ratio,elapsed_time) );
        end
    else
        if Quantization_Quality == 0
            title( sprintf('coarse quantize\nFast DCT\nRestored image with %d coeffs\nCompression ratio %.2f\nTime %f',chosen_number_of_dct_coefficient,Compression_Ratio,elapsed_time) );
        else
            title( sprintf('fine quantize\nFast DCT\nRestored image with %d coeffs\nCompression ratio %.2f\nTime %f',chosen_number_of_dct_coefficient,Compression_Ratio,elapsed_time) );
        end        
    end
    %---------------------------------------------------------------------
    
    
    end    % end of coefficient number loop
    end    % end of DCT type loop
    end    % end of quantization qualoty loop
end    % end of image index loop

end    % end of 'jpeg' function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%         I N N E R   F U N C T I O N   I M P L E M E N T A T I O N
%% -----------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ------------------------------------------------------------------------
% Classic_DCT_Block_8x8 - implementation of a 2 Dimensional DCT
% assumption: input matrix is a square matrix !
% ------------------------------------------------------------------------
function out = Classic_DCT_Block_8x8( in )


% get input matrix size
N = size(in,1);

% build the matrix
n = 0:N-1;
for k = 0:N-1
   if (k>0)
      C(k+1,n+1) = cos(pi*(2*n+1)*k/2/N)/sqrt(N)*sqrt(2);
   else
      C(k+1,n+1) = cos(pi*(2*n+1)*k/2/N)/sqrt(N);
   end   
end
out = C*in*(C');
end    % end of Classic_DCT_Block_8x8 function


% ------------------------------------------------------------------------
% pdip_inv_dct2 - implementation of an inverse 2 Dimensional DCT
% assumption: input matrix is a square matrix !
% ------------------------------------------------------------------------
function out = pdip_inv_dct2( in )

% get input matrix size
N = size(in,1);

% build the matrix
n = 0:N-1;
for k = 0:N-1
   if (k>0)
      C(k+1,n+1) = cos(pi*(2*n+1)*k/2/N)/sqrt(N)*sqrt(2);
   else
      C(k+1,n+1) = cos(pi*(2*n+1)*k/2/N)/sqrt(N);
   end   
end
out = (C')*in*C;
end


% ------------------------------------------------------------------------
% Classic_DCT - perform a block DCT for an image
% ------------------------------------------------------------------------
function transform_image = Classic_DCT( input_image )

transform_image = zeros( size( input_image,1 ),size( input_image,2 ) );