clear all;

% Setup data

image_list = mGetLeftImageFileList;
image_total = size(image_list,1);


step_list = [ 4 8 16 32 64 128 256 ]
results = zeros( 2, length( step_list ) );
for step_index = 1:length(step_list)
  
step_size = step_list( step_index )

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;

Overall_PDF = zeros( 64, 1 );
range = [ 0 0 ];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

draw = 0;


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) );
  
  if draw
    figure( 1 );
    mShowImage( image );
  end
  
  image_size = prod( size( image ) );
  size_list( image_count ) = image_size;
  

  image_mean = round( mean(mean(image)) );

  
  M = 8;

  % DCT Image: start
  tic;
  [ coefficients, image_size_vector ] = mDCT8x8( image - image_mean );
  
 
%  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 ) ) + image_mean;
  
  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;
%  Separate_Rate=( mComputeEntropy( real( RLE_update) )+ mComputeEntropy( imag( RLE_update) ) )*size(RLE_update,2)/image_size
  
  if draw
    figure( 2);
    mShowImage(reconstructed_image);
    drawnow;
  end
  
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%Overall_SNR = 10*log10( Overall_signal / Overall_noise )
Overall_SNR = 10*log10( mean( signal_list ./ noise_list ) )
SNR_Results( image_total+1 ) = Overall_SNR;

Overall_Rate = mean( Rate_Results( 1:image_total ) )
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 ]

results( :,step_index ) = [ Overall_SNR; Overall_Rate ];

end

results = fliplr( results )
xr2=results(1,:);
yr2=results(2,:);