[ 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) Cognitive Systems Group, University of Hamburg, Germany | html generated using doxygen and Python VIGRA 1.3.2 (27 Jan 2005)