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