clear all;

% Setup data

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


step_list = [ 1 2 ]%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 );
weight_scale = step_size / 16

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;
  

  M = 8;

  bit_total = 0;
  
  current_level = image;
  Current_Size = size( current_level );
  
  reconstructed_level = zeros( size( image ) );
  
  level_list = [];
  level_index = [];

  % Make list of ideal levels:
  
  while rem( Current_Size, 2 ) == zeros( 1, 2 )
    level_index = mAddRowsWithZeroPadding( level_index, [ size( level_list, 1 )+1, Current_Size ] );  
    level_list = mAddRowsWithZeroPadding( level_list, current_level );
    current_level = mSubsampleAverage( current_level );
    Current_Size = size( current_level );
  end
  
%  level_index = mAddRowsWithZeroPadding( level_index, [ size( level_list, 1 )+1, Current_Size ] );  
%  level_list = mAddRowsWithZeroPadding( level_list, current_level );

  level_count = size( level_index, 1 );
  level_sum = zeros( size( current_level ) );
  
  Current_Size
  
  while level_count > 0
    upper_level_row = level_index( level_count, 1 );
    upper_level_height = level_index( level_count, 2 );
    upper_level_width = level_index( level_count, 3 );
    upper_level = level_list( upper_level_row:upper_level_row+upper_level_height-1, 1:upper_level_width );
    
    [ reconstructed_level, quantized_coefficients ] = mOptimalDCT8x8ImageAndCoeffs( current_level, weight_scale );
    
    [ PDF PDF_range RLE_count ] = mComputeRLEPDF( quantized_coefficients, 1, M );
    PlogP = log( PDF.^PDF)./log(2);
    Level_Rate = -sum(sum( PlogP' ))*RLE_count/prod( Current_Size )
    Level_Bits = ceil( Level_Rate*prod( Current_Size ) );
    bit_total = bit_total + Level_Bits
    
    level_sum = mSupersampleAverage( level_sum + reconstructed_level );
    if draw
      mShowImage( level_sum ); 
      drawnow;
    end

    current_level = upper_level - level_sum;
    Current_Size = size( current_level )
    
    level_count = level_count - 1;
  end
  
  [ reconstructed_level, quantized_coefficients ] = mOptimalDCT8x8ImageAndCoeffs( current_level, weight_scale );
    
  [ PDF PDF_range RLE_count ] = mComputeRLEPDF( quantized_coefficients, 1, M );
  PlogP = log( PDF.^PDF)./log(2);
  Level_Rate = -sum(sum( PlogP' ))*RLE_count/prod( Current_Size )
  Level_Bits = ceil( Level_Rate*prod( Current_Size ) );
  bit_total = bit_total + Level_Bits
  
  level_sum = level_sum + reconstructed_level;
  if draw
    mShowImage( level_sum ); 
    drawnow;
  end
  
  Rate = bit_total / image_size
    
  signal = mean( mean( image.^2 ) );
  noise  = mean( mean( (image-level_sum).^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, 1, 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;
  
  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 )
xp=results(1,:);
yp=results(2,:);