
clear all, close all;

% Display Stuff
mygray = repmat(linspace(0, 1, 256)',1,3);
fps    = 5;
Nplay  = 1;


% 1) Define Sensor Characteristics
R = 100; C = R; N = 16; 
std_v = 60; % Read noise (in e-)

% 2) Define Exposure Characteristics 
dExpTime     = 10e-3;
ph_10ms_pix  = 30; % photons per 10ms of exposure per pixel


% 3) Definition of the scene
    scene1 = double(imread('lena.tif','tif')); %scene = double(imread('testpat1.tif','tif'));
    imshow(scene,[]),axis image

    % Just a different scene
    AA     = ones(2,2)*255; BB = zeros(2,2); CC = [AA,BB;BB,AA];
    scene2 = repmat(C,128,128);

    % More different scenes 
    scene3 = makePatternREC; %makePatternSMOOTH;
    figure, colormap(mygray), image(255*scene3/max(pat(:))),axis image;

scene = scene1; 
pos0  = [210 200];        
off = 16;            
mask  = zeros(size(scene)); mask( pos0(1)+off : pos0(1)+R-off-1, pos0(2)+off : pos0(2)+C-off-1 ) = 1;
scene = scene.*mask;
figure, imshow(scene,[]),axis image;
originalImage = scene(pos0(1):pos0(1)+R-1, pos0(2):pos0(2)+C-1);
figure, imagesc(originalImage), colormap(mygray), axis image;


% 4) Define Motion
TotalShift = 4;
positions  = createPositionsLin(N, TotalShift);  
            %createPositionsCirc(TimeSteps, TotalAngle, Radius)

% END OF DEFINITIONS
% BEGIN OF MAIN ROUTINES

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the Charge Distribution in the sensor's spacetime, with the following parameters
% Sensor parameters (as described above)
% Exposure parameters
% The scene
% Motion parameters
noiseflag    = 1; adcflag = 1; % If 1, adds noise/does quantization
[Q,convgain] = createChargeDistribution(R,C,N, scene, positions, pos0, dExpTime, noiseflag, adcflag, std_v, ph_10ms_pix);
Qmax         = max(Q(:));
% Play as movie
figure, colormap(mygray)
for tp = 1:N+1
    image(255*Q(:,:,tp)/Qmax),axis image;
    F(tp) = getframe;
end
movie(F,Nplay,fps) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the Differential Images
% from the discrete Charge Distribution
dQ = calculateDifferentialImages(Q);
dQmax = max(dQ(:));
% Play as movie
figure, colormap(mygray)
for tp = 1:N
    imagesc(dQ(:,:,tp)),axis image;
    F2(tp) = getframe;
end
movie(F2,Nplay,fps) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the Motion Vector Field
% using 
% Differential Images
% Noise level lamda used to reject false motion vectors because of noise
% Half block size hbs and search range sr
hbs = 1; sr = 1;
lambda = 2*sqrt(2)*std_v/(2*hbs+1); % Noise suppression factor for motion vector estimation
%MVF = calculateMVFblock(dQ, bs);
%MVF = calculateMVFglob(dQ, bs);
MVF = calculateMVFpix(dQ, hbs, sr, lambda);
% Play as movie
figure
for tp = 1:N-1 
    quiver(MVF(:,:,tp,2), MVF(:,:,tp,1)),  axis ij, axis equal,
    F3(tp) = getframe;
end
movie(F3,Nplay,fps) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the datastructure 'motion' from the MVF
% This datastructure contains for each pixel in frame 1 a history of 
% where it is later at which point in time.
motion = convert(MVF);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the datastructures 'modQ' and 'num' from the charge distribution Q 
% and the array num. 
% 'Qmod' contains a new charge distribution that contains 
% for each pixel the amounts of charge that it created during 
% its path over different pixels of the sensor surface
% The array 'num' keeps track of the numer of pixels per timestep over which each point 
% of the scene has moved on its trajectory. This information is
% needed for correct statistical evaluation in the routine 'blurFreeImConst'.
[modQ,num] = qModifier(Q,motion);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct the blur free image with higher SNR and dynamic range 
% from the information in modQ, and num that was gathered over the 
% spacetime trajectory
im = blurFreeImConst(modQ, num, dExpTime, std_v, convgain);
figure,imagesc(im),colormap(mygray),axis image;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do the same but use a perfect MVF (calculated from the original motion) 
% as opposed to the MVF estimted from the image data
rMVF = realMVF(off, off, positions, [R,C]);
rmotion = convert(rMVF);
[rmodQ, rnum] = qModifier(Q, rmotion);
rim = blurFreeImConst(rmodQ,rnum,dExpTime,std_v,convgain);
figure,imagesc(rim),colormap(mygray),axis image;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate some figures of merit
stdOneCap = std2(Q(:,:,2)-originalImage)
stdBlurred = std2(Q(:,:,16)-originalImage)
stdIm = std2(im - originalImage)
stdrIm = std2(rim - originalImage)







