%% Make a data set
N=4             %Number of particles
sz=128          %Size of images
Photons=500;    %Photons per emitter per frame
sigma=1.3;      %PSF size (pixels)
bg=1;           %background photons per pixel
T=20;           %number of time frames
D=0.1;          %(micron^2/s)
timestep=.05;   %TimeStep (s)
pixelsize=.1;   %pixel size (micron)
Dim=2;          %Number of Dimensions

%This creates a trajectory
displacements=sqrt(2*D*timestep)*randn(T,Dim,N);
x=sz*repmat(rand(1,N),[T 1 1])+squeeze(cumsum(displacements(:,1,:)));
y=sz*repmat(rand(1,N),[T 1 1])+squeeze(cumsum(displacements(:,2,:)));

data=newim([sz sz T]);
for tt=1:T %a simple 'for' loop
    tmp=newim([sz sz]);
    for nn=1:N
        tmp=gaussianblob(tmp,[x(tt,nn) y(tt,nn)],sigma,Photons)+bg;
    end
    data(:,:,tt-1)=tmp;
    tt
end
data=noise(data,'poisson')



%% One method of finding centers using local maxima

Thresh = 10;

%finding local max
unisz=round(sigma*2+1);  %uniform filter size
im_unif=unif(data,[unisz unisz 0],'rectangular')-unif(data,[2*unisz 2*unisz 0],'rectangular');
maxsz=2*round(sigma*2+1)-1;
im_max=(im_unif>=.999*maxf(im_unif,[maxsz maxsz 0],'rectangular'))&(im_unif>Thresh);

a=find(permute(im_max,[2 1 3]));  %pixel number of all positive pixels
z=floor(a/sz/sz);               %calculating z value from pixel number
pnum=mod(a,sz*sz);              %pixel number in each frame
x=mod(pnum,sz);                 %calculating x value from pixel number
y=floor(pnum/sz);               %calculating y value from pixel number

%combine into array
allcds=[x y z]













