% A didactic script for using the NESIF Matlab function
% with synthetic data.
% See the paper "Measuring the noise of digital imaging sensors 
% by stacking raw images affected by vibrations and illumination flickering" 
% SIAM Journal on Imaging Sciences
% vol. 8, no. 1, p. 611-643
% Frederic Sur, Michel Grediac
% 2015

clear all
close all


% use either "synthetic" or "boat" experiment
% (uncomment the corresponding lines)

% % boat
% Image=double(imread('boat.tif'));
% h = fspecial('gaussian',[3 3],.5); 
% Image=imfilter(Image,h);
% Image=Image/max(Image(:));

% synthetic 
Image=zeros(512,512);
for i=1:512,
   Image(:,i)=1+floor(8*i);
end;
for i=0:4,
   Image((60+90*i):(110+90*i),380:470)=288*i+1;
end;
for i=0:4,
   Image((60+90*i):(110+90*i),230:320)=288*(i+5)+1;
end;
for i=0:4,
   Image((60+90*i):(110+90*i),80:170)=288*(i+10)+1;
end;
h = fspecial('gaussian',[3 3],.5); 
Image=imfilter(Image,h);
Image=Image/2^12;


% Data for Section 4.2
% definition of camera and external parameters

cameraparameters=1; % 1, 2, or 3
extparameters=2; % 1, 2, or 3

switch cameraparameters,
   case 1,
      gain=5;
      meand=200.0;
      stdd=15;
      amplitude_electron=600;
   case 2,
      gain=1;
      meand=100.0;
      stdd=8;
      amplitude_electron=3000;
   case 3,
      gain=0.3;
      meand=40;
      stdd=3;
      amplitude_electron=10000;   
end;

switch extparameters,
    case 1,
       stda=0.3;
       stdb=0.2;
       stdg=0.1;
    case 2,
       stda=0.1;
       stdb=0.05;
       stdg=0.01;
    case 3,
       stda=0.01;
       stdb=0.02;
       stdg=0.05;
end;

nr_im=100;  % number of images in the stack
taille_xx=size(Image,1);
taille_yy=size(Image,2);
[I,J]=meshgrid(1:taille_xx,1:taille_yy);
rng('shuffle')
gamma_tab=1+stdg*randn(1,nr_im);  % illumination flickering

% generate the stack of 12-bit images
disp('generating the synthetic dataset...')
drawnow('update')

stack=zeros(taille_xx,taille_yy,nr_im);
for i=1:nr_im
    alpha=stda*randn;
    beta=stdb*randn;
    preimage=amplitude_electron*gamma_tab(i)*interp2(I,J,Image,I+alpha,J+beta);
    stack(:,:,i)=uint16(gain*poissrnd(preimage)+meand+stdd*randn(size(preimage)));
end;
stack=stack(10:end-10,10:end-10,:);
stack(stack>2^12)=2^12;
disp('done')
disp('')
drawnow('update')

output=NESIF(stack,0,12);

disp('estimating camera parameters   (and ground truth GT)')
disp(['gain: ',num2str(output(1)),'   ( GT: ',num2str(gain),' )']) 
disp(['offset: ',num2str(output(2)),'   ( GT: ',num2str(meand),' )']) 
disp(['readout / quantization noise: ',num2str(output(3)),'   ( GT: ',num2str(stdd),' )']) 
disp(['sigma^2-g.mu: ',num2str(output(4)),'   ( GT: ',num2str(stdd^2-gain*meand),' )']) 
disp(' ')
disp('estimating external parameters  (and ground truth)')
disp(['flickering st.d.: ',num2str(output(5)),'   ( GT: ',num2str(stdg),' )']) 
disp(['x-axis vibration st. d.: ',num2str(output(6)),'   ( GT: ',num2str(stda),' )']) 
disp(['y-axis vibration st. d.: ',num2str(output(7)),'   ( GT: ',num2str(stdb),' )']) 
disp(['x-y vibration covariance: ',num2str(output(8)),'   ( GT: ',num2str(0),' )']) 
disp(' ')
