clear all;
%clf reset;


% Setup data

left_image_list = mGetLeftImageFileList;
right_image_list = mGetRightImageFileList;
	   
image_total = 1%size(left_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 );

%step_list = [  8 16 64 256 ];
step_list = [  16];

results = zeros( 4, length( step_list ) );

for step_index = 1:length(step_list)
  step_size = step_list( step_index )


Left_SNR_Results        = ones( image_total+1, 1 )*NaN; 
Left_Rate_Results       = ones( image_total, 1 )*NaN;
Right_SNR_Results       = ones( image_total+1, 1 )*NaN; 
Right_Rate_Results      = ones( image_total, 1 )*NaN;


%weight_scale=.5
draw = 0;


left_RLE_list = [];
right_RLE_list = [];
Total_DCT_Time = 0;

'********************************************************************'


for image_count = 1:image_total

  step_size
  left_image_name = deblank( left_image_list( image_count,: ) )
  left_image = pgmRead( sprintf( '../images/%s', left_image_name) );
  right_image_name = deblank( right_image_list( image_count,: ) )
  right_image = pgmRead( sprintf( '../images/%s', right_image_name) );
  
  if draw
    figure( 1 );
    mShowImage( left_image );
  end
  
  image_size = prod( size( left_image ) );
  size_list( image_count ) = image_size;
  

  M = 8;

  left_bit_total = 0;
  right_bit_total = 0;
  
  left_current_level = left_image;
  right_current_level = right_image;
  Current_Size = size( left_current_level );
  
  reconstructed_left_level = zeros( size( left_image ) );
  reconstructed_right_level = zeros( size( right_image ) );
  
  left_level_list = [];
  right_level_list = [];
  level_index = [];

  % Make list of ideal levels:
  
  while rem( Current_Size, 2 ) == zeros( 1, 2 )
    level_index = mAddRowsWithZeroPadding( level_index, [ size( left_level_list, 1 )+1, Current_Size ] );  
    left_level_list = mAddRowsWithZeroPadding( left_level_list, left_current_level );
    right_level_list = mAddRowsWithZeroPadding( right_level_list, right_current_level );
    left_current_level = mSubsampleAverage( left_current_level );
    right_current_level = mSubsampleAverage( right_current_level );
    Current_Size = size( left_current_level );
  end
  
  level_count = size( level_index, 1 );
  left_level_sum = zeros( size( left_current_level ) );
  right_level_sum = left_level_sum;
  
%%%%  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 );
    left_upper_level = left_level_list( upper_level_row:upper_level_row+upper_level_height-1, 1:upper_level_width );
    right_upper_level = right_level_list( upper_level_row:upper_level_row+upper_level_height-1, 1:upper_level_width );
    
%    [ reconstructed_left_level, quantized_left_coefficients ] = mOptimalDCT8x8ImageAndCoeffs( left_current_level, weight_scale );
    [ reconstructed_left_level, quantized_left_coefficients ] = mSimpleDCT8x8ImageAndCoeffs( left_current_level, step_size );
%%%%   
%    [ i_reconstructed_right_level, i_quantized_right_coefficients ] = mOptimalDCT8x8ImageAndCoeffs( right_current_level, weight_scale );
%    [ j_reconstructed_right_level, j_quantized_right_coefficients ] = mOptimalDCT8x8ImageAndCoeffs( right_current_level-left_current_level, weight_scale );
    [ i_reconstructed_right_level, i_quantized_right_coefficients ] = mSimpleDCT8x8ImageAndCoeffs( right_current_level, step_size );
