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