OpenSIM Matlab code

OpenSIM Matlab code

SIMbasic/ApodizationFunction.m

@@ -0,0 +1,30 @@
+function [FsumA] = ApodizationFunction(Fsum,k2fa,k2fb,k2fc,Kotf,Index)
+
+t = size(Fsum,1);
+
+OTFo_mask = OTFmaskShifted(0.*k2fa,Kotf,t);
+OTFap_mask = OTFmaskShifted(k2fa,Kotf,t);
+OTFam_mask = OTFmaskShifted(-k2fa,Kotf,t);
+OTFbp_mask = OTFmaskShifted(k2fb,Kotf,t);
+OTFbm_mask = OTFmaskShifted(-k2fb,Kotf,t);
+OTFcp_mask = OTFmaskShifted(k2fc,Kotf,t);
+OTFcm_mask = OTFmaskShifted(-k2fc,Kotf,t);
+ApoMask = OTFo_mask.*OTFap_mask.*OTFam_mask.*OTFbp_mask.*OTFbm_mask.*OTFcp_mask.*OTFcm_mask;
+clear OTFoz OTFo_mask OTFap_mask OTFam_mask OTFbp_mask OTFbm_mask OTFcp_mask OTFcm_mask
+%{
+figure;
+imshow(ApoMask,[ ])
+%}
+
+DistApoMask = bwdist(ApoMask);
+maxApoMask = max(max(DistApoMask));
+ApoFunc = double(DistApoMask./maxApoMask).^Index;
+%{
+figure;
+imshow(DistApoMask,[ ])
+figure;
+mesh(double(ApoFunc))
+%}
+
+FsumA = Fsum.*ApoFunc;
+clear ApoMask DistApoMask ApoFunc

SIMbasic/ApproxFreqDuplex.m

