/************************************************************************/
/*                                                                      */
/*               Copyright 1998-2004 by Ullrich Koethe                  */
/*       Cognitive Systems Group, University of Hamburg, Germany        */
/*                                                                      */
/*    This file is part of the VIGRA computer vision library.           */
/*    ( Version 1.3.2, Jan 27 2005 )                                    */
/*    You may use, modify, and distribute this software according       */
/*    to the terms stated in the LICENSE file included in               */
/*    the VIGRA distribution.                                           */
/*                                                                      */
/*    The VIGRA Website is                                              */
/*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
/*    Please direct questions, bug reports, and contributions to        */
/*        koethe@informatik.uni-hamburg.de                              */
/*                                                                      */
/*  THIS SOFTWARE IS PROVIDED AS IS AND WITHOUT ANY EXPRESS OR          */
/*  IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED      */
/*  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */
/*                                                                      */
/************************************************************************/

#ifndef VIGRA_GAUSSIANS_HXX
#define VIGRA_GAUSSIANS_HXX

#include <cmath>
#include "vigra/config.hxx"
#include "vigra/mathutil.hxx"
#include "vigra/array_vector.hxx"

namespace vigra {

#if 0
/** \addtogroup MathFunctions Mathematical Functions
*/
//@{
#endif
/*! The Gaussian function and its derivatives.

    Implemented as a unary functor. Since it supports the <tt>radius()</tt> function
    it can also be used as a kernel in \ref resamplingConvolveImage().

    <b>\#include</b> "<a href="gaussians_8hxx-source.html">vigra/gaussians.hxx</a>"<br>
    Namespace: vigra

    \ingroup MathFunctions
*/
template <class T = double>
class Gaussian
{
  public:
  
        /** the value type if used as a kernel in \ref resamplingConvolveImage().
        */
    typedef T            value_type;  
        /** the functor's argument type
        */
    typedef T            argument_type;  
        /** the functor's result type
        */
    typedef T            result_type; 
    
        /** Create functor for the given standard deviation <tt>sigma</tt> and 
            derivative order <i>n</i>. The functor then realizes the function
             
            \f[ f_{\sigma,n}(x)=\frac{\partial^n}{\partial x^n}
                 \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{x^2}{2\sigma^2}}
            \f]
             
            Precondition:
            \code
            sigma > 0.0
            \endcode
        */
    explicit Gaussian(T sigma = 1.0, unsigned int derivativeOrder = 0)
    : sigma_(sigma),
      sigma2_(-0.5 / sigma / sigma),
      norm_(0.0),
      order_(derivativeOrder),
      hermitePolynomial_(derivativeOrder / 2 + 1)
    {
        vigra_precondition(sigma_ > 0.0,
            "Gaussian::Gaussian(): sigma > 0 required.");
        switch(order_)
        {
            case 1:
            case 2:
                norm_ = -1.0 / (VIGRA_CSTD::sqrt(2.0 * M_PI) * sq(sigma) * sigma);
                break;
            case 3:
                norm_ = 1.0 / (VIGRA_CSTD::sqrt(2.0 * M_PI) * sq(sigma) * sq(sigma) * sigma);
                break;
            default:
                norm_ = 1.0 / VIGRA_CSTD::sqrt(2.0 * M_PI) / sigma;
        }
        calculateHermitePolynomial();
    }

        /** Function (functor) call.
        */
    result_type operator()(argument_type x) const;

        /** Get the standard deviation of the Gaussian.
        */
    value_type sigma() const
        { return sigma_; }
    
        /** Get the derivative order of the Gaussian.
        */
    unsigned int derivativeOrder() const
        { return order_; }
    
        /** Get the required filter radius for a discrete approximation of the Gaussian.
            The radius is given as a multiple of the Gaussian's standard deviation 
            (default: <tt>sigma * (3 + 1/2 * derivativeOrder()</tt> -- the second term
            accounts for the fact that the derivatives of the Gaussian become wider 
            with increasing order). The result is rounded to the next higher integer.
        */
    double radius(double sigmaMultiple = 3.0) const
        { return VIGRA_CSTD::ceil(sigma_ * (sigmaMultiple + 0.5 * derivativeOrder())); }

  private:
    void calculateHermitePolynomial();    
    T horner(T x) const;
    
    T sigma_, sigma2_, norm_;
    unsigned int order_;
    ArrayVector<T> hermitePolynomial_;
};

template <class T>
typename Gaussian<T>::result_type 
Gaussian<T>::operator()(argument_type x) const
{
    T x2 = x * x;
    T g  = norm_ * VIGRA_CSTD::exp(x2 * sigma2_);
    switch(order_)
    {
        case 0:
            return g;
        case 1:
            return x * g;
        case 2:
            return (1.0 - sq(x / sigma_)) * g;
        case 3:
            return (3.0 - sq(x / sigma_)) * x * g;
        default:
            return order_ % 2 == 0 ?
                       g * horner(x2)
                     : x * g * horner(x2);
    }
}

template <class T>
T Gaussian<T>::horner(T x) const
{
    int i = order_ / 2;
    T res = hermitePolynomial_[i];
    for(--i; i >= 0; --i)
        res = x * res + hermitePolynomial_[i];
    return res;
}
    
template <class T>
void Gaussian<T>::calculateHermitePolynomial()
{
    if(order_ == 0)
    {
        hermitePolynomial_[0] = 1.0;
    }
    else if(order_ == 1)
    {
        hermitePolynomial_[0] = -1.0 / sigma_ / sigma_;
    }
    else
    {
        // calculate Hermite polynomial for requested derivative 
        // recursively according to
        //     (0)
        //    h   (x) = 1
        //
        //     (1)
        //    h   (x) = -x / s^2
        //
        //     (n+1)                        (n)           (n-1)
        //    h     (x) = -1 / s^2 * [ x * h   (x) + n * h     (x) ]
        //
        T s2 = -1.0 / sigma_ / sigma_;
        ArrayVector<T> hn(3*order_+3, 0.0);
        typename ArrayVector<T>::iterator hn0 = hn.begin(),
                                          hn1 = hn0 + order_+1,
                                          hn2 = hn1 + order_+1,
                                          ht;
        hn2[0] = 1.0;
        hn1[1] = s2;
        for(unsigned int i = 2; i <= order_; ++i)
        {
            hn0[0] = s2 * (i-1) * hn2[0];
            for(unsigned int j = 1; j <= i; ++j)
                hn0[j] = s2 * (hn1[j-1] + (i-1) * hn2[j]);
            ht = hn2;
            hn2 = hn1;
            hn1 = hn0;
            hn0 = ht;
        }
        // keep only non-zero coefficients of the polynomial
        for(unsigned int i = 0; i < hermitePolynomial_.size(); ++i)
            hermitePolynomial_[i] = order_ % 2 == 0 ?
                                         hn1[2*i]
                                       : hn1[2*i+1];
    }
}


////@}

} // namespace vigra


#endif /* VIGRA_GAUSSIANS_HXX */
