source: trunk/projects/slim-curve/src/main/c/GCI_Phasor.c @ 7607

Revision 7607, 3.4 KB checked in by paulbarber, 9 years ago (diff)

Added single frequency phasor routine for lifetime.
Added Non-negative least squares for spectral unmixing.

Line 
1#include "GCI_Phasor.h"
2#include <math.h>
3
4#ifndef NULL
5#define NULL 0
6#endif
7
8#define PHASOR_ERR_NO_ERROR                         0
9#define PHASOR_ERR_INVALID_DATA                    -1
10#define PHASOR_ERR_INVALID_WINDOW                  -2
11#define PHASOR_ERR_INVALID_MODEL                   -3
12#define PHASOR_ERR_FUNCTIONALITY_NOT_SUPPORTED     -4
13
14
15// See Clayton 2004 or Leray 2008
16// Classic Phasor or Polar approach to FLIM
17
18// u = integral(data(t) * cos(wt)) dt / integral(data(t)) dt
19// v = integral(data(t) * sin(wt)) dt / integral(data(t)) dt
20//
21// tau phi = taup = 1/w * (v/u)
22// tau mod = taum = 1/w * sqrt(1/(u^2 + v^2) - 1)
23// tau average = tau = (taup + taum) / 2
24
25// Z must have been estimated previously so that it can be subtracted from the data here
26// A is estimated by making the photon count in the fit the same as the data
27
28// chisq is calculated for comparison with other methods but is not used as there is no optimisation with this method
29
30
31int    GCI_Phasor(float xincr, float y[], int fit_start, int fit_end,
32                                                          float *Z, float *U, float *V, float *taup, float *taum, float *tau, float *fitted, float *residuals,
33                                                          float *chisq)
34{
35    // Z must contain a bg estimate
36       
37        int   i, ret = PHASOR_ERR_NO_ERROR, nBins;
38        float *data, u, v, A, w, I, Ifit, bg, chisq_local, res, sigma2;
39
40        data = &(y[fit_start]); 
41        nBins = (fit_end - fit_start);
42        bg = *Z;
43
44    if (!data)
45        return (PHASOR_ERR_INVALID_DATA);
46    if (nBins<0)
47        return (PHASOR_ERR_INVALID_WINDOW);
48
49        // rep frequency, lets use the period of the measurement, but we can stay in the units of bins
50        w = 2.0*3.1415926535897932384626433832795028841971/(float)nBins; //2.0*PI/(float)nBins;
51       
52        // integral over data
53        for (i=0, I=0.0; i<nBins; i++) 
54                I += (data[i]-bg);
55
56        // Phasor coords
57    // Take care that values correspond to the centre of the bin, hence i+0.5
58        for (i=0, u=0.0; i<nBins; i++) 
59                u += (data[i]-bg) * cos(w*((float)i+0.5));
60        u /= I;
61
62        for (i=0, v=0.0; i<nBins; i++) 
63                v += (data[i]-bg) * sin(w*((float)i+0.5));
64        v /= I;
65
66        // taus, convert now to real time with xincr
67        *taup = (xincr/w) * (v/u);
68        *taum = (xincr/w) * sqrt(1.0/(u*u + v*v) - 1.0);
69
70        *tau = ((*taup) + (*taum))/2.0;
71
72        *U = u;
73        *V = v;
74
75        /* Now calculate the fitted curve and chi-squared if wanted. */
76        if (fitted == NULL)
77                return 0;
78
79        memset(fitted, 0, fit_end*sizeof(float));
80        memset(residuals, 0, fit_end*sizeof(float));
81       
82        // integral over nominal fit data
83        for (Ifit=0.0, i=fit_start; i<fit_end; i++) 
84                Ifit += exp(-(i-fit_start)*xincr/(*tau));
85
86        // Estimate A
87        A = I / Ifit;
88
89        // Calculate fit
90        for (i=fit_start; i<fit_end; i++)
91                fitted[i] = bg + A * exp(-(i-fit_start)*xincr/(*tau));
92
93        // OK, so now fitted contains our data for the timeslice of interest.
94        // We can calculate a chisq value and plot the graph, along with
95        // the residuals.
96
97        if (residuals == NULL && chisq == NULL)
98                return 0;
99
100        chisq_local = 0.0;
101        for (i=0; i<fit_start; i++) {
102                res = y[i]-fitted[i];
103                if (residuals != NULL)
104                        residuals[i] = res;
105        }
106
107
108//      case NOISE_POISSON_FIT:
109                /* Summation loop over all data */
110                for (i=fit_start ; i<fit_end; i++) {
111                        res = y[i] - fitted[i];
112                        if (residuals != NULL)
113                                residuals[i] = res;
114                        /* don't let variance drop below 1 */
115                        sigma2 = (fitted[i] > 1 ? 1.0/fitted[i] : 1.0);
116                        chisq_local += res * res * sigma2;
117                }
118
119        if (chisq != NULL)
120                *chisq = chisq_local;
121
122    return (ret);
123}
124
Note: See TracBrowser for help on using the repository browser.