#ifndef GAUSS_H
#define GAUSS_H

#include <boost/multi_array.hpp>
#include <cassert>
#include <cmath>
#include <deque>
#include "image.h"
#include "image_conv.h"
#include "image_op.h"

namespace mimas {

/** @defgroup gauss Gaussian blur and Gauss gradient
    Methods for blurring images with a gauss-bell and gauss-gradient.
    The filter-parameter \f$\sigma\f$ can be choosen and the size of the filter
    is computed by choosing an upper bound for the approximation-error.

    The following example demonstrates how to blur an image:
    \include gauss_tool/main.cc

    @author Stuart Meikle (stu@stumeikle.org)
    @author Jan Wedekind (jan@wedesoft.de)    
    @date Fri Apr 07 18:52:00 2006
    @{ */
/** Compute gauss-bell.
    The values of the cells are computed by using differences of the integral
    of the gauss-function:
    \f$\displaystyle\int_{-r}^{+r}{\frac{1}{\sqrt{2\,\pi}\,\sigma}\,e^{-\displaystyle\frac{x^2}{2\,\sigma^2}}\,\mathrm{d}x}\ =\ \mathrm{erf}(\displaystyle\frac{r}{\sqrt{2}\,\sigma})\f$

    The coefficients are normalised afterwards such that the sum of all
    elements of the filter is \f$1\f$.
    @param sigma Standard deviation of gauss-distribution.
    @param maxError Maximum error boundary
    (relative to range of pixel-values).
    @see gaussBlur */
template< typename T >
std::deque< T > gaussBlurFilter( T sigma, T maxError = (T)( 1.0 / 256.0 ) );

/** Blur 2-D array.
    Perform gaussian blur on 2-D array.
    @param x Input array.
    @param sigma Standard deviation.
    @param maxError Maximum error boundary (relative to range of pixel-values).
    @see gaussBlurFilter */
template< typename T >
boost::multi_array< T, 2 > gaussBlur
  ( const boost::const_multi_array_ref< T, 2 > &x, T sigma, T maxError );

/** Blur image.
    Perform gaussian blur on 2-D image.
    @param x Input image.
    @param sigma Standard deviation.
    @param maxError Maximum error boundary (relative to range of pixel-values).
    @see gaussBlurFilter */
template< typename T >
image< T > gaussBlur( const image< T > &x, T sigma,
                      T maxError = (T)( 1.0 / 256.0 ) )
{
  boost::const_multi_array_ref< T, 2 > data
    ( x.rawData(), boost::extents[ x.getHeight() ][ x.getWidth() ] );
  image< T > retVal; retVal.init( x.getWidth(), x.getHeight() );
  boost::multi_array_ref< T, 2 >
    ( retVal.rawData(),
      boost::extents[ retVal.getHeight() ][ retVal.getWidth() ] ) =
    gaussBlur< T >( data, sigma, maxError );
  return retVal;
}

/** Compute gauss-gradient.
    The values of the cells are computed by using differences of the
    integral (the gauss-function):
    \f$\frac{1}{\sqrt{2\,\pi}\,\sigma}\,e^{-\displaystyle\frac{x^2}{2\,\sigma^2}}\,\mathrm{d}x\big\|_r^\infty\f$

    The coefficients are normalised afterwards such that the sum of the square
    of all elements of the filter is \f$1\f$.
    @param sigma Standard deviation of gauss-distribution.
    @param maxError Maximum error boundary
    (relative to range of pixel-values).
    @see gaussBlur */
template< typename T >
std::deque< T > gaussGradientFilter( T sigma,
                                     T maxError = (T)( 1.0 / 256.0 ) );

/** Take x-gradient of 2-D array.
    Compute gauss-gradient of 2-D array in x-direction.
    @param x Input array.
    @param sigma Standard deviation.
    @param maxError Maximum error boundary (relative to range of pixel-values).
    @see gaussGradientFilter */
template< typename T >
boost::multi_array< T, 2 > gaussGradientX
  ( const boost::const_multi_array_ref< T, 2 > &x, T sigma, T maxError );

/** Take x-gradient of 2-D image.
    Compute gauss-gradient of 2-D image in x-direction.
    @param x Input image.
    @param sigma Standard deviation.
    @param maxError Maximum error boundary (relative to range of pixel-values).
    @see gaussGradientFilter */
template< typename T >
image< T > gaussGradientX( const image< T > &x, T sigma,
                           T maxError = (T)( 1.0 / 256.0 ) );

/** Take y-gradient of 2-D array.
    Compute gauss-gradient of 2-D array in y-direction.
    @param x Input array.
    @param sigma Standard deviation.
    @param maxError Maximum error boundary (relative to range of pixel-values).
    @see gaussGradientFilter */
template< typename T >
boost::multi_array< T, 2 > gaussGradientY
  ( const boost::const_multi_array_ref< T, 2 > &x, T sigma, T maxError );

/** Take y-gradient of 2-D image.
    Compute gauss-gradient of 2-D image in y-direction.
    @param x Input image.
    @param sigma Standard deviation.
    @param maxError Maximum error boundary (relative to range of pixel-values).
    @see gaussGradientFilter */
template< typename T >
image< T > gaussGradientY( const image< T > &x, T sigma,
                           T maxError = (T)( 1.0 / 256.0 ) );

/** Square of gradient-norm.
    Compute square of gradient-norm for 2-D image
    @param im Input image.
    @param sigma Standard deviation.
    @param maxError Maximum error boundary (relative to range of pixel-values).
    @see gaussGradientFilter */
template< typename T >
image< T > gaussGradientSqr( const image< T > &im, T sigma,
                             T maxError = (T)( 1.0 / 256.0 ) )
{
  image< T >
    gradX( gaussGradientX( im, sigma, maxError ) ),
    gradY( gaussGradientY( im, sigma, maxError ) );
  return sumSquares( gradX, gradY );
}

/** Gradient-norm.
    Compute gradient-norm for 2-D image
    @param im Input image.
    @param sigma Standard deviation.
    @param maxError Maximum error boundary (relative to range of pixel-values).
    @see gaussGradientFilter */
template< typename T >
image< T > gaussGradientNorm( const image< T > &im, T sigma,
                              T maxError = (T)( 1.0 / 256.0 ) )
{
  return squareRoot( gaussGradientSqr( im, sigma, maxError ) );
}

///@}

};

#include "gauss.tcc"

#endif