load hart_house_soldiers_tower_newer_Finv.txt
Finv=hart_house_soldiers_tower_newer_Finv;
load hart_house_soldiers_tower_newer_c.txt
c=hart_house_soldiers_tower_newer_c;

k=logspace(-3,0,31);  % linear exposure values from -30dB to 0dB
K=log(k);             % log exposure values

%% photoq cement
disp('determining size by loading first frame f1');
f1=loadframe(1);
[M,N]=size(f1);
qsum=zeros(M,N);
Csum=zeros(M,N);
for i=1:31
  disp(sprintf('loading frame %d',i));
  f_i = loadframe(i);
  %% in Matlab may need to convert uint8 to double (not necessary in Octave)
  q=exp(Finv(f_i+1))/k(i); % foreach i, get estimate q=q\hat_i with inverse F
  C=c(f_i+1); % certainty image for running total of certainties
  Csum = Csum + C; % Cumulative sum of Certainty images
  qsum=qsum+q.*C;
end
q=qsum./Csum; % final quantity of light is weighted sum