unzipped

Spectral-FLIM GUI/Spectral-FLIM GUI/Calc_mIRF.m

@@ -0,0 +1,86 @@
+function IRF = Calc_mIRF(head, tcspc)
+
+maxres = max([head.Resolution]);
+Resolution = max([maxres 0.032]);
+Pulse      = 1e9/head.SyncRate;
+
+tau = Resolution.*((1:size(tcspc,2))-0.5)';
+IRF = zeros(size(tcspc));
+nex = 2;
+
+t0 = 1.5;
+w1 = 0.5^2;  % width of excitation peak
+T1 = 0.10;    % time constant of rise term in IRF
+T2 = 0.10;    % time constant of decay term in IRF
+a  = 0.1;     % rel. amplitude of second peak in IRF
+b  = 0.1;     % rel. amplitude of second peak in IRF
+dt = 0.0;     % time shift of second peak in IRF
+
+for PIE = 1:size(tcspc,3)
+
+    p   = [t0     w1     T1   T2   a  b   dt    1 3]';
+    pl  = [t0-1.0 1e-3 1e-4 1e-4   0  0  -0.3   0.5*ones(1, nex)]';
+    pu  = [t0+1.0 1e-1  1    1      1  1   0.5  5*ones(1, nex)]';
+
+    tc = squeeze(sum(tcspc(:,:,PIE),2));
+    [tmp, ord] = sort(tc,'descend');
+
+    ch = 1;
+    ind = ord(ch);
+    y = squeeze(tcspc(ind,:,PIE));
+
+    err = zeros(1,10);
+
+    for casc=1:10
+        [ts, s] = min(err);
+        r0 = p(:, s);
+        for sub=1:10
+            rf = r0.*[2.^(1.1*(rand(size(r0))-0.5)./casc)];  % randomize start values
+            rf = max([rf pl],[],2);
+            rf = min([rf pu],[],2);
+            p(:,sub) = Simplex('TCSPC_Fun',rf,pl,pu,[],[],tau, y);
+            err(sub) = TCSPC_Fun(p(:,sub), tau, y);
+        end
+    end
+
+    err1 = min(err);
+    p1   = mean(p(:,err==err1),2);
+    [tmp, c1, bla1, tmp1] = TCSPC_Fun(p1, tau, y);
+
+    IRF(ind,:,PIE) = IRF_Fun(p1(1:7),tau);
+    IRF(ind,:,PIE) = IRF(ind,:,PIE)./max(IRF(ind,:,PIE));
+
+    para = p1(2:7);
+    p    = [p1(1);  p1(8:end)];
+    pl   = [  0.5; 0.5*ones(nex, 1)];
+    pu   = [  3.0;  10*ones(nex, 1)];
+
+    for ch = 2:size(tcspc,1)
+        ind = ord(ch);
+        y = squeeze(tcspc(ind,:,PIE));
+
+        err = zeros(1,10);
+
+        for casc=1:10
+            [ts, s] = min(err);
+            r0 = p(:, s);
+            for sub=1:10
+                rf = r0.*[2.^(1.05*(rand(size(r0))-0.5)./casc)];  % randomize start values
+                rf = max([rf pl],[],2);
+                rf = min([rf pu],[],2);
+                p(:,sub) = Simplex('TCSPC_Fun',rf,pl,pu,[],[],tau, y, para);
+                err(sub) = TCSPC_Fun(p(:,sub), tau, y, para);
+            end
+        end
+
+        err1 = min(err);
+        p1   = mean(p(:,err==err1),2);
+        [tmp, c1, bla1, tmp1] = TCSPC_Fun(p1, tau, y, para);
+
+        IRF(ind,:,PIE) = IRF_Fun([p1(1); para; p1(2:end)],tau);
+        IRF(ind,:,PIE) = IRF(ind,:,PIE)./max(IRF(ind,:,PIE));
+    end;    
+end
+
+IRF(IRF<1e-4) = 0;
+IRF = IRF./repmat(sum(IRF,2),[1 size(IRF,2) 1]);

Spectral-FLIM GUI/Spectral-FLIM GUI/Convergence.m

