polaroid-pp

Schlieren and contour plot tool
git clone https://git.0xfab.ch/polaroid-pp.git
Log | Files | Refs | Submodules | README | LICENSE

Types.h (4145B)


      1 // File       : Types.h
      2 // Date       : Tue Apr 26 14:49:41 2016
      3 // Author     : Fabian Wermelinger
      4 // Description: Various types
      5 // Copyright 2016 ETH Zurich. All Rights Reserved.
      6 #ifndef TYPES_H_QATFIWCK
      7 #define TYPES_H_QATFIWCK
      8 
      9 #include <cassert>
     10 #include <algorithm>
     11 #include <cmath>
     12 #include "common.h"
     13 
     14 template <int _SX, int _EX, int _SY, int _EY>
     15 class Slice2D
     16 {
     17     int m_width, m_height, m_N;
     18     Real* m_data;
     19 
     20     int m_max3Ddim;
     21     int m_sliceID;
     22     char m_sliceDir;
     23 
     24 public:
     25     Slice2D() : m_width(0), m_height(0), m_N(0), m_data(nullptr) {}
     26     Slice2D(const int width, const int height) : m_width(width), m_height(height), m_N((width+_EX-_SX)*(height+_EY-_SY))
     27     {
     28         m_data = new Real[m_N];
     29         // for (int i = 0; i < m_N; ++i)
     30         //     m_data[i] = 0.0;
     31     }
     32     virtual ~Slice2D() { delete [] m_data; }
     33 
     34     inline Real operator()(const int ix, const int iy) const
     35     {
     36         assert(ix >= _SX); assert(ix < m_width+_EX);
     37         assert(iy >= _SY); assert(iy < m_height+_EY);
     38 
     39         return m_data[(iy-_SY)*m_width + (ix-_SX)];
     40     }
     41 
     42     inline Real& operator()(const int ix, const int iy)
     43     {
     44         assert(ix >= _SX); assert(ix < m_width+_EX);
     45         assert(iy >= _SY); assert(iy < m_height+_EY);
     46 
     47         return m_data[(iy-_SY)*m_width + (ix-_SX)];
     48     }
     49 
     50     inline void resize(const int width, const int height)
     51     {
     52         m_width = width;
     53         m_height = height;
     54         m_N = (width+_EX-_SX)*(height+_EY-_SY);
     55         if (m_data) delete [] m_data;
     56         m_data = new Real[m_N];
     57         // for (int i = 0; i < m_N; ++i)
     58         //     m_data[i] = 0.0;
     59     }
     60 
     61     inline int width() const { return m_width; }
     62     inline int height() const { return m_height; }
     63     inline Real min() const
     64     {
     65         Real f = m_data[0];
     66         for (int i = 1; i < m_N; ++i)
     67             f = std::min(f, m_data[i]);
     68         return f;
     69     }
     70     inline Real max() const
     71     {
     72         Real f = m_data[0];
     73         for (int i = 1; i < m_N; ++i)
     74             f = std::max(f, m_data[i]);
     75         return f;
     76     }
     77     inline void set_max3Ddim(const int d) { m_max3Ddim = d; }
     78     inline int get_max3Ddim() const { return m_max3Ddim; }
     79     inline void set_sliceID(const int d) { m_sliceID = d; }
     80     inline int get_sliceID() const { return m_sliceID; }
     81     inline void set_sliceDir(const char d) { m_sliceDir = d; }
     82     inline char get_sliceDir() const { return m_sliceDir; }
     83 };
     84 
     85 typedef Slice2D<0,0,0,0> Slice;
     86 
     87 
     88 class Slice2D_Statistics
     89 {
     90     double m_mean; // 1st moment
     91     double m_var;  // 2nd moment
     92     double m_skew; // 3rd moment
     93     double m_kurt; // 4th moment
     94     double m_min, m_max;
     95 
     96     void _compute_stat(const Slice& s)
     97     {
     98         m_min = s.min();
     99         m_max = s.max();
    100 
    101         int k = 0;
    102         double mean = 0;
    103         double M2 = 0, M3 = 0, M4 = 0;
    104         for (int h=0; h < s.height(); ++h)
    105             for (int w=0; w < s.width(); ++w)
    106             {
    107                 const int k1 = k;
    108                 ++k;
    109                 const double delta = s(w,h) - mean;
    110                 const double delta_k = delta / k;
    111                 const double delta_k2 = delta_k * delta_k;
    112                 const double term1 = delta * delta_k * k1;
    113                 mean += delta_k;
    114                 M4 += term1 * delta_k2 * (k*k - 3*k + 3) + 6 * delta_k2 * M2 - 4 * delta_k * M3;
    115                 M3 += term1 * delta_k * (k - 2) - 3 * delta_k * M2;
    116                 M2 += term1;
    117             }
    118         assert(k > 1);
    119 
    120         m_mean = mean;
    121         m_var  = M2 / (k - 1);
    122         m_skew = std::sqrt(k) * M3 / std::pow(M2, 1.5);
    123         m_kurt = k * M4 / (M2 * M2) - 3;
    124     }
    125 
    126 public:
    127     Slice2D_Statistics(const Slice& s) :
    128         m_mean(0), m_var(0), m_skew(0), m_kurt(0), m_min(HUGE_VAL), m_max(-HUGE_VAL)
    129     {
    130         _compute_stat(s);
    131     }
    132 
    133     // accessors
    134     double mean() const { return m_mean; }
    135     double var()  const { return m_var; }
    136     double std()  const { return std::sqrt(m_var); }
    137     double skew() const { return m_skew; }
    138     double kurt() const { return m_kurt; }
    139     double min()  const { return m_min; }
    140     double max()  const { return m_max; }
    141 };
    142 
    143 #endif /* TYPES_H_QATFIWCK */