vigra - vigra/gaussians.hxx Source File

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

vigra/gaussians.hxx |

    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> <vigra/gaussians.hxx><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) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.9.0 (Tue Nov 6 2012)