mathExt.h

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////
//              _____________  ______________________    __   __  _____
//             /  ________  |  |   ___   ________   /   |  \ /  \   |
//            |  |       |  |_ |  |_ |  |       /  /    \__  |      |
//            |  |  ___  |  || |  || |  |      /  /        | |      |
//            |  | |   \ |  || |  || |  |     /  /      \__/ \__/ __|__
//            |  | |_@  ||  || |  || |  |    /  /          Institute
//            |  |___/  ||  ||_|  || |  |   /  /_____________________
//             \_______/  \______/ | |__|  /___________________________
//                        |  |__|  |
//                         \______/
//                    University of Utah       
//                           2002
///////////////////////////////////////////////////////////////////////////

//mathExt.h
//Joe Kniss
// first shot at reworking mm.h

#ifndef __MATH_EXTENSIONS_DOT_H
#define __MATH_EXTENSIONS_DOT_H

#include <iostream>
#include <limits>
#include "vec3.h"
#include "quat.h"
#include "mat3.h"

#ifdef max
#undef max
#endif
#ifdef min
#undef min
#endif

namespace gutz {

///////////////////////////////////////////////////////////////////////////
// type shorthand
///////////////////////////////////////////////////////////////////////////
/*typedef unsigned char    uchar;
typedef unsigned char      ubyte;
typedef unsigned short     ushort;
typedef unsigned int       uint;
typedef unsigned long      ulong;
#if __win32
//typedef LARGE_INTEGER      llong;
//typedef ULARGE_INTEGER     ullong;
#else // __real_os
typedef long long          llong;
typedef unsigned long long ullong;
#endif // __win32
*/

///////////////////////////////////////////////////////////////////////////
///@name usefull constants
/// these are often missing from math.h
//@{
///////////////////////////////////////////////////////////////////////////
#ifndef M_E
#define M_E        2.71828182845904523536 
#endif
#ifndef M_LOG2E
#define M_LOG2E    1.44269504088896340736 
#endif
#ifndef M_LOG10E
#define M_LOG10E   0.434294481903251827651 
#endif
#ifndef M_LN2
#define M_LN2      0.693147180559945309417 
#endif
#ifndef M_LN10
#define M_LN10     2.30258509299404568402 
#endif
#ifndef M_PI
#define M_PI       3.14159265358979323846 
#endif
#ifndef M_PI_2
#define M_PI_2     1.57079632679489661923 
#endif
#ifndef M_PI_4
#define M_PI_4     0.785398163397448309616 
#endif
#ifndef M_1_PI
#define M_1_PI     0.318309886183790671538 
#endif
#ifndef M_2_PI
#define M_2_PI     0.636619772367581343076 
#endif
#ifndef M_2_SQRTPI
#define M_2_SQRTPI 1.12837916709551257390 
#endif
#ifndef M_SQRT2
#define M_SQRT2    1.41421356237309504880 
#endif
#ifndef M_SQRT1_2
#define M_SQRT1_2  0.707106781186547524401 
#endif
///@}

///////////////////////////////////////////////////////////////////////////
//@name usefull operations & tests
/// template inlines are better than #defines,
/// - the value being used is not a function if it is not supposed to be \n
///   for instance,  #define blah(x) x + x + x + x  \n
///   if "x" is a function or complex expression it gets evaluated 4 times! \n
///   on the other hand template<class T> inline T blah(T x) { return x+x+x+x; } \n
///   has a well defined return type and behaves as though "x" is value not.
/// - over loading makes sense
/// - return value and usage are obvious
///@{

///////////////////////////////////////////////////////////////////////////
/// clamp [0,1]
template<class T>
inline
T clamp(const T &x) 
{ return (x > 0 ? ( x < 1 ? x : 1) : 0) ; }

/// clamp [c,C]
template<class L, class T, class U>
inline
T clamp(const L &c, const T &x, const U &C) 
{ return (x > c ? ( x < C ? x : C) : c) ;}

///////////////////////////////////////////////////////////////////////////
/// gutz min, just min, but don't want confusion
template<class T>
inline
T g_min(const T &x, const T &y) 
{ return (x > y ? y : x); }

/// gutz max, just max, but don't want confusion
template<class T>
inline
T g_max(const T &x, const T &y) 
{ return (x > y ? x : y); }

///////////////////////////////////////////////////////////////////////////
/// gutz abs, just abs, ...
template<class T>
inline
T g_abs(const T &x) 
{ return (x > 0 ? x : -x); }

///////////////////////////////////////////////////////////////////////////
/// gutz sign, ...
template<class T>
inline
double g_sgn(const T &x) 
{ return (x > 0 ? 1.0 : (x < 0 ? -1.0 : 0.0)); }

///////////////////////////////////////////////////////////////////////////
/// conversions 
#define RAD2DEG(x) (180.0/M_PI*(x))
#define DEG2RAD(x) (M_PI/180.0*(x))
#define SGN(x)     ((x)>0 ? 1.0 : ((x)<0 ? -1.0 : 0))
#define SQ(x)      ((x)*(x))

///////////////////////////////////////////////////////////////////////////
/// affine:
/// given intervals [i,I], [o,O] and a value x which may or may not be
/// inside [i,I], return the value y such that y stands in the same
/// relationship to [o,O] that x does with [i,I].  Or:
/// \code
///    y - o         x - i     
///   -------   =   -------
///    O - o         I - i
/// \endcode
/// It is the callers responsibility to make sure I-i and O-o are 
/// both greater than zero
///
/// you need to make sure the types for i, I, and x match and o and O match
/// the type of o,O will be your return type \n
/// TODO, maybe the divide should happen at "double" precision? I chose
///  float for speed :)
template<class TI,         /// input type 
         class TO>         /// output type
inline
TO g_affine(const TI &i, const TI &x, const TI &I, const TO &o, const TO &O)
{
   return TO((O - o)*(x - i)/float((I - i)) + o);
}

///////////////////////////////////////////////////////////////////////////
/// standard random float range [0,1]
#ifndef DRAND48
#define DRAND48 (rand()/(double)RAND_MAX)
#endif

///////////////////////////////////////////////////////////////////////////
/// generic 2D line intersect
template< class T >
inline
T intersect2D(const T p1x, const T p1y, const T d1x, const T d1y,
              const T p2x, const T p2y, const T d2x, const T d2y)
{
   /// [d1x, d1y] dot Perp([d2x,d2y])
   const T div = d1x * -d2y + d1y * d2x;
   if(!div) /// [d1x,d1y] and [d2x,d2y] are parallel
      return std::numeric_limits<T>::infinity();
   /// t = ( p2 dot(perp(d2)) - p1 dot(perp(d2)) )/ div
   return (p2x * -d2y + p2y * d2x - p1x * -d2y - p1y * d2x)/div;
   
   /// here is the original implementation (easier to understand)
   //const vec2<T> od(-r.d.y, r.d.x); /// ortho (perpendicular) to r.d
   //const T div = d.dot(od);  /// in direction dot r.d-perp
   //if(!div) /// in direction parallel to r.d-perp
   //   return std::numeric_limits<T>::infinity(); 
   //return (r.p.dot(od) - p.dot(od))/div; 
}

///////////////////////////////////////////////////////////////////////////
/// projection (vector-vector)
///  project p1 ONTO p2
template< class T >
inline
vec3<T> project(const vec3<T> &p1, const vec3<T> &p2)
{
   return p2/p2.dot(p2) * p2.dot(p1);
}

///////////////////////////////////////////////////////////////////////////
/// ray/plane intersection.
///  returns infinity if no hit
template< class T >
inline
T intersectRayPlane(const vec3<T> &rpos, const vec3<T> &rdir,
                    const vec3<T> &ppos, const vec3<T> &pnorm)
{
   const T rdDpn = rdir.dot(pnorm);
   if(!rdDpn) return std::numeric_limits<T>::max();
   return ( ppos.dot(pnorm) - rpos.dot(pnorm) )/rdDpn;
}

///////////////////////////////////////////////////////////////////////////
/// ray/sphere intersection.
///   returns infinity if no hit, otherwise returns a vec2 containing the
///   min and max t's for the intersection in that order
template< class T >
inline
vec2<T> intersectRaySphere(const vec3<T> &rpos, const vec3<T> &rdir,
                           const vec3<T> &cent, const T r)
{
   vec3<T> x = rpos - cent;
   T a = rdir.dot(rdir);
   T b = rdir.dot(x) * 2;
   T c = x.dot(x) - r*r;
   T d = b*b - 4*a*c;
   if( d <= 0 ) return vec2<T>(std::numeric_limits<T>::max());
   return vec2<T>((-b-sqrt(d))/(2*a), (-b+sqrt(d))/(2*a));
}

///////////////////////////////////////////////////////////////////////////
/// get the rotation that takes one axis to another,
/// vectors are assumed to be PreNormalized
template< class T >
mat3<T> axis2axisPN(const vec3<T> &v1, const vec3<T> &v2)
{
   vec3<T> axis = -v1.cross(v2);
   float theta = v1.dot(v2);
   if(theta > 0)
   {
      theta = T(M_PI - 2 * asin( (-v1-v2).norm()/2.0 ));
   }
   else
   {
      theta = T(2 * asin( (v1-v2).norm()/2.0 ));
   }
   return mat3<T>(theta,axis);
}
template< class T >
mat3<T> axis2axis(const vec3<T> &v1, const vec3<T> &v2)
{
   return axis2axisPN(v1/v1.norm(), v2/v2.norm());
}

///@}

///////////////////////////////////////////////////////////////
///@name generic policies
///@{

/// deep copies + delete data
template<class T>
struct OwnAllocPolicy {
   /// true = failure, false = success
   bool cpyVec(T *v, const T *c, int num) const
   {
      if(v) delVec(v);
      v = new T[num];
      if(!c) return true;
      for(int i=0; i<num; ++i)
         v[i] = c[i];
      return false;
   }
   void delVec(T *v) const { if(v) delete[] v; }
   void delVal(T *v) const { if(v) delete v; }
};

/// alias copies (shallow) + DON'T delete data
template<class T>
struct WrapAllocPolicy {
   bool cpyVec(T *v, T *c, int num) const
   {
      v = c;
      if(!c) return true;
      return false;
   }
   void delVec(T *v) const {}
   void delVal(T *v) const {}
};

///@}
/////////////////////////////////////////////////////////////////

} //end namespace gutz


#endif

---

Send questions, comments, and bug reports to: