Statistics.h

Go to the documentation of this file.

 
#ifndef Statistics_HPP
#define Statistics_HPP
 
// common header
#include <vector>
#include <stdio.h>
 
#include "Helper.h"
 
#ifdef ARRAYFIRE
 
#define NOMINMAX 
#include <arrayfire.h>
#include <af/util.h>
 
#endif
 
// what about Boost.compute? in Boost 1.61 not yet using newest version of boost...
 
class Statistics
{
public:
    template <typename T>
    static float variance(int start, int offset, const T *in, int len)
    {
        double n = len;
        double xbar = 0;
        double _var = 0;
 
#ifdef ARRAYFIRE
        if (true)
        {
            doInitArrayFire();
            af::array a_in = af::array(len, (in + start + offset));
 
            _var = af::var<double>(a_in, false);
            return (float)_var;
        }
#endif
 
        for (int i4 = start + offset; i4 < len + start; i4++)
        {
            xbar += *(in + i4) / n;
        }
 
        for (int i4 = start + offset; i4 < len + start; i4++)
        {
            _var += (*(in + i4) - xbar) *(*(in + i4) - xbar);
        }
 
        _var = _var / (n - 1.0);
 
        return (float)_var;
    }
 
    // not used?
    template <typename T>
    static void histogramIt(const T  *data_in, int len, int  num_bins, T _min, T  *&data_out)
    {
#ifdef ARRAYFIRE
        if (true)
        {
            doInitArrayFire();
            af::array hist_in = af::array(len, data_in);
            af::array hist_out = histogram(hist_in, num_bins, _min, 65530);
            hist_out.host(data_out);
        }
#endif
    }
 
    template<typename T>
    static void saveHistogramToFile(const T  *out_data, int len) {}
 
    template<typename T>
    static void cdfIt(const T *in_data, int len, int num_bins, float _min, T *&cdf_out) {}
 
    template< typename T, typename F>
    static void minMaxData(const T  *data_out, int width, int height, F &minVal, F &maxVal, int trim = 0)
    {
        int image_size = width * height;
        maxVal = 0;
        minVal = 65530;
        T myMax, myMin;
        myMax = 0;
        myMin = 65535;
 
#ifdef ARRAYFIRE
        af::array a_in = af::array(image_size, data_out);
        unsigned uIdx;
 
        af::min(&myMin, &uIdx, a_in);
 
        // want to use mask to eliminate all of the values before trim, i.e. first 10 rows
        if (trim > 0)
        {
            doSetTrimMaskArray(trim, width, height);
            af::max(&myMax, &uIdx, a_in*trimMaskArray);
        }
        else
        {
            af::max(&myMax, &uIdx, a_in);
        }
#else
 
        // find range of values in out_data
        T val;
        for (int row = 0; row < height; row++)
        {
            for (int col = 0; col < width; col++)
            {
                // ignore first few lines
                //if(!isZeroPad)
                if (col >= trim)
                {
                    val = *(data_out + row * width + col);
                    if (val > myMax)
                    {
                        myMax = val;
                    }
                    else if (val < myMin)
                    {
                        myMin = val;
                    }
                }
            }
        }
#endif
        minVal = static_cast<F>(myMin);
        maxVal = static_cast<F>(myMax);
    }
 
    template< typename T, typename F>
    static void minMaxData(std::vector<T> *data_out, int w, int h, F &minVal, F &maxVal, int trim = 0)
    {
        return minMaxData(data_out->data(), w, h, minVal, maxVal, trim);
    }
 
    template< typename T>
    static void minMaxData(const T *data_in, int len, T &maxVal, T &minVal)
    {
        maxVal = 0;
        minVal = 65535;
#ifdef ARRAYFIRE
        doInitArrayFire();
        af::array a_in = af::array(len, data_in);
        unsigned uIdx;
        T myMax, myMin;
        af::min(&myMin, &uIdx, a_in);
        af::max(&myMax, &uIdx, a_in);
        minVal = myMin;
        maxVal = myMax;
#else
        // find range of values in out_data
        for (int i = 0; i < len; i++)
        {
            T val = *(data_in + i);
            if (val > maxVal)
            {
                maxVal = val;
            }
            else if (val < minVal)
            {
                minVal = val;
            }
        }
#endif
    }
 
