clear all;
clf reset;


% Setup data

image_list = mGetImageFileList;
	   
image_total = 2%size(image_list,1);

size_list = zeros( 1, image_total );
MSE_list = zeros( 1, image_total );
signal_list = zeros( 1, image_total );
noise_list = zeros( 1, image_total );

SNR_Results        = ones( image_total+1, 1 )*NaN; 
Rate_Results       = ones( image_total+1, 1 )*NaN;


step_size=32
RLE_list = [];
Total_DCT_Time = 0;

for image_count = 1:image_total
  
  image_name = deblank( image_list( image_count,: ) )
  image = pgmRead( sprintf( '../images/%s', image_name) );

  figure( 1 );
  mShowImage( image );
  
  image_size = prod( size( image ) );
  size_list( image_count ) = image_size;
  
  M = 16;

  % DCT Image: start
  tic;
  [ coefficients, image_size_vector ] = mDCT8x8( image );
  
 
  DC_coefficients = coefficients( 1:M:image_size_vector( 1 ), 1:M:image_size_vector( 2 ) );
  quantized_coefficients = quant( coefficients, step_size );
  quantized_coefficients( 1:M:image_size_vector( 1 ), 1:M:image_size_vector( 2 ) ) = DC_coefficients;
  
  reconstructed_image = round( mInverseDCT8x8( quantized_coefficients, image_size_vector ) );
  
  Total_DCT_Time = Total_DCT_Time + toc
  % DCT Image: end
  
  signal = mean( mean( image.^2 ) );
  noise  = mean( mean( (image-reconstructed_image).^2 ) );
  
  signal_list( image_count ) = signal;
  noise_list( image_count ) = noise;
  
  SNR = 10*log10( signal/noise )
  SNR_Results( image_count ) = SNR;
  
  [ PDF_update range_update RLE_count RLE_update ] = ComputeRLEPDF( quantized_coefficients, step_size, M );
  
  RLE_list = [ RLE_list RLE_update ];
  
  % x*log(x) = log(x^x)
  PlogP = log( PDF_update.^PDF_update)./log(2);
  Rate = -sum(sum( PlogP' ))*RLE_count/image_size
  Rate_Results( image_count ) = Rate;
  rate2=mCompute2DEntropy(RLE_update )*size(RLE_update,2)/image_size
  Separate_Rate=( mComputeEntropy( real( RLE_update) )+ mComputeEntropy( imag( RLE_update) ) )*size(RLE_update,2)/image_size
  
  figure( 2);
  mShowImage(reconstructed_image);
  drawnow;
  
end


Overall_signal = sum( signal_list.*size_list ) / sum( size_list );
Overall_noise = sum( noise_list.*size_list ) / sum( size_list );

Overall_SNR = 10*log10( Overall_signal / Overall_noise )

SNR_Results( image_total+1 ) = Overall_SNR;

%Overall_PDF = Overall_PDF/sum(size_list);

% x*log(x) = log(x^x)
%PlogP = log( Overall_PDF.^Overall_PDF)./log(2);
%Overall_Rate = -sum(sum( PlogP' ))
Overall_Rate = mCompute2DEntropy( RLE_list )*size(RLE_list,2)/sum( size_list );

Rate_Results( image_total+1 ) = Overall_Rate;

SNR_Results = [ NaN, step_size; [ 1:image_total NaN ]' SNR_Results ]
Rate_Results = [ NaN, step_size; [ 1:image_total NaN ]' Rate_Results ]

Total_DCT_Time