Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
@aahmadza
aahmadza authored Jun 30, 2020
1 parent 366fe6b commit 586f716
Showing 1 changed file with 284 additions and 0 deletions.
284 changes: 284 additions & 0 deletions iPED.m
View file
Original file line number Diff line number Diff line change
@@ -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



0 comments on commit 586f716

Please sign in to comment.