    static void minMaxDataTrim(const unsigned short *data_out, int w, int h, int trim, double &minVal, double &maxVal);
    static void minMaxDataAvg(const unsigned short *data_out, int w, int h, int trim, double &minVal, double &maxVal);
 
    template< typename T, typename F>
    static void avgMatrixToLine(const T* in, int startRow, int stopRow, int w, F *&out)
    {
        double n = (double)(stopRow - startRow);
        for (int row = startRow; row < stopRow; row++)
        {
            for (int col = 0; col < w; col++)
            {
                if (row == startRow)
                {
                    *(out + col) = 0;
                }
                *(out + col) += static_cast<F>(*(in + row * w + col) / n);
            }
        }
    }
 
    static void movingFilterLine(float *&in, int len, int filter);
    static void movingFilterLine(const float *in, int len, int filter, float *&out);
    static void gradientLine(float *&in, int len);
    static void gradientLine(const float *in, int len, float *&out);
    static void filterVector(std::vector<float> *sumOfVar);
 
    template< typename T>
    static T localMaxVal(std::vector<T> *in, int &myMaxIndex)
    {
        T myMax = 0;
        int len = (int)in->size();
        myMaxIndex = len / 2;
        return localMaxVal(in->data(), len, myMaxIndex);
    }
 
    template<typename T>
    static T localMaxVal(const T* in, int len, int &myMaxIndex)
    {
        T myMax = 0;
        myMaxIndex = len / 2;
 
#ifdef ARRAYFIRE
        if (false)
        {
            doInitArrayFire();
            af::array a_in = af::array(len, in);
 
            unsigned uIdx;
            af::max(&myMax, &uIdx, a_in);
            myMaxIndex = uIdx;
        }
#endif
        for (int i = 1; i < len - 1; i++)
        {
            T val = *(in + i);
            T diff1 = val - *(in + i - 1);
            T diff2 = *(in + i + 1) - val;
 
            // look for sign change
            bool a = diff1 / std::abs(diff1) != diff2 / std::abs(diff2);
            if ((val > myMax) && a)
            {
                myMax = val;
                myMaxIndex = i;
            }
            //if(a && t)    smessage(" candidate "+i2S(myMax)+" at "+i2S(myMaxIndex) );
        }
        //smessage(" myMax= "+i2S(myMax)+" at "+i2S(myMaxIndex) );
 
        return myMax;
    }
 
    template<typename T>
    static T maxVal(std::vector<T> *in, int &myMaxIndex)
    {
        T myMax = 0;
        int len = (int)in->size();
 
        float myMax = 0;
        myMaxIndex = len / 2;
 
#ifdef ARRAYFIRE
        doInitArrayFire();
        array a_in = array(len, in);
        unsigned uIdx;
        max<T>(&myMax, &uIdx, a_in);
        myMaxIndex = uIdx;
        return myMax;
#endif
 
        for (int i4 = 0; i4 < len; i4++)
        {
            if ((*in)[i4] > myMax)
            {
                myMax = (*in)[i4];
                myMaxIndex = i4;
            }
        }
 
        return myMax;
    }
 
    static bool isOscilating(float *&in, int len, int posGradLimit);
 
protected:
 
    static bool isVerbose;
    static bool initArrayFire;
    static void doInitArrayFire();
    static void doSetTrimMaskArray(int trimVal, int w, int h);
 
#ifdef ARRAYFIRE
 
    static af::array trimMaskArray;
    static int trim_val;
    static int mask_width;
    static int mask_height;
 
#endif
};
#endif