1 | /* |
---|
2 | This file is part of the SLIM-curve package for exponential curve fitting of spectral lifetime data. |
---|
3 | |
---|
4 | Copyright (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 | |
---|
51 | int 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 | |
---|