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

Revision 7803, 4.4 KB checked in by paulbarber, 8 years ago (diff)

Fixed some memory leaks and some compile warnings.

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