OpenVDB  1.1.0
Math.h
Go to the documentation of this file.
1 
2 //
3 // Copyright (c) 2012-2013 DreamWorks Animation LLC
4 //
5 // All rights reserved. This software is distributed under the
6 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
7 //
8 // Redistributions of source code must retain the above copyright
9 // and license notice and the following restrictions and disclaimer.
10 //
11 // * Neither the name of DreamWorks Animation nor the names of
12 // its contributors may be used to endorse or promote products derived
13 // from this software without specific prior written permission.
14 //
15 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
19 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL,
20 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
21 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
22 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
23 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
25 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 // IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE
27 // LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00.
28 //
30 //
34 
35 #ifndef OPENVDB_MATH_HAS_BEEN_INCLUDED
36 #define OPENVDB_MATH_HAS_BEEN_INCLUDED
37 
38 #include <assert.h>
39 #include <algorithm> //for std::max
40 #include <cmath> //for floor, ceil and sqrt
41 #include <math.h> //for pow, fabs(float,double,long double) etc
42 #include <cstdlib> //for srand, abs(int)
43 #include <limits> //for std::numeric_limits<Type>::max()
44 #include <string>
45 #include <boost/numeric/conversion/conversion_traits.hpp>
46 #include <openvdb/Platform.h>
47 #include <openvdb/version.h>
48 
49 // Compile pragmas
50 
51 // Intel(r) compiler fires remark #1572: floating-point equality and inequality
52 // comparisons are unrealiable when == or != is used with floating point operands.
53 #if defined(__INTEL_COMPILER)
54  #define OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN \
55  _Pragma("warning (push)") \
56  _Pragma("warning (disable:1572)")
57  #define OPENVDB_NO_FP_EQUALITY_WARNING_END \
58  _Pragma("warning (pop)")
59 #else
60  // For GCC, #pragma GCC diagnostic ignored "-Wfloat-equal"
61  // isn't working until gcc 4.2+,
62  // Trying
63  // #pragma GCC system_header
64  // creates other problems, most notably "warning: will never be executed"
65  // in from templates, unsure of how to work around.
66  // If necessary, could use integer based comparisons for equality
67  #define OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN
68  #define OPENVDB_NO_FP_EQUALITY_WARNING_END
69 #endif
70 
71 namespace openvdb {
73 namespace OPENVDB_VERSION_NAME {
74 
79 template<typename T> inline T zeroVal() { return T(0); }
81 template<> inline std::string zeroVal<std::string>() { return ""; }
83 template<> inline bool zeroVal<bool>() { return false; }
84 
85 template<typename T> inline T toleranceValue() { return T(1e-8); }
86 template<> inline float toleranceValue<float>() { return float(1e-6); }
87 
88 
90 
91 
92 
93 inline std::string operator+(const std::string& s, bool) { return s; }
94 inline std::string operator+(const std::string& s, int) { return s; }
95 inline std::string operator+(const std::string& s, float) { return s; }
96 inline std::string operator+(const std::string& s, double) { return s; }
98 
99 
103 template<typename T> inline T negative(const T& val) { return T(-val); }
104 template<> inline bool negative(const bool& val) { return !val; }
106 template<> inline std::string negative(const std::string& val) { return val; }
107 
108 
109 namespace math {
110 
112 
114 inline void randSeed(unsigned int seed)
115 {
116  srand(seed);
117 }
118 
120 inline double randUniform()
121 {
122  return (double)(rand() / (RAND_MAX + 1.0));
123 }
124 
127 {
128  protected:
129  int my_min, my_range;
130  public:
131  RandomInt(unsigned int seed, int min, int max) : my_min(min), my_range(max-min+1) {
132  assert(min<max && "RandomInt: invalid arguments");
133  randSeed(seed);
134  }
135  void setRange(int min, int max) {my_min=min; my_range=max-min+1;}
136  int operator() (void) const {return rand() % my_range + my_min;}
137  int operator() (int min, int max) const {return rand() % (max-min+1) + min;}
138 };
139 
140 
141 // ==========> Clamp/Abs <==================
142 
144 template <typename Type>
145 inline Type Clamp(Type x, Type min, Type max) {
146  assert(min<max);
147  return x > min ? x < max ? x : max : min;
148 }
149 
151 template <class Type>
152 inline Type Clamp01(Type x) {
153  return x > Type(0) ? x < Type(1) ? x : Type(1) : Type(0);
154 }
156 template <class Type>
157 inline bool ClampTest01(Type &x) {
158  if (x>=Type(0) && x<=Type(1)) return false;
159  x = x< Type(0) ? Type(0) : Type(1);
160  return true;
161 }
162 
164 template <class Type>
165 inline Type SmoothUnitStep(Type x, Type min, Type max) {
166  assert(min<max);
167  const Type t = (x-min)/(max-min);
168  return t > 0 ? t < 1 ? (3-2*t)*t*t : Type(1) : Type(0);
169 }
170 
172 inline int32_t Abs(int32_t i)
173 {
174  return abs(i);
175 }
176 
178 inline int64_t Abs(int64_t i)
179 {
180 #ifdef _MSC_VER
181  return (i < int64_t(0) ? -i : i);
182 #else
183  return abs(i);
184 #endif
185 }
186 
188 inline float Abs(float x)
189 {
190  return fabs(x);
191 }
192 
194 inline double Abs(double x)
195 {
196  return fabs(x);
197 }
198 
200 inline long double Abs(long double x)
201 {
202  return fabs(x);
203 }
204 
206 inline uint32_t Abs(uint32_t i)
207 {
208  return i;
209 }
210 
212 inline uint64_t Abs(uint64_t i)
213 {
214  return i;
215 }
217 
218 
219 template<typename Type>
220 inline bool
221 isZero(const Type& x)
222 {
224  return x == zeroVal<Type>();
226 }
227 
228 template<typename Type>
229 inline bool
230 isNegative(const Type& x)
231 {
232  return x < zeroVal<Type>();
233 }
234 
235 template<typename Type>
236 inline bool
237 isApproxEqual(const Type& a, const Type& b)
238 {
239  const Type tolerance = Type(zeroVal<Type>() + toleranceValue<Type>());
240  return !(Abs(a - b) > tolerance);
241 }
242 
243 template<typename Type>
244 inline bool
245 isApproxEqual(const Type& a, const Type& b, const Type& tolerance)
246 {
247  return !(Abs(a - b) > tolerance);
248 }
249 
250 #define OPENVDB_EXACT_IS_APPROX_EQUAL(T) \
251  template<> inline bool isApproxEqual<T>(const T& a, const T& b) { return a == b; } \
252  template<> inline bool isApproxEqual<T>(const T& a, const T& b, const T&) { return a == b; } \
253 
254 
257 
258 
259 template<typename T0, typename T1>
260 inline bool
261 isExactlyEqual(const T0& a, const T1& b)
262 {
264  return a == b;
266 }
267 
268 
269 template<typename Type>
270 inline bool
271 isRelOrApproxEqual(const Type& a, const Type& b, const Type& absTol, const Type& relTol)
272 {
273  // First check to see if we are inside the absolute tolerance
274  // Necessary for numbers close to 0
275  if (!(Abs(a - b) > absTol)) return true;
276 
277  // Next check to see if we are inside the relative tolerance
278  // to handle large numbers that aren't within the abs tolerance
279  // but could be the closest floating point representation
280  double relError;
281  if (Abs(b) > Abs(a)) {
282  relError = Abs((a - b) / b);
283  } else {
284  relError = Abs((a - b) / a);
285  }
286  return (relError <= relTol);
287 }
288 
289 template<>
290 inline bool
291 isRelOrApproxEqual(const bool& a, const bool& b, const bool&, const bool&)
292 {
293  return (a == b);
294 }
295 
296 
298 
299 
300 // Avoid strict aliasing issues by using type punning
301 // http://cellperformance.beyond3d.com/articles/2006/06/understanding-strict-aliasing.html
302 // Using "casting through a union(2)"
303 inline int32_t
304 floatToInt32(const float aFloatValue)
305 {
306  union FloatOrInt32 { float floatValue; int32_t int32Value; };
307  const FloatOrInt32* foi = reinterpret_cast<const FloatOrInt32*>(&aFloatValue);
308  return foi->int32Value;
309 }
310 
311 inline int64_t
312 doubleToInt64(const double aDoubleValue)
313 {
314  union DoubleOrInt64 { double doubleValue; int64_t int64Value; };
315  const DoubleOrInt64* dol = reinterpret_cast<const DoubleOrInt64*>(&aDoubleValue);
316  return dol->int64Value;
317 }
318 
319 
320 // aUnitsInLastPlace is the allowed difference between the least significant digits
321 // of the numbers' floating point representation
322 // Please read refernce paper before trying to use isUlpsEqual
323 // http://www.cygnus-software.com/papers/comparingFloats/comparingFloats.htm
324 inline bool
325 isUlpsEqual(const double aLeft, const double aRight, const int64_t aUnitsInLastPlace)
326 {
327  int64_t longLeft = doubleToInt64(aLeft);
328  // Because of 2's complement, must restore lexicographical order
329  if (longLeft < 0) {
330  longLeft = INT64_C(0x8000000000000000) - longLeft;
331  }
332 
333  int64_t longRight = doubleToInt64(aRight);
334  // Because of 2's complement, must restore lexicographical order
335  if (longRight < 0) {
336  longRight = INT64_C(0x8000000000000000) - longRight;
337  }
338 
339  int64_t difference = labs(longLeft - longRight);
340  return (difference <= aUnitsInLastPlace);
341 }
342 
343 inline bool
344 isUlpsEqual(const float aLeft, const float aRight, const int32_t aUnitsInLastPlace)
345 {
346  int32_t intLeft = floatToInt32(aLeft);
347  // Because of 2's complement, must restore lexicographical order
348  if (intLeft < 0) {
349  intLeft = 0x80000000 - intLeft;
350  }
351 
352  int32_t intRight = floatToInt32(aRight);
353  // Because of 2's complement, must restore lexicographical order
354  if (intRight < 0) {
355  intRight = 0x80000000 - intRight;
356  }
357 
358  int32_t difference = abs(intLeft - intRight);
359  return (difference <= aUnitsInLastPlace);
360 }
361 
362 // ==========> Pow <==================
363 
365 template<typename Type>
366 inline Type Pow2(Type x)
367 {
368  return x*x;
369 }
370 
372 template<typename Type>
373 inline Type Pow3(Type x)
374 {
375  return x*x*x;
376 }
377 
379 template<typename Type>
380 inline Type Pow4(Type x)
381 {
382  return Pow2(Pow2(x));
383 }
384 
386 template <typename Type>
387 Type Pow(Type x, int n)
388 {
389  Type ans = 1;
390  if (n < 0) {
391  n=-n;
392  x=1/x;
393  }
394  for (int i = 0; i < n; i++) ans *= x;
395  return ans;
396 }
397 
399 inline float Pow(float b, float e)
400 {
401  assert( b >= 0.0f && "Pow(float,float): base is negative" );
402  return powf(b,e);
403 }
404 
406 inline double Pow(double b, double e)
407 {
408  assert( b >= 0.0 && "Pow(double,double): base is negative" );
409  return pow(b,e);
410 }
411 
412 // ==========> Max <==================
413 
415 template< typename Type >
416 inline const Type& Max( const Type& a, const Type& b )
417 {
418  return std::max(a,b) ;
419 }
420 
422 template< typename Type >
423 inline const Type& Max( const Type& a, const Type& b, const Type& c )
424 {
425  return std::max( std::max(a,b), c ) ;
426 }
427 
429 template< typename Type >
430 inline const Type& Max( const Type& a, const Type& b, const Type& c, const Type& d )
431 {
432  return std::max( std::max(a,b), std::max(c,d) ) ;
433 }
434 
436 template< typename Type >
437 inline const Type& Max( const Type& a, const Type& b, const Type& c,
438  const Type& d, const Type& e )
439 {
440  return std::max( std::max(a,b), Max(c,d,e) ) ;
441 }
442 
444 template< typename Type >
445 inline const Type& Max( const Type& a, const Type& b, const Type& c,
446  const Type& d, const Type& e, const Type& f )
447 {
448  return std::max( Max(a,b,c), Max(d,e,f) ) ;
449 }
450 
452 template< typename Type >
453 inline const Type& Max( const Type& a, const Type& b, const Type& c, const Type& d,
454  const Type& e, const Type& f, const Type& g )
455 {
456  return std::max( Max(a,b,c,d), Max(e,f,g) ) ;
457 }
458 
460 template< typename Type >
461 inline const Type& Max( const Type& a, const Type& b, const Type& c, const Type& d,
462  const Type& e, const Type& f, const Type& g, const Type& h )
463 {
464  return std::max( Max(a,b,c,d), Max(e,f,g,h) ) ;
465 }
466 
467 // ==========> Min <==================
468 
470 template< typename Type >
471 inline const Type& Min( const Type& a, const Type& b )
472 {
473  return std::min(a,b) ;
474 }
475 
477 template< typename Type >
478 inline const Type& Min( const Type& a, const Type& b, const Type& c )
479 {
480  return std::min( std::min(a,b), c ) ;
481 }
482 
484 template< typename Type >
485 inline const Type& Min( const Type& a, const Type& b, const Type& c, const Type& d )
486 {
487  return std::min( std::min(a,b), std::min(c,d) ) ;
488 }
489 
491 template< typename Type >
492 inline const Type& Min( const Type& a, const Type& b, const Type& c,
493  const Type& d, const Type& e )
494 {
495  return std::min( std::min(a,b), Min(c,d,e) ) ;
496 }
497 
499 template< typename Type >
500 inline const Type& Min( const Type& a, const Type& b, const Type& c,
501  const Type& d, const Type& e, const Type& f )
502 {
503  return std::min( Min(a,b,c), Min(d,e,f) ) ;
504 }
505 
507 template< typename Type >
508 inline const Type& Min( const Type& a, const Type& b, const Type& c, const Type& d,
509  const Type& e, const Type& f, const Type& g )
510 {
511  return std::min( Min(a,b,c,d), Min(e,f,g) ) ;
512 }
513 
515 template< typename Type >
516 inline const Type& Min( const Type& a, const Type& b, const Type& c, const Type& d,
517  const Type& e, const Type& f, const Type& g, const Type& h )
518 {
519  return std::min( Min(a,b,c,d), Min(e,f,g,h) ) ;
520 }
521 
523 template <typename Type>
524 inline int Sign(const Type &x) {return ( (x)<0 ? -1 : (x)==0 ? 0 : 1);}
525 
526 
528 inline float Sqrt(float x) {return sqrtf(x);}
529 inline double Sqrt(double x){return sqrt(x);}
530 inline long double Sqrt(long double x) {return sqrtl(x);}
531 
532 
534 inline int Mod(int i, int j) {return (i%j);};
535 inline float Mod(float x, float y) {return fmodf(x,y);}
536 inline double Mod(double x, double y){return fmod(x,y);}
537 inline long double Mod(long double x, long double y) {return fmodl(x,y);}
538 
540 template <typename Type>
541 inline Type Reminder(Type x, Type y) {return Mod(x,y);}
542 
544 inline float RoundUp(float x) { return ceilf(x); }
545 inline double RoundUp(double x) { return ceil(x); }
546 inline long double RoundUp(long double x) { return ceill(x); }
547 template <typename Type>
548 inline Type RoundUp(Type x, Type base)
549 {
550  Type reminder=Reminder(x,base);
551  return reminder ? x-reminder+base : x;
552 }
553 
554 
556 inline float RoundDown(float x) { return floorf(x); }
557 inline double RoundDown(double x){ return floor(x); }
558 inline long double RoundDown(long double x) { return floorl(x); }
559 template <typename Type>
560 inline Type RoundDown(Type x, Type base)
561 {
562  Type reminder=Reminder(x,base);
563  if (!reminder)
564  return x;
565  else
566  return x-reminder;
567 }
568 
570 template <typename Type>
571 inline Type IntegerPart(Type x)
572 {
573  return (x > 0 ? RoundDown(x) : RoundUp(x));
574 }
575 
577 template <typename Type>
578 inline Type FractionalPart(Type x) { return Mod(x,Type(1)); }
579 
581 inline int Floor(float x) { return (int)RoundDown(x); }
582 inline int Floor(double x) { return (int)RoundDown(x); }
583 inline int Floor(long double x) { return (int)RoundDown(x); }
584 
586 inline int Ceil(float x) { return (int)RoundUp(x); }
587 inline int Ceil(double x) { return (int)RoundUp(x); }
588 inline int Ceil(long double x) { return (int)RoundUp(x); }
589 
591 template <typename Type>
592 inline Type Round(Type x) { return RoundDown(x+0.5); }
593 
595 template <typename Type>
596 inline Type Chop(Type x, Type delta) { return (Abs(x) < delta ? 0 : x); }
597 
599 template <typename Type>
600 inline Type Truncate(Type x,unsigned int digits)
601 {
602  Type tenth=Pow(10,digits);
603  return RoundDown(x*tenth+0.5)/tenth;
604 }
605 
607 template <typename Type>
608 inline Type Inv(Type x)
609 {
610  assert(x);
611  return Type(1)/x;
612 }
613 
614 
615 enum Axis {
616  X_AXIS = 0,
617  Y_AXIS = 1,
618  Z_AXIS = 2
619 };
620 
621 // enum values are consistent with their historical mx analogs.
631 };
632 
633 
634 template <typename S, typename T>
635 struct promote {
636  typedef typename boost::numeric::conversion_traits<S, T>::supertype type;
637 };
638 
639 
640 template<typename T> struct tolerance { static T value() { return 0; } };
641 template<> struct tolerance<float> { static float value() { return 1e-8f; } };
642 template<> struct tolerance<double> { static double value() { return 1e-15; } };
643 
644 } // namespace math
645 } // namespace OPENVDB_VERSION_NAME
646 } // namespace openvdb
647 
648 #endif // OPENVDB_MATH_MATH_HAS_BEEN_INCLUDED
649 
650 // Copyright (c) 2012-2013 DreamWorks Animation LLC
651 // All rights reserved. This software is distributed under the
652 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )