Add files via upload

MATLAB_Open_3DSIM/Open_3DSIM.m

@@ -50,7 +50,7 @@
     NA = 1.4;                   % objective lens NA
     refmed = 1.47;              % refractive index medium
     refcov = 1.512;             % refractive index cover slip
-    refimm = 1.518;             % refractive index immersion medium
+    refimm = 1.515;             % refractive index immersion medium
     exwavelength = 488;         % excitation wavelengths
     emwavelength = 528;         % emission wavelengths
     fwd = 140e3;                % free working distance from objective to cover slip

MATLAB_Open_3DSIM/help_functions/SimOtfProvider2.m

@@ -0,0 +1,89 @@
+function [ otf ] = SimOtfProvider2( dataparams,SIMparams,NA,lambda,numpixelsx,a)
+    param.attStrength = 0;
+    param.imgSize = numpixelsx;
+    param.cyclesPerMicron = 1/(numpixelsx*SIMparams.SIMpixelsize(1)*1e-3);
+    ret.na=NA;
+    ret.lambda=lambda;
+    % ret.cutoff=1000/(0.61*lambda/NA);
+    ret.cutoff=1000/(0.5*lambda/NA);
+
+    ret.imgSize=param.imgSize;
+    ret.cyclesPerMicron=param.cyclesPerMicron;
+    ret.sampleLateral=ceil(ret.cutoff/ret.cyclesPerMicron)+1;
+
+    ret.estimateAValue=a;
+    ret.maxBand=2;
+    ret.attStrength=param.attStrength;
+    ret.attFWHM=1.0;
+    ret.useAttenuation=1;
+
+    ret=fromEstimate(ret);
+
+    ret.otf=zeros(param.imgSize,param.imgSize);
+    ret.otfatt=zeros(param.imgSize,param.imgSize);
+    ret.onlyatt=zeros(param.imgSize,param.imgSize);
+
+    ret.otf=writeOtfVector(ret.otf,ret,1,0,0);
+    otf = ret.otf./max(max(ret.otf));
+    ret.onlyatt=getonlyatt(ret,0,0);
+    ret.otfatt=ret.otf.*ret.onlyatt;
+end
+
+function [va]=valIdealOTF(dist)
+    if dist<0 || dist>1
+        va=0;
+        return;
+    end
+    va=(1/pi)*(2*acos(dist)-sin(2*acos(dist)));
+end
+
+function [va]= valAttenuation(dist,str,fwhm)
+    va=(1-str*(exp(-power(dist,2)/(power(0.5*fwhm,2)))).^1);
+end
+
+function [ret] = fromEstimate(ret)
+    ret.isMultiband=0;
+    ret.isEstimate=1;
+    vals1=zeros(1,ret.sampleLateral);
+    valsAtt=zeros(1,ret.sampleLateral);
+    valsOnlyAtt=zeros(1,ret.sampleLateral);
+
+
+    for I=1:ret.sampleLateral
+        v=abs(I-1)/ret.sampleLateral;
+        r1=valIdealOTF(v)*power(ret.estimateAValue,v);
+        vals1(I)=r1;
+    end
+
+    for I=1:ret.sampleLateral
+        dist=abs(I-1)*ret.cyclesPerMicron;
+        valsOnlyAtt(I)=valAttenuation(dist,ret.attStrength,ret.attFWHM);
+        valsAtt(I)=vals1(I)*valsOnlyAtt(I);
+    end
+
+    ret.vals=vals1;
+
+    ret.valsAtt=valsAtt;
+    ret.valsOnlyAtt=valsOnlyAtt;
+
+end
+
+function [ onlyatt ] = getonlyatt( ret,kx,ky )
+    w=ret.imgSize;
+    h=ret.imgSize;
+    siz=[h w];
+    cnt=siz/2+1;
+    kx=kx+cnt(2);
+    ky=ky+cnt(1);
+    onlyatt=zeros(h,w);
+
+    y=1:h;
+    x=1:w;
+    [x,y]=meshgrid(x,y);
+    rad=hypot(y-ky,x-kx);
+    cycl=rad.*ret.cyclesPerMicron;
+    onlyatt=valAttenuation(cycl,ret.attStrength,ret.attFWHM);
+
+end
+
+

MATLAB_Open_3DSIM/help_functions/do_modulationcheck.m

@@ -1,5 +1,7 @@
 function [MCNR_ims,allmodulations] = do_modulationcheck(allimages_in)
 
