[ VIGRA Homepage | Class Index | Function Index | File Index | Main Page ]

details vigra/gaussians.hxx VIGRA

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.5.0, Dec 07 2006 )                                    */
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 /*        koethe@informatik.uni-hamburg.de          or                  */
00012 /*        vigra@kogs1.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 (koethe@informatik.uni-hamburg.de)
Cognitive Systems Group, University of Hamburg, Germany

html generated using doxygen and Python
VIGRA 1.5.0 (7 Dec 2006)