%    [ j_reconstructed_right_level, j_quantized_right_coefficients ] = mSimpleDCT8x8ImageAndCoeffs( right_current_level-left_current_level, step_size );
%%%%   

    [ left_PDF left_PDF_range left_RLE_count ] = mComputeRLEPDF( quantized_left_coefficients, 1, M );
    left_PlogP = log( left_PDF.^left_PDF)./log(2);
    Left_Level_Rate = -sum(sum( left_PlogP' ))*left_RLE_count/prod( Current_Size )
    Left_Level_Bits = ceil( Left_Level_Rate*prod( Current_Size ) );
    left_bit_total = left_bit_total + Left_Level_Bits;
    
    [ i_right_PDF i_right_PDF_range i_right_RLE_count ] = mComputeRLEPDF( i_quantized_right_coefficients, 1, M );
    i_right_PlogP = log( i_right_PDF.^i_right_PDF)./log(2);
    I_Right_Level_Rate = -sum(sum( i_right_PlogP' ))*i_right_RLE_count/prod( Current_Size );
    
%    [ j_right_PDF j_right_PDF_range j_right_RLE_count ] = mComputeRLEPDF( j_quantized_right_coefficients, 1, M );
%    j_right_PlogP = log( j_right_PDF.^j_right_PDF)./log(2);
%    J_Right_Level_Rate = -sum(sum( j_right_PlogP' ))*j_right_RLE_count/prod( Current_Size );
    
%    [ I_Right_Level_Rate J_Right_Level_Rate ( I_Right_Level_Rate > J_Right_Level_Rate ) ]
    
%%%%%%%%
%    if I_Right_Level_Rate < J_Right_Level_Rate 
      Right_Level_Bits = ceil( I_Right_Level_Rate*prod( Current_Size ) );
%    else
%      Right_Level_Bits = ceil( J_Right_Level_Rate*prod( Current_Size ) );
%    end
    right_bit_total = right_bit_total + Right_Level_Bits
    
    
    left_level_sum = mSupersampleAverage( left_level_sum + reconstructed_left_level );
%%%%   
%%%%%%%%    
%    if I_Right_Level_Rate < J_Right_Level_Rate 
      right_level_sum = mSupersampleAverage( right_level_sum + i_reconstructed_right_level );
%    else
%      right_level_sum = mSupersampleAverage( right_level_sum + reconstructed_left_level + j_reconstructed_right_level );
%    end
%%%%   

    if draw
      mShowImage( left_level_sum ); 
      drawnow;
    end

    left_current_level = left_upper_level - left_level_sum;
    right_current_level = right_upper_level - right_level_sum;
    Current_Size = size( left_current_level );
    
    level_count = level_count - 1;
  end
  
%  [ reconstructed_left_level, quantized_left_coefficients ] = mOptimalDCT8x8ImageAndCoeffs( left_current_level, weight_scale );
  [ reconstructed_left_level, quantized_left_coefficients ] = mSimpleDCT8x8ImageAndCoeffs( left_current_level, step_size );
%%%%   
%  [ i_reconstructed_right_level, i_quantized_right_coefficients ] = mOptimalDCT8x8ImageAndCoeffs( right_current_level, weight_scale );
%  [ j_reconstructed_right_level, j_quantized_right_coefficients ] = mOptimalDCT8x8ImageAndCoeffs( right_current_level-left_current_level, weight_scale );
  [ i_reconstructed_right_level, i_quantized_right_coefficients ] = mSimpleDCT8x8ImageAndCoeffs( right_current_level, step_size );
%  [ j_reconstructed_right_level, j_quantized_right_coefficients ] = mSimpleDCT8x8ImageAndCoeffs( right_current_level-left_current_level, step_size );
%%%%   
    
  [ left_PDF left_PDF_range left_RLE_count ] = mComputeRLEPDF( quantized_left_coefficients, 1, M );
  left_PlogP = log( left_PDF.^left_PDF)./log(2);
  Left_Level_Rate = -sum(sum( left_PlogP' ))*left_RLE_count/prod( Current_Size )
  Left_Level_Bits = ceil( Left_Level_Rate*prod( Current_Size ) );
  left_bit_total = left_bit_total + Left_Level_Bits;
  
%%%%
    [ i_right_PDF i_right_PDF_range i_right_RLE_count ] = mComputeRLEPDF( i_quantized_right_coefficients, 1, M );
    i_right_PlogP = log( i_right_PDF.^i_right_PDF)./log(2);
    I_Right_Level_Rate = -sum(sum( i_right_PlogP' ))*i_right_RLE_count/prod( Current_Size );
    
%    [ j_right_PDF j_right_PDF_range j_right_RLE_count ] = mComputeRLEPDF( j_quantized_right_coefficients, 1, M );
%    j_right_PlogP = log( j_right_PDF.^j_right_PDF)./log(2);
%    J_Right_Level_Rate = -sum(sum( j_right_PlogP' ))*j_right_RLE_count/prod( Current_Size );
    
%    [ I_Right_Level_Rate J_Right_Level_Rate ( I_Right_Level_Rate > J_Right_Level_Rate ) ]
    
%    if I_Right_Level_Rate < J_Right_Level_Rate 
      Right_Level_Bits = ceil( I_Right_Level_Rate*prod( Current_Size ) );
%    else
%      Right_Level_Bits = ceil( J_Right_Level_Rate*prod( Current_Size ) );
%    end
    right_bit_total = right_bit_total + Right_Level_Bits
%%%%    
    

  left_level_sum = left_level_sum + reconstructed_left_level;
%%%%   
%  if I_Right_Level_Rate < J_Right_Level_Rate 
    right_level_sum = right_level_sum + i_reconstructed_right_level;
%  else
%    right_level_sum = right_level_sum + reconstructed_left_level + j_reconstructed_right_level;
%  end
%%%%   

  if draw
    mShowImage( left_level_sum ); 
    drawnow;
  end
  
  Left_Rate = left_bit_total / image_size
  Right_Rate = right_bit_total / image_size
    

  left_signal = mean( mean( left_image.^2 ) );
  left_noise  = mean( mean( (left_image-left_level_sum).^2 ) );
  
  right_signal = mean( mean( right_image.^2 ) );
  right_noise  = mean( mean( (right_image-right_level_sum).^2 ) );
  
  left_signal_list( image_count ) = left_signal;
  left_noise_list( image_count ) = left_noise;
  right_signal_list( image_count ) = right_signal;
  right_noise_list( image_count ) = right_noise;
  
  Left_SNR = 10*log10( left_signal/left_noise )
  Left_SNR_Results( image_count ) = Left_SNR;
  Right_SNR = 10*log10( right_signal/right_noise )
  Right_SNR_Results( image_count ) = Right_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

   Left_Rate_Results( image_count ) = Left_Rate;
   Right_Rate_Results( image_count ) = Right_Rate;
  
end

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

%Left_Overall_signal = sum( left_signal_list.*size_list ) / sum( size_list );
%Left_Overall_noise = sum( left_noise_list.*size_list ) / sum( size_list );
%Right_Overall_signal = sum( right_signal_list.*size_list ) / sum( size_list );
%Right_Overall_noise = sum( right_noise_list.*size_list ) / sum( size_list );
Left_Overall_signal = mean( left_signal_list );
Left_Overall_noise = mean( left_noise_list );
Right_Overall_signal = mean( right_signal_list );
Right_Overall_noise = mean( right_noise_list );

Left_Overall_SNR = 10*log10( Left_Overall_signal / Left_Overall_noise )
Right_Overall_SNR = 10*log10( Right_Overall_signal / Right_Overall_noise )

Left_SNR_Results( image_total+1 ) = Left_Overall_SNR;
Right_SNR_Results( image_total+1 ) = Right_Overall_SNR;

Left_SNR_Results = [ NaN, step_size; [ 1:image_total NaN ]' Left_SNR_Results ]
Right_SNR_Results = [ NaN, step_size; [ 1:image_total NaN ]' Right_SNR_Results ]

results( :, step_index ) = [ Left_Overall_SNR; mean( Left_Rate_Results );
                             Right_Overall_SNR; mean( Right_Rate_Results ) ]
			 

end

results