[ VIGRA Homepage | Class Index | Function Index | File Index | Main Page ]
![]() |
vigra/gaussians.hxx | ![]() |
---|
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 1998-2004 by Ullrich Koethe */ 00004 /* Cognitive Systems Group, University of Hamburg, Germany */ 00005 /* */ 00006 /* This file is part of the VIGRA computer vision library. */ 00007 /* ( Version 1.3.2, Jan 27 2005 ) */ 00008 /* You may use, modify, and distribute this software according */ 00009 /* to the terms stated in the LICENSE file included in */ 00010 /* the VIGRA distribution. */ 00011 /* */ 00012 /* The VIGRA Website is */ 00013 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */ 00014 /* Please direct questions, bug reports, and contributions to */ 00015 /* koethe@informatik.uni-hamburg.de */ 00016 /* */ 00017 /* THIS SOFTWARE IS PROVIDED AS IS AND WITHOUT ANY EXPRESS OR */ 00018 /* IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED */ 00019 /* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */ 00020 /* */ 00021 /************************************************************************/ 00022 00023 #ifndef VIGRA_GAUSSIANS_HXX 00024 #define VIGRA_GAUSSIANS_HXX 00025 00026 #include <cmath> 00027 #include "vigra/config.hxx" 00028 #include "vigra/mathutil.hxx" 00029 #include "vigra/array_vector.hxx" 00030 00031 namespace vigra { 00032 00033 #if 0 00034 /** \addtogroup MathFunctions Mathematical Functions 00035 */ 00036 //@{ 00037 #endif 00038 /*! The Gaussian function and its derivatives. 00039 00040 Implemented as a unary functor. Since it supports the <tt>radius()</tt> function 00041 it can also be used as a kernel in \ref resamplingConvolveImage(). 00042 00043 <b>\#include</b> "<a href="gaussians_8hxx-source.html">vigra/gaussians.hxx</a>"<br> 00044 Namespace: vigra 00045 00046 \ingroup MathFunctions 00047 */ 00048 template <class T = double> 00049 class Gaussian 00050 { 00051 public: 00052 00053 /** the value type if used as a kernel in \ref resamplingConvolveImage(). 00054 */ 00055 typedef T value_type; 00056 /** the functor's argument type 00057 */ 00058 typedef T argument_type; 00059 /** the functor's result type 00060 */ 00061 typedef T result_type; 00062 00063 /** Create functor for the given standard deviation <tt>sigma</tt> and 00064 derivative order <i>n</i>. The functor then realizes the function 00065 00066 \f[ f_{\sigma,n}(x)=\frac{\partial^n}{\partial x^n} 00067 \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{x^2}{2\sigma^2}} 00068 \f] 00069 00070 Precondition: 00071 \code 00072 sigma > 0.0 00073 \endcode 00074 */ 00075 explicit Gaussian(T sigma = 1.0, unsigned int derivativeOrder = 0) 00076 : sigma_(sigma), 00077 sigma2_(-0.5 / sigma / sigma), 00078 norm_(0.0), 00079 order_(derivativeOrder), 00080 hermitePolynomial_(derivativeOrder / 2 + 1) 00081 { 00082 vigra_precondition(sigma_ > 0.0, 00083 "Gaussian::Gaussian(): sigma > 0 required."); 00084 switch(order_) 00085 { 00086 case 1: 00087 case 2: 00088 norm_ = -1.0 / (VIGRA_CSTD::sqrt(2.0 * M_PI) * sq(sigma) * sigma); 00089 break; 00090 case 3: 00091 norm_ = 1.0 / (VIGRA_CSTD::sqrt(2.0 * M_PI) * sq(sigma) * sq(sigma) * sigma); 00092 break; 00093 default: 00094 norm_ = 1.0 / VIGRA_CSTD::sqrt(2.0 * M_PI) / sigma; 00095 } 00096 calculateHermitePolynomial(); 00097 } 00098 00099 /** Function (functor) call. 00100 */ 00101 result_type operator()(argument_type x) const; 00102 00103 /** Get the standard deviation of the Gaussian. 00104 */ 00105 value_type sigma() const 00106 { return sigma_; } 00107 00108 /** Get the derivative order of the Gaussian. 00109 */ 00110 unsigned int derivativeOrder() const 00111 { return order_; } 00112 00113 /** Get the required filter radius for a discrete approximation of the Gaussian. 00114 The radius is given as a multiple of the Gaussian's standard deviation 00115 (default: <tt>sigma * (3 + 1/2 * derivativeOrder()</tt> -- the second term 00116 accounts for the fact that the derivatives of the Gaussian become wider 00117 with increasing order). The result is rounded to the next higher integer. 00118 */ 00119 double radius(double sigmaMultiple = 3.0) const 00120 { return VIGRA_CSTD::ceil(sigma_ * (sigmaMultiple + 0.5 * derivativeOrder())); } 00121 00122 private: 00123 void calculateHermitePolynomial(); 00124 T horner(T x) const; 00125 00126 T sigma_, sigma2_, norm_; 00127 unsigned int order_; 00128 ArrayVector<T> hermitePolynomial_; 00129 }; 00130 00131 template <class T> 00132 typename Gaussian<T>::result_type 00133 Gaussian<T>::operator()(argument_type x) const 00134 { 00135 T x2 = x * x; 00136 T g = norm_ * VIGRA_CSTD::exp(x2 * sigma2_); 00137 switch(order_) 00138 { 00139 case 0: 00140 return g; 00141 case 1: 00142 return x * g; 00143 case 2: 00144 return (1.0 - sq(x / sigma_)) * g; 00145 case 3: 00146 return (3.0 - sq(x / sigma_)) * x * g; 00147 default: 00148 return order_ % 2 == 0 ? 00149 g * horner(x2) 00150 : x * g * horner(x2); 00151 } 00152 } 00153 00154 template <class T> 00155 T Gaussian<T>::horner(T x) const 00156 { 00157 int i = order_ / 2; 00158 T res = hermitePolynomial_[i]; 00159 for(--i; i >= 0; --i) 00160 res = x * res + hermitePolynomial_[i]; 00161 return res; 00162 } 00163 00164 template <class T> 00165 void Gaussian<T>::calculateHermitePolynomial() 00166 { 00167 if(order_ == 0) 00168 { 00169 hermitePolynomial_[0] = 1.0; 00170 } 00171 else if(order_ == 1) 00172 { 00173 hermitePolynomial_[0] = -1.0 / sigma_ / sigma_; 00174 } 00175 else 00176 { 00177 // calculate Hermite polynomial for requested derivative 00178 // recursively according to 00179 // (0) 00180 // h (x) = 1 00181 // 00182 // (1) 00183 // h (x) = -x / s^2 00184 // 00185 // (n+1) (n) (n-1) 00186 // h (x) = -1 / s^2 * [ x * h (x) + n * h (x) ] 00187 // 00188 T s2 = -1.0 / sigma_ / sigma_; 00189 ArrayVector<T> hn(3*order_+3, 0.0); 00190 typename ArrayVector<T>::iterator hn0 = hn.begin(), 00191 hn1 = hn0 + order_+1, 00192 hn2 = hn1 + order_+1, 00193 ht; 00194 hn2[0] = 1.0; 00195 hn1[1] = s2; 00196 for(unsigned int i = 2; i <= order_; ++i) 00197 { 00198 hn0[0] = s2 * (i-1) * hn2[0]; 00199 for(unsigned int j = 1; j <= i; ++j) 00200 hn0[j] = s2 * (hn1[j-1] + (i-1) * hn2[j]); 00201 ht = hn2; 00202 hn2 = hn1; 00203 hn1 = hn0; 00204 hn0 = ht; 00205 } 00206 // keep only non-zero coefficients of the polynomial 00207 for(unsigned int i = 0; i < hermitePolynomial_.size(); ++i) 00208 hermitePolynomial_[i] = order_ % 2 == 0 ? 00209 hn1[2*i] 00210 : hn1[2*i+1]; 00211 } 00212 } 00213 00214 00215 ////@} 00216 00217 } // namespace vigra 00218 00219 00220 #endif /* VIGRA_GAUSSIANS_HXX */
© Ullrich Köthe (koethe@informatik.uni-hamburg.de) |
html generated using doxygen and Python
|