[ 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.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)