|| || || || || || || ||

---

MakeGaussian1d.h

Go to the documentation of this file.

#ifndef Impala_Core_Array_MakeGaussian1d_h
#define Impala_Core_Array_MakeGaussian1d_h

//#include "Data/Basis/CxMath.h"
#include "Core/Array/Arrays.h"

namespace Impala
{
namespace Core
{
namespace Array
{


// recursive Gaussian derivative: 
// deriv(x,n) = -(n-1)/sigma^2 * deriv(x,n-2) - x/sigma^2 * deriv(x,n-1)
// with Gauss = deriv(x,0) = b exp(a x^2)

static inline double
Gauss(const double x, const double sigma)
{
    // normalization to sum=1. is performed later
    return exp(-0.5 * x * x / (sigma*sigma));
}

static inline double
Hermite(const double x, const double H0, const double sigma, const int order)
{
    if (!order)
        return H0;

    double a = 1./(sigma*sigma);
    double ax = a*x;
    double Hnmin1 = H0;
    double Hn = -ax*H0;

    int n;
    for( n = 2; n <= order; n++ ) {
        double Hnmin2 = Hnmin1;
        Hnmin1 = Hn;
        Hn = a*(1-n)*Hnmin2 - ax*Hnmin1;
    }

    return (order%2) ? -Hn : Hn;
}

static inline int
filterWidth(double sigma, int deri, double acc, int maxlen)
{
    double     acc2;
    int        fsize;

    acc2 = 1.0 - (1.0 - acc)/2.0;

    fsize = 2*int(acc*sigma+0.5)+1;

    /*  filter kernel is always odd sized, so if maxlen is even 
     *  subtract one.
     */
    if ( !(maxlen%2) )
        maxlen--; 
    if (fsize > maxlen)
        fsize = maxlen;

    return fsize;
}

static inline double*
makeFilter(double sigma, int deri, double acc, int fsize, int maxfsize)
{
    double* filter = new double[fsize];

    fsize = fsize/2; /* adjust to one side ,don't count the centre */
    double* center = filter + fsize;

    /* calc filter */
    int i;
    double sum = Gauss(0.0, sigma);
    center[0] = sum;

    for (i=1; i<=fsize; i++) {
        double val = Gauss(i, sigma);
        sum += val+val;
        center[i] = val;
        center[-i] = val;
    }

    /* normalize to sum=1.0 */
    for (i=-fsize; i<=fsize; i++)
        center[i] /= sum;

    /* replace by Hermite polynomial of order deri */
    if (deri > 0)
        for (i=-fsize; i<=fsize; i++)
            center[i] = Hermite(i, center[i], sigma, deri);

    fsize = fsize * 2 + 1; /* expand to 2-sided width */
    return filter;
}

inline Array2dScalarReal64*
MakeGaussian1d(double sigma, int deri, double acc, int maxfsize, int fsize = 0)
{
    if (fsize < 1) { /* determine filter size if not specified */
        fsize = filterWidth(sigma, deri, acc, maxfsize);
    }

    double* fData = makeFilter(sigma, deri, acc, fsize, maxfsize);
    return ArrayCreate<Array2dScalarReal64>(fsize, 1, 0, 0, fData);
}

} // namespace Array
} // namespace Core
} // namespace Impala

#endif

---