Skip to content
Permalink
e61cf9a17e
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

iPED/iPED.m

Go to file
 
 
Cannot retrieve contributors at this time
284 lines (235 sloc) 10.4 KB
Raw Blame
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