gfx_math.hh

Go to the documentation of this file.

/*
 * LEGAL:   COPYRIGHT (C) 2004 JIM E. BROOKS
 *          THIS SOURCE CODE IS RELEASED UNDER THE TERMS
 *          OF THE GNU GENERAL PUBLIC LICENSE VERSION 2 (GPL 2).
 *****************************************************************************/

#ifndef GFX_MATH_HH
#define GFX_MATH_HH 1

namespace gfx {


// A NaN never equals anything -- not even another NaN.
// C99 isfinite() is true if neither NaN nor infinity.
//#define IFNAN(N) (!((N)==(N))) // (compiler might optimize this incorrectly)
#if __GNUC__  &&  defined(__USE_ISOC99)
#define IF_NAN(N) (UX(!isfinite(N)))
#else
// Tested on gcc 3.4.3 x86.
//template<typename FP> bool IFNAN( volatile FP& n ) { return !(n==n); }
INLINE bool IF_NAN( volatile fp n ) { return !(n==n); }
#endif


// fpx is used in trig functions since they're used in matrix functions
// which need to maintain precision.

// 180' = PI
// 360' = PI * 2
const fpx PI    = 3.14159265358979323846L;   // L for long double
const fpx PI2   = 6.28318530717958647692L;   // pi * 2.0
const fpx PI_90 = 1.57079632679489661926L;   // pi / 2.0 [equivalent to 90']

// Turn a float type into a distinct type
// to ensure a degree/radian aren't confused.
#if SPEED
DECLARE_DISTINCT_TYPE( Radian, fp )
DECLARE_DISTINCT_TYPE( Degree, fp )
#else
DECLARE_DISTINCT_TYPE( Radian, fpx )
DECLARE_DISTINCT_TYPE( Degree, fpx )
#endif

const Radian RADIAN_180 = gfx::PI;   // 180' in radians (half circle)
const Radian RADIAN_360 = gfx::PI2;  // 360' in radians (full circle)

// Compare degrees/radians while ignoring sign.
#define IF_DEG_GT( DEG1, DEG2 ) (ABS(DEG1) > ABS(DEG2))
#define IF_DEG_LT( DEG1, DEG2 ) (ABS(DEG1) < ABS(DEG2))
#define IF_RAD_GT( RAD1, RAD2 ) (ABS(RAD1) > ABS(RAD2))
#define IF_RAD_LT( RAD1, RAD2 ) (ABS(RAD1) < ABS(RAD2))

template<typename FP>
pair<FP,FP>
SinCos( const Radian rad )
{
#if defined(__GNUC__) && defined(CPU_X86)
    pair<FP,FP> si_co;
    asm("fsincos" : "=t" (si_co.second), "=u" (si_co.first) : "0" (rad));
    return si_co;
#elif defined(__GNUC__) && _GLIBCXX_HAVE_SINCOS
    pair<double,double> si_co;  // must be specific in this case
    sincos( rad, &si_co.first, &si_co.second );  // GNU extension
    return make_pair<FP,FP>( si_co.first, si_co.second );
#else
    pair<FP,FP> si_co;
    si_co.first  = sin( rad );
    si_co.second = cos( rad );
    return si_co;
#endif
}

INLINE Degree
TruncateDeg( Degree deg )
{
    if ( EX(deg < 360.0) )
        return deg;
    else
        return deg - 360.0;
}

INTERN const fpx degRadRatio = (1.0 / 360.0) * (gfx::PI2);
#define DEG2RAD( DEG ) ((DEG) * (degRadRatio))

INLINE Radian
Deg2Rad( Degree deg )
{
    return DEG2RAD( deg );
    //return deg / 360.0 * (2.0 * gfx::PI); // x87 FDIV is slower than FMUL
}

INTERN const fpx radDegRatio = (1.0 / gfx::PI2) * 360.0;
#define RAD2DEG( RAD ) ((RAD) * (radDegRatio))

INLINE Degree
Rad2Deg( Radian rad )
{
    return RAD2DEG( rad );
    //return rad / (2.0 * gfx::PI) * 360.0; // x87 FDIV is slower than FMUL
}

INLINE fpx
Sin( Degree deg )
{
    return sin( Deg2Rad( deg ) );
}

INLINE fpx
Cos( Degree deg )
{
    return cos( Deg2Rad( deg ) );
}


} // namespace gfx

#include "gfx_math_matrix.hh"

