%test acs fit against ground truth data
%
%this script was used to produce
%ground_truth_d1_response_function_with_comparasolve_result.eps
%
%acs fit data determined as follows:
%
%comparasolve on Jsum03.txt (3/3 stops = 1 stop difference, e.g. k=2)
%gives a very good fit resulting in:
%[a c k s]=[2.14107947021 0.34932613936 3.17980174559 0.995241982284]
%  2.14108
%  0.34933
%  3.17980
%  0.99524
%
% with only acs (k fixed at 2)
%[a c k s]=[1.50482368858 0.653010561063            2 1.09382809917]
%  1.50482     log(1.50482*1.50482) ans = 0.81735
%  0.65301      0.65301*0.65301     ans = 0.42642
%  1.09383     1.09383*1.09383      ans = 1.1965
%ans =  4.7693e+05
%done estimating a,c,s
%showit log(a*a), (c*c),(k*k),(s*s)
%
%therefore may want to consider having k=sqrt(2) so k*k will be 2:
%fit seems poor
%bfgs_hess: bailing out of outer loop, max iterations exceeded
%  1.93823
%  0.73812
%  1.14176
%
%[a c k s]=[-7.56184538786 -0.380545718411 1.41421356237            1]
%x =
%  -7.56185    log(a*a) ans = 4.0462
%  -0.38055        c*c  ans = 0.1448
%
% therefore a=4.0462; c=0.1448;

function whatever(a,c,k0)

if nargin == 0
  %[a c k s]=[-7.56184538786 -0.380545718411 1.41421356237            1]
  %x =
  %  -7.56185    log(a*a) ans = 4.0462
  %  -0.38055        c*c  ans = 0.1448
  % therefore:
  a=4.0462; c=0.1448;
  k0=140; % arbitrary unknown shift left-right; 140 is approx.
end%if

load ground_truth_d1_response_function_data_consolidated.txt
%in Octave, the variables are same as filenames;
%in Matlab (a weaker Octave-like equivalent inferior to Octave in various ways)
%filenames are truncated.  Uncomment the appropriate one:
Q=ground_truth_d1_response_functi(:,1);
f=ground_truth_d1_response_functi(:,2);
%Q=ground_truth_d1_response_function_data_consolidated(:,1);
%f=ground_truth_d1_response_function_data_consolidated(:,2);

% sort so we can plot as a continuous function:
[Qsort, ind]=sort(Q);
Q=Qsort;
f=f(ind);

q=(cbrt(2)).^Q;

%f_fit = (q.^a./(q.^a+1)).^c; % ac model
f_fit = (k0*q.^a./(k0*q.^a+1)).^c;
f_fit = 255*f_fit;
s=1; f_fit = min(255,f_fit*s); % no saturation for now

plot(Q,f,'+',Q,f_fit);
axis([min(Q),max(Q),0,255])
xlabel('Q (Quantity of light)','FontSize',20) % Q=log_cbrt{2}(q)
ylabel('f (response byte)','FontSize',20)
title('estimated camera response function with ground truth data','FontSize',14)
grid

figure
% now to plot the certainty, we need uniformly spaced Q rather than random Q:
Q=linspace(min(Q),max(Q),1000); %or try Q=linspace(-50,5,1000);
q=(cbrt(2)).^Q;
f_fit = (k0*q.^a./(k0*q.^a+1)).^c;
f_fit = 255*f_fit;
s=1; f_fit = min(255,f_fit*s); % no saturation for now
%plot(f_fit)
certainty=diff(f_fit);
certainty=certainty/max(certainty);
Certainty=log(certainty);
Certainty=Certainty-min(Certainty);
Certainty=Certainty/max(Certainty);
x=linspace(0,255,999);
%plot(x,certainty,x,Certainty,'-');axis([0 255 0 1.05])
plot(x,certainty,x,Certainty,'--');axis([0 255 0 1.05]) % logscale dashed line
ylabel('relative certainty','FontSize',20) % Q=log_cbrt{2}(q)
xlabel('f (response byte)','FontSize',20)
title('estimated camera certainty function','FontSize',14)
grid

% we could also try making f be linspace (0,255,256):
f=linspace(0,255,256);
q=(k0*f.^(1/c).*(1-f.^(1/c))).^(1/a);
Q=log(q);

