clsm_img_gen.m

function [Img1,Img2]=clsm_img_gen(U,V,W)
%clsm_generator function generates synthetic confocal laser scanning microscopy images (1D/2D/3D)
%Developed/updated by Brian Jun (2016) based on synthetic PIV image gnerator version previously updated by Sam Raben (2011)  
%
%
%Inputs
%Assign uniform velocity U,V,W (pixels/frame)
%%%%%%%%%%%%%%%%%%%%%%%%%%
% physical parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%
    d     = 2.2; %particle diameter (pix)
    Lx     = 32; %image size width x (pix)
    Ly     = 32;  %image size height y (pix)   
    Lz     = 32; %image depth z (pix)
    C     = 0.001;        %density (particles/pix^3) 
    Vout  = 24;          %distance outside image to generate particles (pix)
    distr = 'uniform';    % normal or uniform

    pix_res = 0.31; % pixel resolution (um)
    D = 5.01*10^-6; %Diffusion coefficient of particles (um^2/us)
    T_s = 0.01; %pixel dwell time of the scanner (us)
%%%%%%%%%%%%%%%%%%%%%%%%%%

    if numel(d) ~= 1
        d_ave = mean(d);
        N=round(C*(Ly+2*d_ave)*(Lx+2*d_ave+2*Vout)*(Lz+2*d_ave));
        switch lower(distr)
            case 'normal'
                d = d_ave + 3.*randn(N,1);
            case 'uniform'
                d = rand(N,1).*(d(2)-d(1)) + d(1);
        end
    elseif nargin < 11
        N=round(C*(Ly+2*d)*(Lx+2*d+2*Vout)*(Lz+2*d));
        d = repmat(d,N,1);
    else
        N = numel(X);
        d = repmat(d,N,1);
    end

    if nargin < 11
        X = (Lx-16)*rand(N,1)+8; 
        Y = (Ly-16)*rand(N,1)+8;
        Z = (Lz-16)*rand(N,1)+8;
    end

        N = numel(X);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Generate Image 1                                  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    D_p = D*T_s/(pix_res^2); %Diffusity per pixel scan (pix^2)
    sigma = sqrt(2*D_p); 
    
    Img1 = zeros(Lx+1,Ly+1,Lz+1);
    I1 = zeros(Lx+1,Ly+1,Lz+1);
    
    d_one = d(1);
    
    for t_3 = 1:1:Lz
        for t_2 = 1:1:Ly
            for t_1 = 1:1:Lx
                    x0 = U/(Lx*Ly*Lz) + sigma*randn(size(X,1),size(X,2)); 
                    y0 = V/(Lx*Ly*Lz) + sigma*randn(size(Y,1),size(Y,2)); 
                    z0 = W/(Lx*Ly*Lz) + sigma*randn(size(Z,1),size(Z,2));

                    X = X + x0; 
                    Y = Y + y0; 
                    Z = Z + z0;

                for p=1:N
                    if X(p)<=t_1+4 && X(p)>=t_1-4 
                        for i=floor(X(p)-0.75*d_one):ceil(X(p)+0.75*d_one)
                            for j=floor(Y(p)-0.75*d_one):ceil(Y(p)+0.75*d_one)
                                for k=floor(Z(p)-0.75*d_one):ceil(Z(p)+0.75*d_one)

                                    if i>=0 && i<=Lx && j>=0 && j<=Ly && k>=0 && k<=Lz  && sqrt((i-X(p))^2+(j-Y(p))^2+(k-Z(p))^2)<0.75*d_one 
                                        I1(i+1,j+1,k+1) = I1(i+1,j+1,k+1) + 0.4*exp(-8*((i-X(p))^2+(j-Y(p))^2+(k-Z(p))^2)/d_one ^2); 
                                    end
                                end
                            end
                        end
                    end
                end
                
                Img1(t_1, t_2, t_3) = I1(t_1, t_2, t_3); 
                I1 = zeros(Lx+1,Ly+1,Lz+1);
                
            end
        end
    end

        Iimg = 1;
        Img1 = permute(Img1,[2 1 3])/Iimg;
        
        
%%%%%%%%%%GENERATE THE SECOND IMAGE %%%%%%%%%%%%%%%%%

    Img2 = zeros(Lx+1,Ly+1,Lz+1);
    I1 = zeros(Lx+1,Ly+1,Lz+1);
    parlist = zeros(N,1);
    
    for t_3 = 1:1:Lz
        for t_2 = 1:1:Ly
            for t_1 = 1:1:Lx

                    x0 = U/(Lx*Ly*Lz) + sigma*randn(size(X,1),size(X,2)); 
                    y0 = V/(Lx*Ly*Lz) + sigma*randn(size(Y,1),size(Y,2));
                    z0 = W/(Lx*Ly*Lz) + sigma*randn(size(Z,1),size(Z,2));

                    X = X + x0; 
                    Y = Y + y0; 
                    Z = Z + z0;
                
                for p=1:N
                    if X(p)<=t_1+4 && X(p)>=t_1-4 
                        for i=floor(X(p)-0.75*d_one):ceil(X(p)+0.75*d_one)
                            for j=floor(Y(p)-0.75*d_one):ceil(Y(p)+0.75*d_one)
                                for k=floor(Z(p)-0.75*d_one):ceil(Z(p)+0.75*d_one)
                                    if i>=0 && i<=Lx && j>=0 && j<=Ly && k>=0 && k<=Lz  && sqrt((i-X(p))^2+(j-Y(p))^2+(k-Z(p))^2)<0.75*d_one
                                        parlist(p) = p;
                                        I1(i+1,j+1,k+1) = I1(i+1,j+1,k+1) + 0.4*exp(-8*((i-X(p))^2+(j-Y(p))^2+(k-Z(p))^2)/d_one ^2);
                                    end
                                end
                            end
                        end
                    end
                end 
                Img2(t_1, t_2, t_3) = I1(t_1, t_2, t_3); 
                I1 = zeros(Lx+1,Ly+1,Lz+1);
            end
        end
    end

        Iimg = 1;
        Img2 = permute(Img2,[2 1 3])/Iimg;