[ 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
45namespace 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*/
62template <class T = double>
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 */
116
117 /** Get the standard deviation of the Gaussian.
118 */
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
145template <class T>
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
168template <class T>
169T 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
178template <class T>
179void 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 */
Definition gaussians.hxx:64
Gaussian(T sigma=1.0, unsigned int derivativeOrder=0)
Definition gaussians.hxx:89
T value_type
Definition gaussians.hxx:69
double radius(double sigmaMultiple=3.0) const
Definition gaussians.hxx:133
value_type sigma() const
Definition gaussians.hxx:119
T result_type
Definition gaussians.hxx:75
result_type operator()(argument_type x) const
Definition gaussians.hxx:147
T argument_type
Definition gaussians.hxx:72
unsigned int derivativeOrder() const
Definition gaussians.hxx:124
Class for a single RGB value.
Definition rgbvalue.hxx:128
RGBValue()
Definition rgbvalue.hxx:209

© 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.11.2