@@ -0,0 +1,30 @@
+function [maxK2,Ix,Iy] = ApproxFreqDuplex(FiSMap,Kotf)
+% AIM: approx. illumination frequency vector determination
+% INPUT VARIABLES
+%   FiSMap: FT of raw SIM image
+%   Kotf: OTF cut-off frequency
+% OUTPUT VARIABLES
+%   maxK2: illumination frequency vector (approx)
+%   Ix,Iy: coordinates of illumination frequency peaks
+
+FiSMap = abs(FiSMap);
+
+w = size(FiSMap,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+
+Ro = sqrt( (X-wo).^2 + (Y-wo).^2 );
+Z0 = Ro > round(0.5*Kotf);
+Z1 = X > wo;
+
+FiSMap = FiSMap.*Z0.*Z1;
+dumY = max( FiSMap,[],1 );
+[dummy Iy] = max(dumY);
+dumX = max( FiSMap,[],2 );
+[dummy Ix] = max(dumX);
+
+maxK2 = [Ix-(wo+1) Iy-(wo+1)];
+
+

SIMbasic/IlluminationFreqF.m

@@ -0,0 +1,38 @@
+function [k2fa] = IlluminationFreqF(S1aTnoisy,OTFo)
+% AIM: illumination frequency vector determination
+% INPUT VARIABLES
+%   S1aTnoisy: raw SIM image
+%   OTFo: system OTF
+% OUTPUT VARIABLE
+%   k2fa: illumination frequency vector
+
+w = size(OTFo,1);
+wo = w/2;
+
+% computing PSFe for edge tapering SIM images
+PSFd = real(fftshift( ifft2(fftshift(OTFo.^10)) ));
+PSFd = PSFd/max(max(PSFd));
+PSFd = PSFd/sum(sum(PSFd));
+h = 30;
+PSFe = PSFd(wo-h+1:wo+h,wo-h+1:wo+h);
+
+% edge tapering raw SIM image
+S1aTnoisy_et = edgetaper(S1aTnoisy,PSFe);
+fS1aTnoisy_et = fftshift(fft2(S1aTnoisy_et));
+
+% OTF cut-off freq
+Kotf = OTFedgeF(OTFo);
+
+% Approx illumination frequency vector
+[k2fa,~,~] = ApproxFreqDuplex(fS1aTnoisy_et,Kotf);
+
+fS1aTnoisy = fftshift(fft2(S1aTnoisy));
+% illumination frequency vector determination by optimizing
+% autocorrelation of fS1aTnoisy
+OPT = 1;
+PhaseKai2opt0 = @(k2fa0)PhaseKai2opt(k2fa0,fS1aTnoisy,OTFo,OPT);
+options = optimset('LargeScale','off','Algorithm',...
+	'active-set','MaxFunEvals',500,'MaxIter',500,'Display','notify');
+k2fa0 = k2fa;
+[k2fa,fval] = fminsearch(PhaseKai2opt0,k2fa0,options);
+% k2a = sqrt(k2fa*k2fa')

SIMbasic/IlluminationPhaseF.m

@@ -0,0 +1,14 @@
+function [phaseA1] = IlluminationPhaseF(S1aTnoisy,k2fa)
+% AIM: illumination phase shift determination
+% INPUT VARIABLES
+%   S1aTnoisy: raw SIM image
+%   k2fa: illumination frequency vector
+% OUTPUT VARIABLE
+%   phaseA1: illumination phase shift determined
+
+PatternPhaseOpt0 = @(phaseA0)PatternPhaseOpt(phaseA0,S1aTnoisy,k2fa);
+options = optimset('LargeScale','off','Algorithm',...
+	'active-set','MaxFunEvals',500,'MaxIter',500,'Display','notify');
+phaseA0 = 0;
+[phaseA1,fval] = fminsearch(PatternPhaseOpt0,phaseA0,options);
+% phaseA1*180/pi

SIMbasic/MergingHeptaletsF.m

@@ -0,0 +1,68 @@
+function [Fsum,Fperi,Fcent] = MergingHeptaletsF(fDIo,fDIp,fDIm,fDBo,fDBp,fDBm,...
+    fDCo,fDCp,fDCm,Ma,Mb,Mc,npDIo,npDIp,npDIm,...
+    npDBo,npDBp,npDBm,npDCo,npDCp,npDCm,k2fa,k2fb,k2fc,OBJpara,OTFo)
+
+% AIM: To merge all 9 frequency components into one using generalized
+%       Wiener Filter
+% INPUT VARIABLES
+%   [fDIo fDIp fDIm
+%    fDBo fDBp fDBm
+%    fDCo fDCp fDCm]: nine frequency components
+%   Ma,Mb,Mc: modulation factors for the three illumination orientations
+%   [npDIo npDIp npDIm 
+%    npDBo npDBp npDBm 
+%    npDCo,npDCp,npDCm]: noise powers corresponding to nine frequency
+%                       components
+%   k2fa,k2fb,k2fc: illumination frequency vectors for the three
+%                   illumination orientations
+%   OBJpara: Object spectrum parameters
+%   OTFo: system OTF
+% OUTPUT VARIABLES
+%   Fsum: all nine frequency components merged into one using 
+%           generalised Wiener Filter
+%   Fperi: six off-center frequency components merged into one using 
+%           generalised Wiener Filter
+%   Fcent: averaged of the three central frequency components
+
+
+% obtain signal spectrums corresponding to central and off-center
+% frequency components
+[SIGao,SIGam2,SIGap2] = TripletSNR0(OBJpara,k2fa,OTFo,fDIm);
+[SIGbo,SIGbm2,SIGbp2] = TripletSNR0(OBJpara,k2fb,OTFo,fDBm);
+[SIGco,SIGcm2,SIGcp2] = TripletSNR0(OBJpara,k2fc,OTFo,fDCm);
+
+SIGap2 = Ma*SIGap2;
+SIGam2 = Ma*SIGam2;
+SIGbp2 = Mb*SIGbp2;
+SIGbm2 = Mb*SIGbm2;
+SIGcp2 = Mc*SIGcp2;
+SIGcm2 = Mc*SIGcm2;
+
+%% Generalized Wiener-Filter computation
+SNRao = SIGao.*conj(SIGao)./npDIo;
+SNRap = SIGap2.*conj(SIGap2)./npDIp;
+SNRam = SIGam2.*conj(SIGam2)./npDIm;
+
+SNRbo = SIGbo.*conj(SIGbo)./npDBo;
+SNRbp = SIGbp2.*conj(SIGbp2)./npDBp;
+SNRbm = SIGbm2.*conj(SIGbm2)./npDBm;
+
+SNRco = SIGco.*conj(SIGco)./npDCo;
+SNRcp = SIGcp2.*conj(SIGcp2)./npDCp;
+SNRcm = SIGcm2.*conj(SIGcm2)./npDCm;
+
+ComDeno = 0.01 + ( SNRao + SNRap + SNRam + SNRbo + SNRbp + SNRbm + SNRco + SNRcp + SNRcm );
+Fsum = fDIo.*SNRao + fDIp.*SNRap + fDIm.*SNRam...
+    + fDBo.*SNRbo + fDBp.*SNRbp + fDBm.*SNRbm...
+    + fDCo.*SNRco + fDCp.*SNRcp + fDCm.*SNRcm;
+Fsum = Fsum./ComDeno;
+
+ComPeri = 0.01 + ( SNRap + SNRam + SNRbp + SNRbm + SNRcp + SNRcm );
+Fperi = fDIp.*SNRap + fDIm.*SNRam + fDBp.*SNRbp + fDBm.*SNRbm...
+    + fDCp.*SNRcp + fDCm.*SNRcm;
+Fperi = Fperi./ComPeri;
+
+% averaged central frequency component
+Fcent = (fDIo+fDBo+fDCo)/3;
+
+

SIMbasic/ModulationFactor.m

@@ -0,0 +1,86 @@
+function Mm = ModulationFactor(fDp,k1a,OBJparaA,OTFo)
+% AIM: Determination of modulation factor
+% INPUT VARIABLES
+%   fDp: off-center frequency component
+%   k1a: illumination frequency vector
+%   OBJparaA: Object power parameters
+%   OTFo: system OTF
+% OUTPUT VARIABLE
+%   Mm: modulation factor
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+% magnitude of illumination vector
+k2 = sqrt(k1a*k1a');
+
+% vector along illumination direction
+kv = k1a(2) + 1i*k1a(1); 
+Rp = abs(Cv+kv);
+
+% Object power parameters
+Aobj = OBJparaA(1);
+Bobj = OBJparaA(2);
+
+% Object spectrum
+OBJp = Aobj*(Rp+0).^Bobj;
+
+% illumination vector rounded to nearest pixel
+k3 = -round(k1a);
+
+OBJp(wo+1+k3(1),wo+1+k3(2)) = 0.25*OBJp(wo+2+k3(1),wo+1+k3(2))...
+	+ 0.25*OBJp(wo+1+k3(1),wo+2+k3(2))...
+	+ 0.25*OBJp(wo+0+k3(1),wo+1+k3(2))...
+	+ 0.25*OBJp(wo+1+k3(1),wo+0+k3(2));
+
+% signal spectrum
+SIGap = OBJp.*OTFo;
+
+% figure;
+% mesh(log(abs(OBJp)))
+% figure;
+% mesh(log(abs(fDp)))
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% frequency beyond which NoisePower estimate to be computed
+NoiseFreq = Kotf + 20;
+
+% NoisePower determination
+Zo = Ro>NoiseFreq;
+nNoise = fDp.*Zo;
+NoisePower = sum(sum( nNoise.*conj(nNoise) ))./sum(sum(Zo));
+
+% Noise free object power computation
+Fpower = fDp.*conj(fDp) - NoisePower;
+fDp = sqrt(abs(Fpower));
+
+% frequency range over which signal power matching is done to estimate
+% modulation factor
+Zmask = (Ro > 0.2*k2).*(Ro < 0.8*k2).*(Rp > 0.2*k2);
+
+% least square approximation for modulation factor
+Mm = sum(sum(SIGap.*abs(fDp).*Zmask));
+Mm = Mm./sum(sum(SIGap.^2.*Zmask))
+
+%% for visual verification
+%{
+pp = 3;
+figure;
+hold on
+plot(x-wo,log(abs(fDp(wo+pp,:))),'o-','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(SIGap(wo+pp,:))),'go--','LineWidth',2,'MarkerSize',4)
+plot(x-wo,log(abs(Mm*SIGap(wo+pp,:))),'r*--','LineWidth',2,'MarkerSize',4)
+grid on
+box on
+kk
+%}
+
+
+

SIMbasic/OBJparaOpt.m

@@ -0,0 +1,47 @@
+function Esum = OBJparaOpt(OBJpara0,fDIoTnoisy,OTFo)
+% AIM: Determined Sum of Squared Errors (SSE) between `actual signal power' and
+%   `approximated signal power'
+% INPUT VARIABLES
+%   OBJpara0: [Aobj Bobj], Object power parameters
+%   fDIoTnoisy: FT of central frequency component
+%   OTFo: system OTF
+% OUTPUT VARIABLE
+%   Esum: SSE between `actual signal spectrum' and `approximated signal spectrum'
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+Ro(wo+1,wo+1) = 1; % to avoid nan
+
+% approximated signal power calculation
+Aobj = OBJpara0(1);
+Bobj = OBJpara0(2);
+OBJpower = Aobj*(Ro.^Bobj);
+SIGpower = OBJpower.*OTFo;
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% range of frequency over which SSE is computed
+Zloop = (Ro<0.75*Kotf).*(Ro>0.25*Kotf);
+
+% frequency beyond which NoisePower estimate to be computed
+NoiseFreq = Kotf + 20;
+
+% NoisePower determination
+Zo = Ro>NoiseFreq;
+nNoise = fDIoTnoisy.*Zo;
+NoisePower = sum(sum( nNoise.*conj(nNoise) ))./sum(sum(Zo));
+
+% Noise free object power computation 
+Fpower = fDIoTnoisy.*conj(fDIoTnoisy) - NoisePower;
+fDIoTnoisy = sqrt(abs(Fpower));
+
+% SSE computation
+Error = fDIoTnoisy - SIGpower;
+Esum = sum(sum((Error.^2./Ro).*Zloop));
+

SIMbasic/OBJpowerPara.m

@@ -0,0 +1,48 @@
+function OBJparaA = OBJpowerPara(fDIoTnoisy,OTFo)
+% AIM: determination of object power parameters Aobj and Bobj
+% INPUT VARIABLES
+%   fDIoTnoisy: FT of central frequency component
+%   OTFo: system OTFo
+% OUTPUT VARIABLES
+%   OBJparaA: [Aobj Bobj]
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+%% object power parameters through optimization
+OBJparaOpt0 = @(OBJpara0)OBJparaOpt(OBJpara0,fDIoTnoisy,OTFo);
+options = optimset('LargeScale','off','Algorithm',...
+	'active-set','MaxFunEvals',500,'MaxIter',500,'Display','notify');
+
+% obtaining crude initial guesses for Aobj and Bobj 
+Zm = (Ro>0.3*Kotf).*(Ro<0.4*Kotf);
+Aobj = sum(sum(abs(fDIoTnoisy.*Zm)))./sum(sum(Zm));
+Bobj = -0.5;
+OBJpara0 = [Aobj Bobj];
+
+% optimization step
+[OBJparaA,fval] = fminsearch(OBJparaOpt0,OBJpara0,options);
+
+%% plot for cross-checking the result
+%{
+Aobj = OBJparaA(1);
+Bobj = OBJparaA(2);
+OBJo = Aobj*(Ro.^Bobj);
+SIGp = OBJo.*OTFo;
+pp = 3;
+figure;
+hold on
+plot([0:w-1]-wo,log(abs(fDIoTnoisy(wo+pp,:))),'k--','LineWidth',3,'MarkerSize',6)
+plot([0:w-1]-wo,log(abs(SIGp(wo+pp,:))),'mo-','LineWidth',2,'MarkerSize',6)
+legend('acutal signal power','avg. signal power')
+grid on
+box on
+%}

SIMbasic/OTFdoubling.m

@@ -0,0 +1,17 @@
+function OTF2 = OTFdoubling(OTFo,DoubleMatSize)
+
+%% embeds OTFo in doubled range of frequencies  
+
+w = size(OTFo,1);
+wo = w/2;
+if ( DoubleMatSize>0 )
+    t = 2*w;
+    u = linspace(0,t-1,t);
+    v = linspace(0,t-1,t);
+    OTF2 = zeros(2*w,2*w);
+    OTF2(wo+1:w+wo,wo+1:w+wo) = OTFo;
+    clear OTFo
+else
+    OTF2 = OTFo;
+    clear OTFo
+end