namespace mimas {

template< typename T >
std::deque< T > gaussBlurFilter( T sigma, T maxError )
{
  assert( sigma > 0 );

  std::deque< T > retval;

  // Integral of gauss-curve from -0.5 to +0.5.
  T integral = erf( 0.5 / ( sqrt( 2.0 ) * sigma ) );
  retval.push_back( integral );

  while ( 1.0 - integral > maxError ) {
    // Compute integral of gauss-curve from -newSize2 to +newSize2.
    const T
      newSize2 = 0.5 * ( retval.size() + 2 ),
      newIntegral = erf( newSize2 / ( sqrt( 2.0 ) * sigma ) ),
      value = 0.5 * ( newIntegral - integral );
    retval.push_back( value );
    retval.push_front( value );
    integral = newIntegral;
  };

  // Normalise to maintain power of image.
  const T factor = 1.0 / integral;
  for ( typename std::deque< T >::iterator i = retval.begin();
        i != retval.end(); i++ )
    *i *= factor;

  return retval;
}

template< typename T >
boost::multi_array< T, 2 > gaussBlur
  ( const boost::const_multi_array_ref< T, 2 > &x, T sigma, T maxError )
{
  // Split error-budget up for two operations.
  std::deque< T > filter( gaussBlurFilter< T >( sigma, maxError / 2.0 ) );
  boost::multi_array< T, 2 >
    filterx( boost::extents[ 1            ][ filter.size() ] ),
    filtery( boost::extents[ filter.size()][ 1             ] );
  std::copy( filter.begin(), filter.end(), filterx.data() );
  std::copy( filter.begin(), filter.end(), filtery.data() );
  return correlate( correlate( x, filterx ), filtery );
}

template< typename T >
std::deque< T > gaussGradientFilter( T sigma, T maxError )
{
  assert( sigma > 0 );

  std::deque< T > retval;

  // Integral of gauss-gradient from -0.5 to +0.5 is zero.
  retval.push_back( 0.0 );

  // Set to absolute value of integral of gauss-gradient from 0.5 to infinity
  T integral =
    1.0 / ( sqrt( 2.0 * M_PI ) * sigma ) * exp( -0.125 / ( sigma * sigma ) ),
    squareIntegral = 0.0;

  while ( 2.0 * integral > maxError ) {
    // Compute integral of gauss-curve from -newSize2 to +newSize2.
    const T
      newSize2 = 0.5 * ( retval.size() + 2 ),
      newIntegral = 1.0 / ( sqrt( 2.0 * M_PI ) * sigma ) * exp( -newSize2 * newSize2 / ( 2.0 * sigma * sigma ) ),
      value = integral - newIntegral;
    retval.push_back( value );
    retval.push_front( -value );
    integral = newIntegral;
    squareIntegral += value * value;
  };

  // Normalise to maintain power of image.
  const T factor = 1.0 / squareIntegral;
  for ( typename std::deque< T >::iterator i = retval.begin();
        i != retval.end(); i++ )
    *i *= factor;

  return retval;
}

template< typename T >
boost::multi_array< T, 2 > gaussGradientX
  ( const boost::const_multi_array_ref< T, 2 > &x, T sigma, T maxError )
{
  // Split error-budget up for two operations.
  std::deque< T >
    filter1( gaussGradientFilter< T >( sigma, maxError / 2.0 ) ),
    filter2( gaussBlurFilter< T >( sigma, maxError / 2.0 ) );
  boost::multi_array< T, 2 >
    filterx( boost::extents[ 1             ][ filter1.size() ] ),
    filtery( boost::extents[ filter2.size()][ 1              ] );
  std::copy( filter1.begin(), filter1.end(), filterx.data() );
  std::copy( filter2.begin(), filter2.end(), filtery.data() );
  return correlate( correlate( x, filterx ), filtery );
}

template< typename T >
image< T > gaussGradientX( const image< T > &x, T sigma,
                           T maxError )
{
  // Split error-budget up for two operations.
  std::deque< T >
    filter1( gaussGradientFilter< T >( sigma, maxError / 2.0 ) ),
    filter2( gaussBlurFilter< T >( sigma, maxError / 2.0 ) );
  image< T > filterx, filtery;
  filterx.init( filter1.size(), 1 );
  filtery.init( 1, filter2.size() );
  std::copy( filter1.begin(), filter1.end(), filterx.rawData() );
  std::copy( filter2.begin(), filter2.end(), filtery.rawData() );
  return correlate( correlate( x, filterx ), filtery );
}


template< typename T >
boost::multi_array< T, 2 > gaussGradientY
  ( const boost::const_multi_array_ref< T, 2 > &x, T sigma, T maxError )
{
  // Split error-budget up for two operations.
  std::deque< T >
    filter1( gaussBlurFilter< T >( sigma, maxError / 2.0 ) ),
    filter2( gaussGradientFilter< T >( sigma, maxError / 2.0 ) );
  boost::multi_array< T, 2 >
    filterx( boost::extents[ 1             ][ filter1.size() ] ),
    filtery( boost::extents[ filter2.size()][ 1              ] );
  std::copy( filter1.begin(), filter1.end(), filterx.data() );
  std::copy( filter2.begin(), filter2.end(), filtery.data() );
  return correlate( correlate( x, filterx ), filtery );
}

template< typename T >
image< T > gaussGradientY( const image< T > &x, T sigma,
                           T maxError )
{
  // Split error-budget up for two operations.
  std::deque< T >
    filter1( gaussBlurFilter< T >( sigma, maxError / 2.0 ) ),
    filter2( gaussGradientFilter< T >( sigma, maxError / 2.0 ) );
  image< T > filterx, filtery;
  filterx.init( filter1.size(), 1 );
  filtery.init( 1, filter2.size() );
  std::copy( filter1.begin(), filter1.end(), filterx.rawData() );
  std::copy( filter2.begin(), filter2.end(), filtery.rawData() );
  return correlate( correlate( x, filterx ), filtery );
}

};