vector4d.h

Go to the documentation of this file.

/******************************************************************************
 *
 * MantaFlow fluid solver framework
 * Copyright 2011 Tobias Pfaff, Nils Thuerey 
 *
 * This program is free software, distributed under the terms of the
 * GNU General Public License (GPL) 
 * http://www.gnu.org/licenses
 *
 * 4D vector class
 *
 ******************************************************************************/

#ifndef _VECTOR4D_H
#define _VECTOR4D_H

#include "vectorbase.h"


namespace Manta
{

template<class S>
class Vector4D
{
public:
    inline Vector4D() : x(0),y(0),z(0),t(0) {}
    
    inline Vector4D ( const Vector4D<S> &v ) : x(v.x), y(v.y), z(v.z),t(v.t) {}

    inline Vector4D ( const float * v) : x((S)v[0]), y((S)v[1]), z((S)v[2]), t((S)v[3]) {}

    inline Vector4D ( const double * v) : x((S)v[0]), y((S)v[1]), z((S)v[2]), t((S)v[3]) {}
    
    inline Vector4D ( S v) : x(v), y(v), z(v), t(v) {}
        
    inline Vector4D ( S vx, S vy, S vz, S vw) : x(vx), y(vy), z(vz), t(vw) {}

    // Operators
    
    inline const Vector4D<S>& operator= ( const Vector4D<S>& v ) {
        x = v.x;
        y = v.y;
        z = v.z;
        t = v.t;
        return *this;
    }
    inline const Vector4D<S>& operator= ( S s ) {
        x = y = z = t = s;
        return *this;
    }
    inline const Vector4D<S>& operator+= ( const Vector4D<S>& v ) {
        x += v.x;
        y += v.y;
        z += v.z;
        t += v.t;
        return *this;
    }
    inline const Vector4D<S>& operator+= ( S s ) {
        x += s;
        y += s;
        z += s;
        t += s;
        return *this;
    }
    inline const Vector4D<S>& operator-= ( const Vector4D<S>& v ) {
        x -= v.x;
        y -= v.y;
        z -= v.z;
        t -= v.t;
        return *this;
    }
    inline const Vector4D<S>& operator-= ( S s ) {
        x -= s;
        y -= s;
        z -= s;
        t -= s;
        return *this;
    }
    inline const Vector4D<S>& operator*= ( const Vector4D<S>& v ) {
        x *= v.x;
        y *= v.y;
        z *= v.z;
        t *= v.t;
        return *this;
    }
    inline const Vector4D<S>& operator*= ( S s ) {
        x *= s;
        y *= s;
        z *= s;
        t *= s;
        return *this;
    }
    inline const Vector4D<S>& operator/= ( const Vector4D<S>& v ) {
        x /= v.x;
        y /= v.y;
        z /= v.z;
        t /= v.t;
        return *this;        
    }
    inline const Vector4D<S>& operator/= ( S s ) {        
        x /= s;
        y /= s;
        z /= s;
        t /= s;
        return *this;        
    }
    inline Vector4D<S> operator- () const {
        return Vector4D<S> (-x, -y, -z, -t);        
    }
    
    //inline S min() const { return ( x<y ) ? ( ( x<z ) ? x:z ) : ( ( y<z ) ? y:z ); }
    //inline S max() const { return ( x>y ) ? ( ( x>z ) ? x:z ) : ( ( y>z ) ? y:z ); }
        
    inline bool empty() {
        return x==0 && y==0 && z==0 && t==0;
    }

    inline S& operator[] ( unsigned int i ) { 
        return value[i];
    }
    inline const S& operator[] ( unsigned int i ) const {
        return value[i];
    }

    std::string toString() const;
    
    bool isValid() const; 

    union {
        S value[4];
        struct {
            S x;
            S y;
            S z;
            S t;
        };
        struct {
            S X;
            S Y;
            S Z;
            S T;
        };
    };

