% Follow Cheezum et al. (2001) procedure for simulated particle images
% Generate 1000 frames
% (Script written by SS Rogers (2007))

% first make high res "precise" image - 9nm pixel size
% add circular particle of defined radius
% convolve with PSF
highres_1=zeros(1001);
highres_2=zeros(1001);
%one pix particle - "point source"
highres_1(521,521)=10000;
highres_2(501,501)=10000;

%convolve with point spread function (PSF)
x=(-200:200)*9e-9;
o=ones(size(x));
radius=sqrt((x'*o).^2+(o'*x).^2);
a=2*pi*1.45/520e-9;
PSF=(2*besselj(1,radius*a)./radius).^2;
PSF(201,201)=a^2;
highrestemp_1=conv2(highres_1,PSF);
highrestemp_2=conv2(highres_2,PSF);
highres2_1=sparse([],[],[],3401,3401);
highres2_2=sparse([],[],[],3401,3401);
highres2_1(1001:2401,1001:2401)=highrestemp_1;
highres2_2(1001:2401,1001:2401)=highrestemp_2;
imagemask=1201:2201;

% Create 1000 low resolution frames:
% - shift data by whole number of high res pixels
% - integrate over CCD cells
% - scale signal and background
% - add Poisson noise
% - write image file

%generate a particle track - in whole displacements of the fine pixels
%initial particle position is image centre
shiftx_1=zeros(1,1000);
shifty_1=zeros(1,1000);
shiftx_2=zeros(1,1000);
shifty_2=zeros(1,1000);
    % - use Gaussian distribution for random numbers, mean=0, stdev=3
for t=2:1000    
    shiftx_1(t)=shiftx_1(t-1)+round(random('Normal',0,6,1));
    shifty_1(t)=shifty_1(t-1)+round(random('Normal',0,6,1));
    shiftx_2(t)=shiftx_2(t-1)+round(random('Normal',0,6,1));
    shifty_2(t)=shifty_2(t-1)+round(random('Normal',0,6,1));
end

figure
hold all
plot(shiftx_1,shifty_1,'r')
plot(shiftx_2,shifty_2,'b')
drawnow
hold off

%scale to given signal (and background=10)
SignaltoBackground=30;
Signal=10*SignaltoBackground;
SignaltoNoise=(Signal-10)/sqrt(Signal);
disp(['SignaltoNoise = ' num2str(SignaltoNoise)])

%generate images
mkdir sim_images
for t=1:1000
    highres3_1=highres2_1(imagemask-shiftx_1(t),imagemask-shifty_1(t));
    highres3_2=highres2_2(imagemask-shiftx_2(t),imagemask-shifty_2(t));
    %integrate over blocks (CCD 'cells') of 11x11 pixels - such that unshifted image has peak in the centre of a block
    lowres_1=zeros(91);
    lowres_2=zeros(91);
    for n=1:91
        for m=1:91
            lowres_1(n,m)=mean2(highres3_1((n-1)*11+(1:11),(m-1)*11+(1:11)));
            lowres_2(n,m)=mean2(highres3_2((n-1)*11+(1:11),(m-1)*11+(1:11)));
        end
    end
    
    lowres2_1=lowres_1 / max(max(lowres_1)) * (Signal-10) + 10;
    lowres2_2=lowres_2 / max(max(lowres_2)) * (Signal-10) + 10;
    lowres2 = lowres2_1 + lowres2_2;
    
    %add poisson noise
    %use Image Toolbox function IMNOISE - note: input values are interpreted scaled by 1e12
    lowres3 = imnoise(lowres2/1e12,'poisson')*1e12;

    %write image
    filename=['sim_images/' num2str(t,'%0.4d') '.tif'];
    imwrite(uint16(lowres3),filename,'tif');
    disp(t)
end

%rescale x and y trajectory for CCD cell grid
%exacttrajectoryx=shifty/11+46;
%exacttrajectoryy=shiftx/11+46;