[ 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.1 (3 Dec 2010)