K=1; % in this case, k=2 so K=1
% 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

%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);
% note that 1../((1+exp(-a*(Q+K0))).^c) = k0q^a/(k0q^a+1)
%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(Q,256*fk); hold off
fkk=1../((1+exp(-(a*(Q+2*K+K0)+b))).^c); hold on; plot(Q,256*fkk); hold off
fk3=1../((1+exp(-(a*(Q+3*K+K0)+b))).^c); hold on; plot(Q,256*fk3); hold off
fk4=1../((1+exp(-(a*(Q+4*K+K0)+b))).^c); hold on; plot(Q,256*fk4); hold off
fk5=1../((1+exp(-(a*(Q+5*K+K0)+b))).^c); hold on; plot(Q,256*fk5); hold off
fk6=1../((1+exp(-(a*(Q+6*K+K0)+b))).^c); hold on; plot(Q,256*fk6); hold off
fk7=1../((1+exp(-(a*(Q+7*K+K0)+b))).^c); hold on; plot(Q,256*fk7); hold off
fk8=1../((1+exp(-(a*(Q+8*K+K0)+b))).^c); hold on; plot(Q,256*fk8); hold off
fk9=1../((1+exp(-(a*(Q+9*K+K0)+b))).^c); hold on; plot(Q,256*fk9); hold off
fk10=1../((1+exp(-(a*(Q+10*K+K0)+b))).^c); hold on; plot(Q,256*fk10); hold off
fk11=1../((1+exp(-(a*(Q+11*K+K0)+b))).^c); hold on; plot(Q,256*fk11); hold off
xlabel('Quantity of light, Q = log(q)')
ylabel('Pixel integer') % 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); plot(Q,[c1(1) c1]);
c2=diff(256*fk); hold on; plot(Q,[c2(1) c2]); hold off
c3=diff(256*fkk); hold on; plot(Q,[c3(1) c3]); hold off
c4=diff(256*fk3); hold on; plot(Q,[c4(1) c4]); hold off
c5=diff(256*fk4); hold on; plot(Q,[c5(1) c5]); hold off
c6=diff(256*fk5); hold on; plot(Q,[c6(1) c6]); hold off
c7=diff(256*fk6); hold on; plot(Q,[c7(1) c7]); hold off
c8=diff(256*fk7); hold on; plot(Q,[c8(1) c8]); hold off
c9=diff(256*fk8); hold on; plot(Q,[c9(1) c9]); hold off
c10=diff(256*fk9); hold on; plot(Q,[c10(1) c10]); hold off
c11=diff(256*fk10); hold on; plot(Q,[c11(1) c11]); hold off
xlabel('Quantity of light, Q = log(q)')
ylabel('Certainty') % since F^-1(Q)=1/(1+e^-Q)
title('Certainty functions (derivatives of response curves)')
grid

