[ 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.3.3, Aug 18 2005 ) */ 00008 /* You may use, modify, and distribute this software according */ 00009 /* to the terms stated in the LICENSE file included in */ 00010 /* the VIGRA distribution. */ 00011 /* */ 00012 /* The VIGRA Website is */ 00013 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */ 00014 /* Please direct questions, bug reports, and contributions to */ 00015 /* koethe@informatik.uni-hamburg.de */ 00016 /* */ 00017 /* THIS SOFTWARE IS PROVIDED AS IS AND WITHOUT ANY EXPRESS OR */ 00018 /* IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED */ 00019 /* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */ 00020 /* */ 00021 /************************************************************************/ 00022 00023 #ifndef VIGRA_MATHUTIL_HXX 00024 #define VIGRA_MATHUTIL_HXX 00025 00026 #include <cmath> 00027 #include <cstdlib> 00028 #include "vigra/config.hxx" 00029 #include "vigra/numerictraits.hxx" 00030 00031 /*! \page MathConstants Mathematical Constants 00032 00033 <TT>M_PI, M_SQRT2</TT> 00034 00035 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>" 00036 00037 Since <TT>M_PI</TT> and <TT>M_SQRT2</TT> are not officially standardized, 00038 we provide definitions here for those compilers that don't support them. 00039 00040 \code 00041 #ifndef M_PI 00042 # define M_PI 3.14159265358979323846 00043 #endif 00044 00045 #ifndef M_SQRT2 00046 # define M_SQRT2 1.41421356237309504880 00047 #endif 00048 \endcode 00049 */ 00050 #ifndef M_PI 00051 # define M_PI 3.14159265358979323846 00052 #endif 00053 00054 #ifndef M_SQRT2 00055 # define M_SQRT2 1.41421356237309504880 00056 #endif 00057 00058 namespace vigra { 00059 00060 #ifndef __sun 00061 00062 /** \addtogroup MathFunctions Mathematical Functions 00063 00064 Useful mathematical functions and functors. 00065 */ 00066 //@{ 00067 /*! The error function. 00068 00069 If <tt>erf()</tt> is not provided in the C standard math library (as it should according to the 00070 new C99 standard ?), VIGRA implements <tt>erf()</tt> as an approximation of the error 00071 function 00072 00073 \f[ 00074 \mbox{erf}(x) = \int_0^x e^{-x^2} dx 00075 \f] 00076 00077 according to the formula given in Press et al. "Numerical Recipes". 00078 00079 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00080 Namespace: vigra 00081 */ 00082 template <class T> 00083 double erf(T x) 00084 { 00085 double t = 1.0/(1.0+0.5*VIGRA_CSTD::fabs(x)); 00086 double ans = t*VIGRA_CSTD::exp(-x*x-1.26551223+t*(1.00002368+t*(0.37409196+ 00087 t*(0.09678418+t*(-0.18628806+t*(0.27886807+ 00088 t*(-1.13520398+t*(1.48851587+t*(-0.82215223+ 00089 t*0.17087277))))))))); 00090 if (x >= 0.0) 00091 return 1.0 - ans; 00092 else 00093 return ans - 1.0; 00094 } 00095 00096 #else 00097 00098 using VIGRA_CSTD::erf; 00099 00100 #endif 00101 00102 // import functions into namespace vigra which VIGRA is going to overload 00103 00104 using VIGRA_CSTD::pow; 00105 using VIGRA_CSTD::floor; 00106 using VIGRA_CSTD::ceil; 00107 using std::abs; 00108 00109 #define VIGRA_DEFINE_UNSIGNED_ABS(T) \ 00110 inline T abs(T t) { return t; } 00111 00112 VIGRA_DEFINE_UNSIGNED_ABS(bool) 00113 VIGRA_DEFINE_UNSIGNED_ABS(unsigned char) 00114 VIGRA_DEFINE_UNSIGNED_ABS(unsigned short) 00115 VIGRA_DEFINE_UNSIGNED_ABS(unsigned int) 00116 VIGRA_DEFINE_UNSIGNED_ABS(unsigned long) 00117 00118 #undef VIGRA_DEFINE_UNSIGNED_ABS 00119 00120 /*! The rounding function. 00121 00122 Defined for all floating point types. Rounds towards the nearest integer for both 00123 positive and negative inputs. 00124 00125 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00126 Namespace: vigra 00127 */ 00128 inline float round(float t) 00129 { 00130 return t >= 0.0 00131 ? floor(t + 0.5) 00132 : ceil(t - 0.5); 00133 } 00134 00135 inline double round(double t) 00136 { 00137 return t >= 0.0 00138 ? floor(t + 0.5) 00139 : ceil(t - 0.5); 00140 } 00141 00142 inline long double round(long double t) 00143 { 00144 return t >= 0.0 00145 ? floor(t + 0.5) 00146 : ceil(t - 0.5); 00147 } 00148 00149 /*! The square function. 00150 00151 sq(x) is needed so often that it makes sense to define it as a function. 00152 00153 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00154 Namespace: vigra 00155 */ 00156 template <class T> 00157 inline 00158 typename NumericTraits<T>::Promote sq(T t) 00159 { 00160 return t*t; 00161 } 00162 00163 #ifdef VIGRA_NO_HYPOT 00164 /*! Compute the Euclidean distance (length of the hypothenuse of a right-angled triangle). 00165 00166 The hypot() function returns the sqrt(a*a + b*b). 00167 It is implemented in a way that minimizes round-off error. 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 T hypot(T a, T b) 00174 { 00175 T absa = abs(a), absb = abs(b); 00176 if (absa > absb) 00177 return absa * VIGRA_CSTD::sqrt(1.0 + sq(absb/absa)); 00178 else 00179 return absb == NumericTraits<T>::zero() 00180 ? NumericTraits<T>::zero() 00181 : absb * VIGRA_CSTD::sqrt(1.0 + sq(absa/absb)); 00182 } 00183 00184 #else 00185 00186 using ::hypot; 00187 00188 #endif 00189 00190 /*! The sign function. 00191 00192 Returns 1, 0, or -1 depending on the sign of \a t. 00193 00194 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00195 Namespace: vigra 00196 */ 00197 template <class T> 00198 T sign(T t) 00199 { 00200 return t > NumericTraits<T>::zero() 00201 ? NumericTraits<T>::one() 00202 : t < NumericTraits<T>::zero() 00203 ? -NumericTraits<T>::one() 00204 : NumericTraits<T>::zero(); 00205 } 00206 00207 /*! The binary sign function. 00208 00209 Transfers the sign of \a t2 to \a t1. 00210 00211 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00212 Namespace: vigra 00213 */ 00214 template <class T1, class T2> 00215 T1 sign(T1 t1, T2 t2) 00216 { 00217 return t2 >= NumericTraits<T2>::zero() 00218 ? abs(t1) 00219 : -abs(t1); 00220 } 00221 00222 #define VIGRA_DEFINE_NORM(T) \ 00223 inline NormTraits<T>::SquaredNormType squaredNorm(T t) { return sq(t); } \ 00224 inline NormTraits<T>::NormType norm(T t) { return abs(t); } 00225 00226 VIGRA_DEFINE_NORM(bool) 00227 VIGRA_DEFINE_NORM(signed char) 00228 VIGRA_DEFINE_NORM(unsigned char) 00229 VIGRA_DEFINE_NORM(short) 00230 VIGRA_DEFINE_NORM(unsigned short) 00231 VIGRA_DEFINE_NORM(int) 00232 VIGRA_DEFINE_NORM(unsigned int) 00233 VIGRA_DEFINE_NORM(long) 00234 VIGRA_DEFINE_NORM(unsigned long) 00235 VIGRA_DEFINE_NORM(float) 00236 VIGRA_DEFINE_NORM(double) 00237 VIGRA_DEFINE_NORM(long double) 00238 00239 #undef VIGRA_DEFINE_NORM 00240 00241 template <class T> 00242 inline typename NormTraits<std::complex<T> >::SquaredNormType 00243 squaredNorm(std::complex<T> const & t) 00244 { 00245 return sq(t.real()) + sq(t.imag()); 00246 } 00247 00248 #ifdef DOXYGEN // only for documentation 00249 /*! The squared norm of a numerical object. 00250 00251 For scalar types: equals <tt>vigra::sq(t)</tt><br>. 00252 For vectorial types: equals <tt>vigra::dot(t, t)</tt><br>. 00253 For complex types: equals <tt>vigra::sq(t.real()) + vigra::sq(t.imag())</tt><br>. 00254 For matrix types: results in the squared Frobenius norm (sum of squares of the matrix elements). 00255 */ 00256 NormTraits<T>::SquaredNormType squaredNorm(T const & t); 00257 00258 #endif 00259 00260 /*! The norm of a numerical object. 00261 00262 For scalar types: implemented as <tt>abs(t)</tt><br> 00263 otherwise: implemented as <tt>sqrt(squaredNorm(t))</tt>. 00264 00265 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00266 Namespace: vigra 00267 */ 00268 template <class T> 00269 inline typename NormTraits<T>::NormType 00270 norm(T const & t) 00271 { 00272 return VIGRA_CSTD::sqrt(static_cast<typename NormTraits<T>::NormType>(squaredNorm(t))); 00273 } 00274 00275 namespace detail { 00276 00277 // both f1 and f2 are unsigned here 00278 template<class FPT> 00279 inline 00280 FPT safeFloatDivision( FPT f1, FPT f2 ) 00281 { 00282 return f2 < NumericTraits<FPT>::one() && f1 > f2 * NumericTraits<FPT>::max() 00283 ? NumericTraits<FPT>::max() 00284 : (f2 > NumericTraits<FPT>::one() && f1 < f2 * NumericTraits<FPT>::smallestPositive()) || 00285 f1 == NumericTraits<FPT>::zero() 00286 ? NumericTraits<FPT>::zero() 00287 : f1/f2; 00288 } 00289 00290 } // namespace detail 00291 00292 /*! Tolerance based floating-point comparison. 00293 00294 Check whether two floating point numbers are equal within the given tolerance. 00295 This is useful because floating point numbers that should be equal in theory are 00296 rarely exactly equal in practice. If the tolerance \a epsilon is not given, 00297 twice the machine epsilon is used. 00298 00299 <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br> 00300 Namespace: vigra 00301 */ 00302 template <class T1, class T2> 00303 bool closeAtTolerance(T1 l, T2 r, typename PromoteTraits<T1, T2>::Promote epsilon) 00304 { 00305 typedef typename PromoteTraits<T1, T2>::Promote T; 00306 if(l == 0.0 && r != 0.0) 00307 return VIGRA_CSTD::fabs(r) <= epsilon; 00308 if(l != 0.0 && r == 0.0) 00309 return VIGRA_CSTD::fabs(r) <= epsilon; 00310 T diff = VIGRA_CSTD::fabs( l - r ); 00311 T d1 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( r ) ); 00312 T d2 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( l ) ); 00313 00314 return (d1 <= epsilon && d2 <= epsilon); 00315 } 00316 00317 template <class T1, class T2> 00318 bool closeAtTolerance(T1 l, T2 r) 00319 { 00320 typedef typename PromoteTraits<T1, T2>::Promote T; 00321 return closeAtTolerance(l, r, 2.0 * NumericTraits<T>::epsilon()); 00322 } 00323 00324 //@} 00325 00326 } // namespace vigra 00327 00328 #endif /* VIGRA_MATHUTIL_HXX */
© Ullrich Köthe (koethe@informatik.uni-hamburg.de) |
html generated using doxygen and Python
|