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