OpenVDB  1.2.0
Math.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2013 DreamWorks Animation LLC
4 //
5 // All rights reserved. This software is distributed under the
6 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
7 //
8 // Redistributions of source code must retain the above copyright
9 // and license notice and the following restrictions and disclaimer.
10 //
11 // * Neither the name of DreamWorks Animation nor the names of
12 // its contributors may be used to endorse or promote products derived
13 // from this software without specific prior written permission.
14 //
15 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
19 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL,
20 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
21 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
22 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
23 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
25 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 // IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE
27 // LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00.
28 //
30 //
34 
35 #ifndef OPENVDB_MATH_HAS_BEEN_INCLUDED
36 #define OPENVDB_MATH_HAS_BEEN_INCLUDED
37 
38 #include <assert.h>
39 #include <algorithm> //for std::max()
40 #include <cmath> //for floor(), ceil() and sqrt()
41 #include <math.h> //for pow(), fabs() etc
42 #include <cstdlib> //for srand, abs(int)
43 #include <limits> //for std::numeric_limits<Type>::max()
44 #include <string>
45 #include <boost/numeric/conversion/conversion_traits.hpp>
46 #include <openvdb/Platform.h>
47 #include <openvdb/version.h>
48 
49 
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 #else
61  // For GCC, #pragma GCC diagnostic ignored "-Wfloat-equal"
62  // isn't working until gcc 4.2+,
63  // Trying
64  // #pragma GCC system_header
65  // creates other problems, most notably "warning: will never be executed"
66  // in from templates, unsure of how to work around.
67  // If necessary, could use integer based comparisons for equality
68  #define OPENVDB_NO_FP_EQUALITY_WARNING_BEGIN
69  #define OPENVDB_NO_FP_EQUALITY_WARNING_END
70 #endif
71 
72 namespace openvdb {
74 namespace OPENVDB_VERSION_NAME {
75 
80 template<typename T> inline T zeroVal() { return T(0); }
82 template<> inline std::string zeroVal<std::string>() { return ""; }
84 template<> inline bool zeroVal<bool>() { return false; }
85 
86 
88 template<typename T> inline T toleranceValue() { return T(1e-8); }
90 template<> inline float toleranceValue<float>() { return float(1e-6); }
91 template<typename T> struct tolerance { static T value() { return 0; } };
92 template<> struct tolerance<float> { static float value() { return 1e-8f; } };
93 template<> struct tolerance<double> { static double value() { return 1e-15; } };
95 
96 
98 
99 inline std::string operator+(const std::string& s, bool) { return s; }
102 inline std::string operator+(const std::string& s, int) { return s; }
103 inline std::string operator+(const std::string& s, float) { return s; }
104 inline std::string operator+(const std::string& s, double) { return s; }
106 
107 
109 template<typename T> inline T negative(const T& val) { return T(-val); }
113 template<> inline bool negative(const bool& val) { return !val; }
115 template<> inline std::string negative(const std::string& val) { return val; }
117 
118 
119 namespace math {
120 
121 // ==========> Random Values <==================
122 
124 inline void randSeed(unsigned int seed) { srand(seed); }
125 
126 
128 inline double randUniform() { return (double)(rand() / (RAND_MAX + 1.0)); }
129 
130 
133 {
134 protected:
135  int my_min, my_range;
136 public:
137  RandomInt(unsigned int seed, int min, int max): my_min(min), my_range(max-min+1) {
138  assert(min<max && "RandomInt: invalid arguments");
139  randSeed(seed);
140  }
141  void setRange(int min, int max) { my_min = min; my_range = max-min+1; }
142  int operator()() const { return rand() % my_range + my_min; }
143  int operator()(int min, int max) const { return rand() % (max-min+1) + min; }
144 };
145 
146 
147 // ==========> Clamp <==================
148 
150 template<typename Type>
151 inline Type
152 Clamp(Type x, Type min, Type max)
153 {
154  assert(min<max);
155  return x > min ? x < max ? x : max : min;
156 }
157 
158 
160 template<typename Type>
161 inline Type
162 Clamp01(Type x) { return x > Type(0) ? x < Type(1) ? x : Type(1) : Type(0); }
163 
164 
166 template<typename Type>
167 inline bool
168 ClampTest01(Type &x)
169 {
170  if (x >= Type(0) && x <= Type(1)) return false;
171  x = x < Type(0) ? Type(0) : Type(1);
172  return true;
173 }
174 
175 
178 template<typename Type>
179 inline Type
180 SmoothUnitStep(Type x, Type min, Type max)
181 {
182  assert(min < max);
183  const Type t = (x-min)/(max-min);
184  return t > 0 ? t < 1 ? (3-2*t)*t*t : Type(1) : Type(0);
185 }
186 
187 
188 // ==========> Absolute Value <==================
189 
190 
192 inline int32_t Abs(int32_t i) { return abs(i); }
194 inline int64_t Abs(int64_t i)
195 {
196 #ifdef _MSC_VER
197  return (i < int64_t(0) ? -i : i);
198 #else
199  return abs(i);
200 #endif
201 }
202 inline float Abs(float x) { return fabs(x); }
203 inline double Abs(double x) { return fabs(x); }
204 inline long double Abs(long double x) { return fabs(x); }
205 inline uint32_t Abs(uint32_t i) { return i; }
206 inline uint64_t Abs(uint64_t i) { return i; }
208 
209 
211 
212 
213 // ==========> Value Comparison <==================
214 
215 
217 template<typename Type>
218 inline bool
219 isZero(const Type& x)
220 {
222  return x == zeroVal<Type>();
224 }
225 
226 
229 template<typename Type>
230 inline bool
231 isApproxZero(const Type& x)
232 {
233  const Type tolerance = Type(zeroVal<Type>() + toleranceValue<Type>());
234  return x < tolerance && x > -tolerance;
235 }
236 
238 template<typename Type>
239 inline bool
240 isApproxZero(const Type& x, const Type& tolerance)
241 {
242  return x < tolerance && x > -tolerance;
243 }
244 
245 
246 template<typename Type>
247 inline bool
248 isNegative(const Type& x) { return x < zeroVal<Type>(); }
249 
250 
251 template<typename Type>
252 inline bool
253 isApproxEqual(const Type& a, const Type& b)
254 {
255  const Type tolerance = Type(zeroVal<Type>() + toleranceValue<Type>());
256  return !(Abs(a - b) > tolerance);
257 }
258 
259 template<typename Type>
260 inline bool
261 isApproxEqual(const Type& a, const Type& b, const Type& tolerance)
262 {
263  return !(Abs(a - b) > tolerance);
264 }
265 
266 #define OPENVDB_EXACT_IS_APPROX_EQUAL(T) \
267  template<> inline bool isApproxEqual<T>(const T& a, const T& b) { return a == b; } \
268  template<> inline bool isApproxEqual<T>(const T& a, const T& b, const T&) { return a == b; } \
269 
270 
273 
274 
275 template<typename T0, typename T1>
276 inline bool
277 isExactlyEqual(const T0& a, const T1& b)
278 {
280  return a == b;
282 }
283 
284 
285 template<typename Type>
286 inline bool
287 isRelOrApproxEqual(const Type& a, const Type& b, const Type& absTol, const Type& relTol)
288 {
289  // First check to see if we are inside the absolute tolerance
290  // Necessary for numbers close to 0
291  if (!(Abs(a - b) > absTol)) return true;
292 
293  // Next check to see if we are inside the relative tolerance
294  // to handle large numbers that aren't within the abs tolerance
295  // but could be the closest floating point representation
296  double relError;
297  if (Abs(b) > Abs(a)) {
298  relError = Abs((a - b) / b);
299  } else {
300  relError = Abs((a - b) / a);
301  }
302  return (relError <= relTol);
303 }
304 
305 template<>
306 inline bool
307 isRelOrApproxEqual(const bool& a, const bool& b, const bool&, const bool&)
308 {
309  return (a == b);
310 }
311 
312 
313 // Avoid strict aliasing issues by using type punning
314 // http://cellperformance.beyond3d.com/articles/2006/06/understanding-strict-aliasing.html
315 // Using "casting through a union(2)"
316 inline int32_t
317 floatToInt32(const float aFloatValue)
318 {
319  union FloatOrInt32 { float floatValue; int32_t int32Value; };
320  const FloatOrInt32* foi = reinterpret_cast<const FloatOrInt32*>(&aFloatValue);
321  return foi->int32Value;
322 }
323 
324 
325 inline int64_t
326 doubleToInt64(const double aDoubleValue)
327 {
328  union DoubleOrInt64 { double doubleValue; int64_t int64Value; };
329  const DoubleOrInt64* dol = reinterpret_cast<const DoubleOrInt64*>(&aDoubleValue);
330  return dol->int64Value;
331 }
332 
333 
334 // aUnitsInLastPlace is the allowed difference between the least significant digits
335 // of the numbers' floating point representation
336 // Please read refernce paper before trying to use isUlpsEqual
337 // http://www.cygnus-software.com/papers/comparingFloats/comparingFloats.htm
338 inline bool
339 isUlpsEqual(const double aLeft, const double aRight, const int64_t aUnitsInLastPlace)
340 {
341  int64_t longLeft = doubleToInt64(aLeft);
342  // Because of 2's complement, must restore lexicographical order
343  if (longLeft < 0) {
344  longLeft = INT64_C(0x8000000000000000) - longLeft;
345  }
346 
347  int64_t longRight = doubleToInt64(aRight);
348  // Because of 2's complement, must restore lexicographical order
349  if (longRight < 0) {
350  longRight = INT64_C(0x8000000000000000) - longRight;
351  }
352 
353  int64_t difference = labs(longLeft - longRight);
354  return (difference <= aUnitsInLastPlace);
355 }
356 
357 inline bool
358 isUlpsEqual(const float aLeft, const float aRight, const int32_t aUnitsInLastPlace)
359 {
360  int32_t intLeft = floatToInt32(aLeft);
361  // Because of 2's complement, must restore lexicographical order
362  if (intLeft < 0) {
363  intLeft = 0x80000000 - intLeft;
364  }
365 
366  int32_t intRight = floatToInt32(aRight);
367  // Because of 2's complement, must restore lexicographical order
368  if (intRight < 0) {
369  intRight = 0x80000000 - intRight;
370  }
371 
372  int32_t difference = abs(intLeft - intRight);
373  return (difference <= aUnitsInLastPlace);
374 }
375 
376 
378 
379 
380 // ==========> Pow <==================
381 
383 template<typename Type>
384 inline Type Pow2(Type x) { return x*x; }
385 
387 template<typename Type>
388 inline Type Pow3(Type x) { return x*x*x; }
389 
391 template<typename Type>
392 inline Type Pow4(Type x) { return Pow2(Pow2(x)); }
393 
395 template<typename Type>
396 Type
397 Pow(Type x, int n)
398 {
399  Type ans = 1;
400  if (n < 0) {
401  n = -n;
402  x = Type(1)/x;
403  }
404  while (n--) ans *= x;
405  return ans;
406 }
407 
409 inline float
411 Pow(float b, float e)
412 {
413  assert( b >= 0.0f && "Pow(float,float): base is negative" );
414  return powf(b,e);
415 }
416 
417 inline double
418 Pow(double b, double e)
419 {
420  assert( b >= 0.0 && "Pow(double,double): base is negative" );
421  return pow(b,e);
422 }
424 
425 
426 // ==========> Max <==================
427 
429 template<typename Type>
430 inline const Type&
431 Max(const Type& a, const Type& b)
432 {
433  return std::max(a,b) ;
434 }
435 
437 template<typename Type>
438 inline const Type&
439 Max(const Type& a, const Type& b, const Type& c)
440 {
441  return std::max( std::max(a,b), c ) ;
442 }
443 
445 template<typename Type>
446 inline const Type&
447 Max(const Type& a, const Type& b, const Type& c, const Type& d)
448 {
449  return std::max(std::max(a,b), std::max(c,d));
450 }
451 
453 template<typename Type>
454 inline const Type&
455 Max(const Type& a, const Type& b, const Type& c, const Type& d, const Type& e)
456 {
457  return std::max(std::max(a,b), Max(c,d,e));
458 }
459 
461 template<typename Type>
462 inline const Type&
463 Max(const Type& a, const Type& b, const Type& c, const Type& d, const Type& e, const Type& f)
464 {
465  return std::max(Max(a,b,c), Max(d,e,f));
466 }
467 
469 template<typename Type>
470 inline const Type&
471 Max(const Type& a, const Type& b, const Type& c, const Type& d,
472  const Type& e, const Type& f, const Type& g)
473 {
474  return std::max(Max(a,b,c,d), Max(e,f,g));
475 }
476 
478 template<typename Type>
479 inline const Type&
480 Max(const Type& a, const Type& b, const Type& c, const Type& d,
481  const Type& e, const Type& f, const Type& g, const Type& h)
482 {
483  return std::max(Max(a,b,c,d), Max(e,f,g,h));
484 }
485 
486 
487 // ==========> Min <==================
488 
490 template<typename Type>
491 inline const Type&
492 Min(const Type& a, const Type& b) { return std::min(a, b); }
493 
495 template<typename Type>
496 inline const Type&
497 Min(const Type& a, const Type& b, const Type& c) { return std::min(std::min(a, b), c); }
498 
500 template<typename Type>
501 inline const Type&
502 Min(const Type& a, const Type& b, const Type& c, const Type& d)
503 {
504  return std::min(std::min(a, b), std::min(c, d));
505 }
506 
508 template<typename Type>
509 inline const Type&
510 Min(const Type& a, const Type& b, const Type& c, const Type& d, const Type& e)
511 {
512  return std::min(std::min(a,b), Min(c,d,e));
513 }
514 
516 template<typename Type>
517 inline const Type&
518 Min(const Type& a, const Type& b, const Type& c, const Type& d, const Type& e, const Type& f)
519 {
520  return std::min(Min(a,b,c), Min(d,e,f));
521 }
522 
524 template<typename Type>
525 inline const Type&
526 Min(const Type& a, const Type& b, const Type& c, const Type& d,
527  const Type& e, const Type& f, const Type& g)
528 {
529  return std::min(Min(a,b,c,d), Min(e,f,g));
530 }
531 
533 template<typename Type>
534 inline const Type&
535 Min(const Type& a, const Type& b, const Type& c, const Type& d,
536  const Type& e, const Type& f, const Type& g, const Type& h)
537 {
538  return std::min(Min(a,b,c,d), Min(e,f,g,h));
539 }
540 
541 
543 
544 
546 template <typename Type>
547 inline int Sign(const Type &x) { return (zeroVal<Type>() < x) - (x < zeroVal<Type>()); }
548 
549 
552 template <typename Type>
553 inline bool
554 SignChange(const Type& a, const Type& b)
555 {
556  return ( (a<zeroVal<Type>()) ^ (b<zeroVal<Type>()) );
557 }
558 
559 
562 template <typename Type>
563 inline bool
564 ZeroCrossing(const Type& a, const Type& b)
565 {
566  return a * b <= zeroVal<Type>();
567 }
568 
569 
571 inline float Sqrt(float x) { return sqrtf(x); }
573 inline double Sqrt(double x) { return sqrt(x); }
574 inline long double Sqrt(long double x) { return sqrtl(x); }
576 
577 
579 inline float Cbrt(float x) { return cbrtf(x); }
581 inline double Cbrt(double x) { return cbrtf(x); }
582 inline long double Cbrt(long double x) { return cbrtl(x); }
584 
585 
587 inline int Mod(int x, int y) { return (x % y); };
589 inline float Mod(float x, float y) { return fmodf(x,y); }
590 inline double Mod(double x, double y) { return fmod(x,y); }
591 inline long double Mod(long double x, long double y) { return fmodl(x,y); }
592 template<typename Type> inline Type Remainder(Type x, Type y) { return Mod(x,y); }
594 
595 
597 inline float RoundUp(float x) { return ceilf(x); }
599 inline double RoundUp(double x) { return ceil(x); }
600 inline long double RoundUp(long double x) { return ceill(x); }
602 template<typename Type>
604 inline Type
605 RoundUp(Type x, Type base)
606 {
607  Type remainder = Remainder(x, base);
608  return remainder ? x-remainder+base : x;
609 }
610 
611 
613 inline float RoundDown(float x) { return floorf(x); }
615 inline double RoundDown(double x) { return floor(x); }
616 inline long double RoundDown(long double x) { return floorl(x); }
617 template <typename Type> inline Type Round(Type x) { return RoundDown(x+0.5); }
619 template<typename Type>
621 inline Type
622 RoundDown(Type x, Type base)
623 {
624  Type remainder = Remainder(x, base);
625  return remainder ? x-remainder : x;
626 }
627 
628 
630 template<typename Type>
631 inline Type
632 IntegerPart(Type x)
633 {
634  return (x > 0 ? RoundDown(x) : RoundUp(x));
635 }
636 
638 template<typename Type>
639 inline Type FractionalPart(Type x) { return Mod(x,Type(1)); }
640 
641 
643 inline int Floor(float x) { return (int)RoundDown(x); }
645 inline int Floor(double x) { return (int)RoundDown(x); }
646 inline int Floor(long double x) { return (int)RoundDown(x); }
648 
649 
651 inline int Ceil(float x) { return (int)RoundUp(x); }
653 inline int Ceil(double x) { return (int)RoundUp(x); }
654 inline int Ceil(long double x) { return (int)RoundUp(x); }
656 
657 
659 template<typename Type>
660 inline Type Chop(Type x, Type delta) { return (Abs(x) < delta ? zeroVal<Type>() : x); }
661 
662 
664 template<typename Type>
665 inline Type
666 Truncate(Type x, unsigned int digits)
667 {
668  Type tenth = Pow(10,digits);
669  return RoundDown(x*tenth+0.5)/tenth;
670 }
671 
672 
674 
675 
677 template<typename Type>
678 inline Type
679 Inv(Type x)
680 {
681  assert(x);
682  return Type(1)/x;
683 }
684 
685 
686 enum Axis {
687  X_AXIS = 0,
688  Y_AXIS = 1,
689  Z_AXIS = 2
690 };
691 
692 // enum values are consistent with their historical mx analogs.
702 };
703 
704 
705 template <typename S, typename T>
706 struct promote {
707  typedef typename boost::numeric::conversion_traits<S, T>::supertype type;
708 };
709 
710 
713 template<typename Vec3T>
714 size_t
715 MinIndex(const Vec3T& v)
716 {
717  static const size_t hashTable[8] = { 2, 1, 9, 1, 2, 9, 0, 0 };//9 is a dummy value
718  const size_t hashKey =
719  ((v[0] < v[1]) << 2) + ((v[0] < v[2]) << 1) + (v[1] < v[2]);// ?*4+?*2+?*1
720  return hashTable[hashKey];
721 }
722 
723 
726 template<typename Vec3T>
727 size_t
728 MaxIndex(const Vec3T& v)
729 {
730  static const size_t hashTable[8] = { 2, 1, 9, 1, 2, 9, 0, 0 };//9 is a dummy value
731  const size_t hashKey =
732  ((v[0] > v[1]) << 2) + ((v[0] > v[2]) << 1) + (v[1] > v[2]);// ?*4+?*2+?*1
733  return hashTable[hashKey];
734 }
735 
736 } // namespace math
737 } // namespace OPENVDB_VERSION_NAME
738 } // namespace openvdb
739 
740 #endif // OPENVDB_MATH_MATH_HAS_BEEN_INCLUDED
741 
742 // Copyright (c) 2012-2013 DreamWorks Animation LLC
743 // All rights reserved. This software is distributed under the
744 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )