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