% __________________________________________________________ % %\ % welcome %%\ % %%%\ %%%%%%%%%%%%%%%%%%%%\______________________________________________________ clear all clear all clc %%% input parameters %%\ N=500; px_size= 0.0506; %%%%%%%%%%%%%%% number of frames f=0.2041; %%%%%%%%%%% time lapse (s/frame) cutoff=8; %% please insert only 2^n values TauLimit=int64(N/10)+1; %%%%%%%%%%%%%%% maximum time-lag Filter=0; %% bkg subtraction (0=no; 1=yes) av_toll=10000; %% threshold for bkg subtraction bit=16; % image-type (bit=8, bit=16 ...) %_________________________________________________________________________% % import image-stack \%%%%%%%%%%%%%%%%%% %%% string before indexing \%%%%%%%%%%%%%%%%% name_start='PhogrinGFP_preCHOL.lif_Series001_t'; %%% string after indexing \%%%%%%%%%%%%%%%% name_end='_ch00'; %___________________________________________________________ _____________ %\ | %\ START | %\___________| disp(' importing image-stack') disp(' ') for k=1:N if k<11 inputFileName = sprintf([name_start '00%d' name_end '.tif'],k-1); else if k<101 inputFileName = sprintf([name_start '0%d' name_end '.tif'],k-1); else if k<1001 inputFileName = sprintf([name_start '%d' name_end '.tif'],k-1); else end end end A(:,:,k) = imread(inputFileName); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% starting fitting parameters %%%%%%%%%%% %linear stics [scale D offset] %start2=[1 0.005 1]; %parabolic stics [scale D v offset] start3=[1 0.005 0.01 1]; %linear imsd %startS2=[D_1 omega^2]; startS2=[0.001 1]; %parabolic imsd %startS=[D_1 v_1 omega^2]; startS=[0.001 0.005 1]; %%% set rho_\infty = cont_0*omega+fact*vdrift*TauLimit %%%%%%%%% cont_0=1; fact=0; K=double(TauLimit); P=size(A,1); R=px_size*P; A=double(A); BB=A; off2=mean(mean(mean(A))); K1=K*ones(N,1); K2=K1; for k=1:K K1(k)=k-1; K2(N-k+1)=k-1; end A7=A; if Filter==0 else A5=A-av_toll; A5(find(A5<0))=0; A5=A5+av_toll; A5(find(A5==av_toll))=0; A=A5; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot1 dim=0:px_size:R-px_size; orig_h0 = hist (A7(:,:,:),50); check_H0 = sum(sum(orig_h0,2),1); orig_H0 = sum(orig_h0,2); filt_h0 = hist (A(:,:,:),50); filt_H0 = sum(filt_h0,2); max_bar=max(orig_H0); figure (1) subplot (2,3,1) pcolor (dim,dim,A7(:,:,1)) axis square shading interp colormap default title ('Raw image','fontsize',14) %caxis ([0 255]) xlabel ('x (\mum)','fontsize',14) ylabel ('y (\mum)','fontsize',14) set (gca,'fontsize',14) subplot (2,3,3) bar (2^bit/50*(2:50),orig_H0(2:end)) axis square xlabel ('Intensisty','fontsize',14) ylabel ('Counts','fontsize',14) % set(gca,'YScale','log') %ylim ([0 max_bar*(1.05)]) xlim ([0 2^bit]) set (gca,'fontsize',14) subplot (2,3,2) surf (dim,dim,A7(:,:,1)) axis square shading interp xlim ([0 R]) ylim ([0 R]) %zlim ([-10 260]) box on xlabel ('x (\mum)','fontsize',14) ylabel ('y (\mum)','fontsize',14) zlabel ('Intensisty','fontsize',14) set (gca,'fontsize',14) subplot (2,3,4) pcolor (dim,dim,A(:,:,1)) axis square shading interp title ('Bkg subtraction','fontsize',14) %caxis ([0 255]) xlabel ('x (\mum)','fontsize',14) ylabel ('y (\mum)','fontsize',14) set (gca,'fontsize',14) subplot (2,3,6) bar (2^bit/50*(2:50),filt_H0(2:end)) axis square xlabel ('Intensisty','fontsize',14) ylabel ('Counts','fontsize',14) xlim ([0 2^bit]) set (gca,'fontsize',14) subplot (2,3,5) surf (dim,dim,A(:,:,1)) axis square shading interp xlim ([0 R]) ylim ([0 R]) xlabel ('x (\mum)','fontsize',14) ylabel ('y (\mum)','fontsize',14) zlabel ('Intensisty','fontsize',14) set (gca,'fontsize',14) %zlim ([-10 260]) box on %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% spatial autocorrelation function %%%%%%%%% a = (-R/2:px_size:(R-px_size)/2); b = (-R/2:px_size:(R-px_size)/2); disp (' calculating g(\xi, eta, 0)_n ...') for k=1:N prc = sprintf('%d',k); disp ([' computing g(\xi, eta, 0)_' prc]) scale(k)=mean(mean(A(:,:,k)))*P; %Z(:,:,k)=abs(fftshift(ifft2(fft2(A(:,:,k)).*conj(fft2(A(:,:,k))))))./(scale(k)^2*double(TauLimit)); Z(:,:,k)=(abs(fftshift(ifft2(fft2(A(:,:,k)).*conj(fft2(A(:,:,k))))))-mean(mean(A(:,:,k))))./(scale(k)^2*double(TauLimit)); ZB(:,:,k)=(abs(fftshift(ifft2(fft2(BB(:,:,k)).*conj(fft2(BB(:,:,k))))))-mean(mean(BB(:,:,k))))./(scale(k)^2*double(TauLimit)); end Z_mean=mean(Z,3); a1 = (-P/cutoff*px_size:px_size:(P-1)/cutoff*px_size); b1 = a1; %Z_cut=Z(P/2+1+int64(a1/px_size),P/2+1+int64(b1/px_size),:); %Z1=Z_mean(P/2+1+int64(a1/px_size),P/2+1+int64(b1/px_size)); Z2=Z(P/2+1+int64(a1/px_size),P/2+1+int64(a1/px_size),:); A2=A(P/2+1+int64(a1/px_size),P/2+1+int64(b1/px_size),:); ZB2=ZB(P/2+1+int64(a1/px_size),P/2+1+int64(a1/px_size),:); for i=1:N startGxy=[10 0.85 0.0]; lGxy=[0 0 -10]; uGxy=[inf inf inf]; [X Y]=meshgrid(a1,a1); xy(:,:,1)=X; xy(:,:,2)=Y; options=optimset('display','off'); modelFunGxy=@(pgxy,xy)(pgxy(1)*exp(-((xy(:,:,1)).^2+(xy(:,:,2)).^2)./pgxy(2).^2)+pgxy(3)); %options = optimset('Display','off'); Gxy=lsqcurvefit(modelFunGxy,startGxy,xy,Z2(:,:,i),lGxy,uGxy,options); Gaussxy(:,:,i)=modelFunGxy(Gxy,xy); s_Gxy(i,1)=Gxy(2); o_Gxy(i,1)=Gxy(3); %s_G(1,1)=omega^2; Gxy_x(:,i)= Gaussxy(:,length(a1)/2+1,i); Gxy_y(:,i)= Gaussxy(length(a1)/2+1,:,i); end Z1=mean(Z2,3); for i=1:N Gx(:,i)= Z2(:,length(a1)/2+1,i); Gy(:,i)= Z2(length(a1)/2+1,:,i); end omega2=mean(s_Gxy,1); err_omega2=std(s_Gxy,1); omega=omega2; err_omega=err_omega2; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% temporal correlation function %%%% t_lg=TauLimit; [X,Y,num_frame] = size(A); STCorr = zeros(X,Y,TauLimit); ST_std = zeros(X,Y,TauLimit); SeriesMean = squeeze(mean(mean(A))); imgser_fluct = zeros(X,Y,num_frame); for i = 1:num_frame imgser_fluct(:,:,i) = A(:,:,i) - SeriesMean(i); end for tau = 0 : TauLimit-1 clc disp([' iMSD ']) disp (' input parameters') disp ([' image size (px) = ' sprintf('%d*%d',P,P)]) disp ([' pixel size (\mum) = ' sprintf('%.3f',px_size)]) disp ([' time lapse (s) = ' sprintf('%.3f',f)]) disp ([' roi_stics size (px) = ' sprintf('%d*%d',2*P/cutoff,2*P/cutoff)]) disp ([' lag_time (s) = ' sprintf('%.2f',(TauLimit-1).*f)]); disp ([' n_frames = ' sprintf('%d',N)]) disp (' ') prc = sprintf('%.d',(100*(tau)/(TauLimit-1))); disp (' fitting g(\xi, \eta, 0) done') disp ([' omega^2 (\mum^2) = ' sprintf('%.6f', omega^2) ' +- ' sprintf('%.6f',2*err_omega*omega)]) disp (' ') disp ([' computing g(\xi, \eta, \tau) ' prc '%']) lagcorr = zeros(X,Y,(num_frame-tau)); for pair=1:(num_frame-tau) lagcorr(:,:,pair) = fftshift(real(ifft2(fft2(imgser_fluct(:,:,pair)).*conj(fft2(imgser_fluct(:,:,(pair+tau))))))); lagcorr(:,:,pair) = lagcorr(:,:,pair)/SeriesMean(pair)/SeriesMean(pair+tau); end STCorr(:,:,(tau+1)) = mean(lagcorr,3); ST_std(:,:,(tau+1)) = std(lagcorr,0,3); end STCorr_a=STCorr(P/2+1+int64(a1/px_size),P/2+1+int64(b1/px_size),:); W(:,1)=STCorr(P/2+1,P/2+1,:)/(P^2*double(t_lg)); t(:,1)=double(0:t_lg-1); start2=[1 0.0025 1 0.05 1]; t=f*t; t1=t(2:end); W1=W(2:end); STCorr=STCorr/(P^2*double(t_lg)); STCorr_a=STCorr(P/2+1+int64(a1/px_size),P/2+1+int64(b1/px_size),:); % % % % % % iMSD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:t_lg B=STCorr_a(:,:,i); [value, location] = max(B(:)); [R,C] = ind2sub(size(B),location); RR(i,1)=R; CC(i,1)=C; end o2=0; o1=0; ws = warning('off','all'); for i=1:t_lg Iy(:,i)= STCorr_a(:,CC(i,1),i); Ix(:,i)= STCorr_a(RR(i,1),:,i); modelFunyy = @(py,a1)(py(1).*exp(-(a1-py(4)).^2./py(2).^2)+o2); options = statset('nlinfit'); %options.Weigths=[zeros(T) s0]; options.RobustWgtFun= 'bisquare'; options.Display= 'Off'; startG=[1 0 0 1 0]; startG1=[1 0 0 1 0]; startG12=[1 0 0 1 0 0]; lG=[0 -100 -100 0 0 -10]; uG=[100 100 100 100 100 100]; [X Y]=meshgrid(a1,a1); xy(:,:,1)=X; xy(:,:,2)=Y; options=optimset('display','off'); if i==1 modelFunG=@(pg,xy)(pg(1)*exp(-((xy(:,:,1)-pg(2)).^2+(xy(:,:,2)-pg(3)).^2)./pg(4).^2)+pg(5)); %options = optimset('Display','off'); G=lsqcurvefit(modelFunG,startG,xy,STCorr_a(:,:,i),[],[],options); Gaussb(:,:,1)=modelFunG(G,xy); x_G(i,1)=G(2); y_G(i,1)=G(3); s_G(i,1)=G(4).^2; %s_G(1,1)=omega^2; o_G(i,1)=G(5); O_G=o_G; else modelFunG=@(pg,xy)(pg(1)*exp(-((xy(:,:,1)-pg(2)).^2+(xy(:,:,2)-pg(3)).^2)./pg(4).^2)+O_G); G=lsqcurvefit(modelFunG,startG12,xy,STCorr_a(:,:,i),[lG uG],[],options); Gaussb(:,:,i)=modelFunG([G(1) -G(2) G(3) G(4)],xy); x_G(i,1)=-G(2); y_G(i,1)=G(3); s_G(i,1)=G(4).^2; end Gauss_y(:,i)=Gaussb(:,size(STCorr_a,1)/2+1+int64(y_G(i,1)*px_size),i); Gauss_x(:,i)=Gaussb(size(STCorr_a,1)/2+1+int64(x_G(i,1)*px_size),:,i); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% drift fitting %%%%%%%%%%%%%% startX=[0 0]; modelFunX = @(pX,t)(pX(1).*t+pX(2)); options = statset('nlinfit'); options.RobustWgtFun= 'bisquare'; t_1=t(1:t_lg); [wnlmX, RX,JX,CovX] = nlinfit(t_1,x_G,modelFunX,startX,options); errX=sqrt(diag(CovX)); rsqX = 1-(sum(RX.^2,1) ./ sum((x_G - mean(x_G)).^2,1)); resultsX=[wnlmX' errX]; diff=omega^2/(4*wnlmX(2)); fitFunX=modelFunX(wnlmX,t_1); v_driftX=resultsX(1,1); err_v_driftX=resultsX(1,2); startY=[0 0]; modelFunY = @(pY,t)(pY(1).*t+pY(2)); options = statset('nlinfit'); options.RobustWgtFun= 'bisquare'; t_1=t(1:t_lg); [wnlmY, RY,JY,CovY] = nlinfit(t_1,y_G,modelFunY,startY,options); errY=sqrt(diag(CovY)); rsqY = 1-(sum(RY.^2,1) ./ sum((y_G - mean(y_G)).^2,1)); resultsY=[wnlmY' errY]; diff=omega^2/(4*wnlmY(2)); fitFunY=modelFunX(wnlmY,t_1); v_driftY=resultsY(1,1); err_v_driftY=resultsY(1,2); v_drift1=sqrt(v_driftX^2+v_driftY^2); err_v_drift1=(2*err_v_driftX*v_driftX+2*err_v_driftY*v_driftY)/(2*v_drift1); t_1=t(2:end); v_drift=sqrt(x_G(end)^2/(double(t_lg))^2+y_G(end)^2/(double(t_lg))^2); v_drift_x=x_G(end)/(double(t_lg)); v_drift_y=y_G(end)/(double(t_lg)); if v_driftY>0 ALPHA=acos(v_driftX/v_drift1); else ALPHA=-acos(v_driftX/v_drift1); end Sigma=s_G; Sigma_1=Sigma(2:end); y_th(1:t_lg,1)=px_size; square = [-px_size -px_size; -px_size px_size; px_size px_size; px_size -px_size; -px_size -px_size]; disp(' ') disp ([' v_\phi (\mum/s) = ' sprintf('%.6f', v_drift1) ' +- ' sprintf('%.6f',err_v_drift1)]) disp ([' \phi (rad) = ' sprintf('%.3f', ALPHA)]) disp (' _____________') disp (' iMSD analysis') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fitting Sigma options = statset('nlinfit'); options.Robust= 'Off '; lb=[0 0.00005 0]; ub=[10 0.1 1]; t_2=t(2:end); Sigma2=Sigma(2:end); lb=[0 0 0 ]; ub=[1 1 10 ]; modelFunS = @(pS)(pS(3)+4*pS(1).*t_2+pS(2).^2.*t_2.^2-Sigma2); [SS,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(modelFunS,startS,lb,ub,options); dfe = length(t_2)- 3; [Q,R] = qr(jacobian,0); Rinv = inv(R); se = sqrt(sum(Rinv.*Rinv,2)*resnorm/dfe); rsqS=1-(sum(residual.^2,1) ./ sum((Sigma2 - mean(Sigma2)).^2,1)); fitFunS=modelFunS(SS)+Sigma2; par=zeros(3,1); par(1)=se(1); par(2)=se(2); par(3)=se(3); D_Sigma=SS(1); v_Sigma=SS(2); err_D_Sigma=par(1); err_v_Sigma=par(2); lb=[0.00005 0]; ub=[0.1 10]; modelFunS2 = @(pS2)(pS2(2)+4*pS2(1).*t_2-Sigma2); [SS2,resnorm2,residual2,exitflag,output,lambda,jacobian2]=lsqnonlin(modelFunS2,startS2,lb,ub,options); dfe2 = length(t_2)- 3; [Q2,R2] = qr(jacobian2,0); R2inv = inv(R2); se2 = sqrt(sum(R2inv.*R2inv,2)*resnorm2/dfe2); rsqS2=1-(sum(residual2.^2,1) ./ sum((Sigma2 - mean(Sigma2)).^2,1)); fitFunS2=modelFunS2(SS2)+Sigma2; par2=zeros(2,1); par2(1)=se2(1); par2(2)=se2(2); D_Sigma2=SS2(1); err_D_Sigma2=par2(1); v_tot_Sigma=sqrt(v_Sigma^2+v_drift1^2); err_v_tot_Sigma=1/(2*v_tot_Sigma)*(2*v_drift1*err_v_drift1+2*v_Sigma*err_v_Sigma); A_simm=(v_Sigma/v_tot_Sigma)^2; A_asimm=(v_drift1/v_tot_Sigma)^2; a_simm=sqrt(A_simm); a_asimm=sqrt(A_asimm); factors=[A_asimm A_simm; a_asimm a_simm]; theta1=(0:0.005:2*pi); ftheta1=(2./theta1.*sin(theta1./2)).^2; g=A_asimm; vec=find(g>=ftheta1); beta=theta1(vec(1)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% anomalous diffusion lbc=[10^-4 10^-4 0]; ubc=[Inf Inf 2]; startc=[0 10^-3 0.8]; modelFun_conf = @(pS)(pS(1)+4*pS(2).*t_2.^pS(3)-Sigma2); [SSc,resnormc,residualc,exitflagc,outputc,lambdac,jacobianc]=lsqnonlin(modelFun_conf,startc,lbc,ubc,options); dfec = length(t_2)-3; [Qc,Rc] = qr(jacobianc,0); Rinvc = inv(Rc); sec = sqrt(sum(Rinvc.*Rinvc,2)*resnormc/dfec); rsqS_conf=1-(sum(residualc.^2,1) ./ sum((Sigma2 - mean(Sigma2)).^2,1)); fitFun_conf=modelFun_conf(SSc)+Sigma2; par_c=zeros(3,1); par_c(1)=sec(1); par_c(2)=sec(2); par_c(3)=sec(3); pow=SSc(3); err_pow=par_c(3); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% confined diffusion lbc2=[0 10^-4 0 0]; ubc2=[Inf 10^-1 Inf 8]; startc2=[0 10^-3 0.5 1]; modelFun_conf2 = @(pC)(pC(1)+4*pC(2).*t_2+pC(3)./3.*(1-exp(-t_2./pC(4)))-Sigma2); [SSc2,resnormc2,residualc2,exitflagc2,outputc2,lambdac2,jacobianc2]=lsqnonlin(modelFun_conf2,startc2,lbc2,ubc2,options); dfec2 = length(t_2)-4; [Qc2,Rc2] = qr(jacobianc2,0); Rinvc2 = inv(Rc2); sec2 = sqrt(sum(Rinvc2.*Rinvc2,2)*resnormc2/dfec2); rsqS_conf2=1-(sum(residualc2.^2,1) ./ sum((Sigma2 - mean(Sigma2)).^2,1)); fitFun_conf2=modelFun_conf2(SSc2)+Sigma2; par_c2=zeros(4,1); par_c2(1)=sec2(1); par_c2(2)=sec2(2); par_c2(3)=sec2(3); par_c2(4)=sec2(4); D_M=SSc2(2); err_D_M=par_c2(2); D_m=SSc2(3)./(12.*SSc2(4))+D_M; err_D_m=err_D_M+(1/12)*SSc2(3)/SSc2(4)*(par_c2(3)/SSc2(3)+par_c2(3)/SSc2(4)); Lquadro=SSc2(3); tauc=SSc2(4); disp(' ') disp (' - Brownian diffusion') disp ([' offset (\mum^2) = ' sprintf('%.6f', SS2(2)) ' +- ' sprintf('%.6f',par2(2))]) disp ([' D (\mum^2/s) = ' sprintf('%.6f', D_Sigma2) ' +- ' sprintf('%.6f',err_D_Sigma2)]) disp ([' R^2 = ' sprintf('%.6f', rsqS2)]) disp(' ') disp (' - Anomalous diffusion') disp ([' offset (\mum^2) = ' sprintf('%.6f', SSc(1)) ' +- ' sprintf('%.6f',par_c(1))]) disp ([' \alpha = ' sprintf('%.6f', pow) ' +- ' sprintf('%.6f',err_pow)]) disp ([' R^2 = ' sprintf('%.6f', rsqS_conf)]) disp(' ') if pow<1 disp (' - Confined diffusion') disp ([' offset (\mum^2) = ' sprintf('%.6f', SSc2(1)) ' +- ' sprintf('%.6f',par_c2(1))]) disp ([' D_M (\mum^2/s) = ' sprintf('%.6f', D_M) ' +- ' sprintf('%.6f',err_D_M)]) disp ([' D_m (\mum^2/s) = ' sprintf('%.6f', D_m) ' +- ' sprintf('%.6f',err_D_m)]) disp ([' L^2 (\mum^2) = ' sprintf('%.6f', Lquadro) ' +- ' sprintf('%.6f',par_c2(3))]) disp ([' \tau_c (s) = ' sprintf('%.6f', tauc) ' +- ' sprintf('%.6f',par_c2(4))]) disp ([' R^2 = ' sprintf('%.6f', rsqS_conf2)]) %disp(' ') else disp (' - Brownian diffusion + flow motion') disp ([' offset (\mum^2) = ' sprintf('%.6f', SS(3)) ' +- ' sprintf('%.6f',par(3))]) disp ([' D (\mum^2/s) = ' sprintf('%.6f', D_Sigma) ' +- ' sprintf('%.6f',err_D_Sigma)]) disp ([' v_\sigma (\mum/s) = ' sprintf('%.6f', v_Sigma) ' +- ' sprintf('%.6f',err_v_Sigma)]) disp ([' R^2 = ' sprintf('%.6f', rsqS)]) disp(' ') %disp ([' v_\sigma (\mum/s) = ' sprintf('%.6f', SS(2)) ' +- ' sprintf('%.6f',par(2))]) disp ([' v_tot (\mum/s) = ' sprintf('%.6f', v_tot_Sigma) ' +- ' sprintf('%.6f',err_v_tot_Sigma)]) %disp (' ') %disp ([' {a_\phi}^2 = ' sprintf('%.3f', A_asimm)]) %disp ([' {a_\sigma}^2 = ' sprintf('%.3f', A_simm)]) %disp ('') disp ([' \psi (rad) = ' sprintf('%.3f', beta)]) %disp('') end disp (' ______________________________________________') disp('') disp(' DONE') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% g(\rho, \theta, \tau) %%%%%%%% Cont=STCorr_a(size(STCorr_a,1)/2+1+int64((cont_0*omega+0*v_drift*double(t_lg))/px_size),... (size(STCorr_a,1)/2)+1+int64((cont_0*omega+0*v_drift*double(t_lg))/px_size),end); figure (10) [c,h]=contour(1:length(a1),1:length(a1),STCorr_a(:,:,1),Cont); [d,h1]=contour(a1,a1,STCorr_a(:,:,1),Cont); close (10) cc=c(:,2:end); dd_=d(:,2:end); dd_=dd_'; tht=2*pi/50; theta_=-pi:tht:pi-tht; %rho=[dd_(1,1) dd_(1,2)]; %RHO=sqrt(rho(1)^2+rho(2)^2); RHO=cont_0*omega+fact*v_drift*double(t_lg); dd=[-RHO*cos(theta_') RHO*sin(theta_')]; DD_= int64([-RHO*cos(theta_') RHO*sin(theta_')]/px_size); DD=int64([-RHO*cos(theta_') RHO*sin(theta_')]/px_size)+size(STCorr_a,1)/2+1; theta(1:length(theta_)/2)=-acos(double(DD_(length(theta_)/2+1:end,1))/(RHO/px_size)); theta(length(theta_)/2+1:length(theta_))=acos(double(DD_(1:length(theta_)/2,1))/(RHO/px_size)); theta(1)=theta_(1); for j=2:length(theta_) if theta_(j)>=theta_(j-1) theta(j)=theta_(j); end end cc=int64(cc'); mid=mean(double(cc(:,1))); dc=double(cc-mid); tht=2*pi/length(cc); d_theta=zeros(length(DD),t_lg); t2=t/f; tau=t2+1; for k=1:t_lg for i=1:length(theta); d_theta(i,k)=STCorr_a(DD(i,2),DD(i,1),int64(tau(k))); end end F_theta=d_theta(:,end); Ftheta=F_theta'; lbFt=[-1 -pi 0 0 -1 -pi 0]; ubFt=[1 pi 10 10 1 pi 10]; startFt=[0.1 ALPHA 0.4 0.2 0.1 -ALPHA 0.3]; options.Robust='off'; options.MaxIter=500; ref1=int64(mean(find(int64(theta-ALPHA)==0))); Gthetapm=[Ftheta Ftheta Ftheta]; Gtheta=Gthetapm(length(theta)+ref1-length(theta)/2:length(theta)+ref1+length(theta)/2-1); Htheta(1:length(theta)/2)=Gtheta(length(theta)/2+1:end); Htheta(length(theta)/2+1:length(theta))=Gtheta(1:length(theta)/2); %%%%%%%%%%%%%%%%%%%%%%%%% lGt=[0 -pi 0 -10]; uGt=[10 pi 2*pi 10]; startGt=[0.1 0 0.1 0.1]; modelFunGt = @(pG)(pG(1).*((sinc(pG(3).*(theta-pG(2)))).^2)+pG(4)-Gtheta); [SSGt,resnormGt,residualGt,exitflagGt,outputGt,lambdaGt,jacobianGt]=lsqnonlin(modelFunGt,startGt,lGt,uGt,options); dfeGt = length(theta)- 4; [QGt,RGt] = qr(jacobianGt,0); RGtinv = inv(RGt); seGt = sqrt(sum(RGtinv.*RGtinv,2)*resnormGt/dfeGt); rsqGt=1-(sum(residualGt.^2,2) ./ sum((Gtheta - mean(Gtheta)).^2,2)); fitFunGt=modelFunGt(SSGt)+Gtheta; %%%%%%%%%%%%%%%%%%%%%%%%%%% lHt=[0 -pi 0 -10]; uHt=[10 pi 2*pi 10]; startHt=[0.1 0 0.1 0.1]; modelFunHt = @(pH)(pH(1).*(1-(sinc(pH(3).*((theta+pi)-pH(2)))).^2)+pH(4)-Htheta); [SSHt,resnormHt,residualHt,exitflagHt,outputHt,lambdaHt,jacobianHt]=lsqnonlin(modelFunHt,startHt,lHt,uHt,options); dfeHt = length(theta)- 4; [QHt,RHt] = qr(jacobianHt,0); RHtinv = inv(RHt); seHt = sqrt(sum(RHtinv.*RHtinv,2)*resnormHt/dfeHt); rsqHt=1-(sum(residualHt.^2,2) ./ sum((Htheta - mean(Htheta)).^2,2)); fitFunHt=modelFunHt(SSHt)+Htheta; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% lbFt=[-1 -pi 0 0 -1 -pi 0]; ubFt=[1 pi 1 1 1 pi 1]; startFt=[0.1 ALPHA 0.4 0.2 0.1 pi-ALPHA 0.3]; modelFunFt = @(pF)(pF(1).*((sinc(pF(3).*(theta-pF(2)))).^2)+pF(5).*(1-(sinc(pF(7).*(theta-pF(2)))).^2)+pF(4)-Ftheta); [SSFt,resnormFt,residualFt,exitflagFt,outputFt,lambdaFt,jacobianFt]=lsqnonlin(modelFunFt,startFt,lbFt,ubFt,options); dfeFt = length(theta)- 7; [QFt,RFt] = qr(jacobianFt,0); RFtinv = inv(RFt); seFt = sqrt(sum(RFtinv.*RFtinv,2)*resnormFt/dfeFt); rsqFt=1-(sum(residualFt.^2,2) ./ sum((Ftheta - mean(Ftheta)).^2,2)); fitFunFt=modelFunFt(SSFt)+Ftheta; par2=zeros(2,1); par2(1)=seFt(2); par2(2)=seFt(3); phi_Ft=SSFt(2); psi_Ft=2*SSFt(3)/(pi); err_phi_Ft=par2(1); err_psi_Ft=par2(2); ang_distr=[phi_Ft err_phi_Ft psi_Ft err_psi_Ft rsqFt]; thr1=zeros(length(t))+px_size/(2*sqrt(2)); thr2=zeros(length(t))+px_size/(2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure (2) subplot (2,3,1) pcolor (dim,dim,A(:,:,1)) axis square shading interp title ('image','fontsize',14) xlabel ('x (\mum)','fontsize',14) ylabel ('y (\mum)','fontsize',14) set (gca,'fontsize',14) subplot (2,3,2) pcolor (-a1,a1,STCorr_a(:,:,end)) shading interp colormap jet hold on plot (x_G,y_G,'-b',x_G(1),y_G(1),'black+',x_G(end),y_G(end),'blacko','linewidth',1,'markersize',8) axis square plot (dd(:,1),dd(:,2),'--black','linewidth',1.1) %plot (cc(:,1),cc(:,2)) title ('g -> top view','fontsize',14) xlabel ('\xi (\mum)','fontsize',14) ylabel ('\eta (\mum)','fontsize',14) set (gca,'fontsize',14) hold off subplot (2,3,3) p3=zeros(length(a1),1)+max(max(a1))+px_size/2; p4=zeros(length(a1),1)-min(min(a1))-px_size/2; surf (-a1-0*px_size/2,a1,STCorr_a(:,:,end)); hold on plot3 (a1,p4,Gauss_x,'black','linewidth',0.8) plot3 (p3,a1,Gauss_y,'black','linewidth',0.8); hold off xlim([min(min(a1)) max(max(a1))+px_size/2]) ylim([min(min(a1)) max(max(a1))]) xlabel '\xi' ylabel '\eta' zlim ([min(min(min(STCorr_a))) max(max(max(STCorr_a)))]) grid on axis square zlabel ('g(\xi,\eta,\tau)','fontsize',14) xlabel ('\xi (\mum)','fontsize',14) ylabel ('\eta(\mum)','fontsize',14) title ('g -> 3d view','fontsize',14) set (gca,'fontsize',14) box on subplot (2,3,4) plot (t,x_G,'blackx',t,y_G,'black+','linewidth',1.5) legend ('\xi-shift','\eta-shift','location','southwest') legend ('boxon') hold on plot (t,thr2,'-.r',t,-thr2,'-.r','linewidth',1.3) plot (t,thr1,':r',t,-thr1,':r','linewidth',1.1) hold off xlabel ('\tau (s)','fontsize',14) ylabel ('peak position (\mum)','fontsize',14) xlim ([0-f t_2(end)+f]) set (gca,'fontsize',14) title ('net flow','fontsize',14) axis square subplot (2,3,5) plot(t_2,fitFunS2,'--r','linewidth',1.5) hold on %plot(t_2,fitFunS,'--b','linewidth',1.5) plot(t_2,fitFun_conf,'-c','linewidth',1.5) if pow<1 plot(t_2,fitFun_conf2,'--g','linewidth',1.5) else plot (t_2,fitFunS,'--b','linewidth',1.5) end plot(t_2,Sigma2,'blacko','linewidth',1.5) hold off if pow<1 legend ('Brownian','Anomalous','Confined','location','southeast') else legend ('Brownian','Anomalous','Flow','location','southeast') end legend ('boxon') xlabel ('\tau (s)','fontsize',14) ylabel ('\sigma^2 (\mum^2)','fontsize',14) xlim ([0 t_2(end)+f]) set (gca,'fontsize',14) title ('iMSD curve','fontsize',14) axis square subplot (2,3,6) h=surf(t,theta,d_theta); hold on view([-37.5+90 25]) shading interp xlabel ('\tau (s)','FontSize',14) ylabel ('\theta (rad)','FontSize',14) zlabel ('g(\rho_\infty, \theta, \tau)','fontsize',14) grid on colormap jet hold off ylim ([-pi pi]) xlim ([0 t_2(end)+f]) set (gca,'FontSize',14) set(h','edgecolor','k') axis square title ('g -> circular section','fontsize',14) box on %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure (3) subplot (1,3,1) plot(t_2,Sigma2,'blacko','linewidth',1.5) hold on plot(t_2,fitFunS2,'r','linewidth',1.5) title ('Brownian Diffusion','FontSize',14) legend ('iMSD data',sprintf('R^2=%.5f',rsqS2),'location','northwest') legend ('boxoff') hold off xlabel ('\tau (s)','fontsize',14) ylabel ('\sigma^2 (\mum^2)','fontsize',14) xlim ([0 t_2(end)+f]) set (gca,'fontsize',14) axis square subplot (1,3,2) plot(t_2,Sigma2,'blacko','linewidth',1.5) hold on plot(t_2,fitFun_conf,'c','linewidth',1.5) legend ('iMSD data',sprintf('R^2=%.5f',rsqS_conf),'location','northwest') legend ('boxoff') title (sprintf('Anomalous Diffusion \\alpha = %.3f',pow),'FontSize',14) hold off xlabel ('\tau (s)','fontsize',14) ylabel ('\sigma^2 (\mum^2)','fontsize',14) xlim ([0 t_2(end)+f]) set (gca,'fontsize',14) axis square if pow<1 subplot (1,3,3) plot(t_2,Sigma2,'blacko','linewidth',1.5) hold on plot(t_2,fitFun_conf2,'g','linewidth',1.5) title ('Confined Diffusion','FontSize',14) legend ('iMSD data',sprintf('R^2=%.5f',rsqS_conf2),'location','northwest') legend ('boxoff') hold off xlabel ('\tau (s)','fontsize',14) ylabel ('\sigma^2 (\mum^2)','fontsize',14) xlim ([0 t_2(end)+f]) set (gca,'fontsize',14) axis square else subplot (1,3,3) plot(t_2,Sigma2,'blacko','linewidth',1.5) hold on plot(t_2,fitFunS,'b','linewidth',1.5) legend ('iMSD data',sprintf('R^2=%.5f',rsqS),'location','northwest') legend ('boxoff') title ('Brownian Diffusion + Flow','FontSize',14) hold off xlabel ('\tau (s)','fontsize',14) ylabel ('\sigma^2 (\mum^2)','fontsize',14) xlim ([0 t_2(end)+f]) set (gca,'fontsize',14) axis square end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% END %%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% export %%%% pre=[name_start '.xlsx']; delete (pre); head1=cellstr('t'); head2=cellstr('\xi_0'); head3=cellstr('\eta_0'); head4=cellstr('Sigma^2'); head=[head1 head2 head3 head4]; filename=[name_start '.xlsx']; xlswrite(filename,head,'Foglio1') xlswrite(filename,[t_2 x_G(2:end) y_G(2:end) Sigma_1],'Foglio1','A2') txt1=cellstr('INPUT'); txt2=cellstr('sqrt(image size) (px)'); txt3=cellstr('pixel size (\mum)'); txt4=cellstr('sqrt(fitting-ROI size) (px)'); txt5=cellstr('time lapse (s)'); txt6=cellstr('max lag-time (s)'); txt7=cellstr('number of frames'); text1=cellstr('omega (\mum)'); text2=cellstr('D Brownian (\mum^2/s)'); text3=cellstr('offset (\mum^2)'); text4=cellstr('R^2'); text5=cellstr('\alpha Anomalous'); text6=cellstr('offset (\mum)^2'); text7=cellstr('R^2'); text8=cellstr('D_M (\mum^2/s)'); text9=cellstr('D_m (\mum^2/s)'); text10=cellstr('L^2 (\mum^2)'); text11=cellstr('\tau_c (s)'); text12=cellstr('offset (\mum^2)'); text13=cellstr('R^2'); text14=cellstr('D Br.+Flow (\mum^2/s)'); text15=cellstr('v_\phi (\mum/s)'); text16=cellstr('v_\sigma (\mum/s)'); text17=cellstr('offset (\mum^2)'); text18=cellstr('R^2'); txt0=[txt1;txt2;txt3;txt4;txt5;txt6;txt7;cellstr('bkg-corr')]; nm0=[P;px_size;2*P/cutoff;f;t_2(end);double(N);Filter]; text0=[text1;text15;cellstr('\phi');text2;text3;text4;text5;text6;text7;text8;text9;text10;text11;text12;text13;text14;text16;text17;text18]; head2= [cellstr('OUTPUT') cellstr('mean') cellstr('error')]; nm2=[omega^2;v_drift1;ALPHA;D_Sigma2;SS2(2);rsqS2;pow;SSc(1);rsqS_conf;D_M;D_m;Lquadro;tauc;SSc2(1);rsqS_conf2;D_Sigma;v_Sigma;SS(3);rsqS]; err2=[2*omega*err_omega;err_v_drift1;NaN; err_D_Sigma2;par2(2);NaN;err_pow;par_c(1);NaN;err_D_M;err_D_m;par_c2(3);par_c2(4);par_c2(1);NaN;err_D_Sigma;err_v_Sigma;par(3);NaN]; xlswrite(filename,txt0,'Foglio2','A1') xlswrite(filename,nm0,'Foglio2','B2') xlswrite(filename,head2,'Foglio2','A10') xlswrite(filename, text0,'Foglio2','A11') xlswrite(filename, [nm2 err2],'Foglio2','B11') xlswrite(filename, cellstr(sprintf('x%d',P)),'Foglio2','C2') xlswrite(filename, cellstr(sprintf('x%d',2*P/cutoff)),'Foglio2','C4') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% animation return figure for i=1:t_lg subplot (1,2,2) pcolor (-a1, b1,STCorr_a(:,:,i)) shading interp; axis square T=sprintf('%d',i-1); title (['r(\xi,\eta, ' T ')'],'FontSize',12) %set (xticks, ([min(a1):max(a1)])) xlabel ('\xi (\mum)','FontSize',12) ylabel ('\eta (\mum)','FontSize',12') set(gca,'FontSize',12) subplot (1,2,1) surf (-a1, b1,STCorr_a(:,:,i)) %shading interp; axis square T=sprintf('%d',i-1); zlabel (['r(\xi,\eta, ' T ')'],'FontSize',12) zlim ([0 max(max(max(STCorr)))]) xlabel ('\xi (\mum)','FontSize',12) ylabel ('\eta (\mum)','FontSize',12) %xticks (min(a1):max(a1)) xlim ([min(a1) max(a1)]) ylim ([min(a1) max(a1)]) set(gca,'FontSize',12) box on pause(0.5) %saveas(gcf,sprintf('export_AUTOC%d.png',i-1)) colormap jet end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot (2,3,6) plot (theta,F_theta,':blacko','linewidth',1) xlabel ('\theta (rad)','FontSize',12) ylabel ('r(\rho_\infty, \theta, \tau_{M})','FontSize',12) hold on plot (theta,fitFunFt,'black','linewidth',1.5') hold off set (gca,'FontSize',12) xlim ([-pi-0.05 pi+0.05])