[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

gaussians.hxx |

1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2004 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35
36 #ifndef VIGRA_GAUSSIANS_HXX
37 #define VIGRA_GAUSSIANS_HXX
38
39 #include <cmath>
40 #include "config.hxx"
41 #include "mathutil.hxx"
42 #include "array_vector.hxx"
43 #include "error.hxx"
44
45 namespace vigra {
46
47 #if 0
48 /** \addtogroup MathFunctions Mathematical Functions
49 */
50 //@{
51 #endif
52 /** The Gaussian function and its derivatives.
53
54 Implemented as a unary functor. Since it supports the <tt>radius()</tt> function
55 it can also be used as a kernel in \ref resamplingConvolveImage().
56
57 <b>\#include</b> <vigra/gaussians.hxx><br>
58 Namespace: vigra
59
60 \ingroup MathFunctions
61 */
62 template <class T = double>
63 class Gaussian
64 {
65 public:
66
67 /** the value type if used as a kernel in \ref resamplingConvolveImage().
68 */
69 typedef T value_type;
70 /** the functor's argument type
71 */
72 typedef T argument_type;
73 /** the functor's result type
74 */
75 typedef T result_type;
76
77 /** Create functor for the given standard deviation <tt>sigma</tt> and
78 derivative order <i>n</i>. The functor then realizes the function
79
80 \f[ f_{\sigma,n}(x)=\frac{\partial^n}{\partial x^n}
81 \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{x^2}{2\sigma^2}}
82 \f]
83
84 Precondition:
85 \code
86 sigma > 0.0
87 \endcode
88 */
89 explicit Gaussian(T sigma = 1.0, unsigned int derivativeOrder = 0)
90 : sigma_(sigma),
91 sigma2_(T(-0.5 / sigma / sigma)),
92 norm_(0.0),
93 order_(derivativeOrder),
94 hermitePolynomial_(derivativeOrder / 2 + 1)
95 {
96 vigra_precondition(sigma_ > 0.0,
97 "Gaussian::Gaussian(): sigma > 0 required.");
98 switch(order_)
99 {
100 case 1:
101 case 2:
102 norm_ = T(-1.0 / (VIGRA_CSTD::sqrt(2.0 * M_PI) * sq(sigma) * sigma));
103 break;
104 case 3:
105 norm_ = T(1.0 / (VIGRA_CSTD::sqrt(2.0 * M_PI) * sq(sigma) * sq(sigma) * sigma));
106 break;
107 default:
108 norm_ = T(1.0 / VIGRA_CSTD::sqrt(2.0 * M_PI) / sigma);
109 }
110 calculateHermitePolynomial();
111 }
112
113 /** Function (functor) call.
114 */
115 result_type operator()(argument_type x) const;
116
117 /** Get the standard deviation of the Gaussian.
118 */
119 value_type sigma() const
120 { return sigma_; }
121
122 /** Get the derivative order of the Gaussian.
123 */
124 unsigned int derivativeOrder() const
125 { return order_; }
126
127 /** Get the required filter radius for a discrete approximation of the Gaussian.
128 The radius is given as a multiple of the Gaussian's standard deviation
129 (default: <tt>sigma * (3 + 1/2 * derivativeOrder()</tt> -- the second term
130 accounts for the fact that the derivatives of the Gaussian become wider
131 with increasing order). The result is rounded to the next higher integer.
132 */
133 double radius(double sigmaMultiple = 3.0) const
134 { return VIGRA_CSTD::ceil(sigma_ * (sigmaMultiple + 0.5 * derivativeOrder())); }
135
136 private:
137 void calculateHermitePolynomial();
138 T horner(T x) const;
139
140 T sigma_, sigma2_, norm_;
141 unsigned int order_;
142 ArrayVector<T> hermitePolynomial_;
143 };
144
145 template <class T>
146 typename Gaussian<T>::result_type
147 Gaussian<T>::operator()(argument_type x) const
148 {
149 T x2 = x * x;
150 T g = norm_ * VIGRA_CSTD::exp(x2 * sigma2_);
151 switch(order_)
152 {
153 case 0:
154 return detail::RequiresExplicitCast<result_type>::cast(g);
155 case 1:
156 return detail::RequiresExplicitCast<result_type>::cast(x * g);
157 case 2:
158 return detail::RequiresExplicitCast<result_type>::cast((1.0 - sq(x / sigma_)) * g);
159 case 3:
160 return detail::RequiresExplicitCast<result_type>::cast((3.0 - sq(x / sigma_)) * x * g);
161 default:
162 return order_ % 2 == 0 ?
163 detail::RequiresExplicitCast<result_type>::cast(g * horner(x2))
164 : detail::RequiresExplicitCast<result_type>::cast(x * g * horner(x2));
165 }
166 }
167
168 template <class T>
169 T Gaussian<T>::horner(T x) const
170 {
171 int i = order_ / 2;
172 T res = hermitePolynomial_[i];
173 for(--i; i >= 0; --i)
174 res = x * res + hermitePolynomial_[i];
175 return res;
176 }
177
178 template <class T>
179 void Gaussian<T>::calculateHermitePolynomial()
180 {
181 if(order_ == 0)
182 {
183 hermitePolynomial_[0] = 1.0;
184 }
185 else if(order_ == 1)
186 {
187 hermitePolynomial_[0] = T(-1.0 / sigma_ / sigma_);
188 }
189 else
190 {
191 // calculate Hermite polynomial for requested derivative
192 // recursively according to
193 // (0)
194 // h (x) = 1
195 //
196 // (1)
197 // h (x) = -x / s^2
198 //
199 // (n+1) (n) (n-1)
200 // h (x) = -1 / s^2 * [ x * h (x) + n * h (x) ]
201 //
202 T s2 = T(-1.0 / sigma_ / sigma_);
203 ArrayVector<T> hn(3*order_+3, 0.0);
204 typename ArrayVector<T>::iterator hn0 = hn.begin(),
205 hn1 = hn0 + order_+1,
206 hn2 = hn1 + order_+1,
207 ht;
208 hn2[0] = 1.0;
209 hn1[1] = s2;
210 for(unsigned int i = 2; i <= order_; ++i)
211 {
212 hn0[0] = s2 * (i-1) * hn2[0];
213 for(unsigned int j = 1; j <= i; ++j)
214 hn0[j] = s2 * (hn1[j-1] + (i-1) * hn2[j]);
215 ht = hn2;
216 hn2 = hn1;
217 hn1 = hn0;
218 hn0 = ht;
219 }
220 // keep only non-zero coefficients of the polynomial
221 for(unsigned int i = 0; i < hermitePolynomial_.size(); ++i)
222 hermitePolynomial_[i] = order_ % 2 == 0 ?
223 hn1[2*i]
224 : hn1[2*i+1];
225 }
226 }
227
228
229 ////@}
230
231 } // namespace vigra
232
233
234 #endif /* VIGRA_GAUSSIANS_HXX */
vigra::Gaussian::radius
double radius(double sigmaMultiple=3.0) const
Definition: gaussians.hxx:133
vigra::Gaussian::argument_type
T argument_type
Definition: gaussians.hxx:72
vigra::linalg::exp
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
vigra::ArrayVectorView::begin
const_iterator begin() const
Definition: array_vector.hxx:223
vigra
Definition: accessor.hxx:43
vigra::Gaussian::sigma
value_type sigma() const
Definition: gaussians.hxx:119
vigra::ArrayVector< T >
vigra::Gaussian
Definition: gaussians.hxx:63
vigra::sq
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:344
vigra::Gaussian::Gaussian
Gaussian(T sigma=1.0, unsigned int derivativeOrder=0)
Definition: gaussians.hxx:89
vigra::ArrayVectorView::size
size_type size() const
Definition: array_vector.hxx:330
vigra::Gaussian::result_type
T result_type
Definition: gaussians.hxx:75
vigra::ceil
int ceil(FixedPoint< IntBits, FracBits > v)
rounding up.
Definition: fixedpoint.hxx:675
vigra::Gaussian::value_type
T value_type
Definition: gaussians.hxx:69
vigra::Gaussian::operator()
result_type operator()(argument_type x) const
Definition: gaussians.hxx:147
vigra::sqrt
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616
vigra::Gaussian::derivativeOrder
unsigned int derivativeOrder() const
Definition: gaussians.hxx:124

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.10.0