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

Revision 7781, 4.3 KB checked in by paulbarber, 8 years ago (diff)

Added GPL license to all source code.

Line 
1/*
2This file is part of the SLIM-curve package for exponential curve fitting of spectral lifetime data.
3
4Copyright (c) 2010, 2011, Gray Institute University of Oxford & UW-Madison LOCI.
5
6    This program is free software: you can redistribute it and/or modify
7    it under the terms of the GNU General Public License as published by
8    the Free Software Foundation, either version 3 of the License, or
9    (at your option) any later version.
10
11    This program is distributed in the hope that it will be useful,
12    but WITHOUT ANY WARRANTY; without even the implied warranty of
13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14    GNU General Public License for more details.
15
16    You should have received a copy of the GNU General Public License
17    along with this program.  If not, see <http://www.gnu.org/licenses/>.
18 */
19
20#include "GCI_Phasor.h"
21#include <math.h>
22
23#ifndef NULL
24#define NULL 0
25#endif
26
27#define PHASOR_ERR_NO_ERROR                         0
28#define PHASOR_ERR_INVALID_DATA                    -1
29#define PHASOR_ERR_INVALID_WINDOW                  -2
30#define PHASOR_ERR_INVALID_MODEL                   -3
31#define PHASOR_ERR_FUNCTIONALITY_NOT_SUPPORTED     -4
32
33
34// See Clayton 2004 or Leray 2008
35// Classic Phasor or Polar approach to FLIM
36
37// u = integral(data(t) * cos(wt)) dt / integral(data(t)) dt
38// v = integral(data(t) * sin(wt)) dt / integral(data(t)) dt
39//
40// tau phi = taup = 1/w * (v/u)
41// tau mod = taum = 1/w * sqrt(1/(u^2 + v^2) - 1)
42// tau average = tau = (taup + taum) / 2
43
44// Z must have been estimated previously so that it can be subtracted from the data here
45// A is estimated by making the photon count in the fit the same as the data
46
47// chisq is calculated for comparison with other methods but is not used as there is no optimisation with this method
48
49
50int    GCI_Phasor(float xincr, float y[], int fit_start, int fit_end,
51                                                          float *Z, float *U, float *V, float *taup, float *taum, float *tau, float *fitted, float *residuals,
52                                                          float *chisq)
53{
54    // Z must contain a bg estimate
55       
56        int   i, ret = PHASOR_ERR_NO_ERROR, nBins;
57        float *data, u, v, A, w, I, Ifit, bg, chisq_local, res, sigma2;
58
59        data = &(y[fit_start]); 
60        nBins = (fit_end - fit_start);
61        bg = *Z;
62
63    if (!data)
64        return (PHASOR_ERR_INVALID_DATA);
65    if (nBins<0)
66        return (PHASOR_ERR_INVALID_WINDOW);
67
68        // rep frequency, lets use the period of the measurement, but we can stay in the units of bins
69        w = 2.0*3.1415926535897932384626433832795028841971/(float)nBins; //2.0*PI/(float)nBins;
70       
71        // integral over data
72        for (i=0, I=0.0; i<nBins; i++) 
73                I += (data[i]-bg);
74
75        // Phasor coords
76    // Take care that values correspond to the centre of the bin, hence i+0.5
77        for (i=0, u=0.0; i<nBins; i++) 
78                u += (data[i]-bg) * cos(w*((float)i+0.5));
79        u /= I;
80
81        for (i=0, v=0.0; i<nBins; i++) 
82                v += (data[i]-bg) * sin(w*((float)i+0.5));
83        v /= I;
84
85        // taus, convert now to real time with xincr
86        *taup = (xincr/w) * (v/u);
87        *taum = (xincr/w) * sqrt(1.0/(u*u + v*v) - 1.0);
88
89        *tau = ((*taup) + (*taum))/2.0;
90
91        *U = u;
92        *V = v;
93
94        /* Now calculate the fitted curve and chi-squared if wanted. */
95        if (fitted == NULL)
96                return 0;
97
98        memset(fitted, 0, fit_end*sizeof(float));
99        memset(residuals, 0, fit_end*sizeof(float));
100       
101        // integral over nominal fit data
102        for (Ifit=0.0, i=fit_start; i<fit_end; i++) 
103                Ifit += exp(-(i-fit_start)*xincr/(*tau));
104
105        // Estimate A
106        A = I / Ifit;
107
108        // Calculate fit
109        for (i=fit_start; i<fit_end; i++)
110                fitted[i] = bg + A * exp(-(i-fit_start)*xincr/(*tau));
111
112        // OK, so now fitted contains our data for the timeslice of interest.
113        // We can calculate a chisq value and plot the graph, along with
114        // the residuals.
115
116        if (residuals == NULL && chisq == NULL)
117                return 0;
118
119        chisq_local = 0.0;
120        for (i=0; i<fit_start; i++) {
121                res = y[i]-fitted[i];
122                if (residuals != NULL)
123                        residuals[i] = res;
124        }
125
126
127//      case NOISE_POISSON_FIT:
128                /* Summation loop over all data */
129                for (i=fit_start ; i<fit_end; i++) {
130                        res = y[i] - fitted[i];
131                        if (residuals != NULL)
132                                residuals[i] = res;
133                        /* don't let variance drop below 1 */
134                        sigma2 = (fitted[i] > 1 ? 1.0/fitted[i] : 1.0);
135                        chisq_local += res * res * sigma2;
136                }
137
138        if (chisq != NULL)
139                *chisq = chisq_local;
140
141    return (ret);
142}
143
Note: See TracBrowser for help on using the repository browser.