GenerateGauss.cpp

#include <Base/Math/GenerateGauss.hh>

#include <Base/Math/Vector.hh>
#include <Base/Math/Vector.hh>

using namespace BIAS;
using namespace std;

template <class T>
void GenerateGauss::
Gauss1DSigma(const double &sigma, const double &ratio, BIAS::Vector<T> &result, 
             const bool normalize)
{
  unsigned size;
  Sigma2Size_(sigma, ratio, size);
  Gauss1D(size, sigma, result, normalize);
}


template <class T>
void GenerateGauss::
Gauss1DSize(const unsigned size, const double &ratio, BIAS::Vector<T> &result, 
            double &sigma, const bool normalize)
{
  Size2Sigma_(size, ratio, sigma);
  Gauss1D(size, sigma, result, normalize);
}


template <class T>
void GenerateGauss::
Gauss1D(const unsigned size, const double &sigma, BIAS::Vector<T> &result,
        const bool normalize)
{
  result.newsize(size);
  const double fac = 1.0/(2.0*sigma*sigma);
  // n = 1/sqrt(2+M_PI)/sigma
  const double n = 0.39894228040143267794 / sigma;
  T sum = (T)0;
  int hws = size/2;
  if (size%2==0){
    for (int i=0; i<hws; i++){
      result[hws+i] = result[hws-i-1] = 
        (T)(n*exp(-((double)(i)+0.5)*((double)(i)+0.5)*fac));
      if (normalize)
        sum += ((T)2.0)*result[hws+i];
    }
  } else {
    for (int i=0; i<=hws; i++){
      result[hws+i] = result[hws-i] =
          (T)(n*exp(-(double)(i)*(double)(i)*fac));
      if (normalize)
        sum += ((T)2.0)*result[hws-i];
    }
    if (normalize)
      sum -= result[hws]; // added twice
  }

  if (normalize){
    // normalize filter, so that the sum is 1.0
    result /= sum;
  }
}


template <class T>
void GenerateGauss::
Gauss2DSigma(const double &sigma, const double &ratio,
             BIAS::Matrix<T> &result,const bool normalize)
{
  unsigned size;
  Sigma2Size_(sigma, sqrt(ratio), size);
  Gauss2D(size, sigma, result, normalize);
}


template <class T>
void GenerateGauss::
Gauss2DSize(const unsigned size, const double &ratio,
            BIAS::Matrix<T> &result, double &sigma, const bool normalize)
{
  Size2Sigma_(size, sqrt(ratio), sigma);
  Gauss2D(size, sigma, result, normalize);
}


template <class T>
void GenerateGauss::
Gauss2D(const unsigned size, const double &sigma,
        BIAS::Matrix<T> &result, const bool normalize)
{
  Vector<T> vg;
  Gauss1D(size, sigma, vg, false);
  result = vg.OuterProduct(vg);
  if (normalize){
    const int num = size*size;
    const T *data = result.GetData();
    double sum = 0.0;
    for (int i=0; i<num; i++)
      sum += (double)(data[i]);
    result /= (T)sum;
  }
}


//////////////////////////////////////////////////////////////////////////
// instantiation

#define INST(type)\
  template BIASMathBase_EXPORT void GenerateGauss::Gauss1DSigma(const double &sigma, const double &ratio, BIAS::Vector<type> &result,const bool normalize);\
  template BIASMathBase_EXPORT void GenerateGauss::Gauss1DSize(const unsigned size, const double &ratio, BIAS::Vector<type> &result, double &sigma, const bool normalize); \
  template BIASMathBase_EXPORT void GenerateGauss::Gauss1D(const unsigned size, const double &sigma, BIAS::Vector<type> &result, const bool normalize);\
  template BIASMathBase_EXPORT void GenerateGauss::Gauss2DSigma(const double &sigma, const double &ratio, BIAS::Matrix<type> &result,const bool normalize);\
  template BIASMathBase_EXPORT void GenerateGauss::Gauss2DSize(const unsigned size, const double &ratio, BIAS::Matrix<type> &result, double &sigma, const bool normalize);\
  template BIASMathBase_EXPORT void GenerateGauss::Gauss2D(const unsigned size, const double &sigma, BIAS::Matrix<type> &result, const bool normalize)


INST(float);
INST(double);