#
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) |
---|

Line | |
---|---|

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 |

**Note:**See TracBrowser for help on using the repository browser.