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

| vigra/mathutil.hxx |
--- | --- | ---

    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.4.0, Dec 21 2005 )                                    */
    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_MATHUTIL_HXX
    00039 #define VIGRA_MATHUTIL_HXX
    00040
    00041 #include <cmath>
    00042 #include <cstdlib>
    00043 #include "vigra/config.hxx"
    00044 #include "vigra/sized_int.hxx"
    00045 #include "vigra/numerictraits.hxx"
    00046
    00047 /*! \page MathConstants Mathematical Constants
    00048
    00049     <TT>M_PI, M_SQRT2</TT>
    00050
    00051     <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"
    00052
    00053     Since <TT>M_PI</TT> and <TT>M_SQRT2</TT> are not officially standardized,
    00054     we provide definitions here for those compilers that don't support them.
    00055
    00056     \code
    00057     #ifndef M_PI
    00058     #    define M_PI     3.14159265358979323846
    00059     #endif
    00060
    00061     #ifndef M_SQRT2
    00062     #    define M_SQRT2  1.41421356237309504880
    00063     #endif
    00064     \endcode
    00065 */
    00066 #ifndef M_PI
    00067 #    define M_PI     3.14159265358979323846
    00068 #endif
    00069
    00070 #ifndef M_SQRT2
    00071 #    define M_SQRT2  1.41421356237309504880
    00072 #endif
    00073
    00074 namespace vigra {
    00075
    00076 #ifndef __sun
    00077
    00078 /** \addtogroup MathFunctions Mathematical Functions
    00079
    00080     Useful mathematical functions and functors.
    00081 */
    00082 //@{
    00083     /*! The error function.
    00084
    00085         If <tt>erf()</tt> is not provided in the C standard math library (as it should according to the
    00086         new C99 standard ?), VIGRA implements <tt>erf()</tt> as an approximation of the error
    00087         function
    00088
    00089         \f[
    00090             \mbox{erf}(x) = \int_0^x e^{-x^2} dx
    00091         \f]
    00092
    00093         according to the formula given in Press et al. "Numerical Recipes".
    00094
    00095         <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
    00096         Namespace: vigra
    00097     */
    00098 template <class T>
    00099 double erf(T x)
    00100 {
    00101     double t = 1.0/(1.0+0.5*VIGRA_CSTD::fabs(x));
    00102     double ans = t*VIGRA_CSTD::exp(-x*x-1.26551223+t*(1.00002368+t*(0.37409196+
    00103                                     t*(0.09678418+t*(-0.18628806+t*(0.27886807+
    00104                                     t*(-1.13520398+t*(1.48851587+t*(-0.82215223+
    00105                                     t*0.17087277)))))))));
    00106     if (x >= 0.0)
    00107         return 1.0 - ans;
    00108     else
    00109         return ans - 1.0;
    00110 }
    00111
    00112 #else
    00113
    00114 using VIGRA_CSTD::erf;
    00115
    00116 #endif
    00117
    00118 // import functions into namespace vigra which VIGRA is going to overload
    00119
    00120 using VIGRA_CSTD::pow;
    00121 using VIGRA_CSTD::floor;
    00122 using VIGRA_CSTD::ceil;
    00123 using std::abs;
    00124
    00125 #define VIGRA_DEFINE_UNSIGNED_ABS(T) \
    00126     inline T abs(T t) { return t; }
    00127
    00128 VIGRA_DEFINE_UNSIGNED_ABS(bool)
    00129 VIGRA_DEFINE_UNSIGNED_ABS(unsigned char)
    00130 VIGRA_DEFINE_UNSIGNED_ABS(unsigned short)
    00131 VIGRA_DEFINE_UNSIGNED_ABS(unsigned int)
    00132 VIGRA_DEFINE_UNSIGNED_ABS(unsigned long)
    00133
    00134 #undef VIGRA_DEFINE_UNSIGNED_ABS
    00135
    00136     /*! The rounding function.
    00137
    00138         Defined for all floating point types. Rounds towards the nearest integer for both
    00139         positive and negative inputs.
    00140
    00141         <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
    00142         Namespace: vigra
    00143     */
    00144 inline float round(float t)
    00145 {
    00146     return t >= 0.0
    00147                ? floor(t + 0.5)
    00148                : ceil(t - 0.5);
    00149 }
    00150
    00151 inline double round(double t)
    00152 {
    00153     return t >= 0.0
    00154                ? floor(t + 0.5)
    00155                : ceil(t - 0.5);
    00156 }
    00157
    00158 inline long double round(long double t)
    00159 {
    00160     return t >= 0.0
    00161                ? floor(t + 0.5)
    00162                : ceil(t - 0.5);
    00163 }
    00164
    00165     /*! The square function.
    00166
    00167         sq(x) is needed so often that it makes sense to define it as a function.
    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 inline
    00174 typename NumericTraits<T>::Promote sq(T t)
    00175 {
    00176     return t*t;
    00177 }
    00178
    00179 namespace detail {
    00180
    00181 template <class T>
    00182 class IntSquareRoot
    00183 {
    00184   public:
    00185     static int sqq_table[];
    00186     static UInt32 exec(UInt32 v);
    00187 };
    00188
    00189 template <class T>
    00190 int IntSquareRoot<T>::sqq_table[] = {
    00191            0,  16,  22,  27,  32,  35,  39,  42,  45,  48,  50,  53,  55,  57,
    00192           59,  61,  64,  65,  67,  69,  71,  73,  75,  76,  78,  80,  81,  83,
    00193           84,  86,  87,  89,  90,  91,  93,  94,  96,  97,  98,  99, 101, 102,
    00194          103, 104, 106, 107, 108, 109, 110, 112, 113, 114, 115, 116, 117, 118,
    00195          119, 120, 121, 122, 123, 124, 125, 126, 128, 128, 129, 130, 131, 132,
    00196          133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 144, 145,
    00197          146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155, 156, 157,
    00198          158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168,
    00199          169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178,
    00200          179, 180, 181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188,
    00201          189, 189, 190, 191, 192, 192, 193, 193, 194, 195, 195, 196, 197, 197,
    00202          198, 199, 199, 200, 201, 201, 202, 203, 203, 204, 204, 205, 206, 206,
    00203          207, 208, 208, 209, 209, 210, 211, 211, 212, 212, 213, 214, 214, 215,
    00204          215, 216, 217, 217, 218, 218, 219, 219, 220, 221, 221, 222, 222, 223,
    00205          224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230, 230, 231,
    00206          231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238,
    00207          239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246,
    00208          246, 247, 247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253,
    00209          253, 254, 254, 255
    00210 };
    00211
    00212 template <class T>
    00213 UInt32 IntSquareRoot<T>::exec(UInt32 x)
    00214 {
    00215     unsigned long xn;
    00216     if (x >= 0x10000)
    00217         if (x >= 0x1000000)
    00218             if (x >= 0x10000000)
    00219                 if (x >= 0x40000000) {
    00220                     if (x >= (UInt32)65535*(UInt32)65535)
    00221                         return 65535;
    00222                     xn = sqq_table[x>>24] << 8;
    00223                 } else
    00224                     xn = sqq_table[x>>22] << 7;
    00225             else
    00226                 if (x >= 0x4000000)
    00227                     xn = sqq_table[x>>20] << 6;
    00228                 else
    00229                     xn = sqq_table[x>>18] << 5;
    00230         else {
    00231             if (x >= 0x100000)
    00232                 if (x >= 0x400000)
    00233                     xn = sqq_table[x>>16] << 4;
    00234                 else
    00235                     xn = sqq_table[x>>14] << 3;
    00236             else
    00237                 if (x >= 0x40000)
    00238                     xn = sqq_table[x>>12] << 2;
    00239                 else
    00240                     xn = sqq_table[x>>10] << 1;
    00241
    00242             goto nr1;
    00243         }
    00244     else
    00245         if (x >= 0x100) {
    00246             if (x >= 0x1000)
    00247                 if (x >= 0x4000)
    00248                     xn = (sqq_table[x>>8] >> 0) + 1;
    00249                 else
    00250                     xn = (sqq_table[x>>6] >> 1) + 1;
    00251             else
    00252                 if (x >= 0x400)
    00253                     xn = (sqq_table[x>>4] >> 2) + 1;
    00254                 else
    00255                     xn = (sqq_table[x>>2] >> 3) + 1;
    00256
    00257             goto adj;
    00258         } else
    00259             return sqq_table[x] >> 4;
    00260
    00261     /* Run two iterations of the standard convergence formula */
    00262
    00263     xn = (xn + 1 + x / xn) / 2;
    00264   nr1:
    00265     xn = (xn + 1 + x / xn) / 2;
    00266   adj:
    00267
    00268     if (xn * xn > x) /* Correct rounding if necessary */
    00269         xn--;
    00270
    00271     return xn;
    00272 }
    00273
    00274 } // namespace detail
    00275
    00276 using VIGRA_CSTD::sqrt;
    00277
    00278     /*! Signed integer square root.
    00279
    00280         <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
    00281         Namespace: vigra
    00282     */
    00283 inline Int32 sqrti(Int32 v)
    00284 {
    00285     if(v < 0)
    00286         throw std::domain_error("sqrti(Int32): negative argument.");
    00287     return (Int32)detail::IntSquareRoot<UInt32>::exec((UInt32)v);
    00288 }
    00289
    00290     /*! Unsigned integer square root.
    00291
    00292         <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
    00293         Namespace: vigra
    00294     */
    00295 inline UInt32 sqrti(UInt32 v)
    00296 {
    00297     return detail::IntSquareRoot<UInt32>::exec(v);
    00298 }
    00299
    00300 #ifdef VIGRA_NO_HYPOT
    00301     /*! Compute the Euclidean distance (length of the hypothenuse of a right-angled triangle).
    00302
    00303         The  hypot()  function  returns  the  sqrt(a*a  +  b*b).
    00304         It is implemented in a way that minimizes round-off error.
    00305
    00306         <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
    00307         Namespace: vigra
    00308     */
    00309 template <class T>
    00310 T hypot(T a, T b)
    00311 {
    00312     T absa = abs(a), absb = abs(b);
    00313     if (absa > absb)
    00314         return absa * VIGRA_CSTD::sqrt(1.0 + sq(absb/absa));
    00315     else
    00316         return absb == NumericTraits<T>::zero()
    00317                    ? NumericTraits<T>::zero()
    00318                    : absb * VIGRA_CSTD::sqrt(1.0 + sq(absa/absb));
    00319 }
    00320
    00321 #else
    00322
    00323 using ::hypot;
    00324
    00325 #endif
    00326
    00327
    00328
    00329     /*! The sign function.
    00330
    00331         Returns 1, 0, or -1 depending on the sign of \a t.
    00332
    00333         <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
    00334         Namespace: vigra
    00335     */
    00336 template <class T>
    00337 T sign(T t)
    00338 {
    00339     return t > NumericTraits<T>::zero()
    00340                ? NumericTraits<T>::one()
    00341                : t < NumericTraits<T>::zero()
    00342                     ? -NumericTraits<T>::one()
    00343                     : NumericTraits<T>::zero();
    00344 }
    00345
    00346     /*! The binary sign function.
    00347
    00348         Transfers the sign of \a t2 to \a t1.
    00349
    00350         <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
    00351         Namespace: vigra
    00352     */
    00353 template <class T1, class T2>
    00354 T1 sign(T1 t1, T2 t2)
    00355 {
    00356     return t2 >= NumericTraits<T2>::zero()
    00357                ? abs(t1)
    00358                : -abs(t1);
    00359 }
    00360
    00361 #define VIGRA_DEFINE_NORM(T) \
    00362     inline NormTraits<T>::SquaredNormType squaredNorm(T t) { return sq(t); } \
    00363     inline NormTraits<T>::NormType norm(T t) { return abs(t); }
    00364
    00365 VIGRA_DEFINE_NORM(bool)
    00366 VIGRA_DEFINE_NORM(signed char)
    00367 VIGRA_DEFINE_NORM(unsigned char)
    00368 VIGRA_DEFINE_NORM(short)
    00369 VIGRA_DEFINE_NORM(unsigned short)
    00370 VIGRA_DEFINE_NORM(int)
    00371 VIGRA_DEFINE_NORM(unsigned int)
    00372 VIGRA_DEFINE_NORM(long)
    00373 VIGRA_DEFINE_NORM(unsigned long)
    00374 VIGRA_DEFINE_NORM(float)
    00375 VIGRA_DEFINE_NORM(double)
    00376 VIGRA_DEFINE_NORM(long double)
    00377
    00378 #undef VIGRA_DEFINE_NORM
    00379
    00380 template <class T>
    00381 inline typename NormTraits<std::complex<T> >::SquaredNormType
    00382 squaredNorm(std::complex<T> const & t)
    00383 {
    00384     return sq(t.real()) + sq(t.imag());
    00385 }
    00386
    00387 #ifdef DOXYGEN // only for documentation
    00388     /*! The squared norm of a numerical object.
    00389
    00390         For scalar types: equals <tt>vigra::sq(t)</tt><br>.
    00391         For vectorial types: equals <tt>vigra::dot(t, t)</tt><br>.
    00392         For complex types: equals <tt>vigra::sq(t.real()) + vigra::sq(t.imag())</tt><br>.
    00393         For matrix types: results in the squared Frobenius norm (sum of squares of the matrix elements).
    00394     */
    00395 NormTraits<T>::SquaredNormType squaredNorm(T const & t);
    00396
    00397 #endif
    00398
    00399     /*! The norm of a numerical object.
    00400
    00401         For scalar types: implemented as <tt>abs(t)</tt><br>
    00402         otherwise: implemented as <tt>sqrt(squaredNorm(t))</tt>.
    00403
    00404         <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
    00405         Namespace: vigra
    00406     */
    00407 template <class T>
    00408 inline typename NormTraits<T>::NormType
    00409 norm(T const & t)
    00410 {
    00411     typedef typename NormTraits<T>::SquaredNormType SNT;
    00412     return sqrt(static_cast<typename SquareRootTraits<SNT>::SquareRootArgument>(squaredNorm(t)));
    00413 }
    00414
    00415 namespace detail {
    00416
    00417 // both f1 and f2 are unsigned here
    00418 template<class FPT>
    00419 inline
    00420 FPT safeFloatDivision( FPT f1, FPT f2 )
    00421 {
    00422     return  f2 < NumericTraits<FPT>::one() && f1 > f2 * NumericTraits<FPT>::max()
    00423                 ? NumericTraits<FPT>::max()
    00424                 : (f2 > NumericTraits<FPT>::one() && f1 < f2 * NumericTraits<FPT>::smallestPositive()) ||
    00425                    f1 == NumericTraits<FPT>::zero()
    00426                      ? NumericTraits<FPT>::zero()
    00427                      : f1/f2;
    00428 }
    00429
    00430 } // namespace detail
    00431
    00432     /*! Tolerance based floating-point comparison.
    00433
    00434         Check whether two floating point numbers are equal within the given tolerance.
    00435         This is useful because floating point numbers that should be equal in theory are
    00436         rarely exactly equal in practice. If the tolerance \a epsilon is not given,
    00437         twice the machine epsilon is used.
    00438
    00439         <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
    00440         Namespace: vigra
    00441     */
    00442 template <class T1, class T2>
    00443 bool closeAtTolerance(T1 l, T2 r, typename PromoteTraits<T1, T2>::Promote epsilon)
    00444 {
    00445     typedef typename PromoteTraits<T1, T2>::Promote T;
    00446     if(l == 0.0 && r != 0.0)
    00447         return VIGRA_CSTD::fabs(r) <= epsilon;
    00448     if(l != 0.0 && r == 0.0)
    00449         return VIGRA_CSTD::fabs(r) <= epsilon;
    00450     T diff = VIGRA_CSTD::fabs( l - r );
    00451     T d1   = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( r ) );
    00452     T d2   = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( l ) );
    00453
    00454     return (d1 <= epsilon && d2 <= epsilon);
    00455 }
    00456
    00457 template <class T1, class T2>
    00458 bool closeAtTolerance(T1 l, T2 r)
    00459 {
    00460     typedef typename PromoteTraits<T1, T2>::Promote T;
    00461     return closeAtTolerance(l, r, 2.0 * NumericTraits<T>::epsilon());
    00462 }
    00463
    00464 //@}
    00465
    00466 } // namespace vigra
    00467
    00468 #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.4.0 (21 Dec 2005)