+% copyright Sjoerd Stallinga, TU Delft, 2017-2020
+
 %% define parameter
 [Nx,Ny,numsteps,numfocus,numchannels,numframes,numangles] = size(allimages_in);
 maxorder = (numsteps+1)/2; % this assumes that the # independent image Fourier orders = # phase steps

MATLAB_Open_3DSIM/help_functions/doshift.m

@@ -11,6 +11,10 @@
 qvector(3) = kz*Nz*dataparams.SIMpixelsize(3);
 qvector_z = qvector(3);
 qvector(1:2) = patternpitch;
+% single-layer
+if Nz==1
+    qvector(3)=0;
+end
 
 %% Compute notch-filter
 XSupport = ((1:Nx)-floor(Nx/2)-1); % qx grid
@@ -19,12 +23,17 @@
 [qx,qy,qz] = meshgrid(XSupport,YSupport,ZSupport); 
 notchwidthxy = dataparams.notchwidthxy1; % notch width
 notchdips = dataparams.notchdips1;       % notch depth
-qradsq = ( qx/qvector_xy ).^2+( qy/qvector_xy ).^2+( qz/qvector_z ).^2;
+if Nz ~=1
+    qradsq = ( qx/qvector_xy ).^2+( qy/qvector_xy ).^2+( qz/qvector_z ).^2;
+else
+    qradsq = ( qx/qvector_xy/2 ).^2+( qy/qvector_xy/2 ).^2+( qz/qvector_z/2 ).^2;
+end
 notch = 1 - notchdips*exp(-qradsq/2/notchwidthxy^2);
 Mask = ones(Nx,Ny,Nz);
 Mask(abs(dataparams.OTFem)<0.001)=0;     % notch mask
 notch = notch.*Mask;
 notch = notch./max(max(max(notch)));     % normalize
+OTFtemp(OTFtemp>0.4)=0.4;
 OTFtemp1 = OTFtemp.*notch;
 clear dataparams qx qy qz XSupport YSupport ZSupport notch qradsq
 
@@ -39,13 +48,17 @@
             ftshiftorderims(:,:,2,jNz) = double(shift(squeeze(tempimage(:,:,2,jNz)),[m*qvector(1)/2,m*qvector(2)/2]));
             OTFshift0(:,:,2,jNz) = 0.5*double(shift(squeeze(OTFtemp1(:,:,jNz)),[m*qvector(1)/2,m*qvector(2)/2]));
         end
-        OTFshift0(:,:,2,:) = double(shift(squeeze(OTFshift0(:,:,2,:)),[0,0,qvector(3)])) + double(shift(squeeze(OTFshift0(:,:,2,:)),[0,0,-qvector(3)]));  
+        if Nz~=1
+            OTFshift0(:,:,2,:) = double(shift(squeeze(OTFshift0(:,:,2,:)),[0,0,qvector(3)])) + double(shift(squeeze(OTFshift0(:,:,2,:)),[0,0,-qvector(3)])); 
+        end
     elseif jorder == 3 % -1th frequency
         for jNz = 1:Nz
             ftshiftorderims(:,:,3,jNz) = double(shift(squeeze(tempimage(:,:,3,jNz)),[-m*qvector(1)/2,-m*qvector(2)/2])); 
             OTFshift0(:,:,3,jNz) = 0.5*double(shift(squeeze(OTFtemp1(:,:,jNz)),[-m*qvector(1)/2,-m*qvector(2)/2]));
         end
-        OTFshift0(:,:,3,:) = double(shift(squeeze(OTFshift0(:,:,3,:)),[0,0,qvector(3)])) + double(shift(squeeze(OTFshift0(:,:,3,:)),[0,0,-qvector(3)]));
+        if Nz~=1
+            OTFshift0(:,:,3,:) = double(shift(squeeze(OTFshift0(:,:,3,:)),[0,0,qvector(3)])) + double(shift(squeeze(OTFshift0(:,:,3,:)),[0,0,-qvector(3)]));
+        end
     elseif jorder == 4 % 2th frequency
         for jNz = 1:Nz
             ftshiftorderims(:,:,4,jNz) = double(shift(squeeze(tempimage(:,:,4,jNz)),[m*qvector(1),m*qvector(2)]));

MATLAB_Open_3DSIM/help_functions/get_calibrationOTF.m

@@ -1,4 +1,5 @@
 function [OTF_orders,OTFinc] = get_calibrationOTF(allfilenamesOTFdata,SIMparams,numpixelsx,numpixelsy,rawpixelsize)
+
 % copyright Sjoerd Stallinga, TU Delft, 2017-2020
 
 maxorder = 3;

