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 );
}
};