% mannfamily dataset with range-coord transformed sigmoid
%
% forget the cubic-coordinate transformed sigmoid!!!
%
% richardson function!!! http://envhort.ucdavis.edu/jhl/phasic/3phas.htm
% used in Richardson et al. (1988)
% example use: testrichardson([],5,-5,.2)    ...testrichardson(J,a,b,c)
% example use: testrichardson(J,5,b,.1)
%
% c always less than 1 for most cameras
%
% b does not affect the range range plot
% let a go to a*big and c go to c/big to find slope first, then choose a,c:
% testrichardson(J,a*10,0,c/10,2)
%
% mann family: a = 3.8783; c = 0.1423
%
% Use: J=crosshist(v0,v1)
%      fig_mannfamily_resp

K=2
% original values for noisy (recovered from www page) coon image pair
%a=3.8783
%b=0
%c=0.1423
a=3.8783
b=0
c=0.1563

disp('for now, setting K0 to have 0 at maximal Q value')
K0=1.12

%Q=linspace(-10,10,1000); % Quantity (log quantity) of light on image sensor
Q=linspace(-10,0,1000); % Quantity (log quantity) of light on image sensor

figure(1)

clg
h1=axes('position',[.1 .5 .8 .2]);
%subplot(413)
%ha=axes; set(ha,'fontsize',10)

%the following is NOT the correct form, since should be log log this
%f=1../((1+exp(-(a*(Q)+b))).^c);
 f=1../((1+exp(-(a*(Q+K0)+b))).^c);
%plot(f) % first one is the darker image which is like the slower film
%plot(Q,f,':') % first one is the darker image which is like the slower film
plot(Q,256*f)
%title(sprintf(''))
%fk=1../((1+exp(-(a*(Q+K)+b))).^c);
fk=1../((1+exp(-(a*(Q+K+K0)+b))).^c);
hold on
%plot(fk)
%plot(Q,fk,'--')
plot(Q,256*fk)
hold off
xlabel('Quantity of light, Q')
ylabel('Pixel value') % since F^-1(Q)=1/(1+e^-Q)
title('Response curves recovered from differently exposed images')
axis([-10,0,0,256])
grid

%subplot(414)
h2=axes('position',[.1 .1 .8 .2]);
c1=diff(256*f); c2=diff(256*fk);
plot(Q,[c1(1) c1]); hold on; plot(Q,[c2(1) c2]); hold off
xlabel('Quantity of light, Q')
ylabel('Certainty') % since F^-1(Q)=1/(1+e^-Q)
title('Certainty functions (derivatives of response curves)')
grid

figure
clg
%subplot(222)
%fk=1../((1+exp(-(a*(Q+K)+b))).^c);
%fk=1../((1+exp(-(a*(Q+K+K0)+b))).^c);
plot([0 f],[0 fk])
xlabel('darker image   f(q)')
ylabel('lighter image   g(f) = f(kq)')
title('COMPARAMETRIC PLOT')
axis('square')

%figure;
%clg
%subplot(221)
%tvslog(J.',.0001,'neg')
%axis('xy')
%hold on
%fk=1../((1+exp(-(a*(Q+K)+b))).^c);
%fk=1../((1+exp(-(a*(Q+K+K0)+b))).^c);
%plot(255*[0 f],255*[0 fk])
%xlabel('pixel value, image 1')
%ylabel('pixel value, image 2')
%hold off
%grid

%subplot(223)
%no not right... g=f.^(1/c)./(f.^(1/c)+exp(-a*K));
%x=linspace(0,1,1000)
%still not right... g=x.^(1/c)./(x.^(1/c)+exp(-a*K));
%plot(g)

