[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
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) |
html generated using doxygen and Python
|