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

vigra/mathutil.hxx
00001 /************************************************************************/
00002 /*                                                                      */
00003 /*               Copyright 1998-2005 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_MATHUTIL_HXX
00037 #define VIGRA_MATHUTIL_HXX
00038 
00039 #ifdef _MSC_VER
00040 # pragma warning (disable: 4996) // hypot/_hypot confusion
00041 #endif
00042 
00043 #include <cmath>
00044 #include <cstdlib>
00045 #include "config.hxx"
00046 #include "error.hxx"
00047 #include "tuple.hxx"
00048 #include "sized_int.hxx"
00049 #include "numerictraits.hxx"
00050 
00051 /*! \page MathConstants Mathematical Constants
00052 
00053     <TT>M_PI, M_SQRT2</TT>
00054 
00055     <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>>
00056 
00057     Since <TT>M_PI</TT> and <TT>M_SQRT2</TT> are not officially standardized,
00058     we provide definitions here for those compilers that don't support them.
00059 
00060     \code
00061     #ifndef M_PI
00062     #    define M_PI     3.14159265358979323846
00063     #endif
00064 
00065     #ifndef M_SQRT2
00066     #    define M_SQRT2  1.41421356237309504880
00067     #endif
00068     \endcode
00069 */
00070 #ifndef M_PI
00071 #    define M_PI     3.14159265358979323846
00072 #endif
00073 
00074 #ifndef M_SQRT2
00075 #    define M_SQRT2  1.41421356237309504880
00076 #endif
00077 
00078 namespace vigra {
00079 
00080 /** \addtogroup MathFunctions Mathematical Functions
00081 
00082     Useful mathematical functions and functors.
00083 */
00084 //@{
00085 // import functions into namespace vigra which VIGRA is going to overload
00086 
00087 using VIGRA_CSTD::pow;  
00088 using VIGRA_CSTD::floor;  
00089 using VIGRA_CSTD::ceil;  
00090 
00091 // import abs(float), abs(double), abs(long double) from <cmath>
00092 //    and abs(int), abs(long), abs(long long) from <cstdlib>
00093 using std::abs;  
00094 
00095 // define the missing variants of abs() to avoid 'ambigous overload'
00096 // errors in template functions
00097 #define VIGRA_DEFINE_UNSIGNED_ABS(T) \
00098     inline T abs(T t) { return t; }
00099 
00100 VIGRA_DEFINE_UNSIGNED_ABS(bool)
00101 VIGRA_DEFINE_UNSIGNED_ABS(unsigned char)
00102 VIGRA_DEFINE_UNSIGNED_ABS(unsigned short)
00103 VIGRA_DEFINE_UNSIGNED_ABS(unsigned int)
00104 VIGRA_DEFINE_UNSIGNED_ABS(unsigned long)
00105 VIGRA_DEFINE_UNSIGNED_ABS(unsigned long long)
00106 
00107 #undef VIGRA_DEFINE_UNSIGNED_ABS
00108 
00109 #define VIGRA_DEFINE_MISSING_ABS(T) \
00110     inline T abs(T t) { return t < 0 ? -t : t; }
00111 
00112 VIGRA_DEFINE_MISSING_ABS(signed char)
00113 VIGRA_DEFINE_MISSING_ABS(signed short)
00114 
00115 #undef VIGRA_DEFINE_MISSING_ABS
00116 
00117     /*! The rounding function.
00118 
00119         Defined for all floating point types. Rounds towards the nearest integer 
00120         such that <tt>abs(round(t)) == round(abs(t))</tt> for all <tt>t</tt>.
00121 
00122         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00123         Namespace: vigra
00124     */
00125 inline float round(float t)
00126 {
00127      return t >= 0.0
00128                 ? floor(t + 0.5f)
00129                 : ceil(t - 0.5f);
00130 }
00131 
00132 inline double round(double t)
00133 {
00134      return t >= 0.0
00135                 ? floor(t + 0.5)
00136                 : ceil(t - 0.5);
00137 }
00138 
00139 inline long double round(long double t)
00140 {
00141      return t >= 0.0
00142                 ? floor(t + 0.5)
00143                 : ceil(t - 0.5);
00144 }
00145 
00146 
00147     /*! Round and cast to integer.
00148 
00149         Rounds to the nearest integer like round(), but casts the result to 
00150         <tt>int</tt> (this will be faster and is usually needed anyway).
00151 
00152         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00153         Namespace: vigra
00154     */
00155 inline int roundi(double t)
00156 {
00157      return t >= 0.0
00158                 ? int(t + 0.5)
00159                 : int(t - 0.5);
00160 }
00161 
00162     /*! Round up to the nearest power of 2.
00163 
00164         Efficient algorithm for finding the smallest power of 2 which is not smaller than \a x
00165         (function clp2() from Henry Warren: "Hacker's Delight", Addison-Wesley, 2003,
00166          see http://www.hackersdelight.org/).
00167         If \a x > 2^31, the function will return 0 because integer arithmetic is defined modulo 2^32.
00168 
00169         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00170         Namespace: vigra
00171     */
00172 inline UInt32 ceilPower2(UInt32 x) 
00173 {
00174     if(x == 0) return 0;
00175     
00176     x = x - 1;
00177     x = x | (x >> 1);
00178     x = x | (x >> 2);
00179     x = x | (x >> 4);
00180     x = x | (x >> 8);
00181     x = x | (x >>16);
00182     return x + 1;
00183 } 
00184     
00185     /*! Round down to the nearest power of 2.
00186 
00187         Efficient algorithm for finding the largest power of 2 which is not greater than \a x
00188         (function flp2() from Henry Warren: "Hacker's Delight", Addison-Wesley, 2003,
00189          see http://www.hackersdelight.org/).
00190 
00191         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00192         Namespace: vigra
00193     */
00194 inline UInt32 floorPower2(UInt32 x) 
00195 { 
00196     x = x | (x >> 1);
00197     x = x | (x >> 2);
00198     x = x | (x >> 4);
00199     x = x | (x >> 8);
00200     x = x | (x >>16);
00201     return x - (x >> 1);
00202 }
00203 
00204 namespace detail {
00205 
00206 template <class T>
00207 class IntLog2
00208 {
00209   public:
00210     static Int32 table[64];
00211 };
00212 
00213 template <class T>
00214 Int32 IntLog2<T>::table[64] = {
00215          -1,  0,  -1,  15,  -1,  1,  28,  -1,  16,  -1,  -1,  -1,  2,  21,  
00216          29,  -1,  -1,  -1,  19,  17,  10,  -1,  12,  -1,  -1,  3,  -1,  6,  
00217          -1,  22,  30,  -1,  14,  -1,  27,  -1,  -1,  -1,  20,  -1,  18,  9,  
00218          11,  -1,  5,  -1,  -1,  13,  26,  -1,  -1,  8,  -1,  4,  -1,  25,  
00219          -1,  7,  24,  -1,  23,  -1,  31,  -1};
00220 
00221 } // namespace detail
00222 
00223     /*! Compute the base-2 logarithm of an integer.
00224 
00225         Returns the position of the left-most 1-bit in the given number \a x, or
00226         -1 if \a x == 0. That is,
00227         
00228         \code
00229         assert(k >= 0 && k < 32 && log2i(1 << k) == k);
00230         \endcode
00231         
00232         The function uses Robert Harley's algorithm to determine the number of leading zeros
00233         in \a x (algorithm nlz10() at http://www.hackersdelight.org/). But note that the functions
00234         \ref floorPower2() or \ref ceilPower2() are more efficient and should be preferred when possible.
00235 
00236         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00237         Namespace: vigra
00238     */
00239 inline Int32 log2i(UInt32 x) 
00240 { 
00241     // Propagate leftmost 1-bit to the right.
00242     x = x | (x >> 1);
00243     x = x | (x >> 2);
00244     x = x | (x >> 4);
00245     x = x | (x >> 8);
00246     x = x | (x >>16);
00247     x = x*0x06EB14F9; // Multiplier is 7*255**3. 
00248     return detail::IntLog2<Int32>::table[x >> 26];
00249 }
00250 
00251     /*! The square function.
00252 
00253         <tt>sq(x) = x*x</tt> is needed so often that it makes sense to define it as a function.
00254 
00255         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00256         Namespace: vigra
00257     */
00258 template <class T>
00259 inline 
00260 typename NumericTraits<T>::Promote sq(T t)
00261 {
00262     return t*t;
00263 }
00264 
00265 namespace detail {
00266 
00267 template <class T>
00268 class IntSquareRoot
00269 {
00270   public:
00271     static UInt32 sqq_table[];
00272     static UInt32 exec(UInt32 v);
00273 };
00274 
00275 template <class T>
00276 UInt32 IntSquareRoot<T>::sqq_table[] = {
00277            0,  16,  22,  27,  32,  35,  39,  42,  45,  48,  50,  53,  55,  57,
00278           59,  61,  64,  65,  67,  69,  71,  73,  75,  76,  78,  80,  81,  83,
00279           84,  86,  87,  89,  90,  91,  93,  94,  96,  97,  98,  99, 101, 102,
00280          103, 104, 106, 107, 108, 109, 110, 112, 113, 114, 115, 116, 117, 118,
00281          119, 120, 121, 122, 123, 124, 125, 126, 128, 128, 129, 130, 131, 132,
00282          133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 144, 145,
00283          146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155, 156, 157,
00284          158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168,
00285          169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178,
00286          179, 180, 181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188,
00287          189, 189, 190, 191, 192, 192, 193, 193, 194, 195, 195, 196, 197, 197,
00288          198, 199, 199, 200, 201, 201, 202, 203, 203, 204, 204, 205, 206, 206,
00289          207, 208, 208, 209, 209, 210, 211, 211, 212, 212, 213, 214, 214, 215,
00290          215, 216, 217, 217, 218, 218, 219, 219, 220, 221, 221, 222, 222, 223,
00291          224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230, 230, 231,
00292          231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238,
00293          239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246,
00294          246, 247, 247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253,
00295          253, 254, 254, 255
00296 };
00297 
00298 template <class T>
00299 UInt32 IntSquareRoot<T>::exec(UInt32 x) 
00300 {
00301     UInt32 xn;
00302     if (x >= 0x10000)
00303         if (x >= 0x1000000)
00304             if (x >= 0x10000000)
00305                 if (x >= 0x40000000) {
00306                     if (x >= (UInt32)65535*(UInt32)65535)
00307                         return 65535;
00308                     xn = sqq_table[x>>24] << 8;
00309                 } else
00310                     xn = sqq_table[x>>22] << 7;
00311             else
00312                 if (x >= 0x4000000)
00313                     xn = sqq_table[x>>20] << 6;
00314                 else
00315                     xn = sqq_table[x>>18] << 5;
00316         else {
00317             if (x >= 0x100000)
00318                 if (x >= 0x400000)
00319                     xn = sqq_table[x>>16] << 4;
00320                 else
00321                     xn = sqq_table[x>>14] << 3;
00322             else
00323                 if (x >= 0x40000)
00324                     xn = sqq_table[x>>12] << 2;
00325                 else
00326                     xn = sqq_table[x>>10] << 1;
00327 
00328             goto nr1;
00329         }
00330     else
00331         if (x >= 0x100) {
00332             if (x >= 0x1000)
00333                 if (x >= 0x4000)
00334                     xn = (sqq_table[x>>8] >> 0) + 1;
00335                 else
00336                     xn = (sqq_table[x>>6] >> 1) + 1;
00337             else
00338                 if (x >= 0x400)
00339                     xn = (sqq_table[x>>4] >> 2) + 1;
00340                 else
00341                     xn = (sqq_table[x>>2] >> 3) + 1;
00342 
00343             goto adj;
00344         } else
00345             return sqq_table[x] >> 4;
00346 
00347     /* Run two iterations of the standard convergence formula */
00348 
00349     xn = (xn + 1 + x / xn) / 2;
00350   nr1:
00351     xn = (xn + 1 + x / xn) / 2;
00352   adj:
00353 
00354     if (xn * xn > x) /* Correct rounding if necessary */
00355         xn--;
00356 
00357     return xn;
00358 }
00359 
00360 } // namespace detail
00361 
00362 using VIGRA_CSTD::sqrt;
00363 
00364     /*! Signed integer square root.
00365     
00366         Useful for fast fixed-point computations.
00367 
00368         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00369         Namespace: vigra
00370     */
00371 inline Int32 sqrti(Int32 v)
00372 {
00373     if(v < 0)
00374         throw std::domain_error("sqrti(Int32): negative argument.");
00375     return (Int32)detail::IntSquareRoot<UInt32>::exec((UInt32)v);
00376 }
00377 
00378     /*! Unsigned integer square root.
00379 
00380         Useful for fast fixed-point computations.
00381 
00382         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00383         Namespace: vigra
00384     */
00385 inline UInt32 sqrti(UInt32 v)
00386 {
00387     return detail::IntSquareRoot<UInt32>::exec(v);
00388 }
00389 
00390 #ifdef VIGRA_NO_HYPOT
00391     /*! Compute the Euclidean distance (length of the hypothenuse of a right-angled triangle).
00392 
00393         The  hypot()  function  returns  the  sqrt(a*a  +  b*b).
00394         It is implemented in a way that minimizes round-off error.
00395 
00396         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00397         Namespace: vigra
00398     */
00399 inline double hypot(double a, double b) 
00400 { 
00401     double absa = VIGRA_CSTD::fabs(a), absb = VIGRA_CSTD::fabs(b);
00402     if (absa > absb) 
00403         return absa * VIGRA_CSTD::sqrt(1.0 + sq(absb/absa)); 
00404     else 
00405         return absb == 0.0
00406                    ? 0.0
00407                    : absb * VIGRA_CSTD::sqrt(1.0 + sq(absa/absb)); 
00408 }
00409 
00410 #else
00411 
00412 using ::hypot;
00413 
00414 #endif
00415 
00416     /*! The sign function.
00417 
00418         Returns 1, 0, or -1 depending on the sign of \a t, but with the same type as \a t.
00419 
00420         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00421         Namespace: vigra
00422     */
00423 template <class T>
00424 inline T sign(T t) 
00425 { 
00426     return t > NumericTraits<T>::zero()
00427                ? NumericTraits<T>::one()
00428                : t < NumericTraits<T>::zero()
00429                     ? -NumericTraits<T>::one()
00430                     : NumericTraits<T>::zero();
00431 }
00432 
00433     /*! The integer sign function.
00434 
00435         Returns 1, 0, or -1 depending on the sign of \a t, converted to int.
00436 
00437         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00438         Namespace: vigra
00439     */
00440 template <class T>
00441 inline int signi(T t) 
00442 { 
00443     return t > NumericTraits<T>::zero()
00444                ? 1
00445                : t < NumericTraits<T>::zero()
00446                     ? -1
00447                     : 0;
00448 }
00449 
00450     /*! The binary sign function.
00451 
00452         Transfers the sign of \a t2 to \a t1.
00453 
00454         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00455         Namespace: vigra
00456     */
00457 template <class T1, class T2>
00458 inline T1 sign(T1 t1, T2 t2) 
00459 { 
00460     return t2 >= NumericTraits<T2>::zero()
00461                ? abs(t1)
00462                : -abs(t1);
00463 }
00464 
00465 #define VIGRA_DEFINE_NORM(T) \
00466     inline NormTraits<T>::SquaredNormType squaredNorm(T t) { return sq(t); } \
00467     inline NormTraits<T>::NormType norm(T t) { return abs(t); }
00468 
00469 VIGRA_DEFINE_NORM(bool)
00470 VIGRA_DEFINE_NORM(signed char)
00471 VIGRA_DEFINE_NORM(unsigned char)
00472 VIGRA_DEFINE_NORM(short)
00473 VIGRA_DEFINE_NORM(unsigned short)
00474 VIGRA_DEFINE_NORM(int)
00475 VIGRA_DEFINE_NORM(unsigned int)
00476 VIGRA_DEFINE_NORM(long)
00477 VIGRA_DEFINE_NORM(unsigned long)
00478 VIGRA_DEFINE_NORM(float)
00479 VIGRA_DEFINE_NORM(double)
00480 VIGRA_DEFINE_NORM(long double)
00481 
00482 #undef VIGRA_DEFINE_NORM
00483 
00484 template <class T>
00485 inline typename NormTraits<std::complex<T> >::SquaredNormType
00486 squaredNorm(std::complex<T> const & t)
00487 {
00488     return sq(t.real()) + sq(t.imag());
00489 }
00490 
00491 #ifdef DOXYGEN // only for documentation
00492     /*! The squared norm of a numerical object.
00493 
00494         For scalar types: equals <tt>vigra::sq(t)</tt><br>.
00495         For vectorial types: equals <tt>vigra::dot(t, t)</tt><br>.
00496         For complex types: equals <tt>vigra::sq(t.real()) + vigra::sq(t.imag())</tt><br>.
00497         For matrix types: results in the squared Frobenius norm (sum of squares of the matrix elements).
00498     */
00499 NormTraits<T>::SquaredNormType squaredNorm(T const & t);
00500 
00501 #endif
00502 
00503     /*! The norm of a numerical object.
00504 
00505         For scalar types: implemented as <tt>abs(t)</tt><br>
00506         otherwise: implemented as <tt>sqrt(squaredNorm(t))</tt>.
00507 
00508         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00509         Namespace: vigra
00510     */
00511 template <class T>
00512 inline typename NormTraits<T>::NormType 
00513 norm(T const & t)
00514 {
00515     typedef typename NormTraits<T>::SquaredNormType SNT;
00516     return sqrt(static_cast<typename SquareRootTraits<SNT>::SquareRootArgument>(squaredNorm(t)));
00517 }
00518 
00519     /*! Find the minimum element in a sequence.
00520     
00521         The function returns the iterator refering to the minimum element.
00522         
00523         <b>Required Interface:</b>
00524         
00525         \code
00526         Iterator is a standard forward iterator.
00527         
00528         bool f = *first < NumericTraits<typename std::iterator_traits<Iterator>::value_type>::max();
00529         \endcode
00530 
00531         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00532         Namespace: vigra
00533     */
00534 template <class Iterator>
00535 Iterator argMin(Iterator first, Iterator last)
00536 {
00537     typedef typename std::iterator_traits<Iterator>::value_type Value;
00538     Value vopt = NumericTraits<Value>::max();
00539     Iterator best = last;
00540     for(; first != last; ++first)
00541     {
00542         if(*first < vopt)
00543         {
00544             vopt = *first;
00545             best = first;
00546         }
00547     }
00548     return best;
00549 }
00550 
00551     /*! Find the maximum element in a sequence.
00552     
00553         The function returns the iterator refering to the maximum element.
00554         
00555         <b>Required Interface:</b>
00556         
00557         \code
00558         Iterator is a standard forward iterator.
00559         
00560         bool f = NumericTraits<typename std::iterator_traits<Iterator>::value_type>::min() < *first;
00561         \endcode
00562 
00563         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00564         Namespace: vigra
00565     */
00566 template <class Iterator>
00567 Iterator argMax(Iterator first, Iterator last)
00568 {
00569     typedef typename std::iterator_traits<Iterator>::value_type Value;
00570     Value vopt = NumericTraits<Value>::min();
00571     Iterator best = last;
00572     for(; first != last; ++first)
00573     {
00574         if(vopt < *first)
00575         {
00576             vopt = *first;
00577             best = first;
00578         }
00579     }
00580     return best;
00581 }
00582 
00583     /*! Find the minimum element in a sequence conforming to a condition.
00584     
00585         The function returns the iterator refering to the minimum element,
00586         where only elements conforming to the condition (i.e. where 
00587         <tt>condition(*iterator)</tt> evaluates to <tt>true</tt>) are considered.
00588         If no element conforms to the condition, or the sequence is empty,
00589         the end iterator \a last is returned.
00590         
00591         <b>Required Interface:</b>
00592         
00593         \code
00594         Iterator is a standard forward iterator.
00595         
00596         bool c = condition(*first);
00597         
00598         bool f = *first < NumericTraits<typename std::iterator_traits<Iterator>::value_type>::max();
00599         \endcode
00600 
00601         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00602         Namespace: vigra
00603     */
00604 template <class Iterator, class UnaryFunctor>
00605 Iterator argMinIf(Iterator first, Iterator last, UnaryFunctor condition)
00606 {
00607     typedef typename std::iterator_traits<Iterator>::value_type Value;
00608     Value vopt = NumericTraits<Value>::max();
00609     Iterator best = last;
00610     for(; first != last; ++first)
00611     {
00612         if(condition(*first) && *first < vopt)
00613         {
00614             vopt = *first;
00615             best = first;
00616         }
00617     }
00618     return best;
00619 }
00620 
00621     /*! Find the maximum element in a sequence conforming to a condition.
00622     
00623         The function returns the iterator refering to the maximum element,
00624         where only elements conforming to the condition (i.e. where 
00625         <tt>condition(*iterator)</tt> evaluates to <tt>true</tt>) are considered.
00626         If no element conforms to the condition, or the sequence is empty,
00627         the end iterator \a last is returned.
00628         
00629         <b>Required Interface:</b>
00630         
00631         \code
00632         Iterator is a standard forward iterator.
00633         
00634         bool c = condition(*first);
00635         
00636         bool f = NumericTraits<typename std::iterator_traits<Iterator>::value_type>::min() < *first;
00637         \endcode
00638 
00639         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00640         Namespace: vigra
00641     */
00642 template <class Iterator, class UnaryFunctor>
00643 Iterator argMaxIf(Iterator first, Iterator last, UnaryFunctor condition)
00644 {
00645     typedef typename std::iterator_traits<Iterator>::value_type Value;
00646     Value vopt = NumericTraits<Value>::min();
00647     Iterator best = last;
00648     for(; first != last; ++first)
00649     {
00650         if(condition(*first) && vopt < *first)
00651         {
00652             vopt = *first;
00653             best = first;
00654         }
00655     }
00656     return best;
00657 }
00658 
00659     /*! Compute the eigenvalues of a 2x2 real symmetric matrix.
00660     
00661         This uses the analytical eigenvalue formula 
00662         \f[
00663            \lambda_{1,2} = \frac{1}{2}\left(a_{00} + a_{11} \pm \sqrt{(a_{00} - a_{11})^2 + 4 a_{01}^2}\right)
00664         \f]
00665 
00666         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00667         Namespace: vigra
00668     */
00669 template <class T>
00670 void symmetric2x2Eigenvalues(T a00, T a01, T a11, T * r0, T * r1)
00671 {
00672     double d  = hypot(a00 - a11, 2.0*a01);
00673     *r0 = static_cast<T>(0.5*(a00 + a11 + d));
00674     *r1 = static_cast<T>(0.5*(a00 + a11 - d));
00675     if(*r0 < *r1)
00676         std::swap(*r0, *r1);
00677 }
00678 
00679     /*! Compute the eigenvalues of a 3x3 real symmetric matrix.
00680     
00681         This uses a numerically stable version of the analytical eigenvalue formula according to
00682         <p>
00683         David Eberly: <a href="http://www.geometrictools.com/Documentation/EigenSymmetric3x3.pdf">
00684         <em>"Eigensystems for 3 × 3 Symmetric Matrices (Revisited)"</em></a>, Geometric Tools Documentation, 2006
00685 
00686 
00687         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00688         Namespace: vigra
00689     */
00690 template <class T>
00691 void symmetric3x3Eigenvalues(T a00, T a01, T a02, T a11, T a12, T a22,
00692                              T * r0, T * r1, T * r2)
00693 {
00694     static double inv3 = 1.0 / 3.0, root3 = std::sqrt(3.0);
00695     
00696     double c0 = a00*a11*a22 + 2.0*a01*a02*a12 - a00*a12*a12 - a11*a02*a02 - a22*a01*a01;
00697     double c1 = a00*a11 - a01*a01 + a00*a22 - a02*a02 + a11*a22 - a12*a12;
00698     double c2 = a00 + a11 + a22;
00699     double c2Div3 = c2*inv3;
00700     double aDiv3 = (c1 - c2*c2Div3)*inv3;
00701     if (aDiv3 > 0.0) 
00702         aDiv3 = 0.0;
00703     double mbDiv2 = 0.5*(c0 + c2Div3*(2.0*c2Div3*c2Div3 - c1));
00704     double q = mbDiv2*mbDiv2 + aDiv3*aDiv3*aDiv3;
00705     if (q > 0.0) 
00706         q = 0.0;
00707     double magnitude = std::sqrt(-aDiv3);
00708     double angle = std::atan2(std::sqrt(-q),mbDiv2)*inv3;
00709     double cs = std::cos(angle);
00710     double sn = std::sin(angle);
00711     *r0 = static_cast<T>(c2Div3 + 2.0*magnitude*cs);
00712     *r1 = static_cast<T>(c2Div3 - magnitude*(cs + root3*sn));
00713     *r2 = static_cast<T>(c2Div3 - magnitude*(cs - root3*sn));
00714     if(*r0 < *r1)
00715         std::swap(*r0, *r1);
00716     if(*r0 < *r2)
00717         std::swap(*r0, *r2);
00718     if(*r1 < *r2)
00719         std::swap(*r1, *r2);
00720 }
00721 
00722 namespace detail {
00723 
00724 template <class T>
00725 T ellipticRD(T x, T y, T z)
00726 {
00727     double f = 1.0, s = 0.0, X, Y, Z, m;
00728     for(;;)
00729     {
00730         m = (x + y + 3.0*z) / 5.0;
00731         X = 1.0 - x/m;
00732         Y = 1.0 - y/m;
00733         Z = 1.0 - z/m;
00734         if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01)
00735             break;
00736         double l = VIGRA_CSTD::sqrt(x*y) + VIGRA_CSTD::sqrt(x*z) + VIGRA_CSTD::sqrt(y*z);
00737         s += f / (VIGRA_CSTD::sqrt(z)*(z + l));
00738         f /= 4.0;
00739         x = (x + l)/4.0;
00740         y = (y + l)/4.0;
00741         z = (z + l)/4.0;
00742     }
00743     double a = X*Y;
00744     double b = sq(Z);
00745     double c = a - b;
00746     double d = a - 6.0*b;
00747     double e = d + 2.0*c;
00748     return 3.0*s + f*(1.0+d*(-3.0/14.0+d*9.0/88.0-Z*e*4.5/26.0)
00749                       +Z*(e/6.0+Z*(-c*9.0/22.0+a*Z*3.0/26.0))) / VIGRA_CSTD::pow(m,1.5);
00750 }
00751 
00752 template <class T>
00753 T ellipticRF(T x, T y, T z)
00754 {
00755     double X, Y, Z, m;
00756     for(;;)
00757     {
00758         m = (x + y + z) / 3.0;
00759         X = 1.0 - x/m;
00760         Y = 1.0 - y/m;
00761         Z = 1.0 - z/m;
00762         if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01)
00763             break;
00764         double l = VIGRA_CSTD::sqrt(x*y) + VIGRA_CSTD::sqrt(x*z) + VIGRA_CSTD::sqrt(y*z);
00765         x = (x + l)/4.0;
00766         y = (y + l)/4.0;
00767         z = (z + l)/4.0;
00768     }
00769     double d = X*Y - sq(Z);
00770     double p = X*Y*Z;
00771     return (1.0 - d/10.0 + p/14.0 + sq(d)/24.0 - d*p*3.0/44.0) / VIGRA_CSTD::sqrt(m);
00772 }
00773 
00774 } // namespace detail
00775 
00776     /*! The incomplete elliptic integral of the first kind.
00777 
00778         Computes
00779         
00780         \f[
00781             \mbox{F}(x, k) = \int_0^x \frac{1}{\sqrt{1 - k^2 \sin(t)^2}} dt
00782         \f]
00783         
00784         according to the algorithm given in Press et al. "Numerical Recipes". 
00785 
00786         Note: In some libraries (e.g. Mathematica), the second parameter of the elliptic integral
00787         functions must be k^2 rather than k. Check the documentation when results disagree!
00788 
00789         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00790         Namespace: vigra
00791     */
00792 inline double ellipticIntegralF(double x, double k)
00793 {
00794     double c2 = sq(VIGRA_CSTD::cos(x));
00795     double s = VIGRA_CSTD::sin(x);
00796     return s*detail::ellipticRF(c2, 1.0 - sq(k*s), 1.0);
00797 }
00798 
00799     /*! The incomplete elliptic integral of the second kind.
00800 
00801         Computes
00802         
00803         \f[
00804             \mbox{E}(x, k) = \int_0^x \sqrt{1 - k^2 \sin(t)^2} dt
00805         \f]
00806         
00807         according to the algorithm given in Press et al. "Numerical Recipes". The
00808         complete elliptic integral of the second kind is simply <tt>ellipticIntegralE(M_PI/2, k)</TT>.
00809         
00810         Note: In some libraries (e.g. Mathematica), the second parameter of the elliptic integral
00811         functions must be k^2 rather than k. Check the documentation when results disagree!
00812 
00813         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00814         Namespace: vigra
00815     */
00816 inline double ellipticIntegralE(double x, double k)
00817 {
00818     double c2 = sq(VIGRA_CSTD::cos(x));
00819     double s = VIGRA_CSTD::sin(x);
00820     k = sq(k*s);
00821     return s*(detail::ellipticRF(c2, 1.0-k, 1.0) - k/3.0*detail::ellipticRD(c2, 1.0-k, 1.0));
00822 }
00823 
00824 #if _MSC_VER
00825 
00826 namespace detail {
00827 
00828 template <class T>
00829 double erfImpl(T x)
00830 {
00831     double t = 1.0/(1.0+0.5*VIGRA_CSTD::fabs(x));
00832     double ans = t*VIGRA_CSTD::exp(-x*x-1.26551223+t*(1.00002368+t*(0.37409196+
00833                                     t*(0.09678418+t*(-0.18628806+t*(0.27886807+
00834                                     t*(-1.13520398+t*(1.48851587+t*(-0.82215223+
00835                                     t*0.17087277)))))))));
00836     if (x >= 0.0)
00837         return 1.0 - ans;
00838     else
00839         return ans - 1.0;
00840 }
00841 
00842 } // namespace detail 
00843 
00844     /*! The error function.
00845 
00846         If <tt>erf()</tt> is not provided in the C standard math library (as it should according to the
00847         new C99 standard ?), VIGRA implements <tt>erf()</tt> as an approximation of the error 
00848         function
00849         
00850         \f[
00851             \mbox{erf}(x) = \int_0^x e^{-t^2} dt
00852         \f]
00853         
00854         according to the formula given in Press et al. "Numerical Recipes".
00855 
00856         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00857         Namespace: vigra
00858     */
00859 inline double erf(double x)
00860 {
00861     return detail::erfImpl(x);
00862 }
00863 
00864 #else
00865 
00866 using ::erf;
00867 
00868 #endif
00869 
00870 namespace detail {
00871 
00872 template <class T>
00873 double noncentralChi2CDFApprox(unsigned int degreesOfFreedom, T noncentrality, T arg)
00874 {
00875     double a = degreesOfFreedom + noncentrality;
00876     double b = (a + noncentrality) / sq(a);
00877     double t = (VIGRA_CSTD::pow((double)arg / a, 1.0/3.0) - (1.0 - 2.0 / 9.0 * b)) / VIGRA_CSTD::sqrt(2.0 / 9.0 * b);
00878     return 0.5*(1.0 + erf(t/VIGRA_CSTD::sqrt(2.0)));
00879 }
00880 
00881 template <class T>
00882 void noncentralChi2OneIteration(T arg, T & lans, T & dans, T & pans, unsigned int & j)
00883 {
00884     double tol = -50.0;
00885     if(lans < tol)
00886     {
00887         lans = lans + VIGRA_CSTD::log(arg / j);
00888         dans = VIGRA_CSTD::exp(lans);
00889     }
00890     else
00891     {
00892         dans = dans * arg / j;
00893     }
00894     pans = pans - dans;
00895     j += 2;
00896 }
00897 
00898 template <class T>
00899 std::pair<double, double> noncentralChi2CDF(unsigned int degreesOfFreedom, T noncentrality, T arg, T eps)
00900 {
00901     vigra_precondition(noncentrality >= 0.0 && arg >= 0.0 && eps > 0.0,
00902         "noncentralChi2P(): parameters must be positive.");
00903     if (arg == 0.0 && degreesOfFreedom > 0)
00904         return std::make_pair(0.0, 0.0);
00905 
00906     // Determine initial values
00907     double b1 = 0.5 * noncentrality,
00908            ao = VIGRA_CSTD::exp(-b1),
00909            eps2 = eps / ao,
00910            lnrtpi2 = 0.22579135264473,
00911            probability, density, lans, dans, pans, sum, am, hold;
00912     unsigned int maxit = 500,
00913         i, m;
00914     if(degreesOfFreedom % 2)
00915     {
00916         i = 1;
00917         lans = -0.5 * (arg + VIGRA_CSTD::log(arg)) - lnrtpi2;
00918         dans = VIGRA_CSTD::exp(lans);
00919         pans = erf(VIGRA_CSTD::sqrt(arg/2.0));
00920     }
00921     else
00922     {
00923         i = 2;
00924         lans = -0.5 * arg;
00925         dans = VIGRA_CSTD::exp(lans);
00926         pans = 1.0 - dans;
00927     }
00928     
00929     // Evaluate first term
00930     if(degreesOfFreedom == 0)
00931     {
00932         m = 1;
00933         degreesOfFreedom = 2;
00934         am = b1;
00935         sum = 1.0 / ao - 1.0 - am;
00936         density = am * dans;
00937         probability = 1.0 + am * pans;
00938     }
00939     else
00940     {
00941         m = 0;
00942         degreesOfFreedom = degreesOfFreedom - 1;
00943         am = 1.0;
00944         sum = 1.0 / ao - 1.0;
00945         while(i < degreesOfFreedom)
00946             detail::noncentralChi2OneIteration(arg, lans, dans, pans, i);
00947         degreesOfFreedom = degreesOfFreedom + 1;
00948         density = dans;
00949         probability = pans;
00950     }
00951     // Evaluate successive terms of the expansion
00952     for(++m; m<maxit; ++m)
00953     {
00954         am = b1 * am / m;
00955         detail::noncentralChi2OneIteration(arg, lans, dans, pans, degreesOfFreedom);
00956         sum = sum - am;
00957         density = density + am * dans;
00958         hold = am * pans;
00959         probability = probability + hold;
00960         if((pans * sum < eps2) && (hold < eps2))
00961             break; // converged
00962     }
00963     if(m == maxit)
00964         vigra_fail("noncentralChi2P(): no convergence.");
00965     return std::make_pair(0.5 * ao * density, std::min(1.0, std::max(0.0, ao * probability)));
00966 }
00967 
00968 } // namespace detail
00969 
00970     /*! Chi square distribution. 
00971 
00972         Computes the density of a chi square distribution with \a degreesOfFreedom 
00973         and tolerance \a accuracy at the given argument \a arg
00974         by calling <tt>noncentralChi2(degreesOfFreedom, 0.0, arg, accuracy)</tt>.
00975 
00976         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00977         Namespace: vigra
00978     */
00979 inline double chi2(unsigned int degreesOfFreedom, double arg, double accuracy = 1e-7)
00980 {
00981     return detail::noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy).first;
00982 }
00983 
00984     /*! Cumulative chi square distribution. 
00985 
00986         Computes the cumulative density of a chi square distribution with \a degreesOfFreedom 
00987         and tolerance \a accuracy at the given argument \a arg, i.e. the probability that
00988         a random number drawn from the distribution is below \a arg
00989         by calling <tt>noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy)</tt>.
00990 
00991         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
00992         Namespace: vigra
00993     */
00994 inline double chi2CDF(unsigned int degreesOfFreedom, double arg, double accuracy = 1e-7)
00995 {
00996     return detail::noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy).second;
00997 }
00998 
00999     /*! Non-central chi square distribution. 
01000 
01001         Computes the density of a non-central chi square distribution with \a degreesOfFreedom, 
01002         noncentrality parameter \a noncentrality and tolerance \a accuracy at the given argument 
01003         \a arg. It uses Algorithm AS 231 from Appl. Statist. (1987) Vol.36, No.3 (code ported from 
01004         http://lib.stat.cmu.edu/apstat/231). The algorithm has linear complexity in the number of
01005         degrees of freedom.
01006 
01007         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
01008         Namespace: vigra
01009     */
01010 inline double noncentralChi2(unsigned int degreesOfFreedom, 
01011               double noncentrality, double arg, double accuracy = 1e-7)
01012 {
01013     return detail::noncentralChi2CDF(degreesOfFreedom, noncentrality, arg, accuracy).first;
01014 }
01015 
01016     /*! Cumulative non-central chi square distribution. 
01017 
01018         Computes the cumulative density of a chi square distribution with \a degreesOfFreedom, 
01019         noncentrality parameter \a noncentrality and tolerance \a accuracy at the given argument 
01020         \a arg, i.e. the probability that a random number drawn from the distribution is below \a arg
01021         It uses Algorithm AS 231 from Appl. Statist. (1987) Vol.36, No.3 (code ported from 
01022         http://lib.stat.cmu.edu/apstat/231). The algorithm has linear complexity in the number of
01023         degrees of freedom (see noncentralChi2CDFApprox() for a constant-time algorithm).
01024 
01025         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
01026         Namespace: vigra
01027     */
01028 inline double noncentralChi2CDF(unsigned int degreesOfFreedom, 
01029               double noncentrality, double arg, double accuracy = 1e-7)
01030 {
01031     return detail::noncentralChi2CDF(degreesOfFreedom, noncentrality, arg, accuracy).second;
01032 }
01033 
01034     /*! Cumulative non-central chi square distribution (approximate). 
01035 
01036         Computes approximate values of the cumulative density of a chi square distribution with \a degreesOfFreedom, 
01037         and noncentrality parameter \a noncentrality at the given argument 
01038         \a arg, i.e. the probability that a random number drawn from the distribution is below \a arg
01039         It uses the approximate transform into a normal distribution due to Wilson and Hilferty 
01040         (see Abramovitz, Stegun: "Handbook of Mathematical Functions", formula 26.3.32). 
01041         The algorithm's running time is independent of the inputs, i.e. is should be used
01042         when noncentralChi2CDF() is too slow, and approximate values are sufficient. The accuracy is only 
01043         about 0.1 for few degrees of freedom, but reaches about 0.001 above dof = 5.
01044 
01045         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
01046         Namespace: vigra
01047     */
01048 inline double noncentralChi2CDFApprox(unsigned int degreesOfFreedom, double noncentrality, double arg)
01049 {
01050     return detail::noncentralChi2CDFApprox(degreesOfFreedom, noncentrality, arg);
01051 }
01052 
01053 
01054 
01055 namespace detail {
01056 
01057 // both f1 and f2 are unsigned here
01058 template<class FPT>
01059 inline
01060 FPT safeFloatDivision( FPT f1, FPT f2 )
01061 {
01062     return  f2 < NumericTraits<FPT>::one() && f1 > f2 * NumericTraits<FPT>::max()
01063                 ? NumericTraits<FPT>::max() 
01064                 : (f2 > NumericTraits<FPT>::one() && f1 < f2 * NumericTraits<FPT>::smallestPositive()) || 
01065                    f1 == NumericTraits<FPT>::zero()
01066                      ? NumericTraits<FPT>::zero() 
01067                      : f1/f2;
01068 }
01069 
01070 } // namespace detail
01071     
01072     /*! Tolerance based floating-point comparison.
01073 
01074         Check whether two floating point numbers are equal within the given tolerance.
01075         This is useful because floating point numbers that should be equal in theory are
01076         rarely exactly equal in practice. If the tolerance \a epsilon is not given,
01077         twice the machine epsilon is used.
01078 
01079         <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
01080         Namespace: vigra
01081     */
01082 template <class T1, class T2>
01083 bool 
01084 closeAtTolerance(T1 l, T2 r, typename PromoteTraits<T1, T2>::Promote epsilon)
01085 {
01086     typedef typename PromoteTraits<T1, T2>::Promote T;
01087     if(l == 0.0)
01088         return VIGRA_CSTD::fabs(r) <= epsilon;
01089     if(r == 0.0)
01090         return VIGRA_CSTD::fabs(l) <= epsilon;
01091     T diff = VIGRA_CSTD::fabs( l - r );
01092     T d1   = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( r ) );
01093     T d2   = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( l ) );
01094 
01095     return (d1 <= epsilon && d2 <= epsilon);
01096 }
01097 
01098 template <class T1, class T2>
01099 inline bool closeAtTolerance(T1 l, T2 r)
01100 {
01101     typedef typename PromoteTraits<T1, T2>::Promote T;
01102     return closeAtTolerance(l, r, 2.0 * NumericTraits<T>::epsilon());
01103 }
01104 
01105 //@}
01106 
01107 } // namespace vigra
01108 
01109 #endif /* VIGRA_MATHUTIL_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.7.0 (Thu Aug 25 2011)