@@ -0,0 +1,58 @@
+function [x_out,llh_out,n_iter,conv,b]  = Convergence(I,K)
+tic;
+i_err = (min(K,[],2)>0);
+K = K(i_err,:);
+I = I(:,i_err);
+
+lI = log(I)-1;
+
+n  = size(I,1);
+k  = size(K,2);
+
+hd = repmat(sum(K,1),[n 1]);
+
+X = repmat(sum(I,2),[1 k])./k/5;
+Exitflag = false;
+
+me = zeros(1,20);
+conv=[];
+n_iter=0;
+update = 1.2;
+while ~Exitflag
+
+    rate = update;
+
+    for sl = 1:20
+
+        Z  = X*K';
+
+
+        x  = (X./hd).*((I./Z)*K);
+        dx = (x - X);
+
+        X  = X + rate.*dx;
+
+        me(sl) = 25*max(sum(abs(dx),2));
+    end
+
+    if min(me) < 50
+        update = 1.6;
+    end
+
+    if min(me) < 1
+        Exitflag = true;
+    else
+        Exitflag = false;
+    end
+
+    n_iter=n_iter+20;
+    conv=[conv me];  %#ok<AGROW>
+
+end
+
+m   = Z.*(log(Z)-lI);
+llh = sum(m,2);
+x_out   = x;
+llh_out = llh;
+b=toc;
+end

Spectral-FLIM GUI/Spectral-FLIM GUI/Convol.m

@@ -0,0 +1,20 @@
+function y = Convol(irf, x)
+% convol(irf, x) performs a convolution of the instrumental response 
+% function irf with the decay function x. Periodicity (=length(x)) is assumed.
+
+mm = mean(irf(end-10:end));
+if size(x,1)==1 || size(x,2)==1
+    irf = irf(:);
+    x = x(:);
+end
+p = size(x,1);
+n = numel(irf);
+if p>n
+    irf = [irf; mm*ones(p-n,1)]; 
+else
+    irf = irf(1:p);
+end
+y = real(ifft((fft(irf)*ones(1,size(x,2))).*fft(x)));
+t = rem(rem(0:n-1,p)+p,p)+1;
+y = y(t,:);
+

Spectral-FLIM GUI/Spectral-FLIM GUI/DetectTimeGates.m

