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