%%% -march 18th 2002
%%%  found small errors in code which were corrected
%%%  weighting the forward and reverse unrolls gives
%%%  a good fit through the ground truth data. 
%%%  Currently, the forward and reverse unrolls are 
%%%  shifted on the x-axis (one by 4, the other by -2).
%%%  This is unexplained at the moment.
%%%
%%% -march 15th 2002 
%%%
%%% From the march 14th 2002 program, after the correct
%%% base to the logarithm was applied, errors were 
%%% discovered towards the end of the comparunroll.
%%%
%%% This led to the idea of starting the unrolling process
%%% from the comparagram data. The idea goes as follows:
%%% start from a logical point where the comparagram attains
%%% 255, call this f(2^n) for some unknown n to be determined
%%% later. The x-axis co-ordinate where this occurs will be
%%% f(2^n-1). We then find where the value this occurs on the
%%% comparagram. The value at the x-axis co-ordinate at the 
%%% point this occurs is f(2^n-2). We continue this process
%%% until the comparagram attains 0, at this point we have
%%% a value for n (we know how many steps this took).
%%%
%%% This procedure will be more accurate at the right side of
%%% the comparagram instead of the left. If we align the two
%%% results in the centre and then weight the results appropriately,
%%% the combined result should be accurate throughout.  
clear;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                GENERAL FUNCTIONS              % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% slim_comparagram takes a matrix respresenting a comparagram
% and calculates first moments for each row. This result gives
% a vector representing the essence of the comparagram 

function ResultVector = slim_comparagram(Xorig)

   %% Since we don't want to include outliers, the comparagram is
   %% thresholded in a reasonable way. 

   SUMMEDBIN_THRESH = min(max(Xorig'))-1; 

   count = 0;

   % first slim down the comparagram by finding
   % first moments
   Y=zeros(256,256);
   base_str = ".";
   for i=1:256
      if (mod(i-1,16) == 0)
         str=sprintf("slim_comparagram:slimming%s\r",base_str);
         printf(str);
         fflush(1);         
         base_str=strcat(base_str,".");
      end;
      count =0;
      suma=0;
      for j=1:256
        if(Xorig(i,j) > SUMMEDBIN_THRESH)
           count = count + Xorig(i,j);
           suma=suma+(j-1)*Xorig(i,j);
        end;
      end;
      %vec(i)=round(suma/count);
      ResultVector(i)=suma/count;
      %disp([i vec(i) ]);
   end
   str=sprintf("slim_comparagram:slimming%s done\n",base_str);
   printf(str);
   fflush(1);         
endfunction

% load_comparagram data will load comparagrams in the range
% start_data to end_data. It then comparasums the data and
% saves it in the file Xsummed. If no arguments are given,
% load_comparagram simply loads Xsummed. 

function Xsummed = load_comparagram_data(start_data, end_data)

   %%% Load and comparasum the test data if
   %%% necessary. To reload the data load_raw_data
   %%% must be set to 1. 
   if (nargin == 0)
      load_raw_data = 0;
   else
      load_raw_data = 1;
   end

   if (load_raw_data == 1)
      disp("load_comparagram_data: summing given comparagrams\r");
      fflush(1);
      start_ind = start_data;
      end_ind = end_data;
      Xsummed = zeros(256,256);
      for i=start_ind:end_ind
          str=sprintf("c[r%03d.jpg-r%03d.jpg].u32.rgb.raw",i+3,i);
          fin =fopen(str,'r');
          [X count]=fread(fin, 256*256,"unsigned int",0);
           X=reshape(X,256,256);
           X=X';
           fclose(fin);
           Xsummed = Xsummed +X;
      end
      disp("load_comparagram_data: saving comparasummed data\n");
      fflush(1);
      save Xsummed Xsummed;
   else 
      printf("load_comparagram_data: loading already existing data\r");
      fflush(1);
      load Xsummed;
      disp("load_comparagram_data: comparagram loaded                   ");
      fflush(1);
   end
   size(Xsummed);

endfunction

% bisection_search just looks for a function value (given
% spline lookup data cpgram_values). It assumes the function
% value you are looking for is the end_value. Though this
% seems strange, the function's intended use is for reverse
% comparagram estimation.

function ret_val = bisection_search(cpgram_values,start_value,end_value,tolerence)
   target_value = end_value;
   test_point = (start_value+end_value)/2;
   test_point_value = ppval(cpgram_values,test_point);

   % bisection_search debugging info
   %str=sprintf("bisection_search: searching for %7.7g\n",end_value);  
   %disp(str);
   %fflush(1);
  
   while (abs(test_point_value - target_value) > tolerence)
      if (test_point_value > target_value)
         end_value = test_point;
      else
         start_value =test_point;
      end;
      test_point = (start_value+end_value)/2;
      test_point_value = ppval(cpgram_values,test_point);
 
      % bisection_search debugging info
      % disp(test_point);
      % fflush(1);
   end;   
   % bisection_search debugging info
   % str=sprintf("bisection_search: point located at %7.7g\n",test_point);
   % disp(str);
   % fflush(1);
   ret_val=test_point;

endfunction
       


% create_interpolation_fuction takes a vector of values
% and returns the parameters needed for evaluating a spline
% which interpolates each of the input vector's points. 
%
% Note that the function automatically assumes the first 
% value should be 0 (f(0) = 0), and sets this value
% accordingly. The function also makes inferences as
% to what the one-sided derivative should be at the
% origin. 

function ret_values = create_interpolation_function(vec)

   vec(1)=0.0;
   approx_deriv_at_origin = (vec(8)-vec(1))/8;
   ret_values=csape([0:255],vec,'complete',[approx_deriv_at_origin 0]);

endfunction

% clip_function ensures that a vector has no values greater than max_value

function return_vector = clip_function(in_values, max_value)

   return_vector = in_values;
   saturated = 0;
   %%%% Clip function values at max_value 
   for i=1:length(in_values)
      if(return_vector(i) >max_value || saturated)
         return_vector(i) = max_value;
         saturated = 1;
      end;
   end; 

endfunction

% do_forward_unroll takes spline evaluation values (those
% created with the create_interpolation_function) and unroll
% the underlying function.
 
function unroll_data = do_forward_unroll(cpgram_values,granularity,log_base)

   num_of_unrolled_pts = 1; 
   x_val=granularity;
   flat_spot_detect = 0;% detect a flat spot of the comparagram
   count = 0;
   unrolled_values =0;

   % do the unrolling
   while (ppval(cpgram_values,x_val) <255 && flat_spot_detect ~= 1)
      x_val_new=ppval(cpgram_values,x_val);
      count = count +1;
      unrolled_values(count) = x_val_new;
      if (x_val >200 && abs(x_val_new-  x_val) <.001)
         flat_spot_detect = 1;
         disp('the comparagram contains a saturation point or a flat region');
         disp('the maximum value attained was');
         disp(x_val);
      else
         x_val=x_val_new;
         num_of_unrolled_pts = num_of_unrolled_pts+1;
      end
   end

   % this method will attempt to build co-orinate pairs
   x(1) = 0; y(1) = 0; % want the spline to pass through 0

   for i =2:length(unrolled_values)+1
      x(i) = 2^(i-2);
      y(i) = unrolled_values(i-1);
   end;
   x=x.* (1/x(i)); % Normalize X to be between 0 and 1
 
   x=log(x);            %  match to given log base 
   x=1/log(log_base).*x;
   unroll_data =[x y]; 
endfunction 

% unroll_data unrolls the comparagram starting from the 
% first point close to 255

function unroll_data = do_reverse_unroll(cpgram_values,granularity,log_base)
   % find a reasonable point to start the reversed unroll    
   starting_point = 232.8; % remember to generalize this later
   increment = .001;
   end_point = 255;

   reasonable_all_white_value = 254.5;
   %reasonable_all_white_value = 254.99;
   reasonable_all_black_value = 000.01;

   current_point = starting_point;
   found_point = 0;
 
   disp("do_reverse_unroll: Looking for saturation point");
   fflush(1);
   while (found_point == 0 && current_point < end_point) 
      sat_val = ppval(cpgram_values,current_point);
      if (sat_val >= reasonable_all_white_value)  
         found_point = current_point;
      else
         current_point = current_point + increment;
      end;
   end;
   str= sprintf("do_reverse_unroll: Found saturation point at %7.7g",found_point);
   disp(str); 
 
   % start unrolling, found point is f(2^n)=255
   % this means f(2^n-1) = found_point
   vec(1) = found_point; 
   returned_point = found_point;
   count = 1;
   tolerence = .001;
   %returned_point
   %found_point
   %fflush(1);
  
   while (returned_point > reasonable_all_black_value) 
       returned_point = bisection_search(cpgram_values,0,returned_point,tolerence);
       count = count + 1;
       vec(count) = returned_point;
   end
   
   % unrolling is done, now make sense of the points
   unroll_data(1,1)=0;
   unroll_data(1,2)=0;
   for i=1:length(vec)
      unroll_data(i+1,1) = 2^(i-1); 
      unroll_data(i+1,2) = vec(length(vec)-i+1);
   end
   % add in the last values
    unroll_data(i+2,1) = 2^i;
    unroll_data(i+2,2) = 255; 
   %   unroll_data(i+3,1) = 2^(i+1);
   %   unroll_data(i+3,2) = 255; 
      
   % now just normalize the data and put it into the correct scale 
   dimension = size(unroll_data);
   len = dimension(1,1);
   largest_value =unroll_data(len,1);

   for i=1:len
      unroll_data(i,1)=unroll_data(i,1)/largest_value;
   end

   unroll_data(:,1)=log(unroll_data(:,1));            %  match to given log base 
   unroll_data(:,1)=1/log(log_base).*unroll_data(:,1);

endfunction
 
%%%%%%%%%%%%%%% END GENERAL FUNCTIONS %%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%% MAIN PROGRAM %%%%%%%%%%%%%%%%%%%%%%%%%%%%
Xsummed = load_comparagram_data();
vec = slim_comparagram(Xsummed);
cpgram_values = create_interpolation_function(vec);

% evaluate the created function
y_axis = ppval(cpgram_values,[0:.001:255]);

%%% plot the comparagram
%plot([0:.001:255],y_axis);

%%%% Now deal with the forward and reverse comparagram unrolling
ground_truth_log_base = 2^(1/3);

% do reverse unroll
reverse_unrolled_values = do_reverse_unroll(cpgram_values,.1,ground_truth_log_base);
x1=reverse_unrolled_values(:,1);
x1=x1-2;
y1=reverse_unrolled_values(:,2);
spline_vals1=csape(x1(2:length(x1)), y1(2:length(y1)), 'complete',[0,1]);
y_cam_func1 = ppval(spline_vals1, -40:.1:0);
y_cam_func1 = clip_function(y_cam_func1,255);
len = length(x1);

% do forward unroll
forward_unrolled_values = do_forward_unroll(cpgram_values, .1, ground_truth_log_base);
x=forward_unrolled_values(:,1);
x=x+4;
y=forward_unrolled_values(:,2);
spline_vals=csape(x(2:length(x)), y(2:length(y)), 'complete',[0,0]);
y_cam_func = ppval(spline_vals,-40:.1:0);
y_cam_func = clip_function(y_cam_func,255);

%%%% this particular ground truth data is in log base 2^(1/3)
%%%% load ground truth data 
load -f ground_truth_d1_response_function_data_consolidated.txt
D=ground_truth_d1_response_function_data_consolidated;
[c i]= sort(D(:,1));
E=D(i,:);

%%%% these function need to be adjusted per comparasum for now
x_run =[-40:.1:0];
x1_run =[-40:.1:0];
x2_run =[-40:.1:0];

%%% create the composite unroll
count = 1;
for i=-40:.1:0
   weight1= (i+40)*1/40;
   weight2= 1-weight1;
   y_cam_func_composite(count) = weight2*ppval(spline_vals,i) + weight1*ppval(spline_vals1,i);
   count = count+1;
end    
y_cam_func_composite = clip_function(y_cam_func_composite,255);

%%%%%%%%%%%%%%% Plot the results %%%%%%%%%%%%%%%%%
grid ("on") ;
gset key -13,245;
gset xrange [-37:0];
gset yrange [0:255];
title('Spline Approximations Using Unrolled Data against Ground Truth Data');
xlabel("q");
ylabel("f(q)");

%%%% uncomment these for postscript output
gset terminal postscript 19;
gset output "unroll.ps";
gset pointsize .75;
plot(x(2:length(x)),
     y(2:length(y)),
     '1*;forward unrolled points;',

     x_run,
     y_cam_func,
     '1.;forward unroll interpolated values;',

     x1(2:length(x1)),
     y1(2:length(y1)),
     '1x;reverse unroll interpolated points;',

     x1_run,
     y_cam_func1,
     '2-;reverse unroll interpolated values;',
  
     x2_run,
     y_cam_func_composite,
     '1-;composite of forward and reverse unrolls;',

     E(:,1),
     E(:,2),
     'o;Ground Truth Data;'
);

