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

mathutil.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2005 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 #ifndef VIGRA_MATHUTIL_HXX
37 #define VIGRA_MATHUTIL_HXX
38 
39 #ifdef _MSC_VER
40 # pragma warning (disable: 4996) // hypot/_hypot confusion
41 #endif
42 
43 #include <cmath>
44 #include <cstdlib>
45 #include "config.hxx"
46 #include "error.hxx"
47 #include "tuple.hxx"
48 #include "sized_int.hxx"
49 #include "numerictraits.hxx"
50 
51 /*! \page MathConstants Mathematical Constants
52 
53  <TT>M_PI, M_SQRT2</TT>
54 
55  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>>
56 
57  Since <TT>M_PI</TT> and <TT>M_SQRT2</TT> are not officially standardized,
58  we provide definitions here for those compilers that don't support them.
59 
60  \code
61  #ifndef M_PI
62  # define M_PI 3.14159265358979323846
63  #endif
64 
65  #ifndef M_SQRT2
66  # define M_SQRT2 1.41421356237309504880
67  #endif
68  \endcode
69 */
70 #ifndef M_PI
71 # define M_PI 3.14159265358979323846
72 #endif
73 
74 #ifndef M_SQRT2
75 # define M_SQRT2 1.41421356237309504880
76 #endif
77 
78 namespace vigra {
79 
80 /** \addtogroup MathFunctions Mathematical Functions
81 
82  Useful mathematical functions and functors.
83 */
84 //@{
85 // import functions into namespace vigra which VIGRA is going to overload
86 
87 using VIGRA_CSTD::pow;
88 using VIGRA_CSTD::floor;
89 using VIGRA_CSTD::ceil;
90 
91 // import abs(float), abs(double), abs(long double) from <cmath>
92 // and abs(int), abs(long), abs(long long) from <cstdlib>
93 using std::abs;
94 
95 // define the missing variants of abs() to avoid 'ambigous overload'
96 // errors in template functions
97 #define VIGRA_DEFINE_UNSIGNED_ABS(T) \
98  inline T abs(T t) { return t; }
99 
100 VIGRA_DEFINE_UNSIGNED_ABS(bool)
101 VIGRA_DEFINE_UNSIGNED_ABS(unsigned char)
102 VIGRA_DEFINE_UNSIGNED_ABS(unsigned short)
103 VIGRA_DEFINE_UNSIGNED_ABS(unsigned int)
104 VIGRA_DEFINE_UNSIGNED_ABS(unsigned long)
105 VIGRA_DEFINE_UNSIGNED_ABS(unsigned long long)
106 
107 #undef VIGRA_DEFINE_UNSIGNED_ABS
108 
109 #define VIGRA_DEFINE_MISSING_ABS(T) \
110  inline T abs(T t) { return t < 0 ? -t : t; }
111 
112 VIGRA_DEFINE_MISSING_ABS(signed char)
113 VIGRA_DEFINE_MISSING_ABS(signed short)
114 
115 #undef VIGRA_DEFINE_MISSING_ABS
116 
117  /*! The rounding function.
118 
119  Defined for all floating point types. Rounds towards the nearest integer
120  such that <tt>abs(round(t)) == round(abs(t))</tt> for all <tt>t</tt>.
121 
122  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
123  Namespace: vigra
124  */
125 inline float round(float t)
126 {
127  return t >= 0.0
128  ? floor(t + 0.5f)
129  : ceil(t - 0.5f);
130 }
131 
132 inline double round(double t)
133 {
134  return t >= 0.0
135  ? floor(t + 0.5)
136  : ceil(t - 0.5);
137 }
138 
139 inline long double round(long double t)
140 {
141  return t >= 0.0
142  ? floor(t + 0.5)
143  : ceil(t - 0.5);
144 }
145 
146 
147  /*! Round and cast to integer.
148 
149  Rounds to the nearest integer like round(), but casts the result to
150  <tt>int</tt> (this will be faster and is usually needed anyway).
151 
152  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
153  Namespace: vigra
154  */
155 inline int roundi(double t)
156 {
157  return t >= 0.0
158  ? int(t + 0.5)
159  : int(t - 0.5);
160 }
161 
162  /*! Round up to the nearest power of 2.
163 
164  Efficient algorithm for finding the smallest power of 2 which is not smaller than \a x
165  (function clp2() from Henry Warren: "Hacker's Delight", Addison-Wesley, 2003,
166  see http://www.hackersdelight.org/).
167  If \a x > 2^31, the function will return 0 because integer arithmetic is defined modulo 2^32.
168 
169  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
170  Namespace: vigra
171  */
173 {
174  if(x == 0) return 0;
175 
176  x = x - 1;
177  x = x | (x >> 1);
178  x = x | (x >> 2);
179  x = x | (x >> 4);
180  x = x | (x >> 8);
181  x = x | (x >>16);
182  return x + 1;
183 }
184 
185  /*! Round down to the nearest power of 2.
186 
187  Efficient algorithm for finding the largest power of 2 which is not greater than \a x
188  (function flp2() from Henry Warren: "Hacker's Delight", Addison-Wesley, 2003,
189  see http://www.hackersdelight.org/).
190 
191  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
192  Namespace: vigra
193  */
195 {
196  x = x | (x >> 1);
197  x = x | (x >> 2);
198  x = x | (x >> 4);
199  x = x | (x >> 8);
200  x = x | (x >>16);
201  return x - (x >> 1);
202 }
203 
204 namespace detail {
205 
206 template <class T>
207 class IntLog2
208 {
209  public:
210  static Int32 table[64];
211 };
212 
213 template <class T>
214 Int32 IntLog2<T>::table[64] = {
215  -1, 0, -1, 15, -1, 1, 28, -1, 16, -1, -1, -1, 2, 21,
216  29, -1, -1, -1, 19, 17, 10, -1, 12, -1, -1, 3, -1, 6,
217  -1, 22, 30, -1, 14, -1, 27, -1, -1, -1, 20, -1, 18, 9,
218  11, -1, 5, -1, -1, 13, 26, -1, -1, 8, -1, 4, -1, 25,
219  -1, 7, 24, -1, 23, -1, 31, -1};
220 
221 } // namespace detail
222 
223  /*! Compute the base-2 logarithm of an integer.
224 
225  Returns the position of the left-most 1-bit in the given number \a x, or
226  -1 if \a x == 0. That is,
227 
228  \code
229  assert(k >= 0 && k < 32 && log2i(1 << k) == k);
230  \endcode
231 
232  The function uses Robert Harley's algorithm to determine the number of leading zeros
233  in \a x (algorithm nlz10() at http://www.hackersdelight.org/). But note that the functions
234  \ref floorPower2() or \ref ceilPower2() are more efficient and should be preferred when possible.
235 
236  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
237  Namespace: vigra
238  */
239 inline Int32 log2i(UInt32 x)
240 {
241  // Propagate leftmost 1-bit to the right.
242  x = x | (x >> 1);
243  x = x | (x >> 2);
244  x = x | (x >> 4);
245  x = x | (x >> 8);
246  x = x | (x >>16);
247  x = x*0x06EB14F9; // Multiplier is 7*255**3.
248  return detail::IntLog2<Int32>::table[x >> 26];
249 }
250 
251  /*! The square function.
252 
253  <tt>sq(x) = x*x</tt> is needed so often that it makes sense to define it as a function.
254 
255  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
256  Namespace: vigra
257  */
258 template <class T>
259 inline
260 typename NumericTraits<T>::Promote sq(T t)
261 {
262  return t*t;
263 }
264 
265 namespace detail {
266 
267 template <class T>
268 class IntSquareRoot
269 {
270  public:
271  static UInt32 sqq_table[];
272  static UInt32 exec(UInt32 v);
273 };
274 
275 template <class T>
276 UInt32 IntSquareRoot<T>::sqq_table[] = {
277  0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53, 55, 57,
278  59, 61, 64, 65, 67, 69, 71, 73, 75, 76, 78, 80, 81, 83,
279  84, 86, 87, 89, 90, 91, 93, 94, 96, 97, 98, 99, 101, 102,
280  103, 104, 106, 107, 108, 109, 110, 112, 113, 114, 115, 116, 117, 118,
281  119, 120, 121, 122, 123, 124, 125, 126, 128, 128, 129, 130, 131, 132,
282  133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 144, 145,
283  146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155, 156, 157,
284  158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168,
285  169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178,
286  179, 180, 181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188,
287  189, 189, 190, 191, 192, 192, 193, 193, 194, 195, 195, 196, 197, 197,
288  198, 199, 199, 200, 201, 201, 202, 203, 203, 204, 204, 205, 206, 206,
289  207, 208, 208, 209, 209, 210, 211, 211, 212, 212, 213, 214, 214, 215,
290  215, 216, 217, 217, 218, 218, 219, 219, 220, 221, 221, 222, 222, 223,
291  224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230, 230, 231,
292  231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238,
293  239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246,
294  246, 247, 247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253,
295  253, 254, 254, 255
296 };
297 
298 template <class T>
299 UInt32 IntSquareRoot<T>::exec(UInt32 x)
300 {
301  UInt32 xn;
302  if (x >= 0x10000)
303  if (x >= 0x1000000)
304  if (x >= 0x10000000)
305  if (x >= 0x40000000) {
306  if (x >= (UInt32)65535*(UInt32)65535)
307  return 65535;
308  xn = sqq_table[x>>24] << 8;
309  } else
310  xn = sqq_table[x>>22] << 7;
311  else
312  if (x >= 0x4000000)
313  xn = sqq_table[x>>20] << 6;
314  else
315  xn = sqq_table[x>>18] << 5;
316  else {
317  if (x >= 0x100000)
318  if (x >= 0x400000)
319  xn = sqq_table[x>>16] << 4;
320  else
321  xn = sqq_table[x>>14] << 3;
322  else
323  if (x >= 0x40000)
324  xn = sqq_table[x>>12] << 2;
325  else
326  xn = sqq_table[x>>10] << 1;
327 
328  goto nr1;
329  }
330  else
331  if (x >= 0x100) {
332  if (x >= 0x1000)
333  if (x >= 0x4000)
334  xn = (sqq_table[x>>8] >> 0) + 1;
335  else
336  xn = (sqq_table[x>>6] >> 1) + 1;
337  else
338  if (x >= 0x400)
339  xn = (sqq_table[x>>4] >> 2) + 1;
340  else
341  xn = (sqq_table[x>>2] >> 3) + 1;
342 
343  goto adj;
344  } else
345  return sqq_table[x] >> 4;
346 
347  /* Run two iterations of the standard convergence formula */
348 
349  xn = (xn + 1 + x / xn) / 2;
350  nr1:
351  xn = (xn + 1 + x / xn) / 2;
352  adj:
353 
354  if (xn * xn > x) /* Correct rounding if necessary */
355  xn--;
356 
357  return xn;
358 }
359 
360 } // namespace detail
361 
362 using VIGRA_CSTD::sqrt;
363 
364  /*! Signed integer square root.
365 
366  Useful for fast fixed-point computations.
367 
368  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
369  Namespace: vigra
370  */
371 inline Int32 sqrti(Int32 v)
372 {
373  if(v < 0)
374  throw std::domain_error("sqrti(Int32): negative argument.");
375  return (Int32)detail::IntSquareRoot<UInt32>::exec((UInt32)v);
376 }
377 
378  /*! Unsigned integer square root.
379 
380  Useful for fast fixed-point computations.
381 
382  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
383  Namespace: vigra
384  */
385 inline UInt32 sqrti(UInt32 v)
386 {
387  return detail::IntSquareRoot<UInt32>::exec(v);
388 }
389 
390 #ifdef VIGRA_NO_HYPOT
391  /*! Compute the Euclidean distance (length of the hypothenuse of a right-angled triangle).
392 
393  The hypot() function returns the sqrt(a*a + b*b).
394  It is implemented in a way that minimizes round-off error.
395 
396  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
397  Namespace: vigra
398  */
399 inline double hypot(double a, double b)
400 {
401  double absa = VIGRA_CSTD::fabs(a), absb = VIGRA_CSTD::fabs(b);
402  if (absa > absb)
403  return absa * VIGRA_CSTD::sqrt(1.0 + sq(absb/absa));
404  else
405  return absb == 0.0
406  ? 0.0
407  : absb * VIGRA_CSTD::sqrt(1.0 + sq(absa/absb));
408 }
409 
410 #else
411 
413 
414 #endif
415 
416  /*! The sign function.
417 
418  Returns 1, 0, or -1 depending on the sign of \a t, but with the same type as \a t.
419 
420  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
421  Namespace: vigra
422  */
423 template <class T>
424 inline T sign(T t)
425 {
426  return t > NumericTraits<T>::zero()
427  ? NumericTraits<T>::one()
428  : t < NumericTraits<T>::zero()
429  ? -NumericTraits<T>::one()
430  : NumericTraits<T>::zero();
431 }
432 
433  /*! The integer sign function.
434 
435  Returns 1, 0, or -1 depending on the sign of \a t, converted to int.
436 
437  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
438  Namespace: vigra
439  */
440 template <class T>
441 inline int signi(T t)
442 {
443  return t > NumericTraits<T>::zero()
444  ? 1
445  : t < NumericTraits<T>::zero()
446  ? -1
447  : 0;
448 }
449 
450  /*! The binary sign function.
451 
452  Transfers the sign of \a t2 to \a t1.
453 
454  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
455  Namespace: vigra
456  */
457 template <class T1, class T2>
458 inline T1 sign(T1 t1, T2 t2)
459 {
460  return t2 >= NumericTraits<T2>::zero()
461  ? abs(t1)
462  : -abs(t1);
463 }
464 
465 #define VIGRA_DEFINE_NORM(T) \
466  inline NormTraits<T>::SquaredNormType squaredNorm(T t) { return sq(t); } \
467  inline NormTraits<T>::NormType norm(T t) { return abs(t); }
468 
469 VIGRA_DEFINE_NORM(bool)
470 VIGRA_DEFINE_NORM(signed char)
471 VIGRA_DEFINE_NORM(unsigned char)
472 VIGRA_DEFINE_NORM(short)
473 VIGRA_DEFINE_NORM(unsigned short)
474 VIGRA_DEFINE_NORM(int)
475 VIGRA_DEFINE_NORM(unsigned int)
476 VIGRA_DEFINE_NORM(long)
477 VIGRA_DEFINE_NORM(unsigned long)
478 VIGRA_DEFINE_NORM(float)
479 VIGRA_DEFINE_NORM(double)
480 VIGRA_DEFINE_NORM(long double)
481 
482 #undef VIGRA_DEFINE_NORM
483 
484 template <class T>
485 inline typename NormTraits<std::complex<T> >::SquaredNormType
486 squaredNorm(std::complex<T> const & t)
487 {
488  return sq(t.real()) + sq(t.imag());
489 }
490 
491 #ifdef DOXYGEN // only for documentation
492  /*! The squared norm of a numerical object.
493 
494  For scalar types: equals <tt>vigra::sq(t)</tt><br>.
495  For vectorial types: equals <tt>vigra::dot(t, t)</tt><br>.
496  For complex types: equals <tt>vigra::sq(t.real()) + vigra::sq(t.imag())</tt><br>.
497  For matrix types: results in the squared Frobenius norm (sum of squares of the matrix elements).
498  */
499 NormTraits<T>::SquaredNormType squaredNorm(T const & t);
500 
501 #endif
502 
503  /*! The norm of a numerical object.
504 
505  For scalar types: implemented as <tt>abs(t)</tt><br>
506  otherwise: implemented as <tt>sqrt(squaredNorm(t))</tt>.
507 
508  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
509  Namespace: vigra
510  */
511 template <class T>
512 inline typename NormTraits<T>::NormType
513 norm(T const & t)
514 {
515  typedef typename NormTraits<T>::SquaredNormType SNT;
516  return sqrt(static_cast<typename SquareRootTraits<SNT>::SquareRootArgument>(squaredNorm(t)));
517 }
518 
519  /*! Find the minimum element in a sequence.
520 
521  The function returns the iterator refering to the minimum element.
522 
523  <b>Required Interface:</b>
524 
525  \code
526  Iterator is a standard forward iterator.
527 
528  bool f = *first < NumericTraits<typename std::iterator_traits<Iterator>::value_type>::max();
529  \endcode
530 
531  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
532  Namespace: vigra
533  */
534 template <class Iterator>
535 Iterator argMin(Iterator first, Iterator last)
536 {
537  typedef typename std::iterator_traits<Iterator>::value_type Value;
538  Value vopt = NumericTraits<Value>::max();
539  Iterator best = last;
540  for(; first != last; ++first)
541  {
542  if(*first < vopt)
543  {
544  vopt = *first;
545  best = first;
546  }
547  }
548  return best;
549 }
550 
551  /*! Find the maximum element in a sequence.
552 
553  The function returns the iterator refering to the maximum element.
554 
555  <b>Required Interface:</b>
556 
557  \code
558  Iterator is a standard forward iterator.
559 
560  bool f = NumericTraits<typename std::iterator_traits<Iterator>::value_type>::min() < *first;
561  \endcode
562 
563  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
564  Namespace: vigra
565  */
566 template <class Iterator>
567 Iterator argMax(Iterator first, Iterator last)
568 {
569  typedef typename std::iterator_traits<Iterator>::value_type Value;
570  Value vopt = NumericTraits<Value>::min();
571  Iterator best = last;
572  for(; first != last; ++first)
573  {
574  if(vopt < *first)
575  {
576  vopt = *first;
577  best = first;
578  }
579  }
580  return best;
581 }
582 
583  /*! Find the minimum element in a sequence conforming to a condition.
584 
585  The function returns the iterator refering to the minimum element,
586  where only elements conforming to the condition (i.e. where
587  <tt>condition(*iterator)</tt> evaluates to <tt>true</tt>) are considered.
588  If no element conforms to the condition, or the sequence is empty,
589  the end iterator \a last is returned.
590 
591  <b>Required Interface:</b>
592 
593  \code
594  Iterator is a standard forward iterator.
595 
596  bool c = condition(*first);
597 
598  bool f = *first < NumericTraits<typename std::iterator_traits<Iterator>::value_type>::max();
599  \endcode
600 
601  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
602  Namespace: vigra
603  */
604 template <class Iterator, class UnaryFunctor>
605 Iterator argMinIf(Iterator first, Iterator last, UnaryFunctor condition)
606 {
607  typedef typename std::iterator_traits<Iterator>::value_type Value;
608  Value vopt = NumericTraits<Value>::max();
609  Iterator best = last;
610  for(; first != last; ++first)
611  {
612  if(condition(*first) && *first < vopt)
613  {
614  vopt = *first;
615  best = first;
616  }
617  }
618  return best;
619 }
620 
621  /*! Find the maximum element in a sequence conforming to a condition.
622 
623  The function returns the iterator refering to the maximum element,
624  where only elements conforming to the condition (i.e. where
625  <tt>condition(*iterator)</tt> evaluates to <tt>true</tt>) are considered.
626  If no element conforms to the condition, or the sequence is empty,
627  the end iterator \a last is returned.
628 
629  <b>Required Interface:</b>
630 
631  \code
632  Iterator is a standard forward iterator.
633 
634  bool c = condition(*first);
635 
636  bool f = NumericTraits<typename std::iterator_traits<Iterator>::value_type>::min() < *first;
637  \endcode
638 
639  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
640  Namespace: vigra
641  */
642 template <class Iterator, class UnaryFunctor>
643 Iterator argMaxIf(Iterator first, Iterator last, UnaryFunctor condition)
644 {
645  typedef typename std::iterator_traits<Iterator>::value_type Value;
646  Value vopt = NumericTraits<Value>::min();
647  Iterator best = last;
648  for(; first != last; ++first)
649  {
650  if(condition(*first) && vopt < *first)
651  {
652  vopt = *first;
653  best = first;
654  }
655  }
656  return best;
657 }
658 
659  /*! Compute the eigenvalues of a 2x2 real symmetric matrix.
660 
661  This uses the analytical eigenvalue formula
662  \f[
663  \lambda_{1,2} = \frac{1}{2}\left(a_{00} + a_{11} \pm \sqrt{(a_{00} - a_{11})^2 + 4 a_{01}^2}\right)
664  \f]
665 
666  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
667  Namespace: vigra
668  */
669 template <class T>
670 void symmetric2x2Eigenvalues(T a00, T a01, T a11, T * r0, T * r1)
671 {
672  double d = hypot(a00 - a11, 2.0*a01);
673  *r0 = static_cast<T>(0.5*(a00 + a11 + d));
674  *r1 = static_cast<T>(0.5*(a00 + a11 - d));
675  if(*r0 < *r1)
676  std::swap(*r0, *r1);
677 }
678 
679  /*! Compute the eigenvalues of a 3x3 real symmetric matrix.
680 
681  This uses a numerically stable version of the analytical eigenvalue formula according to
682  <p>
683  David Eberly: <a href="http://www.geometrictools.com/Documentation/EigenSymmetric3x3.pdf">
684  <em>"Eigensystems for 3 × 3 Symmetric Matrices (Revisited)"</em></a>, Geometric Tools Documentation, 2006
685 
686 
687  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
688  Namespace: vigra
689  */
690 template <class T>
691 void symmetric3x3Eigenvalues(T a00, T a01, T a02, T a11, T a12, T a22,
692  T * r0, T * r1, T * r2)
693 {
694  static double inv3 = 1.0 / 3.0, root3 = std::sqrt(3.0);
695 
696  double c0 = a00*a11*a22 + 2.0*a01*a02*a12 - a00*a12*a12 - a11*a02*a02 - a22*a01*a01;
697  double c1 = a00*a11 - a01*a01 + a00*a22 - a02*a02 + a11*a22 - a12*a12;
698  double c2 = a00 + a11 + a22;
699  double c2Div3 = c2*inv3;
700  double aDiv3 = (c1 - c2*c2Div3)*inv3;
701  if (aDiv3 > 0.0)
702  aDiv3 = 0.0;
703  double mbDiv2 = 0.5*(c0 + c2Div3*(2.0*c2Div3*c2Div3 - c1));
704  double q = mbDiv2*mbDiv2 + aDiv3*aDiv3*aDiv3;
705  if (q > 0.0)
706  q = 0.0;
707  double magnitude = std::sqrt(-aDiv3);
708  double angle = std::atan2(std::sqrt(-q),mbDiv2)*inv3;
709  double cs = std::cos(angle);
710  double sn = std::sin(angle);
711  *r0 = static_cast<T>(c2Div3 + 2.0*magnitude*cs);
712  *r1 = static_cast<T>(c2Div3 - magnitude*(cs + root3*sn));
713  *r2 = static_cast<T>(c2Div3 - magnitude*(cs - root3*sn));
714  if(*r0 < *r1)
715  std::swap(*r0, *r1);
716  if(*r0 < *r2)
717  std::swap(*r0, *r2);
718  if(*r1 < *r2)
719  std::swap(*r1, *r2);
720 }
721 
722 namespace detail {
723 
724 template <class T>
725 T ellipticRD(T x, T y, T z)
726 {
727  double f = 1.0, s = 0.0, X, Y, Z, m;
728  for(;;)
729  {
730  m = (x + y + 3.0*z) / 5.0;
731  X = 1.0 - x/m;
732  Y = 1.0 - y/m;
733  Z = 1.0 - z/m;
734  if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01)
735  break;
736  double l = VIGRA_CSTD::sqrt(x*y) + VIGRA_CSTD::sqrt(x*z) + VIGRA_CSTD::sqrt(y*z);
737  s += f / (VIGRA_CSTD::sqrt(z)*(z + l));
738  f /= 4.0;
739  x = (x + l)/4.0;
740  y = (y + l)/4.0;
741  z = (z + l)/4.0;
742  }
743  double a = X*Y;
744  double b = sq(Z);
745  double c = a - b;
746  double d = a - 6.0*b;
747  double e = d + 2.0*c;
748  return 3.0*s + f*(1.0+d*(-3.0/14.0+d*9.0/88.0-Z*e*4.5/26.0)
749  +Z*(e/6.0+Z*(-c*9.0/22.0+a*Z*3.0/26.0))) / VIGRA_CSTD::pow(m,1.5);
750 }
751 
752 template <class T>
753 T ellipticRF(T x, T y, T z)
754 {
755  double X, Y, Z, m;
756  for(;;)
757  {
758  m = (x + y + z) / 3.0;
759  X = 1.0 - x/m;
760  Y = 1.0 - y/m;
761  Z = 1.0 - z/m;
762  if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01)
763  break;
764  double l = VIGRA_CSTD::sqrt(x*y) + VIGRA_CSTD::sqrt(x*z) + VIGRA_CSTD::sqrt(y*z);
765  x = (x + l)/4.0;
766  y = (y + l)/4.0;
767  z = (z + l)/4.0;
768  }
769  double d = X*Y - sq(Z);
770  double p = X*Y*Z;
771  return (1.0 - d/10.0 + p/14.0 + sq(d)/24.0 - d*p*3.0/44.0) / VIGRA_CSTD::sqrt(m);
772 }
773 
774 } // namespace detail
775 
776  /*! The incomplete elliptic integral of the first kind.
777 
778  Computes
779 
780  \f[
781  \mbox{F}(x, k) = \int_0^x \frac{1}{\sqrt{1 - k^2 \sin(t)^2}} dt
782  \f]
783 
784  according to the algorithm given in Press et al. "Numerical Recipes".
785 
786  Note: In some libraries (e.g. Mathematica), the second parameter of the elliptic integral
787  functions must be k^2 rather than k. Check the documentation when results disagree!
788 
789  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
790  Namespace: vigra
791  */
792 inline double ellipticIntegralF(double x, double k)
793 {
794  double c2 = sq(VIGRA_CSTD::cos(x));
795  double s = VIGRA_CSTD::sin(x);
796  return s*detail::ellipticRF(c2, 1.0 - sq(k*s), 1.0);
797 }
798 
799  /*! The incomplete elliptic integral of the second kind.
800 
801  Computes
802 
803  \f[
804  \mbox{E}(x, k) = \int_0^x \sqrt{1 - k^2 \sin(t)^2} dt
805  \f]
806 
807  according to the algorithm given in Press et al. "Numerical Recipes". The
808  complete elliptic integral of the second kind is simply <tt>ellipticIntegralE(M_PI/2, k)</TT>.
809 
810  Note: In some libraries (e.g. Mathematica), the second parameter of the elliptic integral
811  functions must be k^2 rather than k. Check the documentation when results disagree!
812 
813  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
814  Namespace: vigra
815  */
816 inline double ellipticIntegralE(double x, double k)
817 {
818  double c2 = sq(VIGRA_CSTD::cos(x));
819  double s = VIGRA_CSTD::sin(x);
820  k = sq(k*s);
821  return s*(detail::ellipticRF(c2, 1.0-k, 1.0) - k/3.0*detail::ellipticRD(c2, 1.0-k, 1.0));
822 }
823 
824 #if _MSC_VER
825 
826 namespace detail {
827 
828 template <class T>
829 double erfImpl(T x)
830 {
831  double t = 1.0/(1.0+0.5*VIGRA_CSTD::fabs(x));
832  double ans = t*VIGRA_CSTD::exp(-x*x-1.26551223+t*(1.00002368+t*(0.37409196+
833  t*(0.09678418+t*(-0.18628806+t*(0.27886807+
834  t*(-1.13520398+t*(1.48851587+t*(-0.82215223+
835  t*0.17087277)))))))));
836  if (x >= 0.0)
837  return 1.0 - ans;
838  else
839  return ans - 1.0;
840 }
841 
842 } // namespace detail
843 
844  /*! The error function.
845 
846  If <tt>erf()</tt> is not provided in the C standard math library (as it should according to the
847  new C99 standard ?), VIGRA implements <tt>erf()</tt> as an approximation of the error
848  function
849 
850  \f[
851  \mbox{erf}(x) = \int_0^x e^{-t^2} dt
852  \f]
853 
854  according to the formula given in Press et al. "Numerical Recipes".
855 
856  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
857  Namespace: vigra
858  */
859 inline double erf(double x)
860 {
861  return detail::erfImpl(x);
862 }
863 
864 #else
865 
866 using ::erf;
867 
868 #endif
869 
870 namespace detail {
871 
872 template <class T>
873 double noncentralChi2CDFApprox(unsigned int degreesOfFreedom, T noncentrality, T arg)
874 {
875  double a = degreesOfFreedom + noncentrality;
876  double b = (a + noncentrality) / sq(a);
877  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);
878  return 0.5*(1.0 + erf(t/VIGRA_CSTD::sqrt(2.0)));
879 }
880 
881 template <class T>
882 void noncentralChi2OneIteration(T arg, T & lans, T & dans, T & pans, unsigned int & j)
883 {
884  double tol = -50.0;
885  if(lans < tol)
886  {
887  lans = lans + VIGRA_CSTD::log(arg / j);
888  dans = VIGRA_CSTD::exp(lans);
889  }
890  else
891  {
892  dans = dans * arg / j;
893  }
894  pans = pans - dans;
895  j += 2;
896 }
897 
898 template <class T>
899 std::pair<double, double> noncentralChi2CDF(unsigned int degreesOfFreedom, T noncentrality, T arg, T eps)
900 {
901  vigra_precondition(noncentrality >= 0.0 && arg >= 0.0 && eps > 0.0,
902  "noncentralChi2P(): parameters must be positive.");
903  if (arg == 0.0 && degreesOfFreedom > 0)
904  return std::make_pair(0.0, 0.0);
905 
906  // Determine initial values
907  double b1 = 0.5 * noncentrality,
908  ao = VIGRA_CSTD::exp(-b1),
909  eps2 = eps / ao,
910  lnrtpi2 = 0.22579135264473,
911  probability, density, lans, dans, pans, sum, am, hold;
912  unsigned int maxit = 500,
913  i, m;
914  if(degreesOfFreedom % 2)
915  {
916  i = 1;
917  lans = -0.5 * (arg + VIGRA_CSTD::log(arg)) - lnrtpi2;
918  dans = VIGRA_CSTD::exp(lans);
919  pans = erf(VIGRA_CSTD::sqrt(arg/2.0));
920  }
921  else
922  {
923  i = 2;
924  lans = -0.5 * arg;
925  dans = VIGRA_CSTD::exp(lans);
926  pans = 1.0 - dans;
927  }
928 
929  // Evaluate first term
930  if(degreesOfFreedom == 0)
931  {
932  m = 1;
933  degreesOfFreedom = 2;
934  am = b1;
935  sum = 1.0 / ao - 1.0 - am;
936  density = am * dans;
937  probability = 1.0 + am * pans;
938  }
939  else
940  {
941  m = 0;
942  degreesOfFreedom = degreesOfFreedom - 1;
943  am = 1.0;
944  sum = 1.0 / ao - 1.0;
945  while(i < degreesOfFreedom)
946  detail::noncentralChi2OneIteration(arg, lans, dans, pans, i);
947  degreesOfFreedom = degreesOfFreedom + 1;
948  density = dans;
949  probability = pans;
950  }
951  // Evaluate successive terms of the expansion
952  for(++m; m<maxit; ++m)
953  {
954  am = b1 * am / m;
955  detail::noncentralChi2OneIteration(arg, lans, dans, pans, degreesOfFreedom);
956  sum = sum - am;
957  density = density + am * dans;
958  hold = am * pans;
959  probability = probability + hold;
960  if((pans * sum < eps2) && (hold < eps2))
961  break; // converged
962  }
963  if(m == maxit)
964  vigra_fail("noncentralChi2P(): no convergence.");
965  return std::make_pair(0.5 * ao * density, std::min(1.0, std::max(0.0, ao * probability)));
966 }
967 
968 } // namespace detail
969 
970  /*! Chi square distribution.
971 
972  Computes the density of a chi square distribution with \a degreesOfFreedom
973  and tolerance \a accuracy at the given argument \a arg
974  by calling <tt>noncentralChi2(degreesOfFreedom, 0.0, arg, accuracy)</tt>.
975 
976  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
977  Namespace: vigra
978  */
979 inline double chi2(unsigned int degreesOfFreedom, double arg, double accuracy = 1e-7)
980 {
981  return detail::noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy).first;
982 }
983 
984  /*! Cumulative chi square distribution.
985 
986  Computes the cumulative density of a chi square distribution with \a degreesOfFreedom
987  and tolerance \a accuracy at the given argument \a arg, i.e. the probability that
988  a random number drawn from the distribution is below \a arg
989  by calling <tt>noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy)</tt>.
990 
991  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
992  Namespace: vigra
993  */
994 inline double chi2CDF(unsigned int degreesOfFreedom, double arg, double accuracy = 1e-7)
995 {
996  return detail::noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy).second;
997 }
998 
999  /*! Non-central chi square distribution.
1000 
1001  Computes the density of a non-central chi square distribution with \a degreesOfFreedom,
1002  noncentrality parameter \a noncentrality and tolerance \a accuracy at the given argument
1003  \a arg. It uses Algorithm AS 231 from Appl. Statist. (1987) Vol.36, No.3 (code ported from
1004  http://lib.stat.cmu.edu/apstat/231). The algorithm has linear complexity in the number of
1005  degrees of freedom.
1006 
1007  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
1008  Namespace: vigra
1009  */
1010 inline double noncentralChi2(unsigned int degreesOfFreedom,
1011  double noncentrality, double arg, double accuracy = 1e-7)
1012 {
1013  return detail::noncentralChi2CDF(degreesOfFreedom, noncentrality, arg, accuracy).first;
1014 }
1015 
1016  /*! Cumulative non-central chi square distribution.
1017 
1018  Computes the cumulative density of a chi square distribution with \a degreesOfFreedom,
1019  noncentrality parameter \a noncentrality and tolerance \a accuracy at the given argument
1020  \a arg, i.e. the probability that a random number drawn from the distribution is below \a arg
1021  It uses Algorithm AS 231 from Appl. Statist. (1987) Vol.36, No.3 (code ported from
1022  http://lib.stat.cmu.edu/apstat/231). The algorithm has linear complexity in the number of
1023  degrees of freedom (see noncentralChi2CDFApprox() for a constant-time algorithm).
1024 
1025  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
1026  Namespace: vigra
1027  */
1028 inline double noncentralChi2CDF(unsigned int degreesOfFreedom,
1029  double noncentrality, double arg, double accuracy = 1e-7)
1030 {
1031  return detail::noncentralChi2CDF(degreesOfFreedom, noncentrality, arg, accuracy).second;
1032 }
1033 
1034  /*! Cumulative non-central chi square distribution (approximate).
1035 
1036  Computes approximate values of the cumulative density of a chi square distribution with \a degreesOfFreedom,
1037  and noncentrality parameter \a noncentrality at the given argument
1038  \a arg, i.e. the probability that a random number drawn from the distribution is below \a arg
1039  It uses the approximate transform into a normal distribution due to Wilson and Hilferty
1040  (see Abramovitz, Stegun: "Handbook of Mathematical Functions", formula 26.3.32).
1041  The algorithm's running time is independent of the inputs, i.e. is should be used
1042  when noncentralChi2CDF() is too slow, and approximate values are sufficient. The accuracy is only
1043  about 0.1 for few degrees of freedom, but reaches about 0.001 above dof = 5.
1044 
1045  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
1046  Namespace: vigra
1047  */
1048 inline double noncentralChi2CDFApprox(unsigned int degreesOfFreedom, double noncentrality, double arg)
1049 {
1050  return detail::noncentralChi2CDFApprox(degreesOfFreedom, noncentrality, arg);
1051 }
1052 
1053 
1054 
1055 namespace detail {
1056 
1057 // both f1 and f2 are unsigned here
1058 template<class FPT>
1059 inline
1060 FPT safeFloatDivision( FPT f1, FPT f2 )
1061 {
1062  return f2 < NumericTraits<FPT>::one() && f1 > f2 * NumericTraits<FPT>::max()
1063  ? NumericTraits<FPT>::max()
1064  : (f2 > NumericTraits<FPT>::one() && f1 < f2 * NumericTraits<FPT>::smallestPositive()) ||
1065  f1 == NumericTraits<FPT>::zero()
1066  ? NumericTraits<FPT>::zero()
1067  : f1/f2;
1068 }
1069 
1070 } // namespace detail
1071 
1072  /*! Tolerance based floating-point comparison.
1073 
1074  Check whether two floating point numbers are equal within the given tolerance.
1075  This is useful because floating point numbers that should be equal in theory are
1076  rarely exactly equal in practice. If the tolerance \a epsilon is not given,
1077  twice the machine epsilon is used.
1078 
1079  <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br>
1080  Namespace: vigra
1081  */
1082 template <class T1, class T2>
1083 bool
1084 closeAtTolerance(T1 l, T2 r, typename PromoteTraits<T1, T2>::Promote epsilon)
1085 {
1086  typedef typename PromoteTraits<T1, T2>::Promote T;
1087  if(l == 0.0)
1088  return VIGRA_CSTD::fabs(r) <= epsilon;
1089  if(r == 0.0)
1090  return VIGRA_CSTD::fabs(l) <= epsilon;
1091  T diff = VIGRA_CSTD::fabs( l - r );
1092  T d1 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( r ) );
1093  T d2 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( l ) );
1094 
1095  return (d1 <= epsilon && d2 <= epsilon);
1096 }
1097 
1098 template <class T1, class T2>
1099 inline bool closeAtTolerance(T1 l, T2 r)
1100 {
1101  typedef typename PromoteTraits<T1, T2>::Promote T;
1102  return closeAtTolerance(l, r, 2.0 * NumericTraits<T>::epsilon());
1103 }
1104 
1105 //@}
1106 
1107 } // namespace vigra
1108 
1109 #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 (Tue Jul 10 2012)