clear all;
clf reset;


% Setup data

image_list = GetImageFileList;
	   
image_skip = 3;
image_total = size(1:image_skip:size( image_list, 1 ),2 )

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

step_start  = 1;
step_factor = 2;
step_total  = 4;

step_size_list     = ones( 1, step_total )*NaN;
MSE_Results        = ones( image_total+1, step_total )*NaN;
SNR_Results        = ones( image_total+1, step_total )*NaN; 
Rate_Results       = ones( image_total+1, step_total )*NaN;
Impairment_Results = ones( image_total+1, step_total )*NaN;

for step_count = 1:step_total
  
  step_size = 2^(step_start+step_count-1)
  
  step_size_list( step_count ) = step_size;
  
  PDF = zeros( 64, 1 );
  range = [ 0 0 ];

  
  image_count = 0;
  for image_index = 1:image_skip:size( image_list, 1 );
    image_count = image_count+1;
     
    image_name = deblank( image_list( image_index,: ) )
    image = pgmRead( sprintf( '../images/%s', image_name) );
    
    figure( 1 );
    ShowImage( image );
    drawnow;
    
    size_list( image_count ) = prod( size( image ) );
    
    
    [ coefficients, image_size ] = DCT8x8( image );
    
    M = 8;
    
    quantized_coefficients = UniformQuantize( coefficients, step_size );
    
    reconstructed_image = round( InverseDCT8x8( quantized_coefficients, image_size ) );
    
    MSE = mse(image, reconstructed_image )
    MSE_Results( image_count, step_count ) = MSE;
    
    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, step_count ) = SNR;
    
    %  SNR = ComputeSNR( image, reconstructed_image )
    
    [ 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) ) ];
    
    PDF = [ zeros( M^2, range(1)-new_range_min ), PDF, zeros( M^2, new_range_max - range(2) ) ];
    
    PDF = PDF + PDF_update;
    
    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, step_count ) = Rate;
    %max(max(abs(reconstructed_image-image)))
    
    figure( 2);
    ShowImage(reconstructed_image);
    Impairment_Results( image_count, step_count ) = input( 'Enter impairment measure: ' );
    
  end
  

  Overall_MSE = sum( MSE_Results( 1:image_total, step_count )'.*size_list ) / sum( size_list )
  
  MSE_Results( image_total+1, step_count ) = Overall_MSE;
  
  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, step_count ) = Overall_SNR;
  
  PDF = PDF/image_total;
  
%  Overall_PDF_problems = sum( abs(sum(PDF') - ones( 1, M^2 )) > 1e-10 )
  
  % x*log(x) = log(x^x)
  PlogP = log( PDF.^PDF)./log(2);
  
  coefficient_rates = -sum( PlogP' );
  Overall_Rate = mean( coefficient_rates )
  
  Rate_Results( image_total+1, step_count ) = Overall_Rate;

  Impairment_Results( image_count+1, step_count ) = mean( Impairment_Results( 1:image_count, step_count ) );

end
  


MSE_Results = [ NaN, step_size_list; [ 1:image_total NaN ]' MSE_Results ]'
SNR_Results = [ NaN, step_size_list; [ 1:image_total NaN ]' SNR_Results ]'
Rate_Results = [ NaN, step_size_list; [ 1:image_total NaN ]' Rate_Results ]'
Impairment_Results = [ NaN, step_size_list; [ 1:image_total NaN ]' Impairment_Results ]'