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

details vigra/mathutil.hxx VIGRA

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*               Copyright 1998-2005 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.3, Aug 18 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_MATHUTIL_HXX
00024 #define VIGRA_MATHUTIL_HXX
00025 
00026 #include <cmath>
00027 #include <cstdlib>
00028 #include "vigra/config.hxx"
00029 #include "vigra/numerictraits.hxx"
00030 
00031 /*! \page MathConstants Mathematical Constants
00032 
00033     <TT>M_PI, M_SQRT2</TT>
00034 
00035     <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"
00036 
00037     Since <TT>M_PI</TT> and <TT>M_SQRT2</TT> are not officially standardized,
00038     we provide definitions here for those compilers that don't support them.
00039 
00040     \code
00041     #ifndef M_PI
00042     #    define M_PI     3.14159265358979323846
00043     #endif
00044 
00045     #ifndef M_SQRT2
00046     #    define M_SQRT2  1.41421356237309504880
00047     #endif
00048     \endcode
00049 */
00050 #ifndef M_PI
00051 #    define M_PI     3.14159265358979323846
00052 #endif
00053 
00054 #ifndef M_SQRT2
00055 #    define M_SQRT2  1.41421356237309504880
00056 #endif
00057 
00058 namespace vigra {
00059 
00060 #ifndef __sun 
00061 
00062 /** \addtogroup MathFunctions Mathematical Functions
00063 
00064     Useful mathematical functions and functors.
00065 */
00066 //@{
00067     /*! The error function.
00068 
00069         If <tt>erf()</tt> is not provided in the C standard math library (as it should according to the
00070         new C99 standard ?), VIGRA implements <tt>erf()</tt> as an approximation of the error 
00071         function
00072         
00073         \f[
00074             \mbox{erf}(x) = \int_0^x e^{-x^2} dx
00075         \f]
00076         
00077         according to the formula given in Press et al. "Numerical Recipes".
00078 
00079         <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
00080         Namespace: vigra
00081     */
00082 template <class T>
00083 double erf(T x)
00084 {
00085     double t = 1.0/(1.0+0.5*VIGRA_CSTD::fabs(x));
00086     double ans = t*VIGRA_CSTD::exp(-x*x-1.26551223+t*(1.00002368+t*(0.37409196+
00087                                     t*(0.09678418+t*(-0.18628806+t*(0.27886807+
00088                                     t*(-1.13520398+t*(1.48851587+t*(-0.82215223+
00089                                     t*0.17087277)))))))));
00090     if (x >= 0.0)
00091         return 1.0 - ans;
00092     else
00093         return ans - 1.0;
00094 }
00095 
00096 #else
00097 
00098 using VIGRA_CSTD::erf;
00099 
00100 #endif
00101 
00102 // import functions into namespace vigra which VIGRA is going to overload
00103 
00104 using VIGRA_CSTD::pow;  
00105 using VIGRA_CSTD::floor;  
00106 using VIGRA_CSTD::ceil;  
00107 using std::abs;  
00108 
00109 #define VIGRA_DEFINE_UNSIGNED_ABS(T) \
00110     inline T abs(T t) { return t; }
00111 
00112 VIGRA_DEFINE_UNSIGNED_ABS(bool)
00113 VIGRA_DEFINE_UNSIGNED_ABS(unsigned char)
00114 VIGRA_DEFINE_UNSIGNED_ABS(unsigned short)
00115 VIGRA_DEFINE_UNSIGNED_ABS(unsigned int)
00116 VIGRA_DEFINE_UNSIGNED_ABS(unsigned long)
00117 
00118 #undef VIGRA_DEFINE_UNSIGNED_ABS
00119 
00120     /*! The rounding function.
00121 
00122         Defined for all floating point types. Rounds towards the nearest integer for both 
00123         positive and negative inputs.
00124 
00125         <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
00126         Namespace: vigra
00127     */
00128 inline float round(float t)
00129 {
00130     return t >= 0.0
00131                ? floor(t + 0.5)
00132                : ceil(t - 0.5);
00133 }
00134 
00135 inline double round(double t)
00136 {
00137     return t >= 0.0
00138                ? floor(t + 0.5)
00139                : ceil(t - 0.5);
00140 }
00141 
00142 inline long double round(long double t)
00143 {
00144     return t >= 0.0
00145                ? floor(t + 0.5)
00146                : ceil(t - 0.5);
00147 }
00148 
00149     /*! The square function.
00150 
00151         sq(x) is needed so often that it makes sense to define it as a function.
00152 
00153         <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
00154         Namespace: vigra
00155     */
00156 template <class T>
00157 inline 
00158 typename NumericTraits<T>::Promote sq(T t)
00159 {
00160     return t*t;
00161 }
00162 
00163 #ifdef VIGRA_NO_HYPOT
00164     /*! Compute the Euclidean distance (length of the hypothenuse of a right-angled triangle).
00165 
00166         The  hypot()  function  returns  the  sqrt(a*a  +  b*b).
00167         It is implemented in a way that minimizes round-off error.
00168 
00169         <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
00170         Namespace: vigra
00171     */
00172 template <class T>
00173 T hypot(T a, T b) 
00174 { 
00175     T absa = abs(a), absb = abs(b);
00176     if (absa > absb) 
00177         return absa * VIGRA_CSTD::sqrt(1.0 + sq(absb/absa)); 
00178     else 
00179         return absb == NumericTraits<T>::zero()
00180                    ? NumericTraits<T>::zero()
00181                    : absb * VIGRA_CSTD::sqrt(1.0 + sq(absa/absb)); 
00182 }
00183 
00184 #else
00185 
00186 using ::hypot;
00187 
00188 #endif
00189 
00190     /*! The sign function.
00191 
00192         Returns 1, 0, or -1 depending on the sign of \a t.
00193 
00194         <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
00195         Namespace: vigra
00196     */
00197 template <class T>
00198 T sign(T t) 
00199 { 
00200     return t > NumericTraits<T>::zero()
00201                ? NumericTraits<T>::one()
00202                : t < NumericTraits<T>::zero()
00203                     ? -NumericTraits<T>::one()
00204                     : NumericTraits<T>::zero();
00205 }
00206 
00207     /*! The binary sign function.
00208 
00209         Transfers the sign of \a t2 to \a t1.
00210 
00211         <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
00212         Namespace: vigra
00213     */
00214 template <class T1, class T2>
00215 T1 sign(T1 t1, T2 t2) 
00216 { 
00217     return t2 >= NumericTraits<T2>::zero()
00218                ? abs(t1)
00219                : -abs(t1);
00220 }
00221 
00222 #define VIGRA_DEFINE_NORM(T) \
00223     inline NormTraits<T>::SquaredNormType squaredNorm(T t) { return sq(t); } \
00224     inline NormTraits<T>::NormType norm(T t) { return abs(t); }
00225 
00226 VIGRA_DEFINE_NORM(bool)
00227 VIGRA_DEFINE_NORM(signed char)
00228 VIGRA_DEFINE_NORM(unsigned char)
00229 VIGRA_DEFINE_NORM(short)
00230 VIGRA_DEFINE_NORM(unsigned short)
00231 VIGRA_DEFINE_NORM(int)
00232 VIGRA_DEFINE_NORM(unsigned int)
00233 VIGRA_DEFINE_NORM(long)
00234 VIGRA_DEFINE_NORM(unsigned long)
00235 VIGRA_DEFINE_NORM(float)
00236 VIGRA_DEFINE_NORM(double)
00237 VIGRA_DEFINE_NORM(long double)
00238 
00239 #undef VIGRA_DEFINE_NORM
00240 
00241 template <class T>
00242 inline typename NormTraits<std::complex<T> >::SquaredNormType
00243 squaredNorm(std::complex<T> const & t)
00244 {
00245     return sq(t.real()) + sq(t.imag());
00246 }
00247 
00248 #ifdef DOXYGEN // only for documentation
00249     /*! The squared norm of a numerical object.
00250 
00251         For scalar types: equals <tt>vigra::sq(t)</tt><br>.
00252         For vectorial types: equals <tt>vigra::dot(t, t)</tt><br>.
00253         For complex types: equals <tt>vigra::sq(t.real()) + vigra::sq(t.imag())</tt><br>.
00254         For matrix types: results in the squared Frobenius norm (sum of squares of the matrix elements).
00255     */
00256 NormTraits<T>::SquaredNormType squaredNorm(T const & t);
00257 
00258 #endif
00259 
00260     /*! The norm of a numerical object.
00261 
00262         For scalar types: implemented as <tt>abs(t)</tt><br>
00263         otherwise: implemented as <tt>sqrt(squaredNorm(t))</tt>.
00264 
00265         <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
00266         Namespace: vigra
00267     */
00268 template <class T>
00269 inline typename NormTraits<T>::NormType 
00270 norm(T const & t)
00271 {
00272     return VIGRA_CSTD::sqrt(static_cast<typename NormTraits<T>::NormType>(squaredNorm(t)));
00273 }
00274 
00275 namespace detail {
00276 
00277 // both f1 and f2 are unsigned here
00278 template<class FPT>
00279 inline
00280 FPT safeFloatDivision( FPT f1, FPT f2 )
00281 {
00282     return  f2 < NumericTraits<FPT>::one() && f1 > f2 * NumericTraits<FPT>::max()
00283                 ? NumericTraits<FPT>::max() 
00284                 : (f2 > NumericTraits<FPT>::one() && f1 < f2 * NumericTraits<FPT>::smallestPositive()) || 
00285                    f1 == NumericTraits<FPT>::zero()
00286                      ? NumericTraits<FPT>::zero() 
00287                      : f1/f2;
00288 }
00289 
00290 } // namespace detail
00291     
00292     /*! Tolerance based floating-point comparison.
00293 
00294         Check whether two floating point numbers are equal within the given tolerance.
00295         This is useful because floating point numbers that should be equal in theory are
00296         rarely exactly equal in practice. If the tolerance \a epsilon is not given,
00297         twice the machine epsilon is used.
00298 
00299         <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
00300         Namespace: vigra
00301     */
00302 template <class T1, class T2>
00303 bool closeAtTolerance(T1 l, T2 r, typename PromoteTraits<T1, T2>::Promote epsilon)
00304 {
00305     typedef typename PromoteTraits<T1, T2>::Promote T;
00306     if(l == 0.0 && r != 0.0)
00307         return VIGRA_CSTD::fabs(r) <= epsilon;
00308     if(l != 0.0 && r == 0.0)
00309         return VIGRA_CSTD::fabs(r) <= epsilon;
00310     T diff = VIGRA_CSTD::fabs( l - r );
00311     T d1   = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( r ) );
00312     T d2   = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( l ) );
00313 
00314     return (d1 <= epsilon && d2 <= epsilon);
00315 }
00316 
00317 template <class T1, class T2>
00318 bool closeAtTolerance(T1 l, T2 r)
00319 {
00320     typedef typename PromoteTraits<T1, T2>::Promote T;
00321     return closeAtTolerance(l, r, 2.0 * NumericTraits<T>::epsilon());
00322 }
00323 
00324 //@}
00325 
00326 } // namespace vigra
00327 
00328 #endif /* VIGRA_MATHUTIL_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.3 (18 Aug 2005)