/************************************************************************/ /* */ /* Copyright 1998-2004 by Ullrich Koethe */ /* Cognitive Systems Group, University of Hamburg, Germany */ /* */ /* This file is part of the VIGRA computer vision library. */ /* The VIGRA Website is */ /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */ /* Please direct questions, bug reports, and contributions to */ /* koethe@informatik.uni-hamburg.de or */ /* vigra@kogs1.informatik.uni-hamburg.de */ /* */ /* Permission is hereby granted, free of charge, to any person */ /* obtaining a copy of this software and associated documentation */ /* files (the "Software"), to deal in the Software without */ /* restriction, including without limitation the rights to use, */ /* copy, modify, merge, publish, distribute, sublicense, and/or */ /* sell copies of the Software, and to permit persons to whom the */ /* Software is furnished to do so, subject to the following */ /* conditions: */ /* */ /* The above copyright notice and this permission notice shall be */ /* included in all copies or substantial portions of the */ /* Software. */ /* */ /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ /* OTHER DEALINGS IN THE SOFTWARE. */ /* */ /************************************************************************/ #ifndef VIGRA_GAUSSIANS_HXX #define VIGRA_GAUSSIANS_HXX #include #include "vigra/config.hxx" #include "vigra/mathutil.hxx" #include "vigra/array_vector.hxx" #include "vigra/error.hxx" namespace vigra { #if 0 /** \addtogroup MathFunctions Mathematical Functions */ //@{ #endif /*! The Gaussian function and its derivatives. Implemented as a unary functor. Since it supports the radius() function it can also be used as a kernel in \ref resamplingConvolveImage(). \#include "vigra/gaussians.hxx"
Namespace: vigra \ingroup MathFunctions */ template class Gaussian { public: /** the value type if used as a kernel in \ref resamplingConvolveImage(). */ typedef T value_type; /** the functor's argument type */ typedef T argument_type; /** the functor's result type */ typedef T result_type; /** Create functor for the given standard deviation sigma and derivative order n. The functor then realizes the function \f[ f_{\sigma,n}(x)=\frac{\partial^n}{\partial x^n} \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{x^2}{2\sigma^2}} \f] Precondition: \code sigma > 0.0 \endcode */ explicit Gaussian(T sigma = 1.0, unsigned int derivativeOrder = 0) : sigma_(sigma), sigma2_(-0.5 / sigma / sigma), norm_(0.0), order_(derivativeOrder), hermitePolynomial_(derivativeOrder / 2 + 1) { vigra_precondition(sigma_ > 0.0, "Gaussian::Gaussian(): sigma > 0 required."); switch(order_) { case 1: case 2: norm_ = -1.0 / (VIGRA_CSTD::sqrt(2.0 * M_PI) * sq(sigma) * sigma); break; case 3: norm_ = 1.0 / (VIGRA_CSTD::sqrt(2.0 * M_PI) * sq(sigma) * sq(sigma) * sigma); break; default: norm_ = 1.0 / VIGRA_CSTD::sqrt(2.0 * M_PI) / sigma; } calculateHermitePolynomial(); } /** Function (functor) call. */ result_type operator()(argument_type x) const; /** Get the standard deviation of the Gaussian. */ value_type sigma() const { return sigma_; } /** Get the derivative order of the Gaussian. */ unsigned int derivativeOrder() const { return order_; } /** Get the required filter radius for a discrete approximation of the Gaussian. The radius is given as a multiple of the Gaussian's standard deviation (default: sigma * (3 + 1/2 * derivativeOrder() -- the second term accounts for the fact that the derivatives of the Gaussian become wider with increasing order). The result is rounded to the next higher integer. */ double radius(double sigmaMultiple = 3.0) const { return VIGRA_CSTD::ceil(sigma_ * (sigmaMultiple + 0.5 * derivativeOrder())); } private: void calculateHermitePolynomial(); T horner(T x) const; T sigma_, sigma2_, norm_; unsigned int order_; ArrayVector hermitePolynomial_; }; template typename Gaussian::result_type Gaussian::operator()(argument_type x) const { T x2 = x * x; T g = norm_ * VIGRA_CSTD::exp(x2 * sigma2_); switch(order_) { case 0: return g; case 1: return x * g; case 2: return (1.0 - sq(x / sigma_)) * g; case 3: return (3.0 - sq(x / sigma_)) * x * g; default: return order_ % 2 == 0 ? g * horner(x2) : x * g * horner(x2); } } template T Gaussian::horner(T x) const { int i = order_ / 2; T res = hermitePolynomial_[i]; for(--i; i >= 0; --i) res = x * res + hermitePolynomial_[i]; return res; } template void Gaussian::calculateHermitePolynomial() { if(order_ == 0) { hermitePolynomial_[0] = 1.0; } else if(order_ == 1) { hermitePolynomial_[0] = -1.0 / sigma_ / sigma_; } else { // calculate Hermite polynomial for requested derivative // recursively according to // (0) // h (x) = 1 // // (1) // h (x) = -x / s^2 // // (n+1) (n) (n-1) // h (x) = -1 / s^2 * [ x * h (x) + n * h (x) ] // T s2 = -1.0 / sigma_ / sigma_; ArrayVector hn(3*order_+3, 0.0); typename ArrayVector::iterator hn0 = hn.begin(), hn1 = hn0 + order_+1, hn2 = hn1 + order_+1, ht; hn2[0] = 1.0; hn1[1] = s2; for(unsigned int i = 2; i <= order_; ++i) { hn0[0] = s2 * (i-1) * hn2[0]; for(unsigned int j = 1; j <= i; ++j) hn0[j] = s2 * (hn1[j-1] + (i-1) * hn2[j]); ht = hn2; hn2 = hn1; hn1 = hn0; hn0 = ht; } // keep only non-zero coefficients of the polynomial for(unsigned int i = 0; i < hermitePolynomial_.size(); ++i) hermitePolynomial_[i] = order_ % 2 == 0 ? hn1[2*i] : hn1[2*i+1]; } } ////@} } // namespace vigra #endif /* VIGRA_GAUSSIANS_HXX */