namespace gfx {


INLINE fp
Distance( const Vector2& v )
{
    return sqrtf( v.x * v.x
                + v.y * v.y );
}

INLINE fp
Distance( const Vector3& v )
{
    return sqrtf( v.x * v.x
                + v.y * v.y
                + v.z * v.z );
}

template<typename VECTOR3A,typename VECTOR3B>
fp
Distance( const VECTOR3A& v1, const VECTOR3B& v2 )
{
    const fp dx = v2.x - v1.x;
    const fp dy = v2.y - v1.y;
    const fp dz = v2.z - v1.z;
    return sqrtf( dx * dx
                + dy * dy
                + dz * dz );
}

INLINE fpx
DistanceSquared( const Vector3& v1, const Vector3& v2 )
{
    fpx xd = v2.x - v1.x;  // extra precision
    fpx yd = v2.y - v1.y;
    fpx zd = v2.z - v1.z;
    return xd*xd + yd*yd + zd*zd;
}

// Distance^2 between two matrixs.
// Don't pass different kinds of matrixs (eg Eye origin is negative but Dyna is positive).
template<typename MATRIX1, typename MATRIX2>
fpx
DistanceSquaredMatrix( const MATRIX1& m1, const MATRIX2& m2 )
{
    Vector3 v1( m1[Ox], m1[Oy], m1[Oz] );
    Vector3 v2( m2[Ox], m2[Oy], m2[Oz] );
    return DistanceSquared( v1, v2 );
}

INLINE Vector3
Normalize( const Vector3& v )
{
    fp dist = Distance( v );
    return Vector3( v.x / dist,
                    v.y / dist,
                    v.z / dist );
}

INLINE fp
NormalizeCoord( uint axis, const Vector3& v )
{
    fp dist = Distance( v );
         if ( axis == XX ) return v.x / dist;
    else if ( axis == YY ) return v.y / dist;
    else                   return v.z / dist;
}

INLINE Vector3
CrossProduct( const Vector3& v0, // origin of v1 and v2
              const Vector3& v1,
              const Vector3& v2 )
{
    const fp xd1 = v1.x - v0.x;
    const fp yd1 = v1.y - v0.y;
    const fp zd1 = v1.z - v0.z;
    const fp xd2 = v2.x - v0.x;
    const fp yd2 = v2.y - v0.y;
    const fp zd2 = v2.z - v0.z;
    return Vector3( yd1 * zd2 - yd2 * zd1,    // x: y*z
                    zd1 * xd2 - zd2 * xd1,    // y: z*x
                    xd1 * yd2 - xd2 * yd1 );  // z: x*y
}

INLINE fp
DotProduct2( const Vector2& v1, const Vector2& v2 )
{
    return v1.x * v2.x
         + v1.y * v2.y;
}

template<typename VECTOR3A,typename VECTOR3B>
fp
DotProduct3( const VECTOR3A& v1, const VECTOR3B& v2 )
{
    return v1.x * v2.x
         + v1.y * v2.y
         + v1.z * v2.z;
}

// Extended-precision.
template<typename VECTOR3A,typename VECTOR3B>
fpx
DotProduct3_fpx( const VECTOR3A& v1, const VECTOR3B& v2 )
{
    return fpx(v1.x) * fpx(v2.x)
         + fpx(v1.y) * fpx(v2.y)
         + fpx(v1.z) * fpx(v2.z);
}

INLINE Radian
Angle2( const Vector2& u, const Vector2& v )  // Angle2/3() differ in DotProduct2/3()
{
    // Dot product formula indirectly provides the angle.
    // angle = arccos( DotProduct(u,v) / Distance(u) * Distance(v) )
    fp uv = Distance(u) * Distance(v);
    if ( ABS(uv) < 0.1 ) return 0.0;  // prevent NaN and div by 0
    fp theta = DotProduct2(u,v) / uv;
    theta = Clamp( theta, -1.0, 1.0 ); // acos() won't tolerate 1.00000001

    fp rad = acos( theta );
ASSERT(errno == 0);  // if acos() failed
    return rad;
}

template<typename VECTOR>
Radian
Angle3( const VECTOR& u, const VECTOR& v )
{
    // Dot product formula indirectly provides the angle.
    // angle = arccos( DotProduct(u,v) / Distance(u) * Distance(v) )
    fp uv = Distance(u) * Distance(v);
    if ( ABS(uv) < 0.1 ) return 0.0;  // prevent NaN and div by 0
    fp theta = DotProduct3(u,v) / uv;
    theta = Clamp( theta, -1.0, 1.0 ); // acos() won't tolerate 1.00000001

    fp rad = acos( theta );
ASSERT(errno == 0);  // if acos() failed
    return rad;
}

INLINE Vector3
Interpolate( const Vector3& v1, const Vector3& v2, const fp fraction )
{
ASSERT( fraction >= 0.0 && fraction <= 1.0 );
    return Vector3( v1 + ((v2-v1) * fraction) );
}

template<typename T>
T
Midpoint( T x1, T x2 )
{
//  return x1 + (x2 - x1) / 2;
    return (x1 + x2) / 2;  // faster
}

INLINE Vector3
Midpoint( const Vector3& v1, const Vector3& v2 )
{
    return Vector3( Midpoint(v1.x, v2.x),
                    Midpoint(v1.y, v2.y),
                    Midpoint(v1.z, v2.z) );
}


template<typename T>
T
SQUARE( T x )
{
    return x * x;
}

template<typename T>
T
AVG( T x, T y )
{
    return (x + y) / 2;
}

template<typename T>
T
Range( T val, T lo, T hi )
{
    if ( val < lo ) return lo;
    if ( val > hi ) return hi;
    return val;
}

template<typename T>
T
IfInRange( T a0, T a1, T b0, T b1 )
{
    // Prove the reverse: test if out-of-range.
         if ( a0 > b1 )
        return false;
    else if ( a1 < b0 )
        return false;
    else if ( b0 > a1 )
        return false;
    else if ( b1 < a0 )
        return false;
    else
        return true;  // within range
}

INLINE fp
Remap( fp value, fp rangeOld, fp rangeNew )
{
    return value / rangeOld * rangeNew;
}

INLINE fp
Truncate( fp val, fp interval )
{
    float integral;  // float, not fp, as address is passed to modff()
    modff( val/interval, &integral );
    return integral*interval;
}

} // namespace gfx

#endif // GFX_MATH_HH

Palomino 3D Engine documents generated by doxygen 1.5.3 on Fri Nov 23 11:26:10 2007