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