diff --git a/iPED.m b/iPED.m new file mode 100644 index 0000000..3a727da --- /dev/null +++ b/iPED.m @@ -0,0 +1,284 @@ +function [iPED_PDF,converged]=... + iPED(imdir,filename,numdig,filetype,resave,dofilter,N,corstep,gausblur,mbf,maf,ensemf,windowz,stradlemod,imstart,bkga) +% imdir ='E:\bacteria_data\ecoli\pdf_test\'; +% filename='20x_hcb1736_secondaryculture_od0_t'; +% numdig= '%04i' +% filetype= '.tif' +% resave=strcat(imdir,'\results'); +% mkdir(resave); +% N=100; numebr of images +% corstep=1 is the correlation step size +% dofilter=1; Edge Tapering gaussian Filter on the whole image +% gausblur=0; % Gaussian blur for unresolved particles +% mbf=0; %mean subtraction before Gaussian filter +% maf=0; %mean subtraction after Gaussianfilter +% ensemf=1; %ensemble in Fourier +% windowz=64; %window size for gridding +% stradlemod=0; %this is for PIV images +% imstart=1; %starting image number +% bkga % average image of time series to be subtracted from all +% images +tic + +h = waitbar(0,'calculating...','Name','Please wait...',... + 'CreateCancelBtn','setappdata(gcbf,''canceling'',1)'); + +cropp=1; % window gridding +l=0; +filter=ones; +breakflag=0; +kahler=zeros; +pid=zeros; +straddling=stradlemod; % 1 if frame1 can only be correlated to frame2 +convergencecriteria=1; +converged=0; +FFFF=zeros; +eeerrr=zeros; + +for i=corstep:corstep %% correlations steps + % define a whole bunch of parameters, names are self explanatory + max_intensity=0; + SCC_ensembler=zeros; + AC_ensembler=zeros; + SCC_ensemblerr=zeros; + AC_ensemblerr=zeros; + SCC_ensembleFr=zeros; + AC_ensembleFr=zeros; + SCC_ensembleFrr=zeros; + AC_ensembleFrr=zeros; + AC=zeros; + meannoiselevel=zeros; + integrallevel=zeros; + L_infinity=zeros; + + ensgrids=0; + % start the PDF calculation + + for j=imstart:1+straddling:N-i + %Read image 1 + filenum=num2str(j,numdig) + matFilename = strcat(filename,filenum,filetype); + [xr1_0,map]=imread([imdir, matFilename]); + %if colored images are loaded they will be converted to grey scale + %here + if size(xr1_0,3)==3 + xr1_0=rgb2gray(xr1_0); + end + % Let's subtract the background + xr1_0_m=double(xr1_0)-double(bkga); + xr1_0=double(xr1_0); + %load image 2 + filenum=num2str(j+i,numdig); + matFilename = strcat(filename,filenum,filetype); + [xr2_0,map]=imread([imdir, matFilename]); + if size(xr2_0,3)==3 + xr2_0=rgb2gray(xr2_0); + end + xr2_0_m=double(xr2_0)-double(bkga); + xr2_0=double(xr2_0); + + + % for the purpose of image croppings if griding is selceted + + kmax=size(xr1_0_m,2); + + for k=1:windowz:kmax-windowz+1 + for kk=1:windowz:kmax-windowz+1 + % if cropp==1 + xr1_m=imcrop(xr1_0_m,[k kk windowz-1 windowz-1]); + xr2_m=imcrop(xr2_0_m,[k kk windowz-1 windowz-1]); + xr1=imcrop(xr1_0,[k kk windowz-1 windowz-1]); + xr2=imcrop(xr2_0,[k kk windowz-1 windowz-1]); + + + % apply gaussian Blur if selected + if gausblur==1 + xr1_m=imgaussfilt(xr1_m,2); + xr2_m=imgaussfilt(xr2_m,2); + end + + %convert uint8 variables to double precision variables so that + %they can be used for calculations AND subtract the mean if + %selected by the user + + + [X,Y]=meshgrid(1:size(xr1_m,2),1:size(xr1_m,1)); + + if mbf==1 + xr1d_m=double(xr1_m)-mean2(double(xr1_m)); + xr2d_m=double(xr2_m)-mean2(double(xr2_m)); + + else + xr1d_m=double(xr1_m); + xr2d_m=double(xr2_m); + end + + % Now lets apply the Gauss Edge Tapering filter + if dofilter==1 + + PSF = fspecial('gaussian',30,5); + % PSF = fspecial('gaussian',7,1); + x1f = edgetaper(xr1d_m,PSF); + x2f = edgetaper(xr2d_m,PSF); + else + x1f=xr1d_m; + x2f=xr2d_m; + end + + % if user tuned on the mean subtraction after filter + % application + if maf==1 + x1f=x1f-mean2(x1f); + x2f=x2f-mean2(x2f); + end + + + %calculate FT of the ROIs + F_CC=zeros(size(xr1_m)); + G1=fft2(x1f); + G2=fft2(x2f); + + %CC and AC + F_CC=G1.*conj(G2); + % F_AC=sqrt(G1.*conj(G1).*G2.*conj(G2)); + F_AC=0.5*(G1.*conj(G1)+G2.*conj(G2)); + + %ensemble in Spectral domain? + if ensemf==1 + SCC_ensembleFrr=SCC_ensembleFrr+(F_CC); + AC_ensembleFrr=AC_ensembleFrr+(F_AC); + l=l+1; + + else + CC=fftshift(ifft2(F_CC)); + AC=fftshift(ifft2(F_AC)); + l=l+1; + SCC_ensembler=SCC_ensembler+CC; + AC_ensembler=AC_ensembler+AC; + end + %calculate GCC-SCOT and L infinity of the PDF + if ensemf==1 + ensgrids=ensgrids+1; + GCC_ensemble=(fftshift(ifft2(SCC_ensembleFrr./abs(SCC_ensembleFrr)))); + SCC_ensembler=fftshift(ifft2(SCC_ensembleFrr)); + AC_ensembler=fftshift(ifft2(AC_ensembleFrr)); + iPED_PDF=(fftshift(ifft2(fft2(SCC_ensembler)./fft2(AC_ensembler)))); + L_infinity(ensgrids)=norm(iPED_PDF, inf); + if L_infinity(ensgrids)==inf || isnan(L_infinity(ensgrids)) + try + L_infinity(ensgrids)=L_infinity(ensgrids-1); + catch + L_infinity(ensgrids)=1 + end + end + + %convergence check using L_infinity + if convergencecriteria==1 && ensgrids>50 &&converged==0 + % L_infinity(ensgrids); + convx=(1:ensgrids); % Explanatory variable + f = @(F,x) F(1)./sqrt(x)+F(2); + lb = [0,0]; + ub = [inf,inf]; + beta0 = [10,10]; + opts = optimset('Display','off'); + [FFF,resnorm] =lsqcurvefit(f,beta0,convx,L_infinity,lb,ub,opts); + eeerrr(ensgrids)=abs(L_infinity(ensgrids)-FFF(2))/FFF(2)*100; + FFFF(ensgrids)=FFF(2); + + + if mean(eeerrr(ensgrids-49:ensgrids))<10 && std(FFFF(ensgrids-49:ensgrids))<0.1*mean(FFFF(ensgrids-49:ensgrids)) + converged=j + figure + y_predicted=f(FFF,convx); + plot(1:ensgrids,L_infinity(1:ensgrids),'*'); + hold on + plot(convx,y_predicted,'g','LineWidth',3); + plot([ensgrids ensgrids],[0 max(max(L_infinity),max(y_predicted))],'k','LineWidth',3) + axis tight + legend('data','fit','convergence') + % set(gca,'XScale','log') + pltname=strcat(resave,'/timelag',num2str(i),'convergence',num2str(j),'.fig') + saveas(gcf,pltname) + break + end + if converged + breakflag=1 + break + end + end + + + else + %store everything + iPED_PDF=(fftshift(ifft2(fft2(SCC_ensembler)./fft2(AC_ensembler)))); + meannoiselevel(j)=mean(abs(iPED_PDF(:))); + integrallevel(j)=trapz(iPED_PDF(:)); + L_infinity(j)=norm(iPED_PDF, inf); + end + + end + if getappdata(h,'canceling') + breakflag=1; + end + if converged + breakflag=1 + break + end + if getappdata(h,'canceling') + breakflag=1; + break + + end + + end + if getappdata(h,'canceling') + breakflag=1; + break + + end + + if converged + breakflag=1 + break + end + waitbar(j/(N),h,strcat('Image ', {' '}, num2str(j),{' '}, 'out of ',{' '},num2str(N-i),{' '},' .'))%{sprintf('/n')} + end + +% savemeannoise = strcat(resave,'/meannoise',num2str(j+i),'window',num2str(windowz),'_',num2str(i)); +% save(savemeannoise,'meannoiselevel','integrallevel','L_infinity') +% + + if ensemf==1 + + SCC_ensembler=fftshift(ifft2(SCC_ensembleFrr)); + AC_ensembler=fftshift(ifft2(AC_ensembleFrr)); + + end + + SCC_ensemble = SCC_ensembler; + AC_ensemble = AC_ensembler; + iPED_PDF=(fftshift(ifft2(fft2(SCC_ensemble)./fft2(AC_ensemble)))); + + %save everything + savecor = strcat(resave,'/SCC_',num2str(j+i),'window',num2str(windowz),'_',num2str(i)); + save(savecor,'SCC_ensemble') + + savecor = strcat(resave,'/AC_',num2str(j+i),'window',num2str(windowz),'_',num2str(i)); + save(savecor,'AC_ensemble') + + savecor = strcat(resave,'/iPED_PDF_',num2str(j+i),'window',num2str(windowz),'_',num2str(i)); + save(savecor,'iPED_PDF') + + + + if breakflag==1 + delete(h) + break + end +end + +toc +end + + +