MATLAB_Open_3DSIM/help_functions/get_filter2.m

@@ -91,6 +91,17 @@
 Apo = 1-qrad./cutoffmap;
 Apo = Apo + abs(min(min(Apo)));
 
-Filter2 = Apo./(params.OTFshiftfinal+lambdaregul); 
+
+if Nz ==1
+    temp = params.OTFshiftfinal;
+    max_ = max(max(params.OTFshiftfinal));
+    temp(temp<0.5*max_) = 0.4*max_;
+    Filter2=Apo./(temp+lambdaregul); 
+else
+    Filter2 = Apo./(params.OTFshiftfinal+lambdaregul);
+end
 Filter2(isnan(Filter2)) = epsy;
+
+
+
 end

MATLAB_Open_3DSIM/help_functions/get_modelOTF.m

@@ -1,4 +1,7 @@
 function [OTF_orders,OTF_every_slide] = get_modelOTF(dataparams,SIMparams)
+
+% copyright Sjoerd Stallinga, TU Delft, 2017-2020
+
 %% Basic parameters
 Nx = SIMparams.numSIMpixelsx;
 Ny = SIMparams.numSIMpixelsy;

MATLAB_Open_3DSIM/help_functions/get_vectormodelOTF.m

@@ -1,7 +1,5 @@
 function [OTFinc,OTFinc2d_throughfocus] = get_vectormodelOTF(OTFparams)
-
-% copyright Sjoerd Stallinga TUD 2017-2020
-
+
 %% Get pupil matrix
 [~,~,wavevector,wavevectorzmed,~,PupilMatrix] = get_pupil_matrix(OTFparams);
 
@@ -30,4 +28,4 @@
 end
 OTFinc = OTFinc/OTFnorm;
 
-end
+end

MATLAB_Open_3DSIM/help_functions/notchfilter.m

@@ -10,6 +10,10 @@
 qvector(3) = 1.0*double(1.0*kz*Nz*dataparams.SIMpixelsize(3));
 qvector_z = qvector(3);
 qvector(1:2) = double(1.0*patternpitch);
+% single-layer
+if dataparams.Nz==1
+    qvector(3)=0;
+end
 
 %% Calculate notch-filter
 XSupport = ((1:Nx)-floor(Nx/2)-1);         % qx grid
@@ -18,11 +22,17 @@
 [qx,qy,qz] = meshgrid(XSupport,YSupport,ZSupport); 
 notchwidthxy = dataparams.notchwidthxy2;   % notch width
 notchdips = dataparams.notchdips2;         % notch depth
-qradsq = ( qx/qvector_xy ).^2+( qy/qvector_xy ).^2+( qz/qvector_z ).^2;
+if Nz ~=1
+    qradsq = ( qx/qvector_xy ).^2+( qy/qvector_xy ).^2+( qz/qvector_z ).^2;
+else
+    qradsq = ( qx/qvector_xy/2 ).^2+( qy/qvector_xy/2 ).^2+( qz/qvector_z/2 ).^2;
+end
 notch = 1 - notchdips*exp(-qradsq/2/notchwidthxy^2);
 Mask = ones(Nx,Ny,Nz);
 Mask(abs(dataparams.OTFem)<0.01)=0;
-notch = notch.*Mask.*abs(dataparams.OTFem).^(attenuation)./max(max(max(abs(dataparams.OTFem))));
+OTFtemp = dataparams.OTFem;
+OTFtemp(OTFtemp>0.4)=0.4;
+notch = notch.*Mask.*abs(OTFtemp).^(attenuation);
 notch = notch./max(max(max(notch)));
 
 %% shift the separated frequency of 0,¡À1,¡À2 and do the notch operation
@@ -36,14 +46,18 @@
         for jNz = 1:Nz
             anneuation(:,:,jNz) = double(shift(squeeze(notch(:,:,jNz)),[m*qvector(1)/2,m*qvector(2)/2]));
         end
-        anneuation(:,:,:) =double(shift(squeeze(anneuation(:,:,:)),[0,0,qvector(3)])) + double(shift(squeeze(anneuation(:,:,:)),[0,0,-qvector(3)]));
+        if Nz~=1
+            anneuation(:,:,:) =double(shift(squeeze(anneuation(:,:,:)),[0,0,qvector(3)])) + double(shift(squeeze(anneuation(:,:,:)),[0,0,-qvector(3)]));
+        end
         anneuation(anneuation<0) = 0;
         ftshiftorderims(:,:,2,:) = squeeze(tempimage(:,:,2,:)).*anneuation(:,:,:)/2;
     elseif jorder == 3  % -1th frequency
         for jNz = 1:Nz
             anneuation(:,:,jNz) = double(shift(squeeze(notch(:,:,jNz)),[-m*qvector(1)/2,-m*qvector(2)/2]));
         end
