clear all;


% Setup data

image_list = mGetLeftImageFileList;
	   
image_total = 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;


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


step_size=16
draw = 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
  
  size_list( image_count ) = prod( size( image ) );
  
  [ coefficients, image_size ] = mDCT8x8( image );
  
  M = 8;
  
  quantized_coefficients = quant( coefficients, step_size );
  reconstructed_image = round( mInverseDCT8x8( quantized_coefficients, image_size ) );
  
  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 ] = ComputeMxMPDF( quantized_coefficients, step_size, M );
  
  %    PDF_update_problems = sum( abs(sum(PDF_update') - ones( 1, M^2 )) > 1e-10 )
  
  new_range_min = min( [ range(1) range_update(1) ] );
  new_range_max = max( [ range(2) range_update(2) ] );
  
  PDF_update = [ zeros( M^2, range_update(1)-new_range_min ), PDF_update, zeros( M^2, new_range_max - range_update(2) ) ];
  
  Overall_PDF = [ zeros( M^2, range(1)-new_range_min ), Overall_PDF, zeros( M^2, new_range_max - range(2) ) ];
  
  Overall_PDF = Overall_PDF + PDF_update*size_list( image_count );
  
  range = [ new_range_min new_range_max ];
  
  % x*log(x) = log(x^x)
  PlogP = log( PDF_update.^PDF_update)./log(2);
  coefficient_rates = -sum( PlogP' );
  Rate = mean( coefficient_rates )
  Rate_Results( image_count ) = Rate;
  
  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 ]