    // zero element
    static const Vector4D<S> Zero, Invalid;
    
protected:

};

//************************************************************************
// Additional operators
//************************************************************************

template<class S> 
inline Vector4D<S> operator+ ( const Vector4D<S> &v1, const Vector4D<S> &v2 ) {
    return Vector4D<S> ( v1.x+v2.x, v1.y+v2.y, v1.z+v2.z, v1.t+v2.t );
}
template<class S, class S2> 
inline Vector4D<S> operator+ ( const Vector4D<S>& v, S2 s ) {
    return Vector4D<S> ( v.x+s, v.y+s, v.z+s, v.t+s );
}
template<class S, class S2> 
inline Vector4D<S> operator+ ( S2 s, const Vector4D<S>& v ) {
    return Vector4D<S> ( v.x+s, v.y+s, v.z+s, v.t+s );
}

template<class S> 
inline Vector4D<S> operator- ( const Vector4D<S> &v1, const Vector4D<S> &v2 ) {
    return Vector4D<S> ( v1.x-v2.x, v1.y-v2.y, v1.z-v2.z, v1.t-v2.t );
}
template<class S, class S2> 
inline Vector4D<S> operator- ( const Vector4D<S>& v, S2 s ) {
    return Vector4D<S> ( v.x-s, v.y-s, v.z-s, v.t-s );
}
template<class S, class S2>
inline Vector4D<S> operator- ( S2 s, const Vector4D<S>& v ) {
    return Vector4D<S> ( s-v.x, s-v.y, s-v.z, s-v.t );
}

template<class S>
inline Vector4D<S> operator* ( const Vector4D<S> &v1, const Vector4D<S> &v2 ) {
    return Vector4D<S> ( v1.x*v2.x, v1.y*v2.y, v1.z*v2.z, v1.t*v2.t );
}
template<class S, class S2>
inline Vector4D<S> operator* ( const Vector4D<S>& v, S2 s ) {
    return Vector4D<S> ( v.x*s, v.y*s, v.z*s , v.t*s );
}
template<class S, class S2> 
inline Vector4D<S> operator* ( S2 s, const Vector4D<S>& v ) {
    return Vector4D<S> ( s*v.x, s*v.y, s*v.z, s*v.t );
}

template<class S>
inline Vector4D<S> operator/ ( const Vector4D<S> &v1, const Vector4D<S> &v2 ) {
    return Vector4D<S> ( v1.x/v2.x, v1.y/v2.y, v1.z/v2.z, v1.t/v2.t );
}
template<class S, class S2>
inline Vector4D<S> operator/ ( const Vector4D<S>& v, S2 s ) {
    return Vector4D<S> ( v.x/s, v.y/s, v.z/s, v.t/s );
}
template<class S, class S2> 
inline Vector4D<S> operator/ ( S2 s, const Vector4D<S>& v ) {
    return Vector4D<S> ( s/v.x, s/v.y, s/v.z, s/v.t );
}

template<class S>
inline bool operator== (const Vector4D<S>& s1, const Vector4D<S>& s2) {
    return s1.x == s2.x && s1.y == s2.y && s1.z == s2.z && s1.t == s2.t;
}

template<class S>
inline bool operator!= (const Vector4D<S>& s1, const Vector4D<S>& s2) {
    return s1.x != s2.x || s1.y != s2.y || s1.z != s2.z || s1.t != s2.t;
}

//************************************************************************
// External functions
//************************************************************************

template<class S>
inline S dot ( const Vector4D<S> &t, const Vector4D<S> &v ) {
    return t.x*v.x + t.y*v.y + t.z*v.z + t.t*v.t;    
}

/*template<class S>
inline Vector4D<S> cross ( const Vector4D<S> &t, const Vector4D<S> &v ) {
    NYI Vector4D<S> cp (
        ( ( t.y*v.z ) - ( t.z*v.y ) ),
        ( ( t.z*v.x ) - ( t.x*v.z ) ),
        ( ( t.x*v.y ) - ( t.y*v.x ) ) );
    return cp;
}*/


template<class S>
inline S norm ( const Vector4D<S>& v ) {
    S l = v.x*v.x + v.y*v.y + v.z*v.z + v.t*v.t;
    return ( fabs ( l-1. ) < VECTOR_EPSILON*VECTOR_EPSILON ) ? 1. : sqrt ( l );
}

template<class S>
inline S normSquare ( const Vector4D<S>& v ) {
    return v.x*v.x + v.y*v.y + v.z*v.z + v.t*v.t;
}

template<class S>
inline Vector4D<S> getNormalized ( const Vector4D<S>& v ) {
    S l = v.x*v.x + v.y*v.y + v.z*v.z + v.t*v.t;
    if ( fabs ( l-1. ) < VECTOR_EPSILON*VECTOR_EPSILON )
        return v; /* normalized "enough"... */
    else if ( l > VECTOR_EPSILON*VECTOR_EPSILON )
    {
        S fac = 1./sqrt ( l );
        return Vector4D<S> ( v.x*fac, v.y*fac, v.z*fac , v.t*fac );
    }
    else
        return Vector4D<S> ( ( S ) 0 );
}


template<class S>
inline S normalize ( Vector4D<S> &v ) {
    S norm;
    S l = v.x*v.x + v.y*v.y + v.z*v.z + v.t*v.t;
    if ( fabs ( l-1. ) < VECTOR_EPSILON*VECTOR_EPSILON ) {
        norm = 1.;
    } else if ( l > VECTOR_EPSILON*VECTOR_EPSILON ) {
        norm = sqrt ( l );
        v *= 1./norm;        
    } else {
        v = Vector4D<S>::Zero;
        norm = 0.;
    }
    return ( S ) norm;
}


template<class S> std::string Vector4D<S>::toString() const {
    char buf[256];
    snprintf ( buf,256,"[%+4.6f,%+4.6f,%+4.6f,%+4.6f]", ( double ) ( *this ) [0], ( double ) ( *this ) [1], ( double ) ( *this ) [2] , ( double ) ( *this ) [3] );
    // for debugging, optionally increase precision:
    //snprintf ( buf,256,"[%+4.16f,%+4.16f,%+4.16f,%+4.16f]", ( double ) ( *this ) [0], ( double ) ( *this ) [1], ( double ) ( *this ) [2], ( double ) ( *this ) [3] );
    return std::string ( buf );
}

template<> std::string Vector4D<int>::toString() const;


template<class S>
std::ostream& operator<< ( std::ostream& os, const Vector4D<S>& i ) {
    os << i.toString();
    return os;
}

template<class S>
std::istream& operator>> ( std::istream& is, Vector4D<S>& i ) {
    char c;
    char dummy[4];
    is >> c >> i[0] >> dummy >> i[1] >> dummy >> i[2] >> dummy >> i[3] >> c;
    return is;
}

/**************************************************************************/
// Define default vector alias
/**************************************************************************/

typedef Vector4D<Real>  Vec4;

typedef Vector4D<int>  Vec4i;

template<class T> inline Vec4 toVec4 ( T v ) {
    return Vec4 ( v[0],v[1],v[2],v[3] );
}
template<class T> inline Vec4i toVec4i ( T v ) {
    return Vec4i( v[0],v[1],v[2],v[3] );
}


/**************************************************************************/
// Specializations for common math functions
/**************************************************************************/

template<> inline Vec4 clamp<Vec4>(const Vec4& a, const Vec4& b, const Vec4& c) {
    return Vec4 ( clamp(a.x, b.x, c.x),
                  clamp(a.y, b.y, c.y),
                  clamp(a.z, b.z, c.z),
                  clamp(a.t, b.t, c.t) );    
}
template<> inline Vec4 safeDivide<Vec4>(const Vec4 &a, const Vec4& b) { 
    return Vec4(safeDivide(a.x,b.x), safeDivide(a.y,b.y), safeDivide(a.z,b.z), safeDivide(a.t,b.t));
}
template<> inline Vec4 nmod<Vec4>(const Vec4& a, const Vec4& b) {
    return Vec4(nmod(a.x,b.x),nmod(a.y,b.y),nmod(a.z,b.z) ,nmod(a.t,b.t) );
}

    
/**************************************************************************/
// 4d interpolation (note only 4d here, 2d/3d interpolations are in interpol.h)
/**************************************************************************/

#define BUILD_INDEX_4D \
    Real px=pos.x-0.5f, py=pos.y-0.5f, pz=pos.z-0.5f, pt=pos.t-0.5f; \
    int xi = (int)px; \
    int yi = (int)py; \
    int zi = (int)pz; \
    int ti = (int)pt; \
    Real s1 = px-(Real)xi, s0 = 1.-s1; \
    Real t1 = py-(Real)yi, t0 = 1.-t1; \
    Real f1 = pz-(Real)zi, f0 = 1.-f1; \
    Real g1 = pt-(Real)ti, g0 = 1.-g1; \
    /* clamp to border */ \
    if (px < 0.) { xi = 0; s0 = 1.0; s1 = 0.0; } \
    if (py < 0.) { yi = 0; t0 = 1.0; t1 = 0.0; } \
    if (pz < 0.) { zi = 0; f0 = 1.0; f1 = 0.0; } \
    if (pt < 0.) { ti = 0; g0 = 1.0; g1 = 0.0; } \
    if (xi >= size.x-1) { xi = size.x-2; s0 = 0.0; s1 = 1.0; } \
    if (yi >= size.y-1) { yi = size.y-2; t0 = 0.0; t1 = 1.0; } \
    if (zi >= size.z-1) { zi = size.z-2; f0 = 0.0; f1 = 1.0; }  \
    if (ti >= size.t-1) { ti = size.t-2; g0 = 0.0; g1 = 1.0; }  \
    const int sX = 1; \
    const int sY = size.x;    

static inline void checkIndexInterpol4d(const Vec4i& size, int idx) {
    if (idx<0 || idx > size.x * size.y * size.z * size.t) {
        std::ostringstream s;
        s << "Grid interpol4d dim " << size << " : index " << idx << " out of bound ";
        errMsg(s.str());
    }
}
        
template <class T>
inline T interpol4d(const T* data, const Vec4i& size, const IndexInt sZ, const IndexInt sT, const Vec4& pos) {
    BUILD_INDEX_4D
    IndexInt idx = (IndexInt)xi + sY * (IndexInt)yi + sZ * (IndexInt)zi + sT * (IndexInt)ti;
    DEBUG_ONLY(checkIndexInterpol4d(size,idx)); 
    DEBUG_ONLY(checkIndexInterpol4d(size,idx+sX+sY+sZ+sT));
    
    return( ((data[idx]      *t0 + data[idx+sY]      *t1) * s0
           + (data[idx+sX]   *t0 + data[idx+sX+sY]   *t1) * s1) * f0
           +((data[idx+sZ]   *t0 + data[idx+sY+sZ]   *t1) * s0
           + (data[idx+sX+sZ]*t0 + data[idx+sX+sY+sZ]*t1) * s1) * f1 ) * g0 
        +
          ( ((data[idx+sT]      *t0 + data[idx+sT+sY]      *t1) * s0
           + (data[idx+sT+sX]   *t0 + data[idx+sT+sX+sY]   *t1) * s1) * f0
           +((data[idx+sT+sZ]   *t0 + data[idx+sT+sY+sZ]   *t1) * s0
           + (data[idx+sT+sX+sZ]*t0 + data[idx+sT+sX+sY+sZ]*t1) * s1) * f1 ) * g1 ;
}

}; // namespace


#endif