[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
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.6.0, Aug 13 2008 ) */ 00008 /* The VIGRA Website is */ 00009 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */ 00010 /* Please direct questions, bug reports, and contributions to */ 00011 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00012 /* vigra@informatik.uni-hamburg.de */ 00013 /* */ 00014 /* Permission is hereby granted, free of charge, to any person */ 00015 /* obtaining a copy of this software and associated documentation */ 00016 /* files (the "Software"), to deal in the Software without */ 00017 /* restriction, including without limitation the rights to use, */ 00018 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00019 /* sell copies of the Software, and to permit persons to whom the */ 00020 /* Software is furnished to do so, subject to the following */ 00021 /* conditions: */ 00022 /* */ 00023 /* The above copyright notice and this permission notice shall be */ 00024 /* included in all copies or substantial portions of the */ 00025 /* Software. */ 00026 /* */ 00027 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00028 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00029 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00030 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00031 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00032 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00033 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00034 /* OTHER DEALINGS IN THE SOFTWARE. */ 00035 /* */ 00036 /************************************************************************/ 00037 00038 #ifndef VIGRA_GAUSSIANS_HXX 00039 #define VIGRA_GAUSSIANS_HXX 00040 00041 #include <cmath> 00042 #include "config.hxx" 00043 #include "mathutil.hxx" 00044 #include "array_vector.hxx" 00045 #include "error.hxx" 00046 00047 namespace vigra { 00048 00049 #if 0 00050 /** \addtogroup MathFunctions Mathematical Functions 00051 */ 00052 //@{ 00053 #endif 00054 /*! The Gaussian function and its derivatives. 00055 00056 Implemented as a unary functor. Since it supports the <tt>radius()</tt> function 00057 it can also be used as a kernel in \ref resamplingConvolveImage(). 00058 00059 <b>\#include</b> <<a href="gaussians_8hxx-source.html">vigra/gaussians.hxx</a>><br> 00060 Namespace: vigra 00061 00062 \ingroup MathFunctions 00063 */ 00064 template <class T = double> 00065 class Gaussian 00066 { 00067 public: 00068 00069 /** the value type if used as a kernel in \ref resamplingConvolveImage(). 00070 */ 00071 typedef T value_type; 00072 /** the functor's argument type 00073 */ 00074 typedef T argument_type; 00075 /** the functor's result type 00076 */ 00077 typedef T result_type; 00078 00079 /** Create functor for the given standard deviation <tt>sigma</tt> and 00080 derivative order <i>n</i>. The functor then realizes the function 00081 00082 \f[ f_{\sigma,n}(x)=\frac{\partial^n}{\partial x^n} 00083 \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{x^2}{2\sigma^2}} 00084 \f] 00085 00086 Precondition: 00087 \code 00088 sigma > 0.0 00089 \endcode 00090 */ 00091 explicit Gaussian(T sigma = 1.0, unsigned int derivativeOrder = 0) 00092 : sigma_(sigma), 00093 sigma2_(-0.5 / sigma / sigma), 00094 norm_(0.0), 00095 order_(derivativeOrder), 00096 hermitePolynomial_(derivativeOrder / 2 + 1) 00097 { 00098 vigra_precondition(sigma_ > 0.0, 00099 "Gaussian::Gaussian(): sigma > 0 required."); 00100 switch(order_) 00101 { 00102 case 1: 00103 case 2: 00104 norm_ = -1.0 / (VIGRA_CSTD::sqrt(2.0 * M_PI) * sq(sigma) * sigma); 00105 break; 00106 case 3: 00107 norm_ = 1.0 / (VIGRA_CSTD::sqrt(2.0 * M_PI) * sq(sigma) * sq(sigma) * sigma); 00108 break; 00109 default: 00110 norm_ = 1.0 / VIGRA_CSTD::sqrt(2.0 * M_PI) / sigma; 00111 } 00112 calculateHermitePolynomial(); 00113 } 00114 00115 /** Function (functor) call. 00116 */ 00117 result_type operator()(argument_type x) const; 00118 00119 /** Get the standard deviation of the Gaussian. 00120 */ 00121 value_type sigma() const 00122 { return sigma_; } 00123 00124 /** Get the derivative order of the Gaussian. 00125 */ 00126 unsigned int derivativeOrder() const 00127 { return order_; } 00128 00129 /** Get the required filter radius for a discrete approximation of the Gaussian. 00130 The radius is given as a multiple of the Gaussian's standard deviation 00131 (default: <tt>sigma * (3 + 1/2 * derivativeOrder()</tt> -- the second term 00132 accounts for the fact that the derivatives of the Gaussian become wider 00133 with increasing order). The result is rounded to the next higher integer. 00134 */ 00135 double radius(double sigmaMultiple = 3.0) const 00136 { return VIGRA_CSTD::ceil(sigma_ * (sigmaMultiple + 0.5 * derivativeOrder())); } 00137 00138 private: 00139 void calculateHermitePolynomial(); 00140 T horner(T x) const; 00141 00142 T sigma_, sigma2_, norm_; 00143 unsigned int order_; 00144 ArrayVector<T> hermitePolynomial_; 00145 }; 00146 00147 template <class T> 00148 typename Gaussian<T>::result_type 00149 Gaussian<T>::operator()(argument_type x) const 00150 { 00151 T x2 = x * x; 00152 T g = norm_ * VIGRA_CSTD::exp(x2 * sigma2_); 00153 switch(order_) 00154 { 00155 case 0: 00156 return g; 00157 case 1: 00158 return x * g; 00159 case 2: 00160 return (1.0 - sq(x / sigma_)) * g; 00161 case 3: 00162 return (3.0 - sq(x / sigma_)) * x * g; 00163 default: 00164 return order_ % 2 == 0 ? 00165 g * horner(x2) 00166 : x * g * horner(x2); 00167 } 00168 } 00169 00170 template <class T> 00171 T Gaussian<T>::horner(T x) const 00172 { 00173 int i = order_ / 2; 00174 T res = hermitePolynomial_[i]; 00175 for(--i; i >= 0; --i) 00176 res = x * res + hermitePolynomial_[i]; 00177 return res; 00178 } 00179 00180 template <class T> 00181 void Gaussian<T>::calculateHermitePolynomial() 00182 { 00183 if(order_ == 0) 00184 { 00185 hermitePolynomial_[0] = 1.0; 00186 } 00187 else if(order_ == 1) 00188 { 00189 hermitePolynomial_[0] = -1.0 / sigma_ / sigma_; 00190 } 00191 else 00192 { 00193 // calculate Hermite polynomial for requested derivative 00194 // recursively according to 00195 // (0) 00196 // h (x) = 1 00197 // 00198 // (1) 00199 // h (x) = -x / s^2 00200 // 00201 // (n+1) (n) (n-1) 00202 // h (x) = -1 / s^2 * [ x * h (x) + n * h (x) ] 00203 // 00204 T s2 = -1.0 / sigma_ / sigma_; 00205 ArrayVector<T> hn(3*order_+3, 0.0); 00206 typename ArrayVector<T>::iterator hn0 = hn.begin(), 00207 hn1 = hn0 + order_+1, 00208 hn2 = hn1 + order_+1, 00209 ht; 00210 hn2[0] = 1.0; 00211 hn1[1] = s2; 00212 for(unsigned int i = 2; i <= order_; ++i) 00213 { 00214 hn0[0] = s2 * (i-1) * hn2[0]; 00215 for(unsigned int j = 1; j <= i; ++j) 00216 hn0[j] = s2 * (hn1[j-1] + (i-1) * hn2[j]); 00217 ht = hn2; 00218 hn2 = hn1; 00219 hn1 = hn0; 00220 hn0 = ht; 00221 } 00222 // keep only non-zero coefficients of the polynomial 00223 for(unsigned int i = 0; i < hermitePolynomial_.size(); ++i) 00224 hermitePolynomial_[i] = order_ % 2 == 0 ? 00225 hn1[2*i] 00226 : hn1[2*i+1]; 00227 } 00228 } 00229 00230 00231 ////@} 00232 00233 } // namespace vigra 00234 00235 00236 #endif /* VIGRA_GAUSSIANS_HXX */
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|