HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Math.h
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2012-2017 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 //
29 ///////////////////////////////////////////////////////////////////////////
30 //
31 /// @file Math.h
32 /// @brief General-purpose arithmetic and comparison routines, most of which
33 /// accept arbitrary value types (or at least arbitrary numeric value types)
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() etc
42 #include <cstdlib> // for srand(), abs(int)
43 #include <limits> // for std::numeric_limits<Type>::max()
44 #include <string>
45 #include <hboost/numeric/conversion/conversion_traits.hpp>
46 #include <hboost/math/special_functions/cbrt.hpp>
47 #include <hboost/math/special_functions/fpclassify.hpp> // hboost::math::isfinite
48 #include <hboost/random/mersenne_twister.hpp> // for hboost::random::mt19937
49 #include <hboost/random/uniform_01.hpp>
50 #include <hboost/random/uniform_int.hpp>
51 #include <hboost/version.hpp> // for HBOOST_VERSION
52 #include <openvdb/Platform.h>
53 #include <openvdb/version.h>
54 
55 
56 // Compile pragmas
57 
58 #define PRAGMA(x) _Pragma(#x)
59 
60 // Intel(r) compiler fires remark #1572: floating-point equality and inequality
61 // comparisons are unrealiable when == or != is used with floating point operands.
62 #if defined(__INTEL_COMPILER)
63  #define OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN \
64  _Pragma("warning (push)") \
65  _Pragma("warning (disable:1572)")
66  #define OPENVDB_NO_FP_EQUALITY_WARNING_END \
67  _Pragma("warning (pop)")
68 #elif defined(__clang__)
69  #define OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN \
70  PRAGMA(clang diagnostic push) \
71  PRAGMA(clang diagnostic ignored "-Wfloat-equal")
72  #define OPENVDB_NO_FP_EQUALITY_WARNING_END \
73  PRAGMA(clang diagnostic pop)
74 #else
75  // For GCC, #pragma GCC diagnostic ignored "-Wfloat-equal"
76  // isn't working until gcc 4.2+,
77  // Trying
78  // #pragma GCC system_header
79  // creates other problems, most notably "warning: will never be executed"
80  // in from templates, unsure of how to work around.
81  // If necessary, could use integer based comparisons for equality
82  #define OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN
83  #define OPENVDB_NO_FP_EQUALITY_WARNING_END
84 #endif
85 
86 namespace openvdb {
88 namespace OPENVDB_VERSION_NAME {
89 
90 /// @brief Return the value of type T that corresponds to zero.
91 /// @note A zeroVal<T>() specialization must be defined for each @c ValueType T
92 /// that cannot be constructed using the form @c T(0). For example, @c std::string(0)
93 /// treats 0 as @c NULL and throws a @c std::logic_error.
94 template<typename T> inline T zeroVal() { return T(0); }
95 /// Return the @c std::string value that corresponds to zero.
96 template<> inline std::string zeroVal<std::string>() { return ""; }
97 /// Return the @c bool value that corresponds to zero.
98 template<> inline bool zeroVal<bool>() { return false; }
99 
100 /// @todo These won't be needed if we eliminate StringGrids.
101 //@{
102 /// @brief Needed to support the <tt>(zeroVal<ValueType>() + val)</tt> idiom
103 /// when @c ValueType is @c std::string
104 inline std::string operator+(const std::string& s, bool) { return s; }
105 inline std::string operator+(const std::string& s, int) { return s; }
106 inline std::string operator+(const std::string& s, float) { return s; }
107 inline std::string operator+(const std::string& s, double) { return s; }
108 //@}
109 
110 
111 namespace math {
112 
113 /// @brief Return the unary negation of the given value.
114 /// @note A negative<T>() specialization must be defined for each ValueType T
115 /// for which unary negation is not defined.
116 template<typename T> inline T negative(const T& val) { return T(-val); }
117 /// Return the negation of the given boolean.
118 template<> inline bool negative(const bool& val) { return !val; }
119 /// Return the "negation" of the given string.
120 template<> inline std::string negative(const std::string& val) { return val; }
121 
122 
123 //@{
124 /// Tolerance for floating-point comparison
125 template<typename T> struct Tolerance { static T value() { return zeroVal<T>(); } };
126 template<> struct Tolerance<float> { static float value() { return 1e-8f; } };
127 template<> struct Tolerance<double> { static double value() { return 1e-15; } };
128 //@}
129 
130 //@{
131 /// Delta for small floating-point offsets
132 template<typename T> struct Delta { static T value() { return zeroVal<T>(); } };
133 template<> struct Delta<float> { static float value() { return 1e-5f; } };
134 template<> struct Delta<double> { static double value() { return 1e-9; } };
135 //@}
136 
137 
138 // ==========> Random Values <==================
139 
140 /// @brief Simple generator of random numbers over the range [0, 1)
141 /// @details Thread-safe as long as each thread has its own Rand01 instance
142 template<typename FloatType = double, typename EngineType = hboost::mt19937>
143 class Rand01
144 {
145 private:
146  EngineType mEngine;
147  hboost::uniform_01<FloatType> mRand;
148 
149 public:
150  typedef FloatType ValueType;
151 
152  /// @brief Initialize the generator.
153  /// @param engine random number generator
154  Rand01(const EngineType& engine): mEngine(engine) {}
155 
156  /// @brief Initialize the generator.
157  /// @param seed seed value for the random number generator
158  Rand01(unsigned int seed): mEngine(static_cast<typename EngineType::result_type>(seed)) {}
159 
160  /// Set the seed value for the random number generator
161  void setSeed(unsigned int seed)
162  {
163  mEngine.seed(static_cast<typename EngineType::result_type>(seed));
164  }
165 
166  /// Return a const reference to the random number generator.
167  const EngineType& engine() const { return mEngine; }
168 
169  /// Return a uniformly distributed random number in the range [0, 1).
170  FloatType operator()() { return mRand(mEngine); }
171 };
172 
174 
175 
176 /// @brief Simple random integer generator
177 /// @details Thread-safe as long as each thread has its own RandInt instance
178 template<typename IntType = int, typename EngineType = hboost::mt19937>
179 class RandInt
180 {
181 private:
182 #if HBOOST_VERSION >= 104700
183  typedef hboost::random::uniform_int_distribution<IntType> Distr;
184 #else
185  typedef hboost::uniform_int<IntType> Distr;
186 #endif
187  EngineType mEngine;
188  Distr mRand;
189 
190 public:
191  /// @brief Initialize the generator.
192  /// @param engine random number generator
193  /// @param imin,imax generate integers that are uniformly distributed over [imin, imax]
194  RandInt(const EngineType& engine, IntType imin, IntType imax):
195  mEngine(engine),
196  mRand(std::min(imin, imax), std::max(imin, imax))
197  {}
198 
199  /// @brief Initialize the generator.
200  /// @param seed seed value for the random number generator
201  /// @param imin,imax generate integers that are uniformly distributed over [imin, imax]
202  RandInt(unsigned int seed, IntType imin, IntType imax):
203  mEngine(static_cast<typename EngineType::result_type>(seed)),
204  mRand(std::min(imin, imax), std::max(imin, imax))
205  {}
206 
207  /// Change the range over which integers are distributed to [imin, imax].
208  void setRange(IntType imin, IntType imax)
209  {
210  mRand = Distr(std::min(imin, imax), std::max(imin, imax));
211  }
212 
213  /// Set the seed value for the random number generator
214  void setSeed(unsigned int seed)
215  {
216  mEngine.seed(static_cast<typename EngineType::result_type>(seed));
217  }
218 
219  /// Return a const reference to the random number generator.
220  const EngineType& engine() const { return mEngine; }
221 
222  /// Return a randomly-generated integer in the current range.
223  IntType operator()() { return mRand(mEngine); }
224 
225  /// @brief Return a randomly-generated integer in the new range [imin, imax],
226  /// without changing the current range.
227  IntType operator()(IntType imin, IntType imax)
228  {
229  const IntType lo = std::min(imin, imax), hi = std::max(imin, imax);
230 #if HBOOST_VERSION >= 104700
231  return mRand(mEngine, typename Distr::param_type(lo, hi));
232 #else
233  return Distr(lo, hi)(mEngine);
234 #endif
235  }
236 };
237 
239 
240 
241 // ==========> Clamp <==================
242 
243 /// Return @a x clamped to [@a min, @a max]
244 template<typename Type>
245 inline Type
247 {
248  assert( !(min>max) );
249  return x > min ? x < max ? x : max : min;
250 }
251 
252 
253 /// Return @a x clamped to [0, 1]
254 template<typename Type>
255 inline Type
256 Clamp01(Type x) { return x > Type(0) ? x < Type(1) ? x : Type(1) : Type(0); }
257 
258 
259 /// Return @c true if @a x is outside [0,1]
260 template<typename Type>
261 inline bool
263 {
264  if (x >= Type(0) && x <= Type(1)) return false;
265  x = x < Type(0) ? Type(0) : Type(1);
266  return true;
267 }
268 
269 /// @brief Return 0 if @a x < @a 0, 1 if @a x > 1 or else @f$(3-2x)x^2@f$.
270 template<typename Type>
271 inline Type
273 {
274  return x > 0 ? x < 1 ? (3-2*x)*x*x : Type(1) : Type(0);
275 }
276 
277 /// @brief Return 0 if @a x < @a min, 1 if @a x > @a max or else @f$(3-2t)t^2@f$,
278 /// where @f$t = (x-min)/(max-min)@f$.
279 template<typename Type>
280 inline Type
282 {
283  assert(min < max);
284  return SmoothUnitStep((x-min)/(max-min));
285 }
286 
287 
288 // ==========> Absolute Value <==================
289 
290 
291 //@{
292 /// Return the absolute value of the given quantity.
293 inline int32_t Abs(int32_t i) { return abs(i); }
294 inline int64_t Abs(int64_t i)
295 {
296 #ifdef _MSC_VER
297  return (i < int64_t(0) ? -i : i);
298 #else
299  return labs(i);
300 #endif
301 }
302 inline float Abs(float x) { return fabsf(x); }
303 inline double Abs(double x) { return fabs(x); }
304 inline long double Abs(long double x) { return fabsl(x); }
305 inline uint32_t Abs(uint32_t i) { return i; }
306 inline uint64_t Abs(uint64_t i) { return i; }
307 inline bool Abs(bool b) { return b; }
308 // On OSX size_t and uint64_t are different types
309 #if defined(__APPLE__) || defined(MACOSX)
310 inline size_t Abs(size_t i) { return i; }
311 #endif
312 //@}
313 
314 
315 ////////////////////////////////////////
316 
317 
318 // ==========> Value Comparison <==================
319 
320 
321 /// Return @c true if @a x is exactly equal to zero.
322 template<typename Type>
323 inline bool
324 isZero(const Type& x)
325 {
327  return x == zeroVal<Type>();
329 }
330 
331 
332 /// @brief Return @c true if @a x is equal to zero to within
333 /// the default floating-point comparison tolerance.
334 template<typename Type>
335 inline bool
337 {
338  const Type tolerance = Type(zeroVal<Type>() + Tolerance<Type>::value());
339  return !(x > tolerance) && !(x < -tolerance);
340 }
341 
342 /// Return @c true if @a x is equal to zero to within the given tolerance.
343 template<typename Type>
344 inline bool
345 isApproxZero(const Type& x, const Type& tolerance)
346 {
347  return !(x > tolerance) && !(x < -tolerance);
348 }
349 
350 
351 /// Return @c true if @a x is less than zero.
352 template<typename Type>
353 inline bool
354 isNegative(const Type& x) { return x < zeroVal<Type>(); }
355 
356 /// Return @c false, since @c bool values are never less than zero.
357 template<> inline bool isNegative<bool>(const bool&) { return false; }
358 
359 
360 /// Return @c true if @a x is finite.
361 template<typename Type>
362 inline bool
363 isFinite(const Type& x) { return hboost::math::isfinite(x); }
364 
365 
366 /// @brief Return @c true if @a a is equal to @a b to within
367 /// the default floating-point comparison tolerance.
368 template<typename Type>
369 inline bool
370 isApproxEqual(const Type& a, const Type& b)
371 {
372  const Type tolerance = Type(zeroVal<Type>() + Tolerance<Type>::value());
373  return !(Abs(a - b) > tolerance);
374 }
375 
376 
377 /// Return @c true if @a a is equal to @a b to within the given tolerance.
378 template<typename Type>
379 inline bool
380 isApproxEqual(const Type& a, const Type& b, const Type& tolerance)
381 {
382  return !(Abs(a - b) > tolerance);
383 }
384 
385 #define OPENVDB_EXACT_IS_APPROX_EQUAL(T) \
386  template<> inline bool isApproxEqual<T>(const T& a, const T& b) { return a == b; } \
387  template<> inline bool isApproxEqual<T>(const T& a, const T& b, const T&) { return a == b; } \
388  /**/
389 
392 
393 
394 /// @brief Return @c true if @a a is larger than @a b to within
395 /// the given tolerance, i.e., if @a b - @a a < @a tolerance.
396 template<typename Type>
397 inline bool
399 {
400  return (b - a < tolerance);
401 }
402 
403 
404 /// @brief Return @c true if @a a is exactly equal to @a b.
405 template<typename T0, typename T1>
406 inline bool
407 isExactlyEqual(const T0& a, const T1& b)
408 {
410  return a == b;
412 }
413 
414 
415 template<typename Type>
416 inline bool
417 isRelOrApproxEqual(const Type& a, const Type& b, const Type& absTol, const Type& relTol)
418 {
419  // First check to see if we are inside the absolute tolerance
420  // Necessary for numbers close to 0
421  if (!(Abs(a - b) > absTol)) return true;
422 
423  // Next check to see if we are inside the relative tolerance
424  // to handle large numbers that aren't within the abs tolerance
425  // but could be the closest floating point representation
426  double relError;
427  if (Abs(b) > Abs(a)) {
428  relError = Abs((a - b) / b);
429  } else {
430  relError = Abs((a - b) / a);
431  }
432  return (relError <= relTol);
433 }
434 
435 template<>
436 inline bool
437 isRelOrApproxEqual(const bool& a, const bool& b, const bool&, const bool&)
438 {
439  return (a == b);
440 }
441 
442 
443 // Avoid strict aliasing issues by using type punning
444 // http://cellperformance.beyond3d.com/articles/2006/06/understanding-strict-aliasing.html
445 // Using "casting through a union(2)"
446 inline int32_t
447 floatToInt32(const float aFloatValue)
448 {
449  union FloatOrInt32 { float floatValue; int32_t int32Value; };
450  const FloatOrInt32* foi = reinterpret_cast<const FloatOrInt32*>(&aFloatValue);
451  return foi->int32Value;
452 }
453 
454 
455 inline int64_t
456 doubleToInt64(const double aDoubleValue)
457 {
458  union DoubleOrInt64 { double doubleValue; int64_t int64Value; };
459  const DoubleOrInt64* dol = reinterpret_cast<const DoubleOrInt64*>(&aDoubleValue);
460  return dol->int64Value;
461 }
462 
463 
464 // aUnitsInLastPlace is the allowed difference between the least significant digits
465 // of the numbers' floating point representation
466 // Please read the reference paper before trying to use isUlpsEqual
467 // http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
468 inline bool
469 isUlpsEqual(const double aLeft, const double aRight, const int64_t aUnitsInLastPlace)
470 {
471  int64_t longLeft = doubleToInt64(aLeft);
472  // Because of 2's complement, must restore lexicographical order
473  if (longLeft < 0) {
474  longLeft = INT64_C(0x8000000000000000) - longLeft;
475  }
476 
477  int64_t longRight = doubleToInt64(aRight);
478  // Because of 2's complement, must restore lexicographical order
479  if (longRight < 0) {
480  longRight = INT64_C(0x8000000000000000) - longRight;
481  }
482 
483  int64_t difference = labs(longLeft - longRight);
484  return (difference <= aUnitsInLastPlace);
485 }
486 
487 inline bool
488 isUlpsEqual(const float aLeft, const float aRight, const int32_t aUnitsInLastPlace)
489 {
490  int32_t intLeft = floatToInt32(aLeft);
491  // Because of 2's complement, must restore lexicographical order
492  if (intLeft < 0) {
493  intLeft = 0x80000000 - intLeft;
494  }
495 
496  int32_t intRight = floatToInt32(aRight);
497  // Because of 2's complement, must restore lexicographical order
498  if (intRight < 0) {
499  intRight = 0x80000000 - intRight;
500  }
501 
502  int32_t difference = abs(intLeft - intRight);
503  return (difference <= aUnitsInLastPlace);
504 }
505 
506 
507 ////////////////////////////////////////
508 
509 
510 // ==========> Pow <==================
511 
512 /// Return @f$ x^2 @f$.
513 template<typename Type>
514 inline Type Pow2(Type x) { return x*x; }
515 
516 /// Return @f$ x^3 @f$.
517 template<typename Type>
518 inline Type Pow3(Type x) { return x*x*x; }
519 
520 /// Return @f$ x^4 @f$.
521 template<typename Type>
522 inline Type Pow4(Type x) { return Pow2(Pow2(x)); }
523 
524 /// Return @f$ x^n @f$.
525 template<typename Type>
526 Type
527 Pow(Type x, int n)
528 {
529  Type ans = 1;
530  if (n < 0) {
531  n = -n;
532  x = Type(1)/x;
533  }
534  while (n--) ans *= x;
535  return ans;
536 }
537 
538 //@{
539 /// Return @f$ b^e @f$.
540 inline float
541 Pow(float b, float e)
542 {
543  assert( b >= 0.0f && "Pow(float,float): base is negative" );
544  return powf(b,e);
545 }
546 
547 inline double
548 Pow(double b, double e)
549 {
550  assert( b >= 0.0 && "Pow(double,double): base is negative" );
551  return pow(b,e);
552 }
553 //@}
554 
555 
556 // ==========> Max <==================
557 
558 /// Return the maximum of two values
559 template<typename Type>
560 inline const Type&
561 Max(const Type& a, const Type& b)
562 {
563  return std::max(a,b) ;
564 }
565 
566 /// Return the maximum of three values
567 template<typename Type>
568 inline const Type&
569 Max(const Type& a, const Type& b, const Type& c)
570 {
571  return std::max( std::max(a,b), c ) ;
572 }
573 
574 /// Return the maximum of four values
575 template<typename Type>
576 inline const Type&
577 Max(const Type& a, const Type& b, const Type& c, const Type& d)
578 {
579  return std::max(std::max(a,b), std::max(c,d));
580 }
581 
582 /// Return the maximum of five values
583 template<typename Type>
584 inline const Type&
585 Max(const Type& a, const Type& b, const Type& c, const Type& d, const Type& e)
586 {
587  return std::max(std::max(a,b), Max(c,d,e));
588 }
589 
590 /// Return the maximum of six values
591 template<typename Type>
592 inline const Type&
593 Max(const Type& a, const Type& b, const Type& c, const Type& d, const Type& e, const Type& f)
594 {
595  return std::max(Max(a,b,c), Max(d,e,f));
596 }
597 
598 /// Return the maximum of seven values
599 template<typename Type>
600 inline const Type&
601 Max(const Type& a, const Type& b, const Type& c, const Type& d,
602  const Type& e, const Type& f, const Type& g)
603 {
604  return std::max(Max(a,b,c,d), Max(e,f,g));
605 }
606 
607 /// Return the maximum of eight values
608 template<typename Type>
609 inline const Type&
610 Max(const Type& a, const Type& b, const Type& c, const Type& d,
611  const Type& e, const Type& f, const Type& g, const Type& h)
612 {
613  return std::max(Max(a,b,c,d), Max(e,f,g,h));
614 }
615 
616 
617 // ==========> Min <==================
618 
619 /// Return the minimum of two values
620 template<typename Type>
621 inline const Type&
622 Min(const Type& a, const Type& b) { return std::min(a, b); }
623 
624 /// Return the minimum of three values
625 template<typename Type>
626 inline const Type&
627 Min(const Type& a, const Type& b, const Type& c) { return std::min(std::min(a, b), c); }
628 
629 /// Return the minimum of four values
630 template<typename Type>
631 inline const Type&
632 Min(const Type& a, const Type& b, const Type& c, const Type& d)
633 {
634  return std::min(std::min(a, b), std::min(c, d));
635 }
636 
637 /// Return the minimum of five values
638 template<typename Type>
639 inline const Type&
640 Min(const Type& a, const Type& b, const Type& c, const Type& d, const Type& e)
641 {
642  return std::min(std::min(a,b), Min(c,d,e));
643 }
644 
645 /// Return the minimum of six values
646 template<typename Type>
647 inline const Type&
648 Min(const Type& a, const Type& b, const Type& c, const Type& d, const Type& e, const Type& f)
649 {
650  return std::min(Min(a,b,c), Min(d,e,f));
651 }
652 
653 /// Return the minimum of seven values
654 template<typename Type>
655 inline const Type&
656 Min(const Type& a, const Type& b, const Type& c, const Type& d,
657  const Type& e, const Type& f, const Type& g)
658 {
659  return std::min(Min(a,b,c,d), Min(e,f,g));
660 }
661 
662 /// Return the minimum of eight values
663 template<typename Type>
664 inline const Type&
665 Min(const Type& a, const Type& b, const Type& c, const Type& d,
666  const Type& e, const Type& f, const Type& g, const Type& h)
667 {
668  return std::min(Min(a,b,c,d), Min(e,f,g,h));
669 }
670 
671 
672 // ============> Exp <==================
673 
674 /// Return @f$ e^x @f$.
675 template<typename Type>
676 inline Type Exp(const Type& x) { return std::exp(x); }
677 
678 // ============> Sin <==================
679 
680 //@{
681 /// Return @f$ sin(x) @f$.
682 inline float Sin(const float& x) { return sinf(x); }
683 
684 inline double Sin(const double& x) { return sin(x); }
685 //@}
686 
687 // ============> Cos <==================
688 
689 //@{
690 /// Return @f$ cos(x) @f$.
691 inline float Cos(const float& x) { return cosf(x); }
692 
693 inline double Cos(const double& x) { return cos(x); }
694 //@}
695 
696 
697 ////////////////////////////////////////
698 
699 
700 /// Return the sign of the given value as an integer (either -1, 0 or 1).
701 template <typename Type>
702 inline int Sign(const Type &x) { return (zeroVal<Type>() < x) - (x < zeroVal<Type>()); }
703 
704 
705 /// @brief Return @c true if @a a and @a b have different signs.
706 /// @note Zero is considered a positive number.
707 template <typename Type>
708 inline bool
709 SignChange(const Type& a, const Type& b)
710 {
711  return ( (a<zeroVal<Type>()) ^ (b<zeroVal<Type>()) );
712 }
713 
714 
715 /// @brief Return @c true if the interval [@a a, @a b] includes zero,
716 /// i.e., if either @a a or @a b is zero or if they have different signs.
717 template <typename Type>
718 inline bool
719 ZeroCrossing(const Type& a, const Type& b)
720 {
721  return a * b <= zeroVal<Type>();
722 }
723 
724 
725 //@{
726 /// Return the square root of a floating-point value.
727 inline float Sqrt(float x) { return sqrtf(x); }
728 inline double Sqrt(double x) { return sqrt(x); }
729 inline long double Sqrt(long double x) { return sqrtl(x); }
730 //@}
731 
732 
733 //@{
734 /// Return the cube root of a floating-point value.
735 inline float Cbrt(float x) { return hboost::math::cbrt(x); }
736 inline double Cbrt(double x) { return hboost::math::cbrt(x); }
737 inline long double Cbrt(long double x) { return hboost::math::cbrt(x); }
738 //@}
739 
740 
741 //@{
742 /// Return the remainder of @a x / @a y.
743 inline int Mod(int x, int y) { return (x % y); }
744 inline float Mod(float x, float y) { return fmodf(x,y); }
745 inline double Mod(double x, double y) { return fmod(x,y); }
746 inline long double Mod(long double x, long double y) { return fmodl(x,y); }
747 template<typename Type> inline Type Remainder(Type x, Type y) { return Mod(x,y); }
748 //@}
749 
750 
751 //@{
752 /// Return @a x rounded up to the nearest integer.
753 inline float RoundUp(float x) { return ceilf(x); }
754 inline double RoundUp(double x) { return ceil(x); }
755 inline long double RoundUp(long double x) { return ceill(x); }
756 //@}
757 /// Return @a x rounded up to the nearest multiple of @a base.
758 template<typename Type>
759 inline Type
761 {
762  Type remainder = Remainder(x, base);
763  return remainder ? x-remainder+base : x;
764 }
765 
766 
767 //@{
768 /// Return @a x rounded down to the nearest integer.
769 inline float RoundDown(float x) { return floorf(x); }
770 inline double RoundDown(double x) { return floor(x); }
771 inline long double RoundDown(long double x) { return floorl(x); }
772 //@}
773 /// Return @a x rounded down to the nearest multiple of @a base.
774 template<typename Type>
775 inline Type
777 {
778  Type remainder = Remainder(x, base);
779  return remainder ? x-remainder : x;
780 }
781 
782 
783 //@{
784 /// Return @a x rounded to the nearest integer.
785 inline float Round(float x) { return RoundDown(x + 0.5f); }
786 inline double Round(double x) { return RoundDown(x + 0.5); }
787 inline long double Round(long double x) { return RoundDown(x + 0.5l); }
788 //@}
789 
790 
791 /// Return the euclidean remainder of @a x.
792 /// Note unlike % operator this will always return a positive result
793 template<typename Type>
794 inline Type
795 EuclideanRemainder(Type x) { return x - RoundDown(x); }
796 
797 
798 /// Return the integer part of @a x.
799 template<typename Type>
800 inline Type
802 {
803  return (x > 0 ? RoundDown(x) : RoundUp(x));
804 }
805 
806 /// Return the fractional part of @a x.
807 template<typename Type>
808 inline Type
809 FractionalPart(Type x) { return Mod(x,Type(1)); }
810 
811 
812 //@{
813 /// Return the floor of @a x.
814 inline int Floor(float x) { return int(RoundDown(x)); }
815 inline int Floor(double x) { return int(RoundDown(x)); }
816 inline int Floor(long double x) { return int(RoundDown(x)); }
817 //@}
818 
819 
820 //@{
821 /// Return the ceiling of @a x.
822 inline int Ceil(float x) { return int(RoundUp(x)); }
823 inline int Ceil(double x) { return int(RoundUp(x)); }
824 inline int Ceil(long double x) { return int(RoundUp(x)); }
825 //@}
826 
827 
828 /// Return @a x if it is greater or equal in magnitude than @a delta. Otherwise, return zero.
829 template<typename Type>
830 inline Type Chop(Type x, Type delta) { return (Abs(x) < delta ? zeroVal<Type>() : x); }
831 
832 
833 /// Return @a x truncated to the given number of decimal digits.
834 template<typename Type>
835 inline Type
836 Truncate(Type x, unsigned int digits)
837 {
838  Type tenth = Pow(10,digits);
839  return RoundDown(x*tenth+0.5)/tenth;
840 }
841 
842 
843 ////////////////////////////////////////
844 
845 
846 /// Return the inverse of @a x.
847 template<typename Type>
848 inline Type
850 {
851  assert(x);
852  return Type(1)/x;
853 }
854 
855 
856 enum Axis {
857  X_AXIS = 0,
858  Y_AXIS = 1,
859  Z_AXIS = 2
860 };
861 
862 // enum values are consistent with their historical mx analogs.
872 };
873 
874 
875 template <typename S, typename T>
876 struct promote {
877  typedef typename hboost::numeric::conversion_traits<S, T>::supertype type;
878 };
879 
880 
881 /// @brief Return the index [0,1,2] of the smallest value in a 3D vector.
882 /// @note This methods assumes operator[] exists and avoids branching.
883 /// @details If two components of the input vector are equal and smaller than the
884 /// third component, the largest index of the two is always returned.
885 /// If all three vector components are equal the largest index, i.e. 2, is
886 /// returned. In other words the return value corresponds to the largest index
887 /// of the of the smallest vector components.
888 template<typename Vec3T>
889 size_t
890 MinIndex(const Vec3T& v)
891 {
892 #ifndef _MSC_VER // Visual C++ doesn't guarantee thread-safe initialization of local statics
893  static
894 #endif
895  const size_t hashTable[8] = { 2, 1, 9, 1, 2, 9, 0, 0 };//9 is a dummy value
896  const size_t hashKey =
897  ((v[0] < v[1]) << 2) + ((v[0] < v[2]) << 1) + (v[1] < v[2]);// ?*4+?*2+?*1
898  return hashTable[hashKey];
899 }
900 
901 
902 /// @brief Return the index [0,1,2] of the largest value in a 3D vector.
903 /// @note This methods assumes operator[] exists and avoids branching.
904 /// @details If two components of the input vector are equal and larger than the
905 /// third component, the largest index of the two is always returned.
906 /// If all three vector components are equal the largest index, i.e. 2, is
907 /// returned. In other words the return value corresponds to the largest index
908 /// of the largest vector components.
909 template<typename Vec3T>
910 size_t
911 MaxIndex(const Vec3T& v)
912 {
913 #ifndef _MSC_VER // Visual C++ doesn't guarantee thread-safe initialization of local statics
914  static
915 #endif
916  const size_t hashTable[8] = { 2, 1, 9, 1, 2, 9, 0, 0 };//9 is a dummy value
917  const size_t hashKey =
918  ((v[0] > v[1]) << 2) + ((v[0] > v[2]) << 1) + (v[1] > v[2]);// ?*4+?*2+?*1
919  return hashTable[hashKey];
920 }
921 
922 } // namespace math
923 } // namespace OPENVDB_VERSION_NAME
924 } // namespace openvdb
925 
926 #endif // OPENVDB_MATH_MATH_HAS_BEEN_INCLUDED
927 
928 // Copyright (c) 2012-2017 DreamWorks Animation LLC
929 // All rights reserved. This software is distributed under the
930 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
void setSeed(unsigned int seed)
Set the seed value for the random number generator.
Definition: Math.h:161
int32_t floatToInt32(const float aFloatValue)
Definition: Math.h:447
SYS_API double cos(double x)
SYS_API double fmod(double x, double y)
Type Inv(Type x)
Return the inverse of x.
Definition: Math.h:849
Type Truncate(Type x, unsigned int digits)
Return x truncated to the given number of decimal digits.
Definition: Math.h:836
int Ceil(float x)
Return the ceiling of x.
Definition: Math.h:822
Simple generator of random numbers over the range [0, 1)
Definition: Math.h:143
Type Pow(Type x, int n)
Return .
Definition: Math.h:527
Type IntegerPart(Type x)
Return the integer part of x.
Definition: Math.h:801
bool ZeroCrossing(const Type &a, const Type &b)
Return true if the interval [a, b] includes zero, i.e., if either a or b is zero or if they have diff...
Definition: Math.h:719
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:407
Type Pow2(Type x)
Return .
Definition: Math.h:514
const hboost::disable_if_c< VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:128
T negative(const T &val)
Return the unary negation of the given value.
Definition: Math.h:116
size_t MaxIndex(const Vec3T &v)
Return the index [0,1,2] of the largest value in a 3D vector.
Definition: Math.h:911
const GLdouble * v
Definition: glcorearb.h:836
GLsizei const GLchar *const * string
Definition: glcorearb.h:813
Type FractionalPart(Type x)
Return the fractional part of x.
Definition: Math.h:809
GLboolean GLboolean g
Definition: glcorearb.h:1221
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1221
const EngineType & engine() const
Return a const reference to the random number generator.
Definition: Math.h:220
float Cos(const float &x)
Return .
Definition: Math.h:691
SYS_API float powf(float x, float y)
GLint y
Definition: glcorearb.h:102
Rand01(unsigned int seed)
Initialize the generator.
Definition: Math.h:158
bool SignChange(const Type &a, const Type &b)
Return true if a and b have different signs.
Definition: Math.h:709
float RoundDown(float x)
Return x rounded down to the nearest integer.
Definition: Math.h:769
RandInt(const EngineType &engine, IntType imin, IntType imax)
Initialize the generator.
Definition: Math.h:194
SYS_API double pow(double x, double y)
bool zeroVal< bool >()
Return the bool value that corresponds to zero.
Definition: Math.h:98
bool isNegative(const Type &x)
Return true if x is less than zero.
Definition: Math.h:354
int64_t doubleToInt64(const double aDoubleValue)
Definition: Math.h:456
png_uint_32 i
Definition: png.h:2877
Tolerance for floating-point comparison.
Definition: Math.h:125
float RoundUp(float x)
Return x rounded up to the nearest integer.
Definition: Math.h:753
float Cbrt(float x)
Return the cube root of a floating-point value.
Definition: Math.h:735
Simple random integer generator.
Definition: Math.h:179
SYS_API float sinf(float x)
const hboost::disable_if_c< VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:132
bool isNegative< bool >(const bool &)
Return false, since bool values are never less than zero.
Definition: Math.h:357
GLdouble n
Definition: glcorearb.h:2007
GLfloat f
Definition: glcorearb.h:1925
size_t MinIndex(const Vec3T &v)
Return the index [0,1,2] of the smallest value in a 3D vector.
Definition: Math.h:890
IMATH_INTERNAL_NAMESPACE_HEADER_ENTER T abs(T a)
Definition: ImathFun.h:55
Coord Abs(const Coord &xyz)
Definition: Coord.h:254
bool isFinite(const Type &x)
Return true if x is finite.
Definition: Math.h:363
#define OPENVDB_VERSION_NAME
Definition: version.h:43
float Round(float x)
Return x rounded to the nearest integer.
Definition: Math.h:785
Type Pow4(Type x)
Return .
Definition: Math.h:522
Delta for small floating-point offsets.
Definition: Math.h:132
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:727
hboost::numeric::conversion_traits< S, T >::supertype type
Definition: Math.h:877
SYS_API double exp(double x)
Type Pow3(Type x)
Return .
Definition: Math.h:518
bool isRelOrApproxEqual(const Type &a, const Type &b, const Type &absTol, const Type &relTol)
Definition: Math.h:417
#define OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN
Definition: Math.h:82
bool isUlpsEqual(const double aLeft, const double aRight, const int64_t aUnitsInLastPlace)
Definition: Math.h:469
bool isApproxZero(const Type &x)
Return true if x is equal to zero to within the default floating-point comparison tolerance...
Definition: Math.h:336
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1221
bool isApproxEqual(const Type &a, const Type &b)
Return true if a is equal to b to within the default floating-point comparison tolerance.
Definition: Math.h:370
int floor(T x)
Definition: ImathFun.h:150
Type Clamp(Type x, Type min, Type max)
Return x clamped to [min, max].
Definition: Math.h:246
std::string operator+(const std::string &s, bool)
Needed to support the (zeroVal<ValueType>() + val) idiom when ValueType is std::string.
Definition: Math.h:104
#define OPENVDB_EXACT_IS_APPROX_EQUAL(T)
Definition: Math.h:385
RandInt(unsigned int seed, IntType imin, IntType imax)
Initialize the generator.
Definition: Math.h:202
GLfloat GLfloat GLfloat GLfloat h
Definition: glcorearb.h:2001
SYS_API float fmodf(float x, float y)
int Sign(const Type &x)
Return the sign of the given value as an integer (either -1, 0 or 1).
Definition: Math.h:702
typedef int
Definition: png.h:1175
void setSeed(unsigned int seed)
Set the seed value for the random number generator.
Definition: Math.h:214
GLint GLenum GLint x
Definition: glcorearb.h:408
GLuint GLfloat * val
Definition: glcorearb.h:1607
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:622
const EngineType & engine() const
Return a const reference to the random number generator.
Definition: Math.h:167
Type Chop(Type x, Type delta)
Return x if it is greater or equal in magnitude than delta. Otherwise, return zero.
Definition: Math.h:830
int Floor(float x)
Return the floor of x.
Definition: Math.h:814
bool ClampTest01(Type &x)
Return true if x is outside [0,1].
Definition: Math.h:262
IntType operator()()
Return a randomly-generated integer in the current range.
Definition: Math.h:223
float Sin(const float &x)
Return .
Definition: Math.h:682
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
#define const
Definition: zconf.h:214
int ceil(T x)
Definition: ImathFun.h:158
T zeroVal()
Return the value of type T that corresponds to zero.
Definition: Math.h:94
Rand01(const EngineType &engine)
Initialize the generator.
Definition: Math.h:154
Type SmoothUnitStep(Type x)
Return 0 if x < 0, 1 if x > 1 or else .
Definition: Math.h:272
void setRange(IntType imin, IntType imax)
Change the range over which integers are distributed to [imin, imax].
Definition: Math.h:208
Type Remainder(Type x, Type y)
Return the remainder of x / y.
Definition: Math.h:747
SYS_API float cosf(float x)
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:561
SYS_API double sin(double x)
int Mod(int x, int y)
Return the remainder of x / y.
Definition: Math.h:743
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:324
Type Exp(const Type &x)
Return .
Definition: Math.h:676
IntType operator()(IntType imin, IntType imax)
Return a randomly-generated integer in the new range [imin, imax], without changing the current range...
Definition: Math.h:227
Rand01< double, hboost::mt19937 > Random01
Definition: Math.h:173
FloatType operator()()
Return a uniformly distributed random number in the range [0, 1).
Definition: Math.h:170
RandInt< int, hboost::mt19937 > RandomInt
Definition: Math.h:238
Type Clamp01(Type x)
Return x clamped to [0, 1].
Definition: Math.h:256
#define OPENVDB_NO_FP_EQUALITY_WARNING_END
Definition: Math.h:83
bool isApproxLarger(const Type &a, const Type &b, const Type &tolerance)
Return true if a is larger than b to within the given tolerance, i.e., if b - a < tolerance...
Definition: Math.h:398