Rev Author Line No. Line
178 kaklik 1 namespace mimas {
2  
3 template< typename T >
4 std::deque< T > gaussBlurFilter( T sigma, T maxError )
5 {
6 assert( sigma > 0 );
7  
8 std::deque< T > retval;
9  
10 // Integral of gauss-curve from -0.5 to +0.5.
11 T integral = erf( 0.5 / ( sqrt( 2.0 ) * sigma ) );
12 retval.push_back( integral );
13  
14 while ( 1.0 - integral > maxError ) {
15 // Compute integral of gauss-curve from -newSize2 to +newSize2.
16 const T
17 newSize2 = 0.5 * ( retval.size() + 2 ),
18 newIntegral = erf( newSize2 / ( sqrt( 2.0 ) * sigma ) ),
19 value = 0.5 * ( newIntegral - integral );
20 retval.push_back( value );
21 retval.push_front( value );
22 integral = newIntegral;
23 };
24  
25 // Normalise to maintain power of image.
26 const T factor = 1.0 / integral;
27 for ( typename std::deque< T >::iterator i = retval.begin();
28 i != retval.end(); i++ )
29 *i *= factor;
30  
31 return retval;
32 }
33  
34 template< typename T >
35 boost::multi_array< T, 2 > gaussBlur
36 ( const boost::const_multi_array_ref< T, 2 > &x, T sigma, T maxError )
37 {
38 // Split error-budget up for two operations.
39 std::deque< T > filter( gaussBlurFilter< T >( sigma, maxError / 2.0 ) );
40 boost::multi_array< T, 2 >
41 filterx( boost::extents[ 1 ][ filter.size() ] ),
42 filtery( boost::extents[ filter.size()][ 1 ] );
43 std::copy( filter.begin(), filter.end(), filterx.data() );
44 std::copy( filter.begin(), filter.end(), filtery.data() );
45 return correlate( correlate( x, filterx ), filtery );
46 }
47  
48 template< typename T >
49 std::deque< T > gaussGradientFilter( T sigma, T maxError )
50 {
51 assert( sigma > 0 );
52  
53 std::deque< T > retval;
54  
55 // Integral of gauss-gradient from -0.5 to +0.5 is zero.
56 retval.push_back( 0.0 );
57  
58 // Set to absolute value of integral of gauss-gradient from 0.5 to infinity
59 T integral =
60 1.0 / ( sqrt( 2.0 * M_PI ) * sigma ) * exp( -0.125 / ( sigma * sigma ) ),
61 squareIntegral = 0.0;
62  
63 while ( 2.0 * integral > maxError ) {
64 // Compute integral of gauss-curve from -newSize2 to +newSize2.
65 const T
66 newSize2 = 0.5 * ( retval.size() + 2 ),
67 newIntegral = 1.0 / ( sqrt( 2.0 * M_PI ) * sigma ) * exp( -newSize2 * newSize2 / ( 2.0 * sigma * sigma ) ),
68 value = integral - newIntegral;
69 retval.push_back( value );
70 retval.push_front( -value );
71 integral = newIntegral;
72 squareIntegral += value * value;
73 };
74  
75 // Normalise to maintain power of image.
76 const T factor = 1.0 / squareIntegral;
77 for ( typename std::deque< T >::iterator i = retval.begin();
78 i != retval.end(); i++ )
79 *i *= factor;
80  
81 return retval;
82 }
83  
84 template< typename T >
85 boost::multi_array< T, 2 > gaussGradientX
86 ( const boost::const_multi_array_ref< T, 2 > &x, T sigma, T maxError )
87 {
88 // Split error-budget up for two operations.
89 std::deque< T >
90 filter1( gaussGradientFilter< T >( sigma, maxError / 2.0 ) ),
91 filter2( gaussBlurFilter< T >( sigma, maxError / 2.0 ) );
92 boost::multi_array< T, 2 >
93 filterx( boost::extents[ 1 ][ filter1.size() ] ),
94 filtery( boost::extents[ filter2.size()][ 1 ] );
95 std::copy( filter1.begin(), filter1.end(), filterx.data() );
96 std::copy( filter2.begin(), filter2.end(), filtery.data() );
97 return correlate( correlate( x, filterx ), filtery );
98 }
99  
100 template< typename T >
101 image< T > gaussGradientX( const image< T > &x, T sigma,
102 T maxError )
103 {
104 // Split error-budget up for two operations.
105 std::deque< T >
106 filter1( gaussGradientFilter< T >( sigma, maxError / 2.0 ) ),
107 filter2( gaussBlurFilter< T >( sigma, maxError / 2.0 ) );
108 image< T > filterx, filtery;
109 filterx.init( filter1.size(), 1 );
110 filtery.init( 1, filter2.size() );
111 std::copy( filter1.begin(), filter1.end(), filterx.rawData() );
112 std::copy( filter2.begin(), filter2.end(), filtery.rawData() );
113 return correlate( correlate( x, filterx ), filtery );
114 }
115  
116  
117 template< typename T >
118 boost::multi_array< T, 2 > gaussGradientY
119 ( const boost::const_multi_array_ref< T, 2 > &x, T sigma, T maxError )
120 {
121 // Split error-budget up for two operations.
122 std::deque< T >
123 filter1( gaussBlurFilter< T >( sigma, maxError / 2.0 ) ),
124 filter2( gaussGradientFilter< T >( sigma, maxError / 2.0 ) );
125 boost::multi_array< T, 2 >
126 filterx( boost::extents[ 1 ][ filter1.size() ] ),
127 filtery( boost::extents[ filter2.size()][ 1 ] );
128 std::copy( filter1.begin(), filter1.end(), filterx.data() );
129 std::copy( filter2.begin(), filter2.end(), filtery.data() );
130 return correlate( correlate( x, filterx ), filtery );
131 }
132  
133 template< typename T >
134 image< T > gaussGradientY( const image< T > &x, T sigma,
135 T maxError )
136 {
137 // Split error-budget up for two operations.
138 std::deque< T >
139 filter1( gaussBlurFilter< T >( sigma, maxError / 2.0 ) ),
140 filter2( gaussGradientFilter< T >( sigma, maxError / 2.0 ) );
141 image< T > filterx, filtery;
142 filterx.init( filter1.size(), 1 );
143 filtery.init( 1, filter2.size() );
144 std::copy( filter1.begin(), filter1.end(), filterx.rawData() );
145 std::copy( filter2.begin(), filter2.end(), filtery.rawData() );
146 return correlate( correlate( x, filterx ), filtery );
147 }
148  
149 };