1 | #ifndef STATS_H |
---|
2 | #define STATS_H |
---|
3 | |
---|
4 | #include <math.h> |
---|
5 | #ifndef UTILS_H |
---|
6 | #include <Utils/utils.hh> |
---|
7 | #endif |
---|
8 | |
---|
9 | namespace Statistics |
---|
10 | { |
---|
11 | |
---|
12 | // Divide by the following correction factor to convert from |
---|
13 | // MADFM to sigma estimator. |
---|
14 | const float correctionFactor = 0.6744888; |
---|
15 | // Multiply by the following correction factor to convert from |
---|
16 | // trimmedSigma to sigma estimator. |
---|
17 | const double trimToNormal = 1.17036753077; |
---|
18 | |
---|
19 | template <class T> float madfmToSigma(T madfm){ |
---|
20 | return float(madfm)/correctionFactor; |
---|
21 | }; |
---|
22 | |
---|
23 | template <class Type> |
---|
24 | class StatsContainer |
---|
25 | { |
---|
26 | public: |
---|
27 | StatsContainer(){useRobust=true;defined=false;useFDR=false;}; |
---|
28 | virtual ~StatsContainer(){}; |
---|
29 | |
---|
30 | float getMean(){return mean;}; |
---|
31 | void setMean(float f){mean=f;}; |
---|
32 | float getStddev(){return stddev;}; |
---|
33 | void setStddev(float f){stddev=f;}; |
---|
34 | Type getMedian(){return median;}; |
---|
35 | void setMedian(Type f){median=f;}; |
---|
36 | Type getMadfm(){return madfm;}; |
---|
37 | void setMadfm(Type f){madfm=f;}; |
---|
38 | float getThreshold(){return threshold;}; |
---|
39 | void setThreshold(float f){threshold=f;}; |
---|
40 | float getThresholdSNR(){ |
---|
41 | return (threshold - this->getMiddle())/this->getSpread();}; |
---|
42 | void setThresholdSNR(float snr){ |
---|
43 | threshold=this->getMiddle() + snr*this->getSpread();}; |
---|
44 | float getPThreshold(){return pThreshold;}; |
---|
45 | void setPThreshold(float f){pThreshold=f;}; |
---|
46 | bool getRobust(){return useRobust;}; |
---|
47 | void setRobust(bool b){useRobust=b;}; |
---|
48 | bool setUseFDR(){return useFDR;}; |
---|
49 | void setUseFDR(bool b){useFDR=b;}; |
---|
50 | |
---|
51 | float getMiddle(){if(useRobust) return float(median); else return mean;}; |
---|
52 | float getSpread(){ |
---|
53 | if(useRobust) return madfmToSigma(madfm); |
---|
54 | else return stddev; |
---|
55 | }; |
---|
56 | |
---|
57 | float getPValue(float value){ |
---|
58 | float zStat = (value - this->getMiddle()) / this->getSpread(); |
---|
59 | return 0.5 * erfc( zStat / M_SQRT2 ); |
---|
60 | }; |
---|
61 | |
---|
62 | bool isDetection(float value){ |
---|
63 | if(useFDR) return (this->getPValue(value) < this->pThreshold); |
---|
64 | else return (value > this->threshold); |
---|
65 | }; |
---|
66 | |
---|
67 | // Functions to calculate the stats for a given array. |
---|
68 | // The idea here is that there are two options to do the calculations: |
---|
69 | // *The first just uses all the points in the array. If you need to |
---|
70 | // remove BLANK points (or something similar), do this beforehand. |
---|
71 | // *Alternatively, construct a bool array of the same size, showing which |
---|
72 | // points are good, and use the second option. |
---|
73 | void calculate(Type *array, long size); |
---|
74 | void calculate(Type *array, long size, bool *isGood); |
---|
75 | |
---|
76 | private: |
---|
77 | bool defined; // a flag indicating whether the stats are defined. |
---|
78 | |
---|
79 | // basic statistics |
---|
80 | float mean; |
---|
81 | float stddev; |
---|
82 | Type median; |
---|
83 | Type madfm; |
---|
84 | |
---|
85 | float threshold; // a threshold for simple sigma-clipping |
---|
86 | float pThreshold; // a threshold for the FDR case -- the upper limit |
---|
87 | // of P values that detected pixels can have. |
---|
88 | bool useRobust; // whether we use the two robust stats or not |
---|
89 | bool useFDR; // whether the FDR method is used for determining a |
---|
90 | // detection |
---|
91 | |
---|
92 | }; |
---|
93 | |
---|
94 | } |
---|
95 | |
---|
96 | #endif /*STATS_H*/ |
---|