-        anneuation(:,:,:) = double(shift(squeeze(anneuation(:,:,:)),[0,0,qvector(3)])) + double(shift(squeeze(anneuation(:,:,:)),[0,0,-qvector(3)]));
+        if Nz~=1
+            anneuation(:,:,:) = double(shift(squeeze(anneuation(:,:,:)),[0,0,qvector(3)])) + double(shift(squeeze(anneuation(:,:,:)),[0,0,-qvector(3)]));
+        end
         anneuation(anneuation<0) = 0;
         ftshiftorderims(:,:,3,:) = squeeze(tempimage(:,:,3,:)).*anneuation(:,:,:)/2;
     elseif jorder == 4  % 2th frequency

MATLAB_Open_3DSIM/help_functions/om_display.m

@@ -1,5 +1,8 @@
 function [wf, ld_om, cm, ouf, fig] = om_display(ld, cmin, cmax, is_display)
-%%
+
+% copyright Karl Zhanghao, 2019
+
+
 ld_s = mean(ld,3); 
 ld_s = max(ld_s-cmin, 0); 
 ld_s = min(ld_s/cmax, 1);

MATLAB_Open_3DSIM/help_functions/process_data.m

@@ -50,14 +50,6 @@
     Snoisy(:,:,(jangle-1)*5+1:(jangle-1)*5+5) = tempimage;
 end
 
-%%%%%%%%%%%%%%%%%%%%%
-% maxnum = max(max(max(Snoisy)));
-% Snoisy = Snoisy./maxnum;
-% stackfilename = ['Snoisy.tif'];
-% for k = 1:15
-%     imwrite(Snoisy(:,:,k), stackfilename, 'WriteMode','append') % дÈëstackͼÏñ
-% end
-%%%%%%%%%%%%%%%%%%%%%
 %% Get the frequency, phase and modulation depth
 [dataparams,~] = find_illumination_pattern(Snoisy,dataparams);
 freq = dataparams.allpatternpitch;
@@ -98,7 +90,6 @@
 dataparams.allmoduleave = allmodule;
 clear allmodule
 
-%% Get OTFem
 SIMparams.upsampling = [2 2 1];                     % interperation
 numSIMpixelsx = SIMparams.upsampling(1)*numpixelsx; % X pixels of final SIM image
 numSIMpixelsy = SIMparams.upsampling(2)*numpixelsy; % y pixels of final SIM image
@@ -109,11 +100,17 @@
 SIMparams.numSIMpixelsx = numSIMpixelsx;
 SIMparams.numSIMpixelsy = numSIMpixelsy;
 SIMparams.numSIMpixelsz = numSIMpixelsz;
+
+%% Get OTFem
+if dataparams.Nz~=1
 if OTFflag == 1
     [OTFem,~] = get_modelOTF(dataparams,SIMparams);
 elseif OTFflag == 0
     [~,OTFem]= get_calibrationOTF(OTF_name,SIMparams,numpixelsx,numpixelsy,dataparams.rawpixelsize);
 end
+else
+    OTFem = SimOtfProvider2( dataparams,SIMparams,NA,emwavelength,numSIMpixelsx,1);
+end
 dataparams.OTFem = OTFem;
 
 %% edge taper

MATLAB_Open_3DSIM/help_functions/psim_sim3d_recon.m

@@ -1,11 +1,7 @@
 function [ld, psim] = psim_sim3d_recon(raw_sim3d, sim_s, theta, calib)
-%% calib
-% if size(calib{1},1) == 1024
-%     fov = size(raw_sim3d,1);
-%     calib1 = calib{1}; calib2 = calib{2};
-%     calib{1} = calib1((513-fov/2):(512+fov/2),(513-fov/2):(512+fov/2)); 
-%     calib{2} = calib2((513-fov/2):(512+fov/2),(513-fov/2):(512+fov/2));
-% end
+
+% copyright Karl Zhanghao, 2019
+
 %% Calculate frequency components on orientational dimension
 % obtain polarized wide field images
 d1_s = mean(raw_sim3d(:,:,1:5),3); 

MATLAB_Open_3DSIM/output/Readme.txt

@@ -1 +1 @@
-The output images will be saved here
+The output images will be saved here.