[ VIGRA Homepage | Class Index | Function Index | File Index | Main Page ]
![]() |
vigra/mathutil.hxx | ![]() |
---|
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 1998-2005 by Ullrich Koethe */ 00004 /* Cognitive Systems Group, University of Hamburg, Germany */ 00005 /* */ 00006 /* This file is part of the VIGRA computer vision library. */ 00007 /* ( Version 1.4.0, Dec 21 2005 ) */ 00008 /* The VIGRA Website is */ 00009 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */ 00010 /* Please direct questions, bug reports, and contributions to */ 00011 /* koethe@informatik.uni-hamburg.de or */ 00012 /* vigra@kogs1.informatik.uni-hamburg.de */ 00013 /* */ 00014 /* Permission is hereby granted, free of charge, to any person */ 00015 /* obtaining a copy of this software and associated documentation */ 00016 /* files (the "Software"), to deal in the Software without */ 00017 /* restriction, including without limitation the rights to use, */ 00018 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00019 /* sell copies of the Software, and to permit persons to whom the */ 00020 /* Software is furnished to do so, subject to the following */ 00021 /* conditions: */ 00022 /* */ 00023 /* The above copyright notice and this permission notice shall be */ 00024 /* included in all copies or substantial portions of the */ 00025 /* Software. */ 00026 /* */ 00027 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00028 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00029 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00030 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00031 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00032 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00033 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00034 /* OTHER DEALINGS IN THE SOFTWARE. */ 00035 /* */ 00036 /************************************************************************/ 00037 00038 #ifndef VIGRA_MATHUTIL_HXX 00039 #define VIGRA_MATHUTIL_HXX 00040 00041 #include <cmath> 00042 #include <cstdlib> 00043 #include "vigra/config.hxx" 00044 #include "vigra/sized_int.hxx" 00045 #include "vigra/numerictraits.hxx" 00046 00047 /*! \page MathConstants Mathematical Constants 00048 00049 <TT>M_PI, M_SQRT2</TT> 00050 00051 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>" 00052 00053 Since <TT>M_PI</TT> and <TT>M_SQRT2</TT> are not officially standardized, 00054 we provide definitions here for those compilers that don't support them. 00055 00056 \code 00057 #ifndef M_PI 00058 # define M_PI 3.14159265358979323846 00059 #endif 00060 00061 #ifndef M_SQRT2 00062 # define M_SQRT2 1.41421356237309504880 00063 #endif 00064 \endcode 00065 */ 00066 #ifndef M_PI 00067 # define M_PI 3.14159265358979323846 00068 #endif 00069 00070 #ifndef M_SQRT2 00071 # define M_SQRT2 1.41421356237309504880 00072 #endif 00073 00074 namespace vigra { 00075 00076 #ifndef __sun 00077 00078 /** \addtogroup MathFunctions Mathematical Functions 00079 00080 Useful mathematical functions and functors. 00081 */ 00082 //@{ 00083 /*! The error function. 00084 00085 If <tt>erf()</tt> is not provided in the C standard math library (as it should according to the 00086 new C99 standard ?), VIGRA implements <tt>erf()</tt> as an approximation of the error 00087 function 00088 00089 \f[ 00090 \mbox{erf}(x) = \int_0^x e^{-x^2} dx 00091 \f] 00092 00093 according to the formula given in Press et al. "Numerical Recipes". 00094 00095 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00096 Namespace: vigra 00097 */ 00098 template <class T> 00099 double erf(T x) 00100 { 00101 double t = 1.0/(1.0+0.5*VIGRA_CSTD::fabs(x)); 00102 double ans = t*VIGRA_CSTD::exp(-x*x-1.26551223+t*(1.00002368+t*(0.37409196+ 00103 t*(0.09678418+t*(-0.18628806+t*(0.27886807+ 00104 t*(-1.13520398+t*(1.48851587+t*(-0.82215223+ 00105 t*0.17087277))))))))); 00106 if (x >= 0.0) 00107 return 1.0 - ans; 00108 else 00109 return ans - 1.0; 00110 } 00111 00112 #else 00113 00114 using VIGRA_CSTD::erf; 00115 00116 #endif 00117 00118 // import functions into namespace vigra which VIGRA is going to overload 00119 00120 using VIGRA_CSTD::pow; 00121 using VIGRA_CSTD::floor; 00122 using VIGRA_CSTD::ceil; 00123 using std::abs; 00124 00125 #define VIGRA_DEFINE_UNSIGNED_ABS(T) \ 00126 inline T abs(T t) { return t; } 00127 00128 VIGRA_DEFINE_UNSIGNED_ABS(bool) 00129 VIGRA_DEFINE_UNSIGNED_ABS(unsigned char) 00130 VIGRA_DEFINE_UNSIGNED_ABS(unsigned short) 00131 VIGRA_DEFINE_UNSIGNED_ABS(unsigned int) 00132 VIGRA_DEFINE_UNSIGNED_ABS(unsigned long) 00133 00134 #undef VIGRA_DEFINE_UNSIGNED_ABS 00135 00136 /*! The rounding function. 00137 00138 Defined for all floating point types. Rounds towards the nearest integer for both 00139 positive and negative inputs. 00140 00141 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00142 Namespace: vigra 00143 */ 00144 inline float round(float t) 00145 { 00146 return t >= 0.0 00147 ? floor(t + 0.5) 00148 : ceil(t - 0.5); 00149 } 00150 00151 inline double round(double t) 00152 { 00153 return t >= 0.0 00154 ? floor(t + 0.5) 00155 : ceil(t - 0.5); 00156 } 00157 00158 inline long double round(long double t) 00159 { 00160 return t >= 0.0 00161 ? floor(t + 0.5) 00162 : ceil(t - 0.5); 00163 } 00164 00165 /*! The square function. 00166 00167 sq(x) is needed so often that it makes sense to define it as a function. 00168 00169 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00170 Namespace: vigra 00171 */ 00172 template <class T> 00173 inline 00174 typename NumericTraits<T>::Promote sq(T t) 00175 { 00176 return t*t; 00177 } 00178 00179 namespace detail { 00180 00181 template <class T> 00182 class IntSquareRoot 00183 { 00184 public: 00185 static int sqq_table[]; 00186 static UInt32 exec(UInt32 v); 00187 }; 00188 00189 template <class T> 00190 int IntSquareRoot<T>::sqq_table[] = { 00191 0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53, 55, 57, 00192 59, 61, 64, 65, 67, 69, 71, 73, 75, 76, 78, 80, 81, 83, 00193 84, 86, 87, 89, 90, 91, 93, 94, 96, 97, 98, 99, 101, 102, 00194 103, 104, 106, 107, 108, 109, 110, 112, 113, 114, 115, 116, 117, 118, 00195 119, 120, 121, 122, 123, 124, 125, 126, 128, 128, 129, 130, 131, 132, 00196 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 144, 145, 00197 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155, 156, 157, 00198 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168, 00199 169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178, 00200 179, 180, 181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188, 00201 189, 189, 190, 191, 192, 192, 193, 193, 194, 195, 195, 196, 197, 197, 00202 198, 199, 199, 200, 201, 201, 202, 203, 203, 204, 204, 205, 206, 206, 00203 207, 208, 208, 209, 209, 210, 211, 211, 212, 212, 213, 214, 214, 215, 00204 215, 216, 217, 217, 218, 218, 219, 219, 220, 221, 221, 222, 222, 223, 00205 224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230, 230, 231, 00206 231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238, 00207 239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246, 00208 246, 247, 247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253, 00209 253, 254, 254, 255 00210 }; 00211 00212 template <class T> 00213 UInt32 IntSquareRoot<T>::exec(UInt32 x) 00214 { 00215 unsigned long xn; 00216 if (x >= 0x10000) 00217 if (x >= 0x1000000) 00218 if (x >= 0x10000000) 00219 if (x >= 0x40000000) { 00220 if (x >= (UInt32)65535*(UInt32)65535) 00221 return 65535; 00222 xn = sqq_table[x>>24] << 8; 00223 } else 00224 xn = sqq_table[x>>22] << 7; 00225 else 00226 if (x >= 0x4000000) 00227 xn = sqq_table[x>>20] << 6; 00228 else 00229 xn = sqq_table[x>>18] << 5; 00230 else { 00231 if (x >= 0x100000) 00232 if (x >= 0x400000) 00233 xn = sqq_table[x>>16] << 4; 00234 else 00235 xn = sqq_table[x>>14] << 3; 00236 else 00237 if (x >= 0x40000) 00238 xn = sqq_table[x>>12] << 2; 00239 else 00240 xn = sqq_table[x>>10] << 1; 00241 00242 goto nr1; 00243 } 00244 else 00245 if (x >= 0x100) { 00246 if (x >= 0x1000) 00247 if (x >= 0x4000) 00248 xn = (sqq_table[x>>8] >> 0) + 1; 00249 else 00250 xn = (sqq_table[x>>6] >> 1) + 1; 00251 else 00252 if (x >= 0x400) 00253 xn = (sqq_table[x>>4] >> 2) + 1; 00254 else 00255 xn = (sqq_table[x>>2] >> 3) + 1; 00256 00257 goto adj; 00258 } else 00259 return sqq_table[x] >> 4; 00260 00261 /* Run two iterations of the standard convergence formula */ 00262 00263 xn = (xn + 1 + x / xn) / 2; 00264 nr1: 00265 xn = (xn + 1 + x / xn) / 2; 00266 adj: 00267 00268 if (xn * xn > x) /* Correct rounding if necessary */ 00269 xn--; 00270 00271 return xn; 00272 } 00273 00274 } // namespace detail 00275 00276 using VIGRA_CSTD::sqrt; 00277 00278 /*! Signed integer square root. 00279 00280 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00281 Namespace: vigra 00282 */ 00283 inline Int32 sqrti(Int32 v) 00284 { 00285 if(v < 0) 00286 throw std::domain_error("sqrti(Int32): negative argument."); 00287 return (Int32)detail::IntSquareRoot<UInt32>::exec((UInt32)v); 00288 } 00289 00290 /*! Unsigned integer square root. 00291 00292 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00293 Namespace: vigra 00294 */ 00295 inline UInt32 sqrti(UInt32 v) 00296 { 00297 return detail::IntSquareRoot<UInt32>::exec(v); 00298 } 00299 00300 #ifdef VIGRA_NO_HYPOT 00301 /*! Compute the Euclidean distance (length of the hypothenuse of a right-angled triangle). 00302 00303 The hypot() function returns the sqrt(a*a + b*b). 00304 It is implemented in a way that minimizes round-off error. 00305 00306 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00307 Namespace: vigra 00308 */ 00309 template <class T> 00310 T hypot(T a, T b) 00311 { 00312 T absa = abs(a), absb = abs(b); 00313 if (absa > absb) 00314 return absa * VIGRA_CSTD::sqrt(1.0 + sq(absb/absa)); 00315 else 00316 return absb == NumericTraits<T>::zero() 00317 ? NumericTraits<T>::zero() 00318 : absb * VIGRA_CSTD::sqrt(1.0 + sq(absa/absb)); 00319 } 00320 00321 #else 00322 00323 using ::hypot; 00324 00325 #endif 00326 00327 00328 00329 /*! The sign function. 00330 00331 Returns 1, 0, or -1 depending on the sign of \a t. 00332 00333 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00334 Namespace: vigra 00335 */ 00336 template <class T> 00337 T sign(T t) 00338 { 00339 return t > NumericTraits<T>::zero() 00340 ? NumericTraits<T>::one() 00341 : t < NumericTraits<T>::zero() 00342 ? -NumericTraits<T>::one() 00343 : NumericTraits<T>::zero(); 00344 } 00345 00346 /*! The binary sign function. 00347 00348 Transfers the sign of \a t2 to \a t1. 00349 00350 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00351 Namespace: vigra 00352 */ 00353 template <class T1, class T2> 00354 T1 sign(T1 t1, T2 t2) 00355 { 00356 return t2 >= NumericTraits<T2>::zero() 00357 ? abs(t1) 00358 : -abs(t1); 00359 } 00360 00361 #define VIGRA_DEFINE_NORM(T) \ 00362 inline NormTraits<T>::SquaredNormType squaredNorm(T t) { return sq(t); } \ 00363 inline NormTraits<T>::NormType norm(T t) { return abs(t); } 00364 00365 VIGRA_DEFINE_NORM(bool) 00366 VIGRA_DEFINE_NORM(signed char) 00367 VIGRA_DEFINE_NORM(unsigned char) 00368 VIGRA_DEFINE_NORM(short) 00369 VIGRA_DEFINE_NORM(unsigned short) 00370 VIGRA_DEFINE_NORM(int) 00371 VIGRA_DEFINE_NORM(unsigned int) 00372 VIGRA_DEFINE_NORM(long) 00373 VIGRA_DEFINE_NORM(unsigned long) 00374 VIGRA_DEFINE_NORM(float) 00375 VIGRA_DEFINE_NORM(double) 00376 VIGRA_DEFINE_NORM(long double) 00377 00378 #undef VIGRA_DEFINE_NORM 00379 00380 template <class T> 00381 inline typename NormTraits<std::complex<T> >::SquaredNormType 00382 squaredNorm(std::complex<T> const & t) 00383 { 00384 return sq(t.real()) + sq(t.imag()); 00385 } 00386 00387 #ifdef DOXYGEN // only for documentation 00388 /*! The squared norm of a numerical object. 00389 00390 For scalar types: equals <tt>vigra::sq(t)</tt><br>. 00391 For vectorial types: equals <tt>vigra::dot(t, t)</tt><br>. 00392 For complex types: equals <tt>vigra::sq(t.real()) + vigra::sq(t.imag())</tt><br>. 00393 For matrix types: results in the squared Frobenius norm (sum of squares of the matrix elements). 00394 */ 00395 NormTraits<T>::SquaredNormType squaredNorm(T const & t); 00396 00397 #endif 00398 00399 /*! The norm of a numerical object. 00400 00401 For scalar types: implemented as <tt>abs(t)</tt><br> 00402 otherwise: implemented as <tt>sqrt(squaredNorm(t))</tt>. 00403 00404 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00405 Namespace: vigra 00406 */ 00407 template <class T> 00408 inline typename NormTraits<T>::NormType 00409 norm(T const & t) 00410 { 00411 typedef typename NormTraits<T>::SquaredNormType SNT; 00412 return sqrt(static_cast<typename SquareRootTraits<SNT>::SquareRootArgument>(squaredNorm(t))); 00413 } 00414 00415 namespace detail { 00416 00417 // both f1 and f2 are unsigned here 00418 template<class FPT> 00419 inline 00420 FPT safeFloatDivision( FPT f1, FPT f2 ) 00421 { 00422 return f2 < NumericTraits<FPT>::one() && f1 > f2 * NumericTraits<FPT>::max() 00423 ? NumericTraits<FPT>::max() 00424 : (f2 > NumericTraits<FPT>::one() && f1 < f2 * NumericTraits<FPT>::smallestPositive()) || 00425 f1 == NumericTraits<FPT>::zero() 00426 ? NumericTraits<FPT>::zero() 00427 : f1/f2; 00428 } 00429 00430 } // namespace detail 00431 00432 /*! Tolerance based floating-point comparison. 00433 00434 Check whether two floating point numbers are equal within the given tolerance. 00435 This is useful because floating point numbers that should be equal in theory are 00436 rarely exactly equal in practice. If the tolerance \a epsilon is not given, 00437 twice the machine epsilon is used. 00438 00439 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00440 Namespace: vigra 00441 */ 00442 template <class T1, class T2> 00443 bool closeAtTolerance(T1 l, T2 r, typename PromoteTraits<T1, T2>::Promote epsilon) 00444 { 00445 typedef typename PromoteTraits<T1, T2>::Promote T; 00446 if(l == 0.0 && r != 0.0) 00447 return VIGRA_CSTD::fabs(r) <= epsilon; 00448 if(l != 0.0 && r == 0.0) 00449 return VIGRA_CSTD::fabs(r) <= epsilon; 00450 T diff = VIGRA_CSTD::fabs( l - r ); 00451 T d1 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( r ) ); 00452 T d2 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( l ) ); 00453 00454 return (d1 <= epsilon && d2 <= epsilon); 00455 } 00456 00457 template <class T1, class T2> 00458 bool closeAtTolerance(T1 l, T2 r) 00459 { 00460 typedef typename PromoteTraits<T1, T2>::Promote T; 00461 return closeAtTolerance(l, r, 2.0 * NumericTraits<T>::epsilon()); 00462 } 00463 00464 //@} 00465 00466 } // namespace vigra 00467 00468 #endif /* VIGRA_MATHUTIL_HXX */
© Ullrich Köthe (koethe@informatik.uni-hamburg.de) |
html generated using doxygen and Python
|