OpenVDB  10.0.0
Math.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
4 /// @file Math.h
5 /// @brief General-purpose arithmetic and comparison routines, most of which
6 /// accept arbitrary value types (or at least arbitrary numeric value types)
7 
8 #ifndef OPENVDB_MATH_HAS_BEEN_INCLUDED
9 #define OPENVDB_MATH_HAS_BEEN_INCLUDED
10 
11 #include <openvdb/Platform.h>
12 #include <openvdb/version.h>
13 #include <boost/numeric/conversion/conversion_traits.hpp>
14 #include <algorithm> // for std::max()
15 #include <cassert>
16 #include <cmath> // for std::ceil(), std::fabs(), std::pow(), std::sqrt(), etc.
17 #include <cstdlib> // for abs(int)
18 #include <cstring> // for memcpy
19 #include <random>
20 #include <string>
21 #include <type_traits> // for std::is_arithmetic
22 
23 
24 // Compile pragmas
25 
26 // Intel(r) compiler fires remark #1572: floating-point equality and inequality
27 // comparisons are unrealiable when == or != is used with floating point operands.
28 #if defined(__INTEL_COMPILER)
29  #define OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN \
30  _Pragma("warning (push)") \
31  _Pragma("warning (disable:1572)")
32  #define OPENVDB_NO_FP_EQUALITY_WARNING_END \
33  _Pragma("warning (pop)")
34 #elif defined(__clang__)
35  #define OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN \
36  PRAGMA(clang diagnostic push) \
37  PRAGMA(clang diagnostic ignored "-Wfloat-equal")
38  #define OPENVDB_NO_FP_EQUALITY_WARNING_END \
39  PRAGMA(clang diagnostic pop)
40 #else
41  // For GCC, #pragma GCC diagnostic ignored "-Wfloat-equal"
42  // isn't working until gcc 4.2+,
43  // Trying
44  // #pragma GCC system_header
45  // creates other problems, most notably "warning: will never be executed"
46  // in from templates, unsure of how to work around.
47  // If necessary, could use integer based comparisons for equality
48  #define OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN
49  #define OPENVDB_NO_FP_EQUALITY_WARNING_END
50 #endif
51 
52 
53 #ifdef OPENVDB_IS_POD
54 #undef OPENVDB_IS_POD
55 #endif
56 #define OPENVDB_IS_POD(Type) \
57 static_assert(std::is_standard_layout<Type>::value, \
58  #Type" must be a POD type (satisfy StandardLayoutType.)"); \
59 static_assert(std::is_trivial<Type>::value, \
60  #Type" must be a POD type (satisfy TrivialType.)");
61 
62 namespace openvdb {
64 namespace OPENVDB_VERSION_NAME {
65 
66 /// @brief Return the value of type T that corresponds to zero.
67 /// @note A zeroVal<T>() specialization must be defined for each @c ValueType T
68 /// that cannot be constructed using the form @c T(0). For example, @c std::string(0)
69 /// treats 0 as @c nullptr and throws a @c std::logic_error.
70 template<typename T> inline T zeroVal() { return T(0); }
71 /// Return the @c std::string value that corresponds to zero.
72 template<> inline std::string zeroVal<std::string>() { return ""; }
73 /// Return the @c bool value that corresponds to zero.
74 template<> inline bool zeroVal<bool>() { return false; }
75 
76 namespace math {
77 
78 /// @todo These won't be needed if we eliminate StringGrids.
79 //@{
80 /// @brief Needed to support the <tt>(zeroVal<ValueType>() + val)</tt> idiom
81 /// when @c ValueType is @c std::string
82 inline std::string operator+(const std::string& s, bool) { return s; }
83 inline std::string operator+(const std::string& s, int) { return s; }
84 inline std::string operator+(const std::string& s, float) { return s; }
85 inline std::string operator+(const std::string& s, double) { return s; }
86 //@}
87 
88 /// @brief Componentwise adder for POD types.
89 template<typename Type1, typename Type2>
90 inline auto cwiseAdd(const Type1& v, const Type2 s)
91 {
93  return v + s;
95 }
96 
97 /// @brief Componentwise less than for POD types.
98 template<typename Type1, typename Type2>
99 inline bool cwiseLessThan(const Type1& a, const Type2& b)
100 {
102  return a < b;
104 }
105 
106 /// @brief Componentwise greater than for POD types.
107 template<typename Type1, typename Type2>
108 inline bool cwiseGreaterThan(const Type1& a, const Type2& b)
109 {
111  return a > b;
113 }
114 
115 
116 
117 /// @brief Pi constant taken from Boost to match old behaviour
118 /// @note Available in C++20
119 template <typename T> inline constexpr T pi() { return 3.141592653589793238462643383279502884e+00; }
120 template <> inline constexpr float pi() { return 3.141592653589793238462643383279502884e+00F; }
121 template <> inline constexpr double pi() { return 3.141592653589793238462643383279502884e+00; }
122 template <> inline constexpr long double pi() { return 3.141592653589793238462643383279502884e+00L; }
123 
124 
125 /// @brief Return the unary negation of the given value.
126 /// @note A negative<T>() specialization must be defined for each ValueType T
127 /// for which unary negation is not defined.
128 template<typename T> inline T negative(const T& val)
129 {
130 // disable unary minus on unsigned warning
131 #if defined(_MSC_VER)
132 #pragma warning(push)
133 #pragma warning(disable:4146)
134 #endif
135  return T(-val);
136 #if defined(_MSC_VER)
137 #pragma warning(pop)
138 #endif
139 }
140 /// Return the negation of the given boolean.
141 template<> inline bool negative(const bool& val) { return !val; }
142 /// Return the "negation" of the given string.
143 template<> inline std::string negative(const std::string& val) { return val; }
144 
145 
146 //@{
147 /// Tolerance for floating-point comparison
148 template<typename T> struct Tolerance { static T value() { return zeroVal<T>(); } };
149 template<> struct Tolerance<float> { static float value() { return 1e-8f; } };
150 template<> struct Tolerance<double> { static double value() { return 1e-15; } };
151 //@}
152 
153 //@{
154 /// Delta for small floating-point offsets
155 template<typename T> struct Delta { static T value() { return zeroVal<T>(); } };
156 template<> struct Delta<float> { static float value() { return 1e-5f; } };
157 template<> struct Delta<double> { static double value() { return 1e-9; } };
158 //@}
159 
160 
161 // ==========> Random Values <==================
162 
163 /// @brief Simple generator of random numbers over the range [0, 1)
164 /// @details Thread-safe as long as each thread has its own Rand01 instance
165 template<typename FloatType = double, typename EngineType = std::mt19937>
166 class Rand01
167 {
168 private:
169  EngineType mEngine;
170  std::uniform_real_distribution<FloatType> mRand;
171 
172 public:
173  using ValueType = FloatType;
174 
175  /// @brief Initialize the generator.
176  /// @param engine random number generator
177  Rand01(const EngineType& engine): mEngine(engine) {}
178 
179  /// @brief Initialize the generator.
180  /// @param seed seed value for the random number generator
181  Rand01(unsigned int seed): mEngine(static_cast<typename EngineType::result_type>(seed)) {}
182 
183  /// Set the seed value for the random number generator
184  void setSeed(unsigned int seed)
185  {
186  mEngine.seed(static_cast<typename EngineType::result_type>(seed));
187  }
188 
189  /// Return a const reference to the random number generator.
190  const EngineType& engine() const { return mEngine; }
191 
192  /// Return a uniformly distributed random number in the range [0, 1).
193  FloatType operator()() { return mRand(mEngine); }
194 };
195 
197 
198 
199 /// @brief Simple random integer generator
200 /// @details Thread-safe as long as each thread has its own RandInt instance
201 template<typename IntType = int, typename EngineType = std::mt19937>
202 class RandInt
203 {
204 private:
205  using Distr = std::uniform_int_distribution<IntType>;
206  EngineType mEngine;
207  Distr mRand;
208 
209 public:
210  /// @brief Initialize the generator.
211  /// @param engine random number generator
212  /// @param imin,imax generate integers that are uniformly distributed over [imin, imax]
213  RandInt(const EngineType& engine, IntType imin, IntType imax):
214  mEngine(engine),
215  mRand(std::min(imin, imax), std::max(imin, imax))
216  {}
217 
218  /// @brief Initialize the generator.
219  /// @param seed seed value for the random number generator
220  /// @param imin,imax generate integers that are uniformly distributed over [imin, imax]
221  RandInt(unsigned int seed, IntType imin, IntType imax):
222  mEngine(static_cast<typename EngineType::result_type>(seed)),
223  mRand(std::min(imin, imax), std::max(imin, imax))
224  {}
225 
226  /// Change the range over which integers are distributed to [imin, imax].
227  void setRange(IntType imin, IntType imax)
228  {
229  mRand = Distr(std::min(imin, imax), std::max(imin, imax));
230  }
231 
232  /// Set the seed value for the random number generator
233  void setSeed(unsigned int seed)
234  {
235  mEngine.seed(static_cast<typename EngineType::result_type>(seed));
236  }
237 
238  /// Return a const reference to the random number generator.
239  const EngineType& engine() const { return mEngine; }
240 
241  /// Return a randomly-generated integer in the current range.
242  IntType operator()() { return mRand(mEngine); }
243 
244  /// @brief Return a randomly-generated integer in the new range [imin, imax],
245  /// without changing the current range.
246  IntType operator()(IntType imin, IntType imax)
247  {
248  const IntType lo = std::min(imin, imax), hi = std::max(imin, imax);
249  return mRand(mEngine, typename Distr::param_type(lo, hi));
250  }
251 };
252 
254 
255 
256 // ==========> Clamp <==================
257 
258 /// Return @a x clamped to [@a min, @a max]
259 template<typename Type>
260 inline Type
261 Clamp(Type x, Type min, Type max)
262 {
263  assert( !(min>max) );
264  return x > min ? x < max ? x : max : min;
265 }
266 
267 
268 /// Return @a x clamped to [0, 1]
269 template<typename Type>
270 inline Type
271 Clamp01(Type x) { return x > Type(0) ? x < Type(1) ? x : Type(1) : Type(0); }
272 
273 
274 /// Return @c true if @a x is outside [0,1]
275 template<typename Type>
276 inline bool
277 ClampTest01(Type &x)
278 {
279  if (x >= Type(0) && x <= Type(1)) return false;
280  x = x < Type(0) ? Type(0) : Type(1);
281  return true;
282 }
283 
284 /// @brief Return 0 if @a x < @a 0, 1 if @a x > 1 or else (3 &minus; 2 @a x) @a x&sup2;.
285 template<typename Type>
286 inline Type
288 {
289  return x > 0 ? x < 1 ? (3-2*x)*x*x : Type(1) : Type(0);
290 }
291 
292 /// @brief Return 0 if @a x < @a min, 1 if @a x > @a max or else (3 &minus; 2 @a t) @a t&sup2;,
293 /// where @a t = (@a x &minus; @a min)/(@a max &minus; @a min).
294 template<typename Type>
295 inline Type
296 SmoothUnitStep(Type x, Type min, Type max)
297 {
298  assert(min < max);
299  return SmoothUnitStep((x-min)/(max-min));
300 }
301 
302 
303 // ==========> Absolute Value <==================
304 
305 
306 //@{
307 /// Return the absolute value of the given quantity.
308 inline int32_t Abs(int32_t i) { return std::abs(i); }
309 inline int64_t Abs(int64_t i)
310 {
311  static_assert(sizeof(decltype(std::abs(i))) == sizeof(int64_t),
312  "std::abs(int64) broken");
313  return std::abs(i);
314 }
315 inline float Abs(float x) { return std::fabs(x); }
316 inline double Abs(double x) { return std::fabs(x); }
317 inline long double Abs(long double x) { return std::fabs(x); }
318 inline uint32_t Abs(uint32_t i) { return i; }
319 inline uint64_t Abs(uint64_t i) { return i; }
320 inline bool Abs(bool b) { return b; }
321 // On systems like macOS and FreeBSD, size_t and uint64_t are different types
322 template <typename T>
323 inline typename std::enable_if<std::is_same<T, size_t>::value, T>::type
324 Abs(T i) { return i; }
325 //@}
326 
327 
328 ////////////////////////////////////////
329 
330 
331 // ==========> Value Comparison <==================
332 
333 
334 /// Return @c true if @a x is exactly equal to zero.
335 template<typename Type>
336 inline bool
337 isZero(const Type& x)
338 {
340  return x == zeroVal<Type>();
342 }
343 
344 
345 /// @brief Return @c true if @a x is equal to zero to within
346 /// the default floating-point comparison tolerance.
347 template<typename Type>
348 inline bool
349 isApproxZero(const Type& x)
350 {
351  const Type tolerance = Type(zeroVal<Type>() + Tolerance<Type>::value());
352  return !(x > tolerance) && !(x < -tolerance);
353 }
354 
355 /// Return @c true if @a x is equal to zero to within the given tolerance.
356 template<typename Type>
357 inline bool
358 isApproxZero(const Type& x, const Type& tolerance)
359 {
360  return !(x > tolerance) && !(x < -tolerance);
361 }
362 
363 
364 /// Return @c true if @a x is less than zero.
365 template<typename Type>
366 inline bool
367 isNegative(const Type& x) { return x < zeroVal<Type>(); }
368 
369 // Return false, since bool values are never less than zero.
370 template<> inline bool isNegative<bool>(const bool&) { return false; }
371 
372 
373 /// Return @c true if @a x is finite.
374 inline bool
375 isFinite(const float x) { return std::isfinite(x); }
376 
377 /// Return @c true if @a x is finite.
379 inline bool
380 isFinite(const Type& x) { return std::isfinite(static_cast<double>(x)); }
381 
382 
383 /// Return @c true if @a x is an infinity value (either positive infinity or negative infinity).
384 inline bool
385 isInfinite(const float x) { return std::isinf(x); }
386 
387 /// Return @c true if @a x is an infinity value (either positive infinity or negative infinity).
389 inline bool
390 isInfinite(const Type& x) { return std::isinf(static_cast<double>(x)); }
391 
392 
393 /// Return @c true if @a x is a NaN (Not-A-Number) value.
394 inline bool
395 isNan(const float x) { return std::isnan(x); }
396 
397 /// Return @c true if @a x is a NaN (Not-A-Number) value.
399 inline bool
400 isNan(const Type& x) { return std::isnan(static_cast<double>(x)); }
401 
402 
403 /// Return @c true if @a a is equal to @a b to within the given tolerance.
404 template<typename Type>
405 inline bool
406 isApproxEqual(const Type& a, const Type& b, const Type& tolerance)
407 {
408  return !cwiseGreaterThan(Abs(a - b), tolerance);
409 }
410 
411 /// @brief Return @c true if @a a is equal to @a b to within
412 /// the default floating-point comparison tolerance.
413 template<typename Type>
414 inline bool
415 isApproxEqual(const Type& a, const Type& b)
416 {
417  const Type tolerance = Type(zeroVal<Type>() + Tolerance<Type>::value());
418  return isApproxEqual(a, b, tolerance);
419 }
420 
421 #define OPENVDB_EXACT_IS_APPROX_EQUAL(T) \
422  template<> inline bool isApproxEqual<T>(const T& a, const T& b) { return a == b; } \
423  template<> inline bool isApproxEqual<T>(const T& a, const T& b, const T&) { return a == b; } \
424  /**/
425 
428 
429 
430 /// @brief Return @c true if @a a is larger than @a b to within
431 /// the given tolerance, i.e., if @a b - @a a < @a tolerance.
432 template<typename Type>
433 inline bool
434 isApproxLarger(const Type& a, const Type& b, const Type& tolerance)
435 {
436  return (b - a < tolerance);
437 }
438 
439 
440 /// @brief Return @c true if @a a is exactly equal to @a b.
441 template<typename T0, typename T1>
442 inline bool
443 isExactlyEqual(const T0& a, const T1& b)
444 {
446  return a == b;
448 }
449 
450 
451 template<typename Type>
452 inline bool
453 isRelOrApproxEqual(const Type& a, const Type& b, const Type& absTol, const Type& relTol)
454 {
455  // First check to see if we are inside the absolute tolerance
456  // Necessary for numbers close to 0
457  if (!(Abs(a - b) > absTol)) return true;
458 
459  // Next check to see if we are inside the relative tolerance
460  // to handle large numbers that aren't within the abs tolerance
461  // but could be the closest floating point representation
462  double relError;
463  if (Abs(b) > Abs(a)) {
464  relError = Abs((a - b) / b);
465  } else {
466  relError = Abs((a - b) / a);
467  }
468  return (relError <= relTol);
469 }
470 
471 template<>
472 inline bool
473 isRelOrApproxEqual(const bool& a, const bool& b, const bool&, const bool&)
474 {
475  return (a == b);
476 }
477 
478 inline int32_t
479 floatToInt32(const float f)
480 {
481  // switch to std:bit_cast in C++20
482  static_assert(sizeof(int32_t) == sizeof f, "`float` has an unexpected size.");
483  int32_t ret;
484  std::memcpy(&ret, &f, sizeof(int32_t));
485  return ret;
486 }
487 
488 inline int64_t
489 doubleToInt64(const double d)
490 {
491  // switch to std:bit_cast in C++20
492  static_assert(sizeof(int64_t) == sizeof d, "`double` has an unexpected size.");
493  int64_t ret;
494  std::memcpy(&ret, &d, sizeof(int64_t));
495  return ret;
496 }
497 
498 // aUnitsInLastPlace is the allowed difference between the least significant digits
499 // of the numbers' floating point representation
500 // Please read the reference paper before trying to use isUlpsEqual
501 // http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
502 inline bool
503 isUlpsEqual(const double aLeft, const double aRight, const int64_t aUnitsInLastPlace)
504 {
505  int64_t longLeft = doubleToInt64(aLeft);
506  // Because of 2's complement, must restore lexicographical order
507  if (longLeft < 0) {
508  longLeft = INT64_C(0x8000000000000000) - longLeft;
509  }
510 
511  int64_t longRight = doubleToInt64(aRight);
512  // Because of 2's complement, must restore lexicographical order
513  if (longRight < 0) {
514  longRight = INT64_C(0x8000000000000000) - longRight;
515  }
516 
517  int64_t difference = Abs(longLeft - longRight);
518  return (difference <= aUnitsInLastPlace);
519 }
520 
521 inline bool
522 isUlpsEqual(const float aLeft, const float aRight, const int32_t aUnitsInLastPlace)
523 {
524  int32_t intLeft = floatToInt32(aLeft);
525  // Because of 2's complement, must restore lexicographical order
526  if (intLeft < 0) {
527  intLeft = 0x80000000 - intLeft;
528  }
529 
530  int32_t intRight = floatToInt32(aRight);
531  // Because of 2's complement, must restore lexicographical order
532  if (intRight < 0) {
533  intRight = 0x80000000 - intRight;
534  }
535 
536  int32_t difference = Abs(intLeft - intRight);
537  return (difference <= aUnitsInLastPlace);
538 }
539 
540 
541 ////////////////////////////////////////
542 
543 
544 // ==========> Pow <==================
545 
546 /// Return @a x<sup>2</sup>.
547 template<typename Type>
548 inline Type Pow2(Type x) { return x*x; }
549 
550 /// Return @a x<sup>3</sup>.
551 template<typename Type>
552 inline Type Pow3(Type x) { return x*x*x; }
553 
554 /// Return @a x<sup>4</sup>.
555 template<typename Type>
556 inline Type Pow4(Type x) { return Pow2(Pow2(x)); }
557 
558 /// Return @a x<sup>n</sup>.
559 template<typename Type>
560 Type
561 Pow(Type x, int n)
562 {
563  Type ans = 1;
564  if (n < 0) {
565  n = -n;
566  x = Type(1)/x;
567  }
568  while (n--) ans *= x;
569  return ans;
570 }
571 
572 //@{
573 /// Return @a b<sup>e</sup>.
574 inline float
575 Pow(float b, float e)
576 {
577  assert( b >= 0.0f && "Pow(float,float): base is negative" );
578  return powf(b,e);
579 }
580 
581 inline double
582 Pow(double b, double e)
583 {
584  assert( b >= 0.0 && "Pow(double,double): base is negative" );
585  return std::pow(b,e);
586 }
587 //@}
588 
589 
590 // ==========> Max <==================
591 
592 /// Return the maximum of two values
593 template<typename Type>
594 inline const Type&
595 Max(const Type& a, const Type& b)
596 {
597  return std::max(a,b);
598 }
599 
600 /// Return the maximum of three values
601 template<typename Type>
602 inline const Type&
603 Max(const Type& a, const Type& b, const Type& c)
604 {
605  return std::max(std::max(a,b), c);
606 }
607 
608 /// Return the maximum of four values
609 template<typename Type>
610 inline const Type&
611 Max(const Type& a, const Type& b, const Type& c, const Type& d)
612 {
613  return std::max(std::max(a,b), std::max(c,d));
614 }
615 
616 /// Return the maximum of five values
617 template<typename Type>
618 inline const Type&
619 Max(const Type& a, const Type& b, const Type& c, const Type& d, const Type& e)
620 {
621  return std::max(std::max(a,b), Max(c,d,e));
622 }
623 
624 /// Return the maximum of six values
625 template<typename Type>
626 inline const Type&
627 Max(const Type& a, const Type& b, const Type& c, const Type& d, const Type& e, const Type& f)
628 {
629  return std::max(Max(a,b,c), Max(d,e,f));
630 }
631 
632 /// Return the maximum of seven values
633 template<typename Type>
634 inline const Type&
635 Max(const Type& a, const Type& b, const Type& c, const Type& d,
636  const Type& e, const Type& f, const Type& g)
637 {
638  return std::max(Max(a,b,c,d), Max(e,f,g));
639 }
640 
641 /// Return the maximum of eight values
642 template<typename Type>
643 inline const Type&
644 Max(const Type& a, const Type& b, const Type& c, const Type& d,
645  const Type& e, const Type& f, const Type& g, const Type& h)
646 {
647  return std::max(Max(a,b,c,d), Max(e,f,g,h));
648 }
649 
650 
651 // ==========> Min <==================
652 
653 /// Return the minimum of two values
654 template<typename Type>
655 inline const Type&
656 Min(const Type& a, const Type& b) { return std::min(a, b); }
657 
658 /// Return the minimum of three values
659 template<typename Type>
660 inline const Type&
661 Min(const Type& a, const Type& b, const Type& c) { return std::min(std::min(a, b), c); }
662 
663 /// Return the minimum of four values
664 template<typename Type>
665 inline const Type&
666 Min(const Type& a, const Type& b, const Type& c, const Type& d)
667 {
668  return std::min(std::min(a, b), std::min(c, d));
669 }
670 
671 /// Return the minimum of five values
672 template<typename Type>
673 inline const Type&
674 Min(const Type& a, const Type& b, const Type& c, const Type& d, const Type& e)
675 {
676  return std::min(std::min(a,b), Min(c,d,e));
677 }
678 
679 /// Return the minimum of six values
680 template<typename Type>
681 inline const Type&
682 Min(const Type& a, const Type& b, const Type& c, const Type& d, const Type& e, const Type& f)
683 {
684  return std::min(Min(a,b,c), Min(d,e,f));
685 }
686 
687 /// Return the minimum of seven values
688 template<typename Type>
689 inline const Type&
690 Min(const Type& a, const Type& b, const Type& c, const Type& d,
691  const Type& e, const Type& f, const Type& g)
692 {
693  return std::min(Min(a,b,c,d), Min(e,f,g));
694 }
695 
696 /// Return the minimum of eight values
697 template<typename Type>
698 inline const Type&
699 Min(const Type& a, const Type& b, const Type& c, const Type& d,
700  const Type& e, const Type& f, const Type& g, const Type& h)
701 {
702  return std::min(Min(a,b,c,d), Min(e,f,g,h));
703 }
704 
705 
706 // ============> Exp <==================
707 
708 /// Return @a e<sup>x</sup>.
709 template<typename Type>
710 inline Type Exp(const Type& x) { return std::exp(x); }
711 
712 // ============> Sin <==================
713 
714 //@{
715 /// Return sin @a x.
716 inline float Sin(const float& x) { return std::sin(x); }
717 
718 inline double Sin(const double& x) { return std::sin(x); }
719 //@}
720 
721 // ============> Cos <==================
722 
723 //@{
724 /// Return cos @a x.
725 inline float Cos(const float& x) { return std::cos(x); }
726 
727 inline double Cos(const double& x) { return std::cos(x); }
728 //@}
729 
730 
731 ////////////////////////////////////////
732 
733 
734 /// Return the sign of the given value as an integer (either -1, 0 or 1).
735 template <typename Type>
736 inline int Sign(const Type &x) { return (zeroVal<Type>() < x) - (x < zeroVal<Type>()); }
737 
738 
739 /// @brief Return @c true if @a a and @a b have different signs.
740 /// @note Zero is considered a positive number.
741 template <typename Type>
742 inline bool
743 SignChange(const Type& a, const Type& b)
744 {
745  return ( (a<zeroVal<Type>()) ^ (b<zeroVal<Type>()) );
746 }
747 
748 
749 /// @brief Return @c true if the interval [@a a, @a b] includes zero,
750 /// i.e., if either @a a or @a b is zero or if they have different signs.
751 template <typename Type>
752 inline bool
753 ZeroCrossing(const Type& a, const Type& b)
754 {
755  return a * b <= zeroVal<Type>();
756 }
757 
758 
759 //@{
760 /// Return the square root of a floating-point value.
761 inline float Sqrt(float x) { return std::sqrt(x); }
762 inline double Sqrt(double x) { return std::sqrt(x); }
763 inline long double Sqrt(long double x) { return std::sqrt(x); }
764 //@}
765 
766 
767 //@{
768 /// Return the cube root of a floating-point value.
769 inline float Cbrt(float x) { return std::cbrt(x); }
770 inline double Cbrt(double x) { return std::cbrt(x); }
771 inline long double Cbrt(long double x) { return std::cbrt(x); }
772 //@}
773 
774 
775 //@{
776 /// Return the remainder of @a x / @a y.
777 inline int Mod(int x, int y) { return (x % y); }
778 inline float Mod(float x, float y) { return std::fmod(x, y); }
779 inline double Mod(double x, double y) { return std::fmod(x, y); }
780 inline long double Mod(long double x, long double y) { return std::fmod(x, y); }
781 template<typename Type> inline Type Remainder(Type x, Type y) { return Mod(x, y); }
782 //@}
783 
784 
785 //@{
786 /// Return @a x rounded up to the nearest integer.
787 inline float RoundUp(float x) { return std::ceil(x); }
788 inline double RoundUp(double x) { return std::ceil(x); }
789 inline long double RoundUp(long double x) { return std::ceil(x); }
790 //@}
791 /// Return @a x rounded up to the nearest multiple of @a base.
792 template<typename Type>
793 inline Type
794 RoundUp(Type x, Type base)
795 {
796  Type remainder = Remainder(x, base);
797  return remainder ? x-remainder+base : x;
798 }
799 
800 
801 //@{
802 /// Return @a x rounded down to the nearest integer.
803 inline float RoundDown(float x) { return std::floor(x); }
804 inline double RoundDown(double x) { return std::floor(x); }
805 inline long double RoundDown(long double x) { return std::floor(x); }
806 //@}
807 /// Return @a x rounded down to the nearest multiple of @a base.
808 template<typename Type>
809 inline Type
810 RoundDown(Type x, Type base)
811 {
812  Type remainder = Remainder(x, base);
813  return remainder ? x-remainder : x;
814 }
815 
816 
817 //@{
818 /// Return @a x rounded to the nearest integer.
819 inline float Round(float x) { return RoundDown(x + 0.5f); }
820 inline double Round(double x) { return RoundDown(x + 0.5); }
821 inline long double Round(long double x) { return RoundDown(x + 0.5l); }
822 //@}
823 
824 
825 /// Return the euclidean remainder of @a x.
826 /// Note unlike % operator this will always return a positive result
827 template<typename Type>
828 inline Type
829 EuclideanRemainder(Type x) { return x - RoundDown(x); }
830 
831 
832 /// Return the integer part of @a x.
833 template<typename Type>
834 inline Type
835 IntegerPart(Type x)
836 {
837  return (x > 0 ? RoundDown(x) : RoundUp(x));
838 }
839 
840 /// Return the fractional part of @a x.
841 template<typename Type>
842 inline Type
843 FractionalPart(Type x) { return Mod(x,Type(1)); }
844 
845 
846 //@{
847 /// Return the floor of @a x.
848 inline int Floor(float x) { return int(RoundDown(x)); }
849 inline int Floor(double x) { return int(RoundDown(x)); }
850 inline int Floor(long double x) { return int(RoundDown(x)); }
851 //@}
852 
853 
854 //@{
855 /// Return the ceiling of @a x.
856 inline int Ceil(float x) { return int(RoundUp(x)); }
857 inline int Ceil(double x) { return int(RoundUp(x)); }
858 inline int Ceil(long double x) { return int(RoundUp(x)); }
859 //@}
860 
861 
862 /// Return @a x if it is greater or equal in magnitude than @a delta. Otherwise, return zero.
863 template<typename Type>
864 inline Type Chop(Type x, Type delta) { return (Abs(x) < delta ? zeroVal<Type>() : x); }
865 
866 
867 /// Return @a x truncated to the given number of decimal digits.
868 template<typename Type>
869 inline Type
870 Truncate(Type x, unsigned int digits)
871 {
872  Type tenth = static_cast<Type>(Pow(size_t(10), digits));
873  return RoundDown(x*tenth+0.5)/tenth;
874 }
875 
876 ////////////////////////////////////////
877 
878 
879 /// @brief 8-bit integer values print to std::ostreams as characters.
880 /// Cast them so that they print as integers instead.
881 template<typename T>
882 inline auto PrintCast(const T& val) -> typename std::enable_if<!std::is_same<T, int8_t>::value
883  && !std::is_same<T, uint8_t>::value, const T&>::type { return val; }
884 inline int32_t PrintCast(int8_t val) { return int32_t(val); }
885 inline uint32_t PrintCast(uint8_t val) { return uint32_t(val); }
886 
887 
888 ////////////////////////////////////////
889 
890 
891 /// Return the inverse of @a x.
892 template<typename Type>
893 inline Type
894 Inv(Type x)
895 {
896  assert(x);
897  return Type(1)/x;
898 }
899 
900 
901 enum Axis {
902  X_AXIS = 0,
903  Y_AXIS = 1,
904  Z_AXIS = 2
905 };
906 
907 // enum values are consistent with their historical mx analogs.
917 };
918 
919 
920 template <typename S, typename T>
921 struct promote {
922  using type = typename boost::numeric::conversion_traits<S, T>::supertype;
923 };
924 
925 /// @brief Return the index [0,1,2] of the smallest value in a 3D vector.
926 /// @note This methods assumes operator[] exists.
927 /// @details The return value corresponds to the largest index of the of
928 /// the smallest vector components.
929 template<typename Vec3T>
930 size_t
931 MinIndex(const Vec3T& v)
932 {
933  size_t r = 0;
934  for (size_t i = 1; i < 3; ++i) {
935  // largest index (backwards compatibility)
936  if (v[i] <= v[r]) r = i;
937  }
938  return r;
939 }
940 
941 /// @brief Return the index [0,1,2] of the largest value in a 3D vector.
942 /// @note This methods assumes operator[] exists.
943 /// @details The return value corresponds to the largest index of the of
944 /// the largest vector components.
945 template<typename Vec3T>
946 size_t
947 MaxIndex(const Vec3T& v)
948 {
949  size_t r = 0;
950  for (size_t i = 1; i < 3; ++i) {
951  // largest index (backwards compatibility)
952  if (v[i] >= v[r]) r = i;
953  }
954  return r;
955 }
956 
957 } // namespace math
958 } // namespace OPENVDB_VERSION_NAME
959 } // namespace openvdb
960 
961 #endif // OPENVDB_MATH_MATH_HAS_BEEN_INCLUDED
Simple random integer generator.
Definition: Math.h:202
#define OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN
Bracket code with OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN/_END, to inhibit warnings about type conve...
Definition: Platform.h:204
void setSeed(unsigned int seed)
Set the seed value for the random number generator.
Definition: Math.h:184
Type Clamp01(Type x)
Return x clamped to [0, 1].
Definition: Math.h:271
Type Truncate(Type x, unsigned int digits)
Return x truncated to the given number of decimal digits.
Definition: Math.h:870
int64_t doubleToInt64(const double d)
Definition: Math.h:489
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:106
static float value()
Definition: Math.h:149
bool isRelOrApproxEqual(const bool &a, const bool &b, const bool &, const bool &)
Definition: Math.h:473
Type IntegerPart(Type x)
Return the integer part of x.
Definition: Math.h:835
Type Exp(const Type &x)
Return ex.
Definition: Math.h:710
Definition: Math.h:902
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:753
void setRange(IntType imin, IntType imax)
Change the range over which integers are distributed to [imin, imax].
Definition: Math.h:227
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:246
const EngineType & engine() const
Return a const reference to the random number generator.
Definition: Math.h:190
size_t MaxIndex(const Vec3T &v)
Return the index [0,1,2] of the largest value in a 3D vector.
Definition: Math.h:947
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:110
bool cwiseLessThan(const Type1 &a, const Type2 &b)
Componentwise less than for POD types.
Definition: Math.h:99
Rand01(unsigned int seed)
Initialize the generator.
Definition: Math.h:181
long double Cbrt(long double x)
Return the cube root of a floating-point value.
Definition: Math.h:771
RandInt(unsigned int seed, IntType imin, IntType imax)
Initialize the generator.
Definition: Math.h:221
Definition: Coord.h:587
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:443
Type Pow4(Type x)
Return x4.
Definition: Math.h:556
auto cwiseAdd(const Type1 &v, const Type2 s)
Componentwise adder for POD types.
Definition: Math.h:90
static double value()
Definition: Math.h:157
double Sin(const double &x)
Return sin x.
Definition: Math.h:718
#define OPENVDB_NO_TYPE_CONVERSION_WARNING_END
Definition: Platform.h:205
Type RoundDown(Type x, Type base)
Return x rounded down to the nearest multiple of base.
Definition: Math.h:810
Rand01(const EngineType &engine)
Initialize the generator.
Definition: Math.h:177
Definition: Math.h:921
bool isUlpsEqual(const float aLeft, const float aRight, const int32_t aUnitsInLastPlace)
Definition: Math.h:522
bool isFinite(const Type &x)
Return true if x is finite.
Definition: Math.h:380
FloatType operator()()
Return a uniformly distributed random number in the range [0, 1).
Definition: Math.h:193
const Type & Min(const Type &a, const Type &b, const Type &c, const Type &d, const Type &e, const Type &f, const Type &g, const Type &h)
Return the minimum of eight values.
Definition: Math.h:699
double Cos(const double &x)
Return cos x.
Definition: Math.h:727
Type Pow2(Type x)
Return x2.
Definition: Math.h:548
Type Clamp(Type x, Type min, Type max)
Return x clamped to [min, max].
Definition: Math.h:261
double Pow(double b, double e)
Return be.
Definition: Math.h:582
bool cwiseGreaterThan(const Type1 &a, const Type2 &b)
Componentwise greater than for POD types.
Definition: Math.h:108
void setSeed(unsigned int seed)
Set the seed value for the random number generator.
Definition: Math.h:233
bool SignChange(const Type &a, const Type &b)
Return true if a and b have different signs.
Definition: Math.h:743
size_t MinIndex(const Vec3T &v)
Return the index [0,1,2] of the smallest value in a 3D vector.
Definition: Math.h:931
bool ClampTest01(Type &x)
Return true if x is outside [0,1].
Definition: Math.h:277
Definition: Exceptions.h:13
Simple generator of random numbers over the range [0, 1)
Definition: Math.h:166
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:337
long double Round(long double x)
Return x rounded to the nearest integer.
Definition: Math.h:821
ValueT value
Definition: GridBuilder.h:1290
int Sign(const Type &x)
Return the sign of the given value as an integer (either -1, 0 or 1).
Definition: Math.h:736
long double Sqrt(long double x)
Return the square root of a floating-point value.
Definition: Math.h:763
#define OPENVDB_NO_FP_EQUALITY_WARNING_END
Definition: Math.h:49
std::enable_if< std::is_same< T, size_t >::value, T >::type Abs(T i)
Return the absolute value of the given quantity.
Definition: Math.h:324
Definition: Math.h:903
bool isInfinite(const Type &x)
Return true if x is an infinity value (either positive infinity or negative infinity).
Definition: Math.h:390
Axis
Definition: Math.h:901
bool isNegative< bool >(const bool &)
Definition: Math.h:370
Tolerance for floating-point comparison.
Definition: Math.h:148
Type Remainder(Type x, Type y)
Return the remainder of x / y.
Definition: Math.h:781
constexpr long double pi()
Definition: Math.h:122
#define OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN
Definition: Math.h:48
const Type & Max(const Type &a, const Type &b, const Type &c, const Type &d, const Type &e, const Type &f, const Type &g, const Type &h)
Return the maximum of eight values.
Definition: Math.h:644
Type SmoothUnitStep(Type x, Type min, Type max)
Return 0 if x < min, 1 if x > max or else (3 − 2 t) t², where t = (x − min)/(max − min)...
Definition: Math.h:296
int Ceil(long double x)
Return the ceiling of x.
Definition: Math.h:858
const EngineType & engine() const
Return a const reference to the random number generator.
Definition: Math.h:239
static double value()
Definition: Math.h:150
typename boost::numeric::conversion_traits< S, T >::supertype type
Definition: Math.h:922
IntType operator()()
Return a randomly-generated integer in the current range.
Definition: Math.h:242
RandInt(const EngineType &engine, IntType imin, IntType imax)
Initialize the generator.
Definition: Math.h:213
int32_t floatToInt32(const float f)
Definition: Math.h:479
int Floor(long double x)
Return the floor of x.
Definition: Math.h:850
Type EuclideanRemainder(Type x)
Definition: Math.h:829
Type Inv(Type x)
Return the inverse of x.
Definition: Math.h:894
std::string negative(const std::string &val)
Return the "negation" of the given string.
Definition: Math.h:143
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:415
static T value()
Definition: Math.h:155
Delta for small floating-point offsets.
Definition: Math.h:155
std::string operator+(const std::string &s, double)
Needed to support the (zeroVal<ValueType>() + val) idiom when ValueType is std::string.
Definition: Math.h:85
static T value()
Definition: Math.h:148
Type FractionalPart(Type x)
Return the fractional part of x.
Definition: Math.h:843
long double Mod(long double x, long double y)
Return the remainder of x / y.
Definition: Math.h:780
static float value()
Definition: Math.h:156
Type Chop(Type x, Type delta)
Return x if it is greater or equal in magnitude than delta. Otherwise, return zero.
Definition: Math.h:864
uint32_t PrintCast(uint8_t val)
Definition: Math.h:885
bool isNan(const Type &x)
Return true if x is a NaN (Not-A-Number) value.
Definition: Math.h:400
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
T zeroVal()
Return the value of type T that corresponds to zero.
Definition: Math.h:70
bool zeroVal< bool >()
Return the bool value that corresponds to zero.
Definition: Math.h:74
#define OPENVDB_EXACT_IS_APPROX_EQUAL(T)
Definition: Math.h:421
Definition: Math.h:904
Type Pow3(Type x)
Return x3.
Definition: Math.h:552
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:434
bool isNegative(const Type &x)
Return true if x is less than zero.
Definition: Math.h:367
RotationOrder
Definition: Math.h:908
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:212
bool isApproxZero(const Type &x, const Type &tolerance)
Return true if x is equal to zero to within the given tolerance.
Definition: Math.h:358
Type RoundUp(Type x, Type base)
Return x rounded up to the nearest multiple of base.
Definition: Math.h:794