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