I am sharing this simple JPEG encoder on code by the following copyright holder. Procured the file from matlab file exchange site. I have gone through the source of JPEG encoder and modified code so that encoding and decoding is done for one file.
(Please scroll through this copyright notification for code with explanation)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Copyright (c) 2009, Muhammad Anwarul Azim
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
File Name: jpeg.m
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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 ) );
for m = 0:15
for n = 0:15
transform_image( m*8+[1:8],n*8+[1:8] ) = Classic_DCT_Block_8x8( input_image( m*8+[1:8],n*8+[1:8] ) );
end
end
end
% ------------------------------------------------------------------------
% image_8x8_block_flowgraph_forward_dct - perform a block Flowgraph forward DCT for an image
% ------------------------------------------------------------------------
function transform_image = image_8x8_block_flowgraph_forward_dct( input_image )
transform_image = zeros( size( input_image,1 ),size( input_image,2 ) );
for m = 0:15
for n = 0:15
transform_image( m*8+[1:8],n*8+[1:8] ) = flowgraph_forward_dct( input_image( m*8+[1:8],n*8+[1:8] ) );
end
end
end
% ------------------------------------------------------------------------
% image_8x8_block_inv_dct - perform a block inverse DCT for an image
% ------------------------------------------------------------------------
function restored_image = image_8x8_block_inv_dct( transform_image )
restored_image = zeros( size( transform_image,1 ),size( transform_image,2 ) );
for m = 0:15
for n = 0:15
restored_image( m*8+[1:8],n*8+[1:8] ) = pdip_inv_dct2( transform_image( m*8+[1:8],n*8+[1:8] ) );
end
end
end
% ------------------------------------------------------------------------
% image_8x8_block_flowgraph_inverse_dct - perform a block Flowgraph inverse DCT for an image
% ------------------------------------------------------------------------
function restored_image = image_8x8_block_flowgraph_inverse_dct( transform_image )
restored_image = zeros( size( transform_image,1 ),size( transform_image,2 ) );
for m = 0:15
for n = 0:15
restored_image( m*8+[1:8],n*8+[1:8] ) = flowgraph_inverse_dct( transform_image( m*8+[1:8],n*8+[1:8] ) );
end
end
end
% ------------------------------------------------------------------------
% FLOWGRAPH forward dct (Chen,Fralick and Smith)
% ------------------------------------------------------------------------
function [DCT_8x8] = flowgraph_forward_dct(in_8x8)
% constant cosine values will be used for both forward & inverse flowgraph DCT
c1=0.980785;
c2=0.923880;
c3=0.831470;
c4=0.707107;
c5=0.555570;
c6=0.382683;
c7=0.195090;
%---------------------------row calculation FDCT--------------------------
for row_number=1:8
%sample image value initialization from input matrix
f0=in_8x8(row_number,1);
f1=in_8x8(row_number,2);
f2=in_8x8(row_number,3);
f3=in_8x8(row_number,4);
f4=in_8x8(row_number,5);
f5=in_8x8(row_number,6);
f6=in_8x8(row_number,7);
f7=in_8x8(row_number,8);
%first stage of FLOWGRAPH (Chen,Fralick and Smith)
i0=f0+f7;
i1=f1+f6;
i2=f2+f5;
i3=f3+f4;
i4=f3-f4;
i5=f2-f5;
i6=f1-f6;
i7=f0-f7;
%second stage of FLOWGRAPH (Chen,Fralick and Smith)
j0=i0+i3;
j1=i1+i2;
j2=i1-i2;
j3=i0-i3;
j4=i4;
j5=(i6-i5)*c4;
j6=(i6+i5)*c4;
j7=i7;
%third stage of FLOWGRAPH (Chen,Fralick and Smith)
k0=(j0+j1)*c4;
k1=(j0-j1)*c4;
k2=(j2*c6)+(j3*c2);
k3=(j3*c6)-(j2*c2);
k4=j4+j5;
k5=j4-j5;
k6=j7-j6;
k7=j7+j6;
%fourth stage of FLOWGRAPH; 1-dimensional DCT coefficients
F0=k0/2;
F1=(k4*c7+k7*c1)/2;
F2=k2/2;
F3=(k6*c3-k5*c5)/2;
F4=k1/2;
F5=(k5*c3+k6*c5)/2;
F6=k3/2;
F7=(k7*c7-k4*c1)/2;
%DCT coefficient assignment
One_D_DCT_Row_8x8(row_number,1)=F0;
One_D_DCT_Row_8x8(row_number,2)=F1;
One_D_DCT_Row_8x8(row_number,3)=F2;
One_D_DCT_Row_8x8(row_number,4)=F3;
One_D_DCT_Row_8x8(row_number,5)=F4;
One_D_DCT_Row_8x8(row_number,6)=F5;
One_D_DCT_Row_8x8(row_number,7)=F6;
One_D_DCT_Row_8x8(row_number,8)=F7;
end %end of row calculations
%---------------------------end: row calculation FDCT---------------------
%--------------------------- column calculation FDCT----------------------
for column_number=1:8 %start of column calculation
%sample image value initialization
f0=One_D_DCT_Row_8x8(1,column_number);
f1=One_D_DCT_Row_8x8(2,column_number);
f2=One_D_DCT_Row_8x8(3,column_number);
f3=One_D_DCT_Row_8x8(4,column_number);
f4=One_D_DCT_Row_8x8(5,column_number);
f5=One_D_DCT_Row_8x8(6,column_number);
f6=One_D_DCT_Row_8x8(7,column_number);
f7=One_D_DCT_Row_8x8(8,column_number);
%first stage of FLOWGRAPH (Chen,Fralick and Smith)
i0=f0+f7;
i1=f1+f6;
i2=f2+f5;
i3=f3+f4;
i4=f3-f4;
i5=f2-f5;
i6=f1-f6;
i7=f0-f7;
%second stage of FLOWGRAPH (Chen,Fralick and Smith)
j0=i0+i3;
j1=i1+i2;
j2=i1-i2;
j3=i0-i3;
j4=i4;
j5=(i6-i5)*c4;
j6=(i6+i5)*c4;
j7=i7;
%third stage of FLOWGRAPH (Chen,Fralick and Smith)
k0=(j0+j1)*c4;
k1=(j0-j1)*c4;
k2=(j2*c6)+(j3*c2);
k3=(j3*c6)-(j2*c2);
k4=j4+j5;
k5=j4-j5;
k6=j7-j6;
k7=j7+j6;
%fourth stage of FLOWGRAPH; Desired DCT coefficients
F0=k0/2;
F1=(k4*c7+k7*c1)/2;
F2=k2/2;
F3=(k6*c3-k5*c5)/2;
F4=k1/2;
F5=(k5*c3+k6*c5)/2;
F6=k3/2;
F7=(k7*c7-k4*c1)/2;
%DCT coefficient assignment
DCT_8x8(1,column_number)=F0;
DCT_8x8(2,column_number)=F1;
DCT_8x8(3,column_number)=F2;
DCT_8x8(4,column_number)=F3;
DCT_8x8(5,column_number)=F4;
DCT_8x8(6,column_number)=F5;
DCT_8x8(7,column_number)=F6;
DCT_8x8(8,column_number)=F7;
end %end of column calculations
%---------------------------end: column calculation FDCT------------------
end % end of function flowgraph_forward_dct
% ------------------------------------------------------------------------
% FLOWGRAPH Inverse dct (Chen,Fralick and Smith)
% ------------------------------------------------------------------------
function [out_8x8] = flowgraph_inverse_dct(DCT_8x8)
% constant cosine values will be used for both forward & inverse flowgraph DCT
c1=0.980785;
c2=0.923880;
c3=0.831470;
c4=0.707107;
c5=0.555570;
c6=0.382683;
c7=0.195090;
%---------------------------row calculation Inverse DCT-------------------
for row_number=1:8
%DCT coefficient initialization
F0=DCT_8x8(row_number,1);
F1=DCT_8x8(row_number,2);
F2=DCT_8x8(row_number,3);
F3=DCT_8x8(row_number,4);
F4=DCT_8x8(row_number,5);
F5=DCT_8x8(row_number,6);
F6=DCT_8x8(row_number,7);
F7=DCT_8x8(row_number,8);
% first stage of FLOWGRAPH (Chen,Fralick and Smith)
k0=F0/2;
k1=F4/2;
k2=F2/2;
k3=F6/2;
k4=(F1/2*c7-F7/2*c1);
k5=(F5/2*c3-F3/2*c5);
k6=F5/2*c5+F3/2*c3;
k7=F1/2*c1+F7/2*c7;
% second stage of FLOWGRAPH (Chen,Fralick and Smith)
j0=(k0+k1)*c4;
j1=(k0-k1)*c4;
j2=(k2*c6-k3*c2);
j3=k2*c2+k3*c6;
j4=k4+k5;
j5=(k4-k5);
j6=(k7-k6);
j7=k7+k6;
% third stage of FLOWGRAPH (Chen,Fralick and Smith)
i0=j0+j3;
i1=j1+j2;
i2=(j1-j2);
i3=(j0-j3);
i4=j4;
i5=(j6-j5)*c4;
i6=(j5+j6)*c4;
i7=j7;
% fourth stage of FLOWGRAPH (Chen,Fralick and Smith)
f0=i0+i7;
f1=i1+i6;
f2=i2+i5;
f3=i3+i4;
f4=(i3-i4);
f5=(i2-i5);
f6=(i1-i6);
f7=(i0-i7);
%1 dimensional sample image vale assignment only after row calculations
One_D_IDCT_Row_8x8(row_number,1)=f0;
One_D_IDCT_Row_8x8(row_number,2)=f1;
One_D_IDCT_Row_8x8(row_number,3)=f2;
One_D_IDCT_Row_8x8(row_number,4)=f3;
One_D_IDCT_Row_8x8(row_number,5)=f4;
One_D_IDCT_Row_8x8(row_number,6)=f5;
One_D_IDCT_Row_8x8(row_number,7)=f6;
One_D_IDCT_Row_8x8(row_number,8)=f7;
end
%---------------------------end: row calculation Inverse DCT--------------
%---------------------------column calculation Inverse DCT----------------
for column_number=1:8
%DCT coefficient initialization
F0=One_D_IDCT_Row_8x8(1,column_number);
F1=One_D_IDCT_Row_8x8(2,column_number);
F2=One_D_IDCT_Row_8x8(3,column_number);
F3=One_D_IDCT_Row_8x8(4,column_number);
F4=One_D_IDCT_Row_8x8(5,column_number);
F5=One_D_IDCT_Row_8x8(6,column_number);
F6=One_D_IDCT_Row_8x8(7,column_number);
F7=One_D_IDCT_Row_8x8(8,column_number);
% first stage of FLOWGRAPH (Chen,Fralick and Smith)
k0=F0/2;
k1=F4/2;
k2=F2/2;
k3=F6/2;
k4=(F1/2*c7-F7/2*c1);
k5=(F5/2*c3-F3/2*c5);
k6=F5/2*c5+F3/2*c3;
k7=F1/2*c1+F7/2*c7;
% second stage of FLOWGRAPH (Chen,Fralick and Smith)
j0=(k0+k1)*c4;
j1=(k0-k1)*c4;
j2=(k2*c6-k3*c2);
j3=k2*c2+k3*c6;
j4=k4+k5;
j5=(k4-k5);
j6=(k7-k6);
j7=k7+k6;
% third stage of FLOWGRAPH (Chen,Fralick and Smith)
i0=j0+j3;
i1=j1+j2;
i2=(j1-j2);
i3=(j0-j3);
i4=j4;
i5=(j6-j5)*c4;
i6=(j5+j6)*c4;
i7=j7;
% fourth stage of FLOWGRAPH (Chen,Fralick and Smith)
f0=i0+i7;
f1=i1+i6;
f2=i2+i5;
f3=i3+i4;
f4=(i3-i4);
f5=(i2-i5);
f6=(i1-i6);
f7=(i0-i7);
% Desired sample image values assignment only after 2 dimensional inverse transformation
out_8x8(1,column_number)=f0;
out_8x8(2,column_number)=f1;
out_8x8(3,column_number)=f2;
out_8x8(4,column_number)=f3;
out_8x8(5,column_number)=f4;
out_8x8(6,column_number)=f5;
out_8x8(7,column_number)=f6;
out_8x8(8,column_number)=f7;
end
%---------------------------end: column calculation Inverse DCT-----------
end % end of function flowgraph_inverse_dct
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Example Input image is: (You can download it and convert it to .tiff of size 100x100)
The output after Compression and Decompression is:
(you could have also got this code from github)
(Please scroll through this copyright notification for code with explanation)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Copyright (c) 2009, Muhammad Anwarul Azim
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution
* Neither the name of the Saitama University nor the names
of its contributors may be used to endorse or promote products derived
from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
File Name: jpeg.m
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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 ) );
for m = 0:15
for n = 0:15
transform_image( m*8+[1:8],n*8+[1:8] ) = Classic_DCT_Block_8x8( input_image( m*8+[1:8],n*8+[1:8] ) );
end
end
end
% ------------------------------------------------------------------------
% image_8x8_block_flowgraph_forward_dct - perform a block Flowgraph forward DCT for an image
% ------------------------------------------------------------------------
function transform_image = image_8x8_block_flowgraph_forward_dct( input_image )
transform_image = zeros( size( input_image,1 ),size( input_image,2 ) );
for m = 0:15
for n = 0:15
transform_image( m*8+[1:8],n*8+[1:8] ) = flowgraph_forward_dct( input_image( m*8+[1:8],n*8+[1:8] ) );
end
end
end
% ------------------------------------------------------------------------
% image_8x8_block_inv_dct - perform a block inverse DCT for an image
% ------------------------------------------------------------------------
function restored_image = image_8x8_block_inv_dct( transform_image )
restored_image = zeros( size( transform_image,1 ),size( transform_image,2 ) );
for m = 0:15
for n = 0:15
restored_image( m*8+[1:8],n*8+[1:8] ) = pdip_inv_dct2( transform_image( m*8+[1:8],n*8+[1:8] ) );
end
end
end
% ------------------------------------------------------------------------
% image_8x8_block_flowgraph_inverse_dct - perform a block Flowgraph inverse DCT for an image
% ------------------------------------------------------------------------
function restored_image = image_8x8_block_flowgraph_inverse_dct( transform_image )
restored_image = zeros( size( transform_image,1 ),size( transform_image,2 ) );
for m = 0:15
for n = 0:15
restored_image( m*8+[1:8],n*8+[1:8] ) = flowgraph_inverse_dct( transform_image( m*8+[1:8],n*8+[1:8] ) );
end
end
end
% ------------------------------------------------------------------------
% FLOWGRAPH forward dct (Chen,Fralick and Smith)
% ------------------------------------------------------------------------
function [DCT_8x8] = flowgraph_forward_dct(in_8x8)
% constant cosine values will be used for both forward & inverse flowgraph DCT
c1=0.980785;
c2=0.923880;
c3=0.831470;
c4=0.707107;
c5=0.555570;
c6=0.382683;
c7=0.195090;
%---------------------------row calculation FDCT--------------------------
for row_number=1:8
%sample image value initialization from input matrix
f0=in_8x8(row_number,1);
f1=in_8x8(row_number,2);
f2=in_8x8(row_number,3);
f3=in_8x8(row_number,4);
f4=in_8x8(row_number,5);
f5=in_8x8(row_number,6);
f6=in_8x8(row_number,7);
f7=in_8x8(row_number,8);
%first stage of FLOWGRAPH (Chen,Fralick and Smith)
i0=f0+f7;
i1=f1+f6;
i2=f2+f5;
i3=f3+f4;
i4=f3-f4;
i5=f2-f5;
i6=f1-f6;
i7=f0-f7;
%second stage of FLOWGRAPH (Chen,Fralick and Smith)
j0=i0+i3;
j1=i1+i2;
j2=i1-i2;
j3=i0-i3;
j4=i4;
j5=(i6-i5)*c4;
j6=(i6+i5)*c4;
j7=i7;
%third stage of FLOWGRAPH (Chen,Fralick and Smith)
k0=(j0+j1)*c4;
k1=(j0-j1)*c4;
k2=(j2*c6)+(j3*c2);
k3=(j3*c6)-(j2*c2);
k4=j4+j5;
k5=j4-j5;
k6=j7-j6;
k7=j7+j6;
%fourth stage of FLOWGRAPH; 1-dimensional DCT coefficients
F0=k0/2;
F1=(k4*c7+k7*c1)/2;
F2=k2/2;
F3=(k6*c3-k5*c5)/2;
F4=k1/2;
F5=(k5*c3+k6*c5)/2;
F6=k3/2;
F7=(k7*c7-k4*c1)/2;
%DCT coefficient assignment
One_D_DCT_Row_8x8(row_number,1)=F0;
One_D_DCT_Row_8x8(row_number,2)=F1;
One_D_DCT_Row_8x8(row_number,3)=F2;
One_D_DCT_Row_8x8(row_number,4)=F3;
One_D_DCT_Row_8x8(row_number,5)=F4;
One_D_DCT_Row_8x8(row_number,6)=F5;
One_D_DCT_Row_8x8(row_number,7)=F6;
One_D_DCT_Row_8x8(row_number,8)=F7;
end %end of row calculations
%---------------------------end: row calculation FDCT---------------------
%--------------------------- column calculation FDCT----------------------
for column_number=1:8 %start of column calculation
%sample image value initialization
f0=One_D_DCT_Row_8x8(1,column_number);
f1=One_D_DCT_Row_8x8(2,column_number);
f2=One_D_DCT_Row_8x8(3,column_number);
f3=One_D_DCT_Row_8x8(4,column_number);
f4=One_D_DCT_Row_8x8(5,column_number);
f5=One_D_DCT_Row_8x8(6,column_number);
f6=One_D_DCT_Row_8x8(7,column_number);
f7=One_D_DCT_Row_8x8(8,column_number);
%first stage of FLOWGRAPH (Chen,Fralick and Smith)
i0=f0+f7;
i1=f1+f6;
i2=f2+f5;
i3=f3+f4;
i4=f3-f4;
i5=f2-f5;
i6=f1-f6;
i7=f0-f7;
%second stage of FLOWGRAPH (Chen,Fralick and Smith)
j0=i0+i3;
j1=i1+i2;
j2=i1-i2;
j3=i0-i3;
j4=i4;
j5=(i6-i5)*c4;
j6=(i6+i5)*c4;
j7=i7;
%third stage of FLOWGRAPH (Chen,Fralick and Smith)
k0=(j0+j1)*c4;
k1=(j0-j1)*c4;
k2=(j2*c6)+(j3*c2);
k3=(j3*c6)-(j2*c2);
k4=j4+j5;
k5=j4-j5;
k6=j7-j6;
k7=j7+j6;
%fourth stage of FLOWGRAPH; Desired DCT coefficients
F0=k0/2;
F1=(k4*c7+k7*c1)/2;
F2=k2/2;
F3=(k6*c3-k5*c5)/2;
F4=k1/2;
F5=(k5*c3+k6*c5)/2;
F6=k3/2;
F7=(k7*c7-k4*c1)/2;
%DCT coefficient assignment
DCT_8x8(1,column_number)=F0;
DCT_8x8(2,column_number)=F1;
DCT_8x8(3,column_number)=F2;
DCT_8x8(4,column_number)=F3;
DCT_8x8(5,column_number)=F4;
DCT_8x8(6,column_number)=F5;
DCT_8x8(7,column_number)=F6;
DCT_8x8(8,column_number)=F7;
end %end of column calculations
%---------------------------end: column calculation FDCT------------------
end % end of function flowgraph_forward_dct
% ------------------------------------------------------------------------
% FLOWGRAPH Inverse dct (Chen,Fralick and Smith)
% ------------------------------------------------------------------------
function [out_8x8] = flowgraph_inverse_dct(DCT_8x8)
% constant cosine values will be used for both forward & inverse flowgraph DCT
c1=0.980785;
c2=0.923880;
c3=0.831470;
c4=0.707107;
c5=0.555570;
c6=0.382683;
c7=0.195090;
%---------------------------row calculation Inverse DCT-------------------
for row_number=1:8
%DCT coefficient initialization
F0=DCT_8x8(row_number,1);
F1=DCT_8x8(row_number,2);
F2=DCT_8x8(row_number,3);
F3=DCT_8x8(row_number,4);
F4=DCT_8x8(row_number,5);
F5=DCT_8x8(row_number,6);
F6=DCT_8x8(row_number,7);
F7=DCT_8x8(row_number,8);
% first stage of FLOWGRAPH (Chen,Fralick and Smith)
k0=F0/2;
k1=F4/2;
k2=F2/2;
k3=F6/2;
k4=(F1/2*c7-F7/2*c1);
k5=(F5/2*c3-F3/2*c5);
k6=F5/2*c5+F3/2*c3;
k7=F1/2*c1+F7/2*c7;
% second stage of FLOWGRAPH (Chen,Fralick and Smith)
j0=(k0+k1)*c4;
j1=(k0-k1)*c4;
j2=(k2*c6-k3*c2);
j3=k2*c2+k3*c6;
j4=k4+k5;
j5=(k4-k5);
j6=(k7-k6);
j7=k7+k6;
% third stage of FLOWGRAPH (Chen,Fralick and Smith)
i0=j0+j3;
i1=j1+j2;
i2=(j1-j2);
i3=(j0-j3);
i4=j4;
i5=(j6-j5)*c4;
i6=(j5+j6)*c4;
i7=j7;
% fourth stage of FLOWGRAPH (Chen,Fralick and Smith)
f0=i0+i7;
f1=i1+i6;
f2=i2+i5;
f3=i3+i4;
f4=(i3-i4);
f5=(i2-i5);
f6=(i1-i6);
f7=(i0-i7);
%1 dimensional sample image vale assignment only after row calculations
One_D_IDCT_Row_8x8(row_number,1)=f0;
One_D_IDCT_Row_8x8(row_number,2)=f1;
One_D_IDCT_Row_8x8(row_number,3)=f2;
One_D_IDCT_Row_8x8(row_number,4)=f3;
One_D_IDCT_Row_8x8(row_number,5)=f4;
One_D_IDCT_Row_8x8(row_number,6)=f5;
One_D_IDCT_Row_8x8(row_number,7)=f6;
One_D_IDCT_Row_8x8(row_number,8)=f7;
end
%---------------------------end: row calculation Inverse DCT--------------
%---------------------------column calculation Inverse DCT----------------
for column_number=1:8
%DCT coefficient initialization
F0=One_D_IDCT_Row_8x8(1,column_number);
F1=One_D_IDCT_Row_8x8(2,column_number);
F2=One_D_IDCT_Row_8x8(3,column_number);
F3=One_D_IDCT_Row_8x8(4,column_number);
F4=One_D_IDCT_Row_8x8(5,column_number);
F5=One_D_IDCT_Row_8x8(6,column_number);
F6=One_D_IDCT_Row_8x8(7,column_number);
F7=One_D_IDCT_Row_8x8(8,column_number);
% first stage of FLOWGRAPH (Chen,Fralick and Smith)
k0=F0/2;
k1=F4/2;
k2=F2/2;
k3=F6/2;
k4=(F1/2*c7-F7/2*c1);
k5=(F5/2*c3-F3/2*c5);
k6=F5/2*c5+F3/2*c3;
k7=F1/2*c1+F7/2*c7;
% second stage of FLOWGRAPH (Chen,Fralick and Smith)
j0=(k0+k1)*c4;
j1=(k0-k1)*c4;
j2=(k2*c6-k3*c2);
j3=k2*c2+k3*c6;
j4=k4+k5;
j5=(k4-k5);
j6=(k7-k6);
j7=k7+k6;
% third stage of FLOWGRAPH (Chen,Fralick and Smith)
i0=j0+j3;
i1=j1+j2;
i2=(j1-j2);
i3=(j0-j3);
i4=j4;
i5=(j6-j5)*c4;
i6=(j5+j6)*c4;
i7=j7;
% fourth stage of FLOWGRAPH (Chen,Fralick and Smith)
f0=i0+i7;
f1=i1+i6;
f2=i2+i5;
f3=i3+i4;
f4=(i3-i4);
f5=(i2-i5);
f6=(i1-i6);
f7=(i0-i7);
% Desired sample image values assignment only after 2 dimensional inverse transformation
out_8x8(1,column_number)=f0;
out_8x8(2,column_number)=f1;
out_8x8(3,column_number)=f2;
out_8x8(4,column_number)=f3;
out_8x8(5,column_number)=f4;
out_8x8(6,column_number)=f5;
out_8x8(7,column_number)=f6;
out_8x8(8,column_number)=f7;
end
%---------------------------end: column calculation Inverse DCT-----------
end % end of function flowgraph_inverse_dct
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Example Input image is: (You can download it and convert it to .tiff of size 100x100)
The output after Compression and Decompression is:
(you could have also got this code from github)
No comments:
Post a Comment