@@ -0,0 +1,161 @@
+function [timegate, Ngate, tcspc] = DetectTimeGates(tcspcdata, cnum, Resolution, pic)
+%. 
+% 
+%  DetectTimeGates divides the TCSPC histogramm into time-windows of the
+%  same length in order to separate the data from different excitation
+%  pulses.
+%
+%  usage: [timegate, Ngate, tcspc] = DetectTimeGates(tcspcdata, cnum, Resolution, pic)
+%
+%    tcspcdata  : TCSPC histogram of a PIE experiment
+%    cnum       : number of excitation pulses (default: 1)
+%    Resolution : bin-size of TSCPC data in ns (default: 0.002 ns)
+%    pic        : flag to display aligned data
+%
+%  DetectTimeGates returns the following data:
+%
+%    timegate   : array of the size [cnum*dnum, 4] containing the
+%                 respective beginning and end of each time-window. 
+%
+%                 The order of the list is 
+%
+%                 timegate(  1   ,:): time-window for pulse  1   in detection channel 1
+%                 timegate(  1   ,:): time-window for pulse  1   in detection channel 2
+%            ...  timegate(  n   ,:): time-window for pulse  1   in detection channel n
+%                 timegate( n+1  ,:): time-window for pulse  2   in detection channel 1 
+%                 timegate( n+2  ,:): time-window for pulse  2   in detection channel 2
+%            ...  timegate( 2*n  ,:): time-window for pulse  2   in detection channel n
+%            ...  timegate(cnum*n,:): time-window for pulse cnum in detection channel n
+%
+%                 timegate(:, 1)    : begin of time-window 
+%                 timegate(:, 2)    : end of time-window 
+%                 timegate(:, 3)    : if > 0: continuation of time-window
+%                 timegate(:, 4)    : if > 0: end of time-window 
+%
+%    Ngate      : Number of bins in each time-window
+%    tcscpc     : Sorted and aligned TCSPC histogram in the same order as the timegates
+%
+% .........................................................................................
+
+if (nargin < 4)||isempty(pic)
+    pic = 0;
+end
+pic = (pic>0);
+
+if (nargin < 3)||isempty(Resolution)
+    Resolution = 0.002;
+end
+
+if (nargin < 2)||isempty(cnum)
+    cnum = 1;
+end
+
+NChannels = size(tcspcdata,1);
+dnum      = size(tcspcdata,2);
+timegate  = zeros(dnum*cnum,4); 
+
+%% Determine maxima of all detectors
+
+tmpdata = zeros(size(tcspcdata));
+
+for k = 1:dnum
+    tmpdata(:,k) = smooth(tcspcdata(:,k),ceil(0.1./Resolution));
+end
+
+irf_w = floor(1.5/Resolution);
+win   = floor(NChannels/cnum);     % Max #channels / pulse;
+
+im = zeros(dnum, cnum);
+m  = zeros(dnum, cnum);
+bg = zeros(dnum, cnum);
+
+[tpm,i1] = max(tmpdata); %#ok<ASGLU>
+off      = floor(mod(i1,win)-irf_w);
+
+for l = 1:cnum
+    for k = 1:dnum
+        ind  = off(k)+(l-1)*win+(1:win-2*irf_w);
+        ind(ind<1)         = ind(ind<1) + NChannels;
+        ind(ind>NChannels) = ind(ind>NChannels) - NChannels;
+        [m(k,l),im(k,l)] = max(tmpdata(ind,k));
+        im(k,l) = ind(im(k,l));
+        bg(k,l) = mean(tmpdata(ind(2*irf_w:end),k));
+    end
+end
+
+ind = (m<3*bg)|(m<10);
+
+for l = 1:cnum
+    for k = 1:dnum
+        if ind(k,l) == 1
+            mi = find(im(k,:) == i1(k));
+            im(k,l) = im(k,mi)+floor((l-mi)*NChannels/cnum);
+        end
+    end
+end
+
+%% Determine time-gates for photon assignment
+if cnum>1
+    Ngate = zeros(1,cnum);
+    for k = 1:cnum
+        p = k-1;
+        if p < 1
+            p = cnum;
+        end
+        Ngate(k) = min([abs(im(:,k)-im(:,p)); NChannels-abs(im(:,p)-im(:,k))]);
+    end
+    Ngate = min(Ngate);
+else
+    Ngate = NChannels;
+end
+
+%% Define final windows for each pulse and detector
+
+g = zeros(cnum*dnum, 4);
+
+for k = 1:dnum
+    for l = 1:cnum
+        d = im(k,l)-irf_w;
+        if d > 0
+            if d <= (NChannels - Ngate)+1
+                g(k+(l-1)*dnum,1:2) = [d d+Ngate-1];
+            else
+                g(k+(l-1)*dnum,:) = [d NChannels 1 d+Ngate-NChannels-1];
+            end;
+        else
+            d = d + NChannels;
+            if d > NChannels - Ngate +1
+                g(k+(l-1)*dnum,:) = [d NChannels 1 d+Ngate-NChannels-1];
+            else
+                g(k+(l-1)*dnum,1:2) = [d d+Ngate-1];
+            end;
+        end;
+    end
+end
+
+%% Sort windows to according excitation pulse
+
+channels = repmat(1:dnum,[1 cnum]);
+
+for n = 1:cnum*dnum
+    mch = mod(mean(g(n,1):(g(n,2)+g(n,4))),NChannels);
+    pie = floor(mch/(NChannels/(cnum-0.5)));
+    timegate(channels(n)+pie*dnum,:) = g(n,:);
+end
+
+tcspc = zeros(Ngate,cnum*dnum);
+
+for n = 1:cnum*dnum
+    if (timegate(n,3)==0)
+        tmp =  tcspcdata(timegate(n,1):timegate(n,2),channels(n));
+    else
+        tmp = [tcspcdata(timegate(n,1):timegate(n,2),channels(n)); tcspcdata(timegate(n,3):timegate(n,4),channels(n))];
+    end
+    tcspc(:,n) = tmp;
+end
+
+if pic
+    figure;
+    semilogy((1:Ngate), tcspc);
+end
+