RenderTypes.h

Go to the documentation of this file.

/*
 * Eidolon Biomedical Framework
 * Copyright (C) 2016-8 Eric Kerfoot, King's College London, all rights reserved
 * 
 * This file is part of Eidolon.
 *
 * Eidolon is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * Eidolon is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License along
 * with this program (LICENSE.txt).  If not, see <http://www.gnu.org/licenses/>
 */

#ifndef RENDERTYPES_H
#define RENDERTYPES_H

#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <cerrno>
#include <string>
#include <fstream>
#include <sstream>
#include <iostream>
#include <cstring>
#include <cmath>
#include <vector>
#include <map>
#include <algorithm>
#include <utility>
#include <limits>
#include <stdexcept>

#ifdef WIN32
  typedef void* _WId;
  #define WIN32_LEAN_AND_MEAN
  #include <windows.h>
  #include <tchar.h>
  #ifdef RENDER_EXPORT
  #  define DLLEXPORT __declspec( dllexport )
  #else
  #  define DLLEXPORT __declspec( dllimport )
  #endif
  #define PLATFORM_ID "Windows"
  #pragma warning(disable:4244) // disable double->float conversion complaint
  #pragma warning(disable:4250) // disable inheritance via dominance complaint
  #pragma warning(disable:4996) // disable strerror warning
  #pragma warning(disable:4100) // disable unreferenced formal parameter warning
#elif defined(__APPLE__)
  #include <unistd.h>
  #include <pthread.h>
  #include <sys/mman.h>
  #include <sys/stat.h>
  #include <sys/time.h>
  #include <sys/posix_shm.h>
  #include <fcntl.h>
  #define MAXSHMNAMLEN PSHMNAMLEN
  #define DLLEXPORT
  #define PLATFORM_ID "OSX"
  typedef long _WId;
#else
  #include <unistd.h>
  #include <pthread.h>
  #include <sys/mman.h>
  #include <sys/stat.h>
  #include <fcntl.h>
  #include <limits.h>
  #define MAXSHMNAMLEN NAME_MAX
  #define DLLEXPORT
  #define PLATFORM_ID "Linux"
  typedef unsigned long _WId;
#endif

#define dPI 3.141592653589793238462
#define fPI 3.141592f
#define fEPSILON 1.0e-10f
#define dEPSILON 1.0e-10

#define _HASH(h,v,s) (((h)<<(s)|((h)>>((sizeof(h)<<3)-(s))))^(v))

#define DBGOUT(c) do { std::cout << c << std::endl; std::cout.flush(); } while(0)

#define SAFE_DELETE(p) do { if((p)!=NULL){ delete (p); (p)=NULL; } } while(0)

namespace RenderTypes {

// various string names
static const char* platformID=PLATFORM_ID;
static const char* parentPIDVar="PARENTPID";
static const char* RenderParamGroup="RenderParam";

// platform-independent basic type definitions
typedef int i32;
typedef long long i64;
typedef unsigned char u8;
typedef unsigned int u32;
typedef unsigned long long u64;

// specific type definitions which are meaningful for data structures and the rendering systems
typedef u32 sval; // size value, fixed at 32bits even on 64bit platforms
typedef double real; // real value data type for internal code and file formats
typedef u32 rgba; // 32bit color data type
typedef u32 indexval; // index value data type for internal code and file formats, may differ from sval in future if larger indices are needed

static const real realInf=std::numeric_limits<real>::infinity();

template<typename T> T _min(const T& a, const T& b) { return a<b ? a : b; }
template<typename T> T _max(const T& a, const T& b) { return a>b ? a : b; }

template<typename T> T clamp(const T& val, const T& minval, const T& maxval)
{
        if(val>maxval)
                return maxval;
        if(val<minval)
                return minval;
        
        return val;
}

template<typename T> T lerpXi(const T& val, const T& minv, const T& maxv)
{
        return minv==maxv ? val : (val-minv)/(maxv-minv);
}

template<typename V,typename T> T lerp(const V& val, const T& v1, const T& v2)
{
        return v1+(v2-v1)*val;
}

template<typename T> int compT( const T& t1, const T& t2)
{
        if(t2<t1)
                return 1;
        if(t1<t2)
                return -1;

        return 0;
}

template<typename T> int compV(const void* t1, const void* t2)
{
        return compT<T>(*((T*)t1),*((T*)t2));
}

template<typename T> int sortTupleFirstCB(const void* v1, const void* v2)
{
        return compT(((T*)v1)->first,((T*)v2)->first);
}

template<typename T> int sortTupleSecondCB(const void* v1, const void* v2)
{
        return compT(((T*)v1)->second,((T*)v2)->second);
}

template<typename T> int sortTupleThirdCB(const void* v1, const void* v2)
{
        return compT(((T*)v1)->third,((T*)v2)->third);
}

template<typename T> int sortTupleFourthCB(const void* v1, const void* v2)
{
        return compT(((T*)v1)->fourth,((T*)v2)->fourth);
}

inline bool equalsEpsilon(real v1, real v2)
{
        return fabs(v1-v2)<=dEPSILON;
}

inline bool isNan(real v)
{
        volatile real vv=v; // must be volatile to prevent optimizations
        return vv!=vv; // NaN is the only value not equal to itself
}

inline real frand()
{
        return real(rand())/real(RAND_MAX);
}

inline real fround(real r)
{
        return floor(0.5+r);
}

inline std::string getPIDStr()
{
#ifdef WIN32
        DWORD self=GetCurrentProcessId();
#else
        pid_t self=getpid();
#endif

        std::ostringstream out;
        out << self;
        return out.str();
}

inline std::string getPPIDStr()
{
#ifdef WIN32
        DWORD self=0;
#else
        pid_t self=getppid();
#endif
        std::ostringstream out;
        out << self;
        return out.str();
}

inline bool isParentProc()
{
        const char* parentpid=getenv(parentPIDVar);
        return parentpid==NULL || std::string(parentpid)==getPIDStr();
}

template<typename T> T swapEndianN(T t)
{
        union {T v; u8 b[sizeof(T)]; } src,dst;
        src.v=t;
        for(size_t x=0;x<sizeof(T);x++)
                dst.b[x]=src.b[sizeof(T)-x-1];
        return dst.v;
}

template<typename T> T swapEndian32(T t)
{
        union {T v; u8 b[4]; } src,dst;
        src.v=t;
        dst.b[0]=src.b[3];
        dst.b[1]=src.b[2];
        dst.b[2]=src.b[1];
        dst.b[3]=src.b[0];
        return dst.v;
}

template<typename T> T swapEndian64(T t)
{
        union {T v; u8 b[8]; } src,dst;
        src.v=t;
        dst.b[0]=src.b[7];
        dst.b[1]=src.b[6];
        dst.b[2]=src.b[5];
        dst.b[3]=src.b[4];
        dst.b[4]=src.b[3];
        dst.b[5]=src.b[2];
        dst.b[6]=src.b[1];
        dst.b[7]=src.b[0];
        return dst.v;
}

template<typename F,typename S,typename T>
class triple
{
public:
                F first;
                S second;
                T third;

                triple() : first(),second(),third() {}
                triple(const F& first,const S& second, const T& third) : first(first),second(second), third(third) {}
                triple(const triple<F,S,T>& t) : first(t.first),second(t.second),third(t.third) {}
};

template<typename F,typename S,typename T,typename U>
class quadruple
{
public:
                F first;
                S second;
                T third;
                U fourth;

                quadruple() : first(),second(),third(), fourth() {}
                quadruple(const F& first,const S& second, const T& third, const U& fourth) : first(first),second(second), third(third),fourth(fourth) {}
                quadruple(const quadruple<F,S,T,U>& t) : first(t.first),second(t.second),third(t.third),fourth(t.fourth) {}
};

typedef std::pair<indexval,indexval> indexpair;
typedef std::pair<real,real> realpair;
typedef triple<real,real,real> realtriple;
typedef std::pair<indexval,realtriple> indextriple;
typedef triple<sval,sval,real> intersect;

template<typename T>
void bswap(T& a, T& b)
{
        u8 tmp,*aa=(u8*)&a,*bb=(u8*)&b;
        for(size_t i=0;i<sizeof(T);i++){
                tmp=aa[i];
                aa[i]=bb[i];
                bb[i]=tmp;
        }
}

class TimingObject
{
public:
        clock_t start,stop;
        double delta;
        bool doPrint;
        bool entered;
        std::string label;

        TimingObject(const std::string& label, bool doPrint=true) : delta(0),doPrint(doPrint), entered(false), label(label)
        { 
#ifdef WIN32
                this->label=label.substr(label.find_last_of("::")+1);
#endif

                start=clock(); 
                if(doPrint){
                        std::cout << this->label << std::endl;
                        std::cout.flush();
                }
        }

        ~TimingObject(){ stopTiming(); }

        void stopTiming()
        {
                stop=clock();
                delta=double(stop-start)/CLOCKS_PER_SEC;
                if(doPrint){
                        std::cout << label << " dT (s) = " << delta << std::endl;
                        std::cout.flush();
                }
        }

        bool loopOnce()
        {
                if(entered)
                        return false;

                entered=true;
                return true;
        }
};

#define TIMING TimingObject __functimer(__FUNCTION__)
#define TIMINGBLOCK(label) for(TimingObject __blockobj(label);__blockobj.loopOnce();)
 
#ifdef WIN32
  #define MutexType CRITICAL_SECTION
  #define trylock_mutex(m) (TryEnterCriticalSection(m)!=0)
  #define lock_mutex EnterCriticalSection
  #define unlock_mutex LeaveCriticalSection
  #define destroy_mutex DeleteCriticalSection
#else
  #define MutexType pthread_mutex_t
  #define trylock_mutex(m) (pthread_mutex_trylock(m)==0)
  #define lock_mutex pthread_mutex_lock
  #define unlock_mutex pthread_mutex_unlock
  #define destroy_mutex pthread_mutex_destroy
#endif

class Mutex
{
        MutexType _mutex;

public:
        class Locker
        {
        protected:
                Mutex *parent;
                bool runblock;

        public:
                Locker(Mutex* p, real timeout=0.0) : parent(p) { runblock=parent->lock(timeout); }

                ~Locker() { if(parent)parent->release(); }

                bool loopOnce()
                {
                        if(!runblock)
                                return false;

                        runblock=false;
                        return true;
                }
        };

        Mutex()
        {
#ifdef WIN32
                InitializeCriticalSection(&_mutex);
#else
                pthread_mutexattr_t attrs;
                pthread_mutexattr_init(&attrs);
                pthread_mutexattr_settype(&attrs,PTHREAD_MUTEX_RECURSIVE);
                pthread_mutex_init(&_mutex,&attrs);
#endif
        }

        ~Mutex() { destroy_mutex(&_mutex); }

        bool lock(real timeout=0.0)
        {
                if(timeout>0){
                        clock_t start=clock();
                        real delta=0;
                        bool result=false;

                        while(!result && delta<timeout){
                                result=trylock_mutex(&_mutex);
                                delta=real(clock()-start)/CLOCKS_PER_SEC;
                        }

                        return result;
                }
                else{
                        lock_mutex(&_mutex);
                        return true;
                }
        }
        
        void release() { unlock_mutex(&_mutex); }
};

#define critical(m) for(Mutex::Locker __locker__=Mutex::Locker(m);__locker__.loopOnce();)
#define trylock(m,timeout) for(Mutex::Locker __locker__=Mutex::Locker(m,timeout);__locker__.loopOnce();)


enum FigureType
{
        FT_LINELIST      =  0, // list of line segments
        FT_POINTLIST     =  1, // list of discrete points
        FT_TRILIST       =  2, // list of triangles
        FT_TRISTRIP      =  3, // strip of triangles 
        FT_BB_POINT      =  4, // BBT_POINT billboard points
        FT_BB_FIXED_PAR  =  5, // BBT_ORIENTED_SELF billboard oriented along an axis of rotation
        FT_BB_FIXED_PERP =  6, // BBT_PERPENDICULAR_SELF billboard oriented along a normal of rotation
        FT_GLYPH         =  7, // list of points represented by glyph meshes
//      FT_INTERPFIGURE  =  8, // figure which interpolates between two given other figures
        FT_RIBBON        =  8, // list of lines represented by ribbons
        FT_TEXVOLUME     =  9, // 3D texture volume
        FT_TEXT          = 10  // text billboard 
};

enum BlendMode
{
        BM_ALPHA,
        BM_COLOR,
        BM_ADD,
        BM_MOD,
        BM_REPLACE
};

enum TextureFormat
{
        TF_RGB24,     // 24-bit RGB pixels PF_R8G8B8
        TF_RGBA32,    // 32-bit RGBA pixels PF_R8G8B8A8
        TF_ARGB32,    // 32-bit ARGB pixels PF_A8R8G8B8
        TF_LUM8,      // 8-bit greyscale, no alpha PF_L8
        TF_LUM16,     // 16-bit greyscale, no alpha PF_L16
        TF_ALPHA8,    // 8-bit alpha mask, no greyscale PF_A8
        TF_ALPHALUM8, // 4-bit alpha, 4-bit greyscale PF_A4L4
        TF_UNKNOWN    // Any other format from the renderer that's not used/understood
};

enum ProgramType
{
        PT_VERTEX  =0,
        PT_FRAGMENT=1,
        PT_GEOMETRY=2
};

enum HAlignType
{
        H_LEFT,
        H_RIGHT,
        H_CENTER
};

enum VAlignType
{
        V_TOP,
        V_BOTTOM,
        V_CENTER
};



#ifdef WIN32
  std::string formatLastErrorMsg(); // Windows error reporting helper
#endif

// OS X shared memory cleanup stuff

//extern std::vector<std::string> sharednamelist; // list of shared names to delete on exit
//void addShared(const std::string& name);
//void cleanupShared();


void initSharedDir(const std::string& path);
std::string getSharedDir();
void addShared(const std::string& name);
void unlinkShared(const std::string& name);

/*****************************************************************************************************************************/
/* Math Objects */
/*****************************************************************************************************************************/

class color
{
        float _r,_g,_b,_a;

public:
        static rgba interpolate(real val,rgba left, rgba right)
        {
                u32 bf = u32(val*255);
                u32 af = 255 - bf;
                
                u32 al = (left & 0x00ff00ff);
                u32 ah = (left & 0xff00ff00) >> 8;
                u32 bl = (right & 0x00ff00ff);
                u32 bh = (right & 0xff00ff00) >> 8;
                
                u32 ml = (al * af + bl * bf);
                u32 mh = (ah * af + bh * bf);
                
                return (mh & 0xff00ff00) | ((ml & 0xff00ff00) >> 8);
        }
        
        color(float r=1.0f, float g=1.0f, float b=1.0f, float a=1.0f) : _r(r),_g(g),_b(b),_a(a)
        {}

        color(const color & c) : _r(c.r()),_g(c.g()),_b(c.b()),_a(c.a())
        {}
        
        color(const rgba& c) : _r(u8(c>>24)/255.0f), _g(u8(c>>16)/255.0f), _b(u8(c>>8)/255.0f), _a(u8(c)/255.0f)
        {}

        float r() const { return _r; }
        float g() const { return _g; }
        float b() const { return _b; }
        float a() const { return _a; }

        float r(float val) { _r=val; return _r; }
        float g(float val) { _g=val; return _g; }
        float b(float val) { _b=val; return _b; }
        float a(float val) { _a=val; return _a; }
        
        void setBuff(float *v) const { v[0]=_r; v[1]=_g; v[2]=_b; v[3]=_a; }

        rgba toRGBA() const
        {
                rgba result=u8(_r*255);
                result=(result<<8)|u8(_g*255);
                result=(result<<8)|u8(_b*255);
                result=(result<<8)|u8(_a*255);
                return result;
        }

        color interpolate(real val,const color& col) const
        {
                if(val>=1)
                        return col;

                if(val<=0)
                        return *this;

                real val1=1.0f-val;
                return color(_r*val1+col.r()*val,_g*val1+col.g()*val,_b*val1+col.b()*val,_a*val1+col.a()*val);
        }

        color unitClamp()
        {
                return color(clamp(_r,0.0f,1.0f),clamp(_g,0.0f,1.0f),clamp(_b,0.0f,1.0f),clamp(_a,0.0f,1.0f));
        }

        bool operator == (const color & c) const { return equalsEpsilon(r(),c.r()) && equalsEpsilon(g(),c.g()) && equalsEpsilon(b(),c.b()) && equalsEpsilon(a(),c.a()); }
        bool operator != (const color & c) const { return !((*this)==c); }
        
        color operator * (const color & c) const { return color(_r*c.r(),_g*c.g(),_b*c.b(),_a*c.a()); }
        color operator * (real r) const { return color(_r*r,_g*r,_b*r,_a*r); }

        color operator + (const color & c) const { return color(_r+c.r(),_g+c.g(),_b+c.b(),_a+c.a()); }
        color operator + (real r) const { return color(_r+r,_g+r,_b+r,_a+r); }  

        color operator - (const color & c) const { return color(_r-c.r(),_g-c.g(),_b-c.b(),_a-c.a()); }
        color operator - (real r) const { return color(_r-r,_g-r,_b-r,_a-r); }

        bool operator < (const color &c) const { return (_r-dEPSILON)<c.r() && (_g-dEPSILON)<c.g() && (_b-dEPSILON)<c.b() && (_a-dEPSILON)<c.a(); }
        bool operator > (const color &c) const { return (_r+dEPSILON)>c.r() && (_g+dEPSILON)>c.g() && (_b+dEPSILON)>c.b() && (_a+dEPSILON)>c.a(); }

        friend std::ostream& operator << (std::ostream &out, const color &c)
        {
                return out << "color(" << c.r() << ", " << c.g() << ", " << c.b() << ", " << c.a() << ")";
        }
};

class vec3
{
        real _x,_y,_z;

public:
        vec3(real val=0) : _x(val),_y(val),_z(val) {}
        vec3(real x, real y, real z=0) : _x(x),_y(y),_z(z) {}
        vec3(const vec3 & v) : _x(v.x()),_y(v.y()),_z(v.z()) {}

        real x() const { return _x; }
        real y() const { return _y; }
        real z() const { return _z; }

        real x(real v) { _x=v; return _x; }
        real y(real v) { _y=v; return _y; }
        real z(real v) { _z=v; return _z; }
        
        void setBuff(float *v) const { v[0]=_x; v[1]=_y; v[2]=_z; }

        vec3 operator + (const vec3& v) const { return vec3(_x+v.x(),_y+v.y(),_z+v.z()); }
        vec3 operator - (const vec3& v) const { return vec3(_x-v.x(),_y-v.y(),_z-v.z()); }
        vec3 operator * (const vec3& v) const { return vec3(_x*v.x(),_y*v.y(),_z*v.z()); }
        vec3 operator / (const vec3& v) const { return vec3(_x/v.x(),_y/v.y(),_z/v.z()); }
        vec3 operator + (real v) const { return vec3(_x+v,_y+v,_z+v); }
        vec3 operator - (real v) const { return vec3(_x-v,_y-v,_z-v); }
        vec3 operator * (real v) const { return vec3(_x*v,_y*v,_z*v); }
        vec3 operator / (real v) const { return vec3(_x/v,_y/v,_z/v); }
        vec3 operator - () const { return vec3(-_x,-_y,-_z); }
        
        vec3 abs() const { return vec3(fabs(_x),fabs(_y),fabs(_z)); }
        vec3 inv() const { return vec3(_x!=0 ? 1/_x : 0,_y!=0 ? 1/_y : 0,_z!=0 ? 1/_z : 0); }
        vec3 sign() const { return vec3(_x>=0 ? 1 : -1,_y>=0 ? 1 : -1,_z>=0 ? 1 : -1); }
        vec3 cross(const vec3& v) const { return vec3(_y * v.z() - _z * v.y(), _z * v.x() - _x * v.z(), _x * v.y() - _y * v.x()); }
        real dot(const vec3& v) const { return _x*v.x()+_y*v.y()+_z*v.z(); }
        real len() const { return sqrt(_x*_x+_y*_y+_z*_z); }
        real lenSq() const { return _x*_x+_y*_y+_z*_z; }
        vec3 norm() const { real l=len(); return l==0.0 ? vec3() : (*this)*(1.0/l); }
        real distTo(const vec3 & v) const { return ((*this)-v).len(); }
        real distToSq(const vec3 & v) const { return ((*this)-v).lenSq(); }
        vec3 clamp(const vec3& v1,const vec3& v2) const { return vec3(RenderTypes::clamp(_x,v1.x(),v2.x()),RenderTypes::clamp(_y,v1.y(),v2.y()),RenderTypes::clamp(_z,v1.z(),v2.z())); }
        
        void setMinVals(const vec3 &v) { _x=_min(_x,v.x()); _y=_min(_y,v.y()); _z=_min(_z,v.z()); }
        void setMaxVals(const vec3 &v) { _x=_max(_x,v.x()); _y=_max(_y,v.y()); _z=_max(_z,v.z()); }

        void normThis() { real l=len(); if(l>0){ _x/=l;_y/=l;_y/=l;} }

        vec3 toPolar() const { real l=len(); return l==0.0 ? vec3() : vec3(atan2(_y,_x),acos(_z/l),l); }

        vec3 toCylindrical() const { return vec3(atan2(_y,_x),_z,sqrt(_y*_y+_x*_x)); }
        
        vec3 fromPolar() const { return vec3(cos(_x)*sin(_y)*_z,sin(_y)*sin(_x)*_z,cos(_y)*_z); }

        vec3 fromCylindrical() const { return vec3(cos(_x)*_z,sin(_x)*_z,_y); }

        bool isZero() const { return equalsEpsilon(_x,0) && equalsEpsilon(_y,0) && equalsEpsilon(_z,0); }

        bool inAABB(const vec3& minv, const vec3& maxv) const
        {
                //return (_x+dEPSILON)>=minv.x() && (_y+dEPSILON)>=minv.y() && (_z+dEPSILON)>=minv.z()
                //      && (_x-dEPSILON)<=maxv.x() && (_y-dEPSILON)<=maxv.y() && (_z-dEPSILON)<=maxv.z();
                return *this>minv && *this<maxv;
        }

        bool inOBB(const vec3& center, const vec3& hx, const vec3& hy, const vec3& hz) const
        {
                vec3 diff=(*this)-center;

                // perform 3 plane distance checks, return true if the distance from plane (center,hx) is less than hx.len(), etc
                return fabs(hx.dot(diff))<=hx.lenSq() && fabs(hy.dot(diff))<=hy.lenSq() && fabs(hz.dot(diff))<=hz.lenSq();
        } 
        
        bool inSphere(const vec3& center,real radius) const { return distToSq(center)<=(radius*radius+dEPSILON); }

        bool onPlane(const vec3& planept, const vec3& planenorm) const { return equalsEpsilon(planeDist(planept,planenorm),0); }
        
        bool isInUnitCube(real margin=0.0) const { return _x>=-margin && _x<=(1.0+margin) && _y>=-margin && _y<=(1.0+margin) && _z>=-margin && _z<=(1.0+margin); }

        bool isParallel(const vec3 &other) const
        {
                return cross(other).isZero();
        }

        bool operator == (const vec3 & v) const { return equalsEpsilon(_x,v.x()) && equalsEpsilon(_y,v.y()) && equalsEpsilon(_z,v.z()); }
        
        bool operator != (const vec3 &v) const { return !equalsEpsilon(_x,v.x()) || !equalsEpsilon(_y,v.y()) || !equalsEpsilon(_z,v.z()); }
        
        bool operator < (const vec3 &v) const { return (_x-dEPSILON)<v.x() && (_y-dEPSILON)<v.y() && (_z-dEPSILON)<v.z(); }
        bool operator > (const vec3 &v) const { return (_x+dEPSILON)>v.x() && (_y+dEPSILON)>v.y() && (_z+dEPSILON)>v.z(); }
        
        int cmp(const vec3 &v) const
        {
                // the order of these statements defines the sorting order where Z is used to sort first and X last
                if(_z<v.z()) return -1;
                if(_z>v.z()) return  1;
                if(_y<v.y()) return -1;
                if(_y>v.y()) return  1;
                if(_x<v.x()) return -1;
                if(_x>v.x()) return  1;
                return 0;
        }

        real angleTo(const vec3 &v) const
        {
                real l=sqrt(lenSq() * v.lenSq());

                if(l<dEPSILON)
                        return 0.0;

                real vl=dot(v)/l;

                if(vl>=(1.0-dEPSILON))
                        return 0.0;
                
                if(vl<=(-1.0+dEPSILON))
                        return dPI;

                return acos(vl);
        }

        vec3 planeNorm(const vec3& v2, const vec3& v3) const { return (v2-*this).cross(v3-*this).norm(); }
        
        vec3 planeNorm(const vec3& v2, const vec3& v3, const vec3& farv) const
        {
                vec3 norm=planeNorm(v2,v3);
                return norm.angleTo(farv-*this)>=(dPI*0.5) ? norm : -norm;
        }

        real planeDist(const vec3& planept, const vec3& planenorm) const { return planenorm.dot((*this)-planept); }

        vec3 planeProject(const vec3& planept, const vec3& planenorm) const { return (*this)-(planenorm*planeDist(planept,planenorm)); }
        
        int planeOrder(const vec3& planenorm,const vec3& v1,const vec3& v2) const
        {
                real order=(v1-*this).cross(v2-*this).dot(planenorm);
                if(order>0)
                        return 1;
                if(order<0)
                        return -1;
                return 0;
        }

        real triArea(const vec3& b, const vec3& c) const
        {
                vec3 bb=b-*this;
                vec3 cc=c-*this;
                
                return bb.len()*cc.len()*sin(bb.angleTo(cc))*0.5;
        }

        real lineDist(vec3 p1,vec3 p2) const
        {
                vec3 p=p2-p1;
                real pl=p.len();
                if(pl<dEPSILON) // p1==p2 so there's no line
                        return -1;

                if(planeDist(p1,p)<0 || planeDist(p2,-p)<0) // this is outside the cylinder area
                        return -1;

                return p.cross(p1-*this).len()/pl;
        }
        
        vec3 lerp(real val,const vec3& v) const
        {
                return vec3(_x+(v.x()-_x)*val,_y+(v.y()-_y)*val,_z+(v.z()-_z)*val);
        }

        i32 hash() const
        {
                union { real f; i64 i;} x,y,z;
                x.f=_x;
                y.f=_y;
                z.f=_z;
                
                i64 hash=_HASH(x.i,_HASH(y.i,z.i,13),14);
                return i32(hash>>32)^i32(hash);
        }
        
        friend std::ostream& operator << (std::ostream &out, const vec3 &v)
        {
                return out << "vec3(" << v.x() << ", " << v.y() << ", " << v.z() << ")";
        }

        static int compX(const void* v1, const void* v2)
        {
                return compT<real>(((const vec3*)v1)->x(),((const vec3*)v2)->x());
        }

        static int compY(const void* v1, const void* v2)
        {
                return compT<real>(((const vec3*)v1)->y(),((const vec3*)v2)->y());
        }

        static int compZ(const void* v1, const void* v2)
        {
                return compT<real>(((const vec3*)v1)->z(),((const vec3*)v2)->z());
        }
        
        static vec3 posInfinity() { return vec3(realInf); }
        static vec3 negInfinity() { return vec3(-realInf); }

        static vec3 X() { return vec3(1,0,0); }

        static vec3 Y() { return vec3(0,1,0); }

        static vec3 Z() { return vec3(0,0,1); }
};

class mat4
{
public:
        union {
                struct{ 
                        real 
                        m00,m01,m02,m03,
                        m10,m11,m12,m13,
                        m20,m21,m22,m23,
                        m30,m31,m32,m33;
                };
                real m[4][4];
        };

        mat4() { clear(); }
        mat4(const real* mat) { memcpy(m,mat,sizeof(real)*16); }
        mat4(real m00,real m01,real m02,real m03, real m10,real m11,real m12,real m13, real m20,real m21,real m22,real m23, real m30,real m31,real m32,real m33) : 
                m00(m00),m01(m01),m02(m02),m03(m03), m10(m10),m11(m11),m12(m12),m13(m13), m20(m20),m21(m21),m22(m22),m23(m23), m30(m30),m31(m31),m32(m32),m33(m33) {}

        void clear() { memset(m,0,sizeof(real)*16); }
        void ident() { clear(); m00=m11=m22=m33=1.0; }
        
        real* getPointer() const { return (real*)m; } 

        vec3 operator * (const vec3& v) const
        {
                vec3 vv= vec3(m00*v.x() + m01*v.y() + m02*v.z() + m03,m10*v.x() + m11*v.y() + m12*v.z() + m13,m20*v.x() + m21*v.y() + m22*v.z() + m23);
                real d=m30*v.x() + m31*v.y() + m32*v.z() + m33;
                return d==0 ? vec3() : (vv/d);
        }

        friend vec3 operator * (const vec3 &v,const mat4& m)
        {
                return m*v;
        }

        mat4 operator * (const mat4& m) const
        {
                return mat4(
                        m.m00*m00 + m.m10*m01 + m.m20*m02 + m.m30*m03,
                        m.m01*m00 + m.m11*m01 + m.m21*m02 + m.m31*m03,
                        m.m02*m00 + m.m12*m01 + m.m22*m02 + m.m32*m03,
                        m.m03*m00 + m.m13*m01 + m.m23*m02 + m.m33*m03,
                        m.m00*m10 + m.m10*m11 + m.m20*m12 + m.m30*m13,
                        m.m01*m10 + m.m11*m11 + m.m21*m12 + m.m31*m13,
                        m.m02*m10 + m.m12*m11 + m.m22*m12 + m.m32*m13,
                        m.m03*m10 + m.m13*m11 + m.m23*m12 + m.m33*m13,
                        m.m00*m20 + m.m10*m21 + m.m20*m22 + m.m30*m23,
                        m.m01*m20 + m.m11*m21 + m.m21*m22 + m.m31*m23,
                        m.m02*m20 + m.m12*m21 + m.m22*m22 + m.m32*m23,
                        m.m03*m20 + m.m13*m21 + m.m23*m22 + m.m33*m23,
                        m.m00*m30 + m.m10*m31 + m.m20*m32 + m.m30*m33,
                        m.m01*m30 + m.m11*m31 + m.m21*m32 + m.m31*m33,
                        m.m02*m30 + m.m12*m31 + m.m22*m32 + m.m32*m33,
                        m.m03*m30 + m.m13*m31 + m.m23*m32 + m.m33*m33
                );
        }
        
        real determinant() const
        {
                //return m00*m11*m22*m33 - m00*m11*m23*m32 - m00*m12*m21*m33 + m00*m12*m23*m31 + 
                //      m00*m13*m21*m32 - m00*m13*m22*m31 - m01*m10*m22*m33 + m01*m10*m23*m32 + 
                //      m01*m12*m20*m33 - m01*m12*m23*m30 - m01*m13*m20*m32 + m01*m13*m22*m30 + 
                //      m02*m10*m21*m33 - m02*m10*m23*m31 - m02*m11*m20*m33 + m02*m11*m23*m30 + 
                //      m02*m13*m20*m31 - m02*m13*m21*m30 - m03*m10*m21*m32 + m03*m10*m22*m31 + 
                //      m03*m11*m20*m32 - m03*m11*m22*m30 - m03*m12*m20*m31 + m03*m12*m21*m30;
                
                real x0 = m00*m11, x1 = m22*m33, x2 = m00*m12, x3 = m23*m31, x4 = m00*m13, x5 = m21*m32, x6 = m01*m10, 
                        x7 = m23*m32, x8 = m01*m12, x9 = m20*m33, x10 = m01*m13, x11 = m22*m30, x12 = m02*m10, x13 = m21*m33, 
                        x14 = m02*m11, x15 = m23*m30, x16 = m02*m13, x17 = m20*m31, x18 = m03*m10, x19 = m22*m31, 
                        x20 = m03*m11, x21 = m20*m32, x22 = m03*m12, x23 = m21*m30;
                        
                return x0*x1 - x0*x7 - x1*x6 + x10*x11 - x10*x21 - x11*x20 + x12*x13 - x12*x3 - x13*x2 + x14*x15 - x14*x9 - x15*x8 + 
                        x16*x17 - x16*x23 - x17*x22 + x18*x19 - x18*x5 - x19*x4 + x2*x3 + x20*x21 + x22*x23 + x4*x5 + x6*x7 + x8*x9;
        }
        
        mat4 inverse() const
        {
                real s0 =  m00*m11 - m01*m10;
                real s1 =  m00*m12 - m02*m10;
                real s2 =  m00*m13 - m03*m10;
                real s3 =  m01*m12 - m02*m11;
                real s4 =  m01*m13 - m03*m11;
                real s5 =  m02*m13 - m03*m12;
                real c5 =  m22*m33 - m23*m32;
                real c4 =  m21*m33 - m23*m31;
                real c3 =  m21*m32 - m22*m31;
                real c2 =  m20*m33 - m23*m30;
                real c1 =  m20*m32 - m22*m30;
                real c0 =  m20*m31 - m21*m30;
                 
                real invdet = 1.0 / (s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0);
                 
                real b00 = ( c3*m13 - c4*m12 + c5*m11 ) * invdet;
                real b01 = ( -c3*m03 + c4*m02 - c5*m01 ) * invdet;
                real b02 = ( m31*s5 - m32*s4 + m33*s3 ) * invdet;
                real b03 = ( -m21*s5 + m22*s4 - m23*s3 ) * invdet;
                real b10 = ( -c1*m13 + c2*m12 - c5*m10 ) * invdet;
                real b11 = ( c1*m03 - c2*m02 + c5*m00 ) * invdet;
                real b12 = ( -m30*s5 + m32*s2 - m33*s1 ) * invdet;
                real b13 = ( m20*s5 - m22*s2 + m23*s1 ) * invdet;
                real b20 = ( c0*m13 - c2*m11 + c4*m10 ) * invdet;
                real b21 = ( -c0*m03 + c2*m01 - c4*m00 ) * invdet;
                real b22 = ( m30*s4 - m31*s2 + m33*s0 ) * invdet;
                real b23 = ( -m20*s4 + m21*s2 - m23*s0 ) * invdet;
                real b30 = ( -c0*m12 + c1*m11 - c3*m10 ) * invdet;
                real b31 = ( c0*m02 - c1*m01 + c3*m00 ) * invdet;
                real b32 = ( -m30*s3 + m31*s1 - m32*s0 ) * invdet;
                real b33 = ( m20*s3 - m21*s1 + m22*s0 ) * invdet;
        
                return mat4(b00,b01,b02,b03,b10,b11,b12,b13,b20,b21,b22,b23,b30,b31,b32,b33);
        }
};

class rotator
{
protected:
        real _w,_x,_y,_z;

public:
        rotator() : _w(1),_x(0),_y(0),_z(0)
        {}

        rotator(const rotator &r)
        {
                set(r);
        }

        rotator(const vec3& axis, real rads)
        {
                setAxis(axis,rads);
        }

                
        rotator(real x, real y, real z, real w)
        {
                set(x,y,z,w);
        }
        
        rotator(const vec3& from, const vec3& to)
        {
                if(from==to)
                        set(0,0,0,1);
                //else if(from.isParallel(to)){
                //      vec3 axis=from.isParallel(vec3(1,0,0)) ? vec3(0,1,0) : vec3(1,0,0);
                //      setAxis(axis.cross(from),dPI);
                //}
                else
                        setAxis(from.cross(to),from.angleTo(to));
        }

        rotator(real yaw, real pitch, real roll)
        {
                real c1 = cos(0.5*roll); // was yaw in the original (X=right, Y=up, Z=towards) axes definition
                real s1 = sin(0.5*roll);
                real c2 = cos(0.5*yaw); // was pitch
                real s2 = sin(0.5*yaw);
                real c3 = cos(0.5*pitch); // was roll
                real s3 = sin(0.5*pitch);
                real c1c2 = c1*c2;
                real s1s2 = s1*s2;
                real c1s2 = c1*s2;
                real s1c2 = s1*c2;
                
                _w =c1c2*c3 - s1s2*s3;
                _x =c1c2*s3 + s1s2*c3;
                _y =s1c2*c3 + c1s2*s3;
                _z =c1s2*c3 - s1c2*s3;

                // TODO: the above defines ordering of rotations as pitch-yaw-roll, is this correct? This instead perhaps:
                //set(rotator(vec3::Z(),yaw)*rotator(vec3::X(),pitch)*rotator(vec3::Y(),roll));
        }
        
        rotator(real m00,real m01,real m02,real m10,real m11,real m12,real m20,real m21,real m22)
        {
                real tr = m00 + m11 + m22;
                
                if (tr > 0) { 
                        real S = sqrt(tr+1.0) * 2; // S=4*qw 
                        _w = 0.25 * S;
                        _x = (m21 - m12) / S;
                        _y = (m02 - m20) / S; 
                        _z = (m10 - m01) / S; 
                } else if(m00 > m11 && m00 > m22) { 
                        real S = sqrt(1.0 + m00 - m11 - m22) * 2; // S=4*qx 
                        _w = (m21 - m12) / S;
                        _x = 0.25 * S;
                        _y = (m01 + m10) / S; 
                        _z = (m02 + m20) / S; 
                } else if (m11 > m22) { 
                        real S = sqrt(1.0 + m11 - m00 - m22) * 2; // S=4*qy
                        _w = (m02 - m20) / S;
                        _x = (m01 + m10) / S; 
                        _y = 0.25 * S;
                        _z = (m12 + m21) / S; 
                } else { 
                        real S = sqrt(1.0 + m22 - m00 - m11) * 2; // S=4*qz
                        _w = (m10 - m01) / S;
                        _x = (m02 + m20) / S;
                        _y = (m12 + m21) / S;
                        _z = 0.25 * S;
                }
        }
        
        rotator(vec3 row1, vec3 col1, vec3 row2, vec3 col2)
        {
                vec3 norm1=col1.cross(row1).norm(); // first plane normal
                vec3 norm2=col2.cross(row2).norm(); // second plane normal                
                rotator rot;
                
                // define rot as the rotation aligning the second plane to the first
                if(norm1==-norm2)
                        rot=rotator(row1,dPI); // flip along the row vector
                else
                        rot=rotator(norm2,norm1);
                
                // combine rot with a rotation which aligns row vectors by rotating in the first plane, then assign values to this object
                //set(rotator(norm1,row1.angleTo(rot*row2))*rot); // y u no work?
                set(rotator(rot*row2,row1)*rot);
        }

        void setAxis(const vec3& axis, real rads)
        {
                if(!equalsEpsilon(rads,0.0) && !axis.isZero()){
                        vec3 na=axis.norm();
                        real srads = (real)sin(rads / 2.0);
                        set(na.x()*srads,na.y()*srads,na.z()*srads,(real)cos(rads / 2.0));
                }
                else
                        set(0,0,0,1);
        }
        
        void set(const rotator& r)
        {
                set(r._x,r._y,r._z,r._w);
        }

        void set(real ry,real rz, real rw)
        {
                set(sqrt(1.0-(ry*ry+rz*rz+rw*rw)),ry,rz,rw);
        }

        void set(real rx,real ry, real rz,real rw)
        {
                _x = rx;
                _y = ry;
                _z = rz;
                _w = rw;
        }

        real w() const { return _w; }
        real x() const { return _x; }
        real y() const { return _y; }
        real z() const { return _z; }

        real getPitch() const
        {
                real test=_x*_y+_z*_w;
                if(test>(0.5-dEPSILON))
                        return 0;
                
                if(test<(-0.5+dEPSILON))
                        return 0;
                
                return atan2(2*_x*_w-2*_y*_z , 1 - 2*_x*_x - 2*_z*_z);
        }

        real getYaw() const
        {       
                real test=_x*_y+_z*_w;
                if(test>(0.5-dEPSILON))
                        return dPI*0.5;
                
                if(test<(-0.5+dEPSILON))
                        return dPI*-0.5;
                
                return asin(2*test); 
        }

        real getRoll() const
        {
                real test=_x*_y+_z*_w;
                if(test>(0.5-dEPSILON))
                        return 2*atan2(_x,_w);
                
                if(test<(-0.5+dEPSILON))
                        return -2*atan2(_x,_w);
                
                return atan2(2*_y*_w-2*_x*_z , 1 - 2*_y*_y - 2*_z*_z);
        }

        vec3 operator * (const vec3 &v) const
        {
                vec3 axis(_x,_y,_z);
                vec3 vc=axis.cross(v);
                vec3 vcc=axis.cross(vc);
                return (vc*(2.0*_w))+(vcc*2.0)+v;
        }
        
        friend vec3 operator * (const vec3 &v,const rotator& r)
        {
                return r*v;
        }

        vec3 operator / (const vec3& v) const
        {
                return inverse()*v;
        }

        friend vec3 operator / (const vec3 &v,const rotator& r)
        {
                return r/v;
        }

        rotator operator * (const rotator & r) const
        {
                rotator rr;
                rr.set( _w * r._x + _x * r._w + _y * r._z - _z * r._y,
                        _w * r._y + _y * r._w + _z * r._x - _x * r._z,
                        _w * r._z + _z * r._w + _x * r._y - _y * r._x,
                        _w * r._w - _x * r._x - _y * r._y - _z * r._z);
                return rr;
        }
        
        rotator operator * (real r) const 
        {
                rotator rr;
                rr.set(_x*r,_y*r,_z*r,_w*r);
                return rr;
        }

        rotator operator + (const rotator & r) const
        {
                rotator rr;
                rr.set(_x+r._x,_y+r._y,_z+r._z,_w+r._w);
                return rr;
        }
        
        rotator operator - (const rotator & r) const
        {
                rotator rr;
                rr.set(_x-r._x,_y-r._y,_z-r._z,_w-r._w);
                return rr;
        }
        
        rotator operator - () const
        {
                rotator rr;
                rr.set(-_x,-_y,-_z,-_w);
                return rr;
        }

        bool operator == (const rotator & v) const
        {
                if(equalsEpsilon(_w,v._w) && equalsEpsilon(_x,v._x) && equalsEpsilon(_y,v._y) && equalsEpsilon(_z,v._z))
                        return true;

                if(equalsEpsilon(-_w,v._w) && equalsEpsilon(-_x,v._x) && equalsEpsilon(-_y,v._y) && equalsEpsilon(-_z,v._z))
                        return true;

                return false;
        }
        
        bool operator != (const rotator &v) const { return !((*this)==v); }

        rotator conjugate() const
        {
                rotator rr;
                rr.set(-_x,-_y,-_z,_w);
                return rr;
        }

        real len() const
        {
                return sqrt(_x*_x+_y*_y+_z*_z+_w*_w);
        }
        
        real dot(const rotator &r) const
        {
                return _x*r._x + _y*r._y + _z*r._z + _w*r._w;
        }

        rotator norm() const
        {
                rotator rr(*this);
                rr.normThis();
                return rr;
        }

        void normThis()
        {
                real n=len();
                if(n!=0.0){
                        _x/=n;
                        _y/=n;
                        _z/=n;
                        _w/=n;
                }
        }

        rotator inverse() const
        {
                rotator rr=conjugate();
                rr.normThis();
                return rr;
        }
        
        rotator interpolate(real val,const rotator& r) const
        {
                if(val>=1)
                        return r;

                if(val<=0)
                        return *this;

                rotator rr;
                rr.set(*this + ((dot(r)<0.0 ? -r : r) - *this)*val);
                rr.normThis();
                
                return rr;
        }

        void toMatrix(real* mat) const
        {
                rotator r=this->norm();
                real x2=r.x()*r.x(),y2=r.y()*r.y(),z2=r.z()*r.z(),
                         xy=r.x()*r.y(),xz=r.x()*r.z(),yz=r.y()*r.z(),
                         wz=r.w()*r.z(),wx=r.w()*r.x(),wy=r.w()*r.y();

                mat[ 0] = 1.0f - 2.0f * ( y2 + z2 ); 
                mat[ 1] = 2.0f * (xy - wz);
                mat[ 2] = 2.0f * (xz + wy);
                mat[ 3] = 0.0f;  
                mat[ 4] = 2.0f * (xy + wz);  
                mat[ 5] = 1.0f - 2.0f * (x2 + z2); 
                mat[ 6] = 2.0f * (yz - wx );  
                mat[ 7] = 0.0f;  
                mat[ 8] = 2.0f * (xz - wy);
                mat[ 9] = 2.0f * (yz + wx);
                mat[10] = 1.0f - 2.0f * (x2 + y2);  
                mat[11] = 0.0f;  
                mat[12] = 0;  
                mat[13] = 0;  
                mat[14] = 0;  
                mat[15] = 1.0f;
        }

        mat4 toMatrix() const
        {
                mat4 m;
                toMatrix((real*)m.m);
                return m;
        }

        i32 hash() const
        {
                union { real f; i64 i;} x,y,z,w;
                x.f=_x;
                y.f=_y;
                z.f=_z;
                w.f=_w;
                
                i64 hash=_HASH(x.i,_HASH(y.i,_HASH(z.i,w.i,12),13),14);
                return i32(hash>>32)^i32(hash);
        }
        
        friend std::ostream& operator << (std::ostream &out,const rotator &r)
        {
                return out << "rotator(" << r.x() << ", " << r.y() << ", " << r.z() << ", " << r.w() << ")";
        }
};

/*****************************************************************************************************************************/
/* Exceptions */
/*****************************************************************************************************************************/

class IndexException : public std::exception
{
protected:
        std::string name;
        std::string msg;
        size_t val, maxval;

        void setMsg()
        {
                std::ostringstream out;
                out << "Bad value " << val << " for index '" << name.c_str() << "' (0 <= " << name.c_str() << " < " << maxval << ")";
                msg=out.str();
        }

public:
        IndexException(std::string name, size_t val, size_t maxval) : name(name),val(val),maxval(maxval)
        {
                setMsg();
        }

        IndexException(const IndexException& ind) : name(ind.name),val(ind.val),maxval(ind.maxval)
        {
                setMsg();
        }

        virtual ~IndexException() throw() {}

        virtual const char* what() const throw() { return msg.c_str(); }
};

class MemException : public std::exception
{
protected:
        std::string msg;
public:
        MemException(const std::string msg) : msg(msg){}

        MemException(const std::string m,int err)
        {
                std::ostringstream out;
                out << m << ": errno: " << strerror(errno);
                msg=out.str();
        }

        MemException(const MemException & op) : msg(op.msg) {}

        virtual ~MemException() throw() {}

        virtual const char* what() const throw() { return msg.c_str(); }
};

class RenderException : public std::exception
{
protected:
        std::string msg;
public:
        RenderException(const std::string msg) : msg(msg){}

        RenderException(const std::string msg, const char* file, int line){
                std::ostringstream out;
                out << file << ":" << line << ":" << msg;
                this->msg=out.str();
        }

        RenderException(const RenderException & op) : msg(op.msg) {}

        virtual ~RenderException() throw() {}

        virtual const char* what() const throw() { return msg.c_str(); }
};

class ValueException : public std::exception
{
protected:
        std::string msg;
public:
        ValueException(const std::string valuename, const std::string msg, const char* file=NULL, int line=-1) 
        {
                std::ostringstream out;
                if(file)
                        out << file << ":" << line << ": ";
                out << "Bad value for " << valuename  << ": " << msg;
                this->msg=out.str();

        }

        ValueException(const ValueException & op) : msg(op.msg) {}

        virtual ~ValueException() throw() {}

        virtual const char* what() const throw() { return msg.c_str(); }
};

inline void checkNull(const std::string valuename, const void* val, const char* file=NULL, int line=-1) throw(ValueException)
{
        if(val==NULL)
                throw ValueException(valuename,"Must not be null",file,line);
}

#define CHECK_NULL(val) checkNull("val",(const void*)val,__FILE__,__LINE__)

/*****************************************************************************************************************************/
/* Data Structure Objects */
/*****************************************************************************************************************************/

void readBinaryFileToBuff(const char* filename,size_t offset,void* dest,size_t len) throw(MemException);

void storeBufftoBinaryFile(const char* filename,void* src,size_t len,int* header, size_t headerlen) throw(MemException);

class MetaType
{
protected:
        //typedef std::pair<std::string,std::string> strpair;
        typedef std::map<std::string,std::string> metamap;
        typedef metamap::iterator iter;
        typedef metamap::const_iterator citer;

        metamap _meta;
public:
        
        bool hasMetaKey(const char* key) const { return _meta.find(key)!=_meta.end(); }

        std::vector<std::string> getMetaKeys() const
        {
                std::vector<std::string> result;
                for(citer i=_meta.begin();i!=_meta.end();i++)
                        result.push_back((*i).first);

                return result;
        }

        std::string meta() const
        {
                std::ostringstream out;
                for(citer c=_meta.begin();c!=_meta.end();c++)
                        out << (*c).first << " = " << (*c).second << std::endl;

                return out.str();
        }

        const char* meta(const char* key) const
        {
                return hasMetaKey(key) ? _meta.find(key)->second.c_str() : "";
        }

        void meta(const char* key, const char* val)
        {
                if(!key || !val)
                        return;
                
                std::string fkey=key,fval=val;
                
                // filter out | characters
                fkey.erase (std::remove(fkey.begin(), fkey.end(), '|'), fkey.end());
                fval.erase (std::remove(fval.begin(), fval.end(), '|'), fval.end());

                _meta[fkey]=fval;
        }
        
        void meta(const MetaType* m)
        {
                for(citer i=m->_meta.begin();i!=m->_meta.end();i++)
                        meta((*i).first.c_str(),(*i).second.c_str());
        }
        
        std::string serializeMeta() const
        {
                std::ostringstream out;
                for(citer i=_meta.begin();i!=_meta.end();i++)
                        out << (*i).first << "||" << (*i).second << "||";
                
                return out.str();
        }
        
        void deserializeMeta(const std::string &s)
        {
                char *buf=new char[s.length()+1];
                strcpy(buf,s.c_str());
                
                char* p=strtok(buf,"|");
                
                while(p!=NULL){
                        std::string key=p;
                        p=strtok(NULL,"|");
                        std::string val=p;
                        _meta[key]=val;
                        p=strtok(NULL,"|");
                }
                
                delete buf;
        }
};

// Metaprogramming types encapsulating the 4 operators +-/*
template<typename R,typename LH, typename RH> struct AddOp { static inline R op(const LH& lh, const RH& rh) { return lh+rh; } };
template<typename R,typename LH, typename RH> struct SubOp { static inline R op(const LH& lh, const RH& rh) { return lh-rh; } };
template<typename R,typename LH, typename RH> struct DivOp { static inline R op(const LH& lh, const RH& rh) { return lh/rh; } };
template<typename R,typename LH, typename RH> struct MulOp { static inline R op(const LH& lh, const RH& rh) { return lh*rh; } };


// Metaprogramming types for encapsulating the endian swap functions
template<typename T> struct SwapEndian { static T swap(T t) { return swapEndianN(t); } };
template<> struct SwapEndian<real> { static real swap(real t) { return swapEndian64(t); } };
template<> struct SwapEndian<indexval> { static indexval swap(indexval t) { return swapEndian32(t); } };
template<> struct SwapEndian<vec3> { static vec3 swap(vec3 t) { return vec3(swapEndian64(t.x()),swapEndian64(t.y()),swapEndian64(t.z())); } };
template<> struct SwapEndian<color> { static color swap(color t) { return color(swapEndian32(t.r()),swapEndian32(t.g()),swapEndian32(t.b()),swapEndian32(t.a())); } };

template<typename T> class Matrix : public MetaType
{
protected:
        std::string _name; // matrix name, should be unique application-wide
        std::string _type; // type of matrix, may be an empty string if types don't make sense for what's stored
        std::string _sharedname; // name of the share memory segment `data' refers to

        T* data; // refers to locally-allocated segments of mmap-addressed shared memory segments
        sval _n_actual; // actual length of allocated data, >=_n
        sval _n,_m; // _n rows, _m columns

        bool _isShared; // true if the memory is shared, false if locally allocated memory
        
#ifdef WIN32
        HANDLE mapFile;
#endif

public:
        
        Matrix(const char* name,sval n, sval m=1,bool isShared=false)  throw(MemException) :
                        _name(name), _type(""),_sharedname(""),data(0),_n_actual(0),_n(n),_m(m),_isShared(false)
        {
                checkDimension("m",m);
                setShared(isShared);
        }

        Matrix(const char* name,const char* type,sval n, sval m=1,bool isShared=false)  throw(MemException) :
                        _name(name), _type(type),_sharedname(""),data(0),_n_actual(0),_n(n),_m(m),_isShared(false)
        {
                checkDimension("m",m);
                setShared(isShared);
        }

        Matrix(const char* name,const char* type,const char* sharedname,const char* serialmeta,sval n, sval m) throw(MemException)  :
                        _name(name), _type(type),_sharedname(sharedname),data(0),_n_actual(n),_n(n),_m(m),_isShared(true)
        {
                checkDimension("n",n);
                checkDimension("m",m);
                deserializeMeta(serialmeta);
                data=createShared();
        }

        Matrix(const char* name,const char* type,const T* array,sval n, sval m,bool isShared=false)  throw(MemException) :
                _name(name), _type(type),_sharedname(""),data(0),_n_actual(0),_n(n),_m(m),_isShared(false)
        {
                checkDimension("n",n);
                checkDimension("m",m);
                setShared(isShared);
                memcpy(data,array,memSize());
        }

        virtual ~Matrix()
        {
                clear();
        }

        T* dataPtr() const { return data; }

        Matrix<T>* clone(const char* newname=NULL, bool isShared=false) const
        {
                Matrix<T>* m=new Matrix<T>(newname ? newname : _name.c_str(),_type.c_str(),_n,_m,isShared);
                memcpy(m->data,data,memSize());
                m->meta(this);
                return m;
        }

        const char* getName() const { return _name.c_str(); }
        const char* getSharedName() const { return _sharedname.c_str(); }
        const char* getType() const { return _type.c_str(); }

        void setName(const char* name) { _name=name; }
        void setType(const char* type) { _type=type; }

        bool isShared() const { return _isShared; }

        void setShared(bool val) throw(MemException)
        {
                if(data && val==_isShared) // do nothing if the shared state to set is the current state and the matrix is allocated
                        return;

                if(val){
                        _n_actual=_n;
                        sval size=memSize();
                        
                        if(size==0)
                                throw MemException("Cannot make empty matrix shared");
                        
                        T* shared=createShared();

                        if(data){
                                memcpy(shared,data,size);
                                delete[] data;
                        }
                        else
                                memset(shared,0,size);

                        data=shared;
                }
                else{
                        T* olddata=data;
                        resize();
                        closeShared(olddata);
                }

                _isShared=val;
        }

        void clear() throw(MemException)
        {
                if(data){
                        if(_isShared){
                                closeShared(data);
                                unlinkShared(_sharedname);
                        }
                        else if(data!=NULL)
                                delete[] data;
                }

                _n_actual=0;
                _n=0;
                _isShared=false;
                data=NULL;
        }

        sval n() const { return _n; }
        
        sval m() const { return _m; }

        sval memSize() const { return sval(sizeof(T)*_n*_m); }

        void fill(const T& t)
        {
                T* d=data;
                for(sval n=0;n<_n*_m;n++)
                        *d++=t;
        }

        template<typename R> void copyFrom(const Matrix<R>* r)
        {
                sval minsize=_min(memSize(),r->memSize());
                if(minsize>0)
                        memcpy(data,r->dataPtr(),minsize);
        }

        Matrix<T>* subMatrix(const char* name,sval n, sval m=1,sval noff=0,sval moff=0,bool isShared=false) const throw(MemException)
        {
                checkDimension("n",n);
                checkDimension("m",m);
                
                if(n>_n || m>_m)
                        throw MemException("Submatrix dimensions may not exceed matrix dimensions");

                if((n+noff)>_n || (m+moff)>_m)
                        throw MemException("Submatrix dimensions plus offsets may not exceed matrix dimensions");

                Matrix<T>* mat=new Matrix<T>(name,_type.c_str(),n,m,isShared);

                for(sval nn=0;nn<n;nn++)
                        for(sval mm=0;mm<m;mm++)
                                mat->at(nn,mm)=at(nn+noff,mm+moff);

                return mat;
        }

        Matrix<T>* reshape(const char* name,sval n, sval m,bool isShared=false) const throw(MemException)
        {
                checkDimension("n",n);
                checkDimension("m",m);
                
                Matrix<T>* mat=new Matrix<T>(name,_type.c_str(),n,m,isShared);
                memcpy(mat->dataPtr(),data,_min(mat->memSize(),memSize()));
                return mat;
        }
        
        // These methods define mathematical operators for single value or matrix right hand sides
        
        template<typename Ctx>
        void applyFunc(T (*op)(Ctx,const T&,sval,sval),Ctx ctx,sval minrow=0,sval mincol=0,sval maxrow=-1,sval maxcol=-1)
        {
                maxcol=_min(_m,maxcol);
                maxrow=_min(_n,maxrow);
                
                for(sval n=minrow;n<maxrow;n++)
                        for(sval m=mincol;m<maxcol;m++)
                                at(n,m)=op(ctx,at(n,m),n,m);
        }
        
        template<typename R, typename OpType>
        void scalarop(const R& r,sval minrow=0,sval mincol=0,sval maxrow=-1,sval maxcol=-1)
        {
                maxcol=_min(_m,maxcol);
                maxrow=_min(_n,maxrow);
                
                for(sval n=minrow;n<maxrow;n++)
                        for(sval m=mincol;m<maxcol;m++)
                                at(n,m)=OpType::op(at(n,m),r);
        }

        template<typename R, typename OpType>
        void matop(const Matrix<R>& mat,sval minrow=0,sval mincol=0,sval maxrow=-1,sval maxcol=-1)
        {
                maxcol=_min(_min(mat.m(),_m),maxcol);
                maxrow=_min(_min(mat.n(),_n),maxrow);
                
                for(sval n=minrow;n<maxrow;n++)
                        for(sval m=mincol;m<maxcol;m++)
                                at(n,m)=OpType::op(at(n,m),mat(n,m));
        }
        
        template<typename R> void add(const R& r,sval minrow=0,sval mincol=0,sval maxrow=-1,sval maxcol=-1)
        {
                return scalarop<R,AddOp<T,T,R> >(r,minrow,mincol,maxrow,maxcol);
        }

        template<typename R> void sub(const R& r,sval minrow=0,sval mincol=0,sval maxrow=-1,sval maxcol=-1)
        {
                return scalarop<R,SubOp<T,T,R> >(r,minrow,mincol,maxrow,maxcol);
        }

        template<typename R> void mul(const R& r,sval minrow=0,sval mincol=0,sval maxrow=-1,sval maxcol=-1)
        {
                return scalarop<R,MulOp<T,T,R> >(r,minrow,mincol,maxrow,maxcol);
        }

        template<typename R> void div(const R& r,sval minrow=0,sval mincol=0,sval maxrow=-1,sval maxcol=-1)
        {
                return scalarop<R,DivOp<T,T,R> >(r,minrow,mincol,maxrow,maxcol);
        }
        
        template<typename R> void addm(const Matrix<R>& mat,sval minrow=0,sval mincol=0,sval maxrow=-1,sval maxcol=-1)
        {
                return matop<R,AddOp<T,T,R> >(mat,minrow,mincol,maxrow,maxcol);
        }
        
        template<typename R> void subm(const Matrix<R>& mat,sval minrow=0,sval mincol=0,sval maxrow=-1,sval maxcol=-1)
        {
                return matop<R,SubOp<T,T,R> >(mat,minrow,mincol,maxrow,maxcol);
        }
        
        template<typename R> void mulm(const Matrix<R>& mat,sval minrow=0,sval mincol=0,sval maxrow=-1,sval maxcol=-1)
        {
                return matop<R,MulOp<T,T,R> >(mat,minrow,mincol,maxrow,maxcol);
        }
        
        template<typename R> void divm(const Matrix<R>& mat,sval minrow=0,sval mincol=0,sval maxrow=-1,sval maxcol=-1)
        {
                return matop<R,DivOp<T,T,R> >(mat,minrow,mincol,maxrow,maxcol);
        }

        void reorderColumns(const sval *orderinds) throw(IndexException)
        {
                T* buff=new T[_m];

                for(sval j=0;j<_m;j++)
                        checkIndex("sortinds",orderinds[j],_m);

                for(sval i=0;i<_n;i++){
                        for(sval j=0;j<_m;j++)
                                buff[j]=at(i,j);

                        for(sval j=0;j<_m;j++)
                                at(i,j)=buff[orderinds[j]];
                }

                delete buff;
        }

        void swapEndian()
        {
                sval len=_n*_m;
                for(sval i=0;i<len;i++)
                        data[i]=SwapEndian<T>::swap(data[i]);
        }

        T& at(sval n, sval m=0) const { return data[m+(_m*n)]; }
        const T& atc(sval n, sval m=0) const { return data[m+(_m*n)]; }
        void ats(sval n, sval m, const T& t) { data[m+(_m*n)]=t; }
        
        T& operator () (sval n, sval m=0) const { return data[m+(_m*n)]; }
        
        T& operator [] (sval n) const { return data[n]; }

        T getAt(sval n, sval m=0) const throw(IndexException) { return data[getIndex(n,m)]; }

        void setAt(const T& t, sval n, sval m=0) throw(IndexException) { data[getIndex(n,m)]=t; }

        void setN(sval _newn) throw(MemException) 
        { 
                checkNotShared(); 
                _n=_newn; 
                resize(); 
        }

        void setM(sval _newm) throw(MemException)
        {
                checkDimension("_newm",_newm);
                checkNotShared();
                _newm=_max<sval>(1,_newm);
                if(_newm>(_m*_n))
                        throw MemException("New m value larger than matrix size");

                _n=sval((_n*_m)/_newm);
                _m=_newm;
        }

        void addRows(sval num) throw(MemException) { setN(_n+num); }

        void reserveRows(sval num) throw(MemException) { checkNotShared(); resize(num); }

        void append(const Matrix<T> &t) throw(MemException)
        {
                if(t.m()!=_m)
                        throw MemException("Column dimensions of `this' and `t' do not match");
                
                sval oldn=_n;
                addRows(t.n());
                memcpy(data+(_m*oldn),t.data,t.memSize());
        }

        void append(const T& t,sval m=0) throw(MemException)
        {
                addRows(1);
                at(_n-1,m)=t;
        }

        void removeRow(sval n) throw(MemException,IndexException)
        {
                checkNotShared();
                checkIndex("n",n,_n);

                if(n<(_n-1))
                        std::memmove(&data[n*_m],&data[(n+1)*_m],(_n-n-1)*_m*sizeof(T));

                setN(_n-1);
        }

        void readBinaryFile(const char* filename,size_t offset) throw(MemException)
        {
                readBinaryFileToBuff(filename,offset,data,memSize());
        }

        void readTextFile(const char* filename,sval numHeaders) throw(MemException)
        {
                setN(0);
                readTextFileMatrix(filename,numHeaders,this);
        }
        
        void storeBinaryFile(const char* filename, int* header, sval headerlen) throw(MemException)
        {
                storeBufftoBinaryFile(filename,data,memSize(),header,headerlen);
        }

        indexpair indexOf(const T& t,sval aftern=0,sval afterm=0) const
        {
                sval numelems=_n*_m;
                sval nm=afterm+aftern*_m;

                while(nm<numelems && data[nm]!=t)
                        nm++;

                return indexpair(nm/_m,nm%_m);
        }

protected:
        
        inline sval getIndex(sval n, sval m) const throw(IndexException)
        {
                checkIndex("n",n,_n);
                checkIndex("m",m,_m);

                return m+(_m*n);
        }

        inline void checkIndex(const char* name, sval val, sval maxval) const throw(IndexException)
        {
                if(val>=maxval)
                        throw IndexException(name,val,maxval);
        }

        inline void checkNotShared() const throw(MemException)
        {
                if(_isShared)
                        throw MemException("Operation may only be performed on non-shared matrices");
        }
        
        inline void checkDimension(const char* name, sval dim) const throw(MemException)
        {
                if(dim==0){
                        std::ostringstream out;
                        out << "Matrix dimension " << name << " must be non-zero";
                        throw MemException(out.str());
                }
        }

        void resize(sval reserveNum=0)
        {
                if((_n+reserveNum)<=_n_actual && !_isShared) // do nothing if we have enough rows to meet demand and we're not shared
                        return;

                // arbitrary min allocation of 1000 rows and 50% reserve, or just _n if this is the first allocation
                sval new_n_actual=data==NULL ? _n : _max<sval>(1000ul,((_n*3)/2)+reserveNum);

                if(new_n_actual<_n_actual && !_isShared)
                        return;

                T* olddata=data;
                data=new T[new_n_actual*_m];

                memset(data,0,new_n_actual*_m*sizeof(T)); // always do this to ensure rows past _n_actual are clean

                if(olddata){
                        memcpy(data,olddata,_n_actual*_m*sizeof(T));
                        if(!_isShared)
                                delete[] olddata;
                }
                
                _n_actual=new_n_actual;
        }

        void chooseSharedName(int counter=0)
        {
                std::ostringstream out;
#ifdef WIN32
                out << "Local\\" ;

                if(counter>0)
                        out << std::hex << counter << std::dec;

                out << GetCurrentProcessId() << _name;
                _sharedname=out.str();
#else

// OSX shm_open requires names to be shorter than normal
#ifdef __APPLE__
                if(counter>0)
                        out << std::hex << counter << std::dec;

                out << getpid() << _name;
#else
                out << "__viz__" << getppid() << "_" << getpid();
                if(counter>0)
                        out << "_" << std::hex << counter << std::dec;

                out << "_" << _name;
#endif // __APPLE__

                _sharedname=out.str();
                
                for(sval i=0;i<_sharedname.size();i++)
                        if(_sharedname[i]=='/')
                                _sharedname[i]='_';

                if(_sharedname.size()>=MAXSHMNAMLEN)
                        _sharedname[MAXSHMNAMLEN-1]='\0';
#endif // WIN32
        }

        T* createShared() throw(MemException)
        {
                sval size=memSize();
                T* ptr;
                bool isCreator=_sharedname=="";

                if(isCreator)
                        chooseSharedName(); // start with the default shared name

#ifdef WIN32
                std::ostringstream out;

                createFileMapping();
                
                // attempt to choose a unique shared name only when creating a segment
                for(int c=1;c<100000 && isCreator && mapFile && GetLastError()==ERROR_ALREADY_EXISTS;c++){
                        CloseHandle(mapFile);
                        chooseSharedName(c);
                        createFileMapping();
                }

                if(!mapFile){
                        out << "Unable to open shared memory handle to " << _sharedname << ": " << formatLastErrorMsg();
                        throw MemException(out.str());
                }

                ptr=(T*)MapViewOfFile(mapFile,FILE_MAP_ALL_ACCESS,0,0,size);

                if(!ptr && isCreator && GetLastError()==ERROR_ACCESS_DENIED){
                        for(int c=1;c<100000 && mapFile && !ptr && (GetLastError()==ERROR_ACCESS_DENIED || GetLastError()==ERROR_ALREADY_EXISTS);c++){
                                CloseHandle(mapFile);
                                chooseSharedName(c);
                                createFileMapping();

                                if(mapFile)
                                        ptr=(T*)MapViewOfFile(mapFile,FILE_MAP_ALL_ACCESS,0,0,size);
                        }
                }

                if(!ptr){
                        CloseHandle(mapFile);
                        out << "Unable to map view of memory file " << _sharedname << ": " << formatLastErrorMsg();
                        throw MemException(out.str());
                }

#else // Linux/OSX

                int shm_fd=shm_open(_sharedname.c_str(), O_CREAT | O_RDWR | (isCreator ? O_EXCL : 0), S_IRUSR | S_IWUSR);

                // attempt to choose a unique shared name only when creating a segment
                for(int c=1;c<100000 && isCreator && shm_fd==-1 && errno==EEXIST;c++){
                        chooseSharedName(c);
                        shm_fd=shm_open(_sharedname.c_str(), O_CREAT | O_RDWR | O_EXCL, S_IRUSR | S_IWUSR);
                }

                if(shm_fd==-1)
                        throw MemException(std::string("Unable to open shared memory descriptor, filename:")+_sharedname,errno);

                if(isCreator && ftruncate(shm_fd, size) == -1)
                        throw MemException(std::string("Unable to extend shared memory section, filename:")+_sharedname,errno);

                ptr = (T*)mmap(0, size, PROT_READ | PROT_WRITE, MAP_SHARED, shm_fd, 0);

                if(!ptr || ptr==MAP_FAILED){
                        if(isCreator)
                                shm_unlink(_sharedname.c_str());

                        throw MemException("Unable to mmap shared memory",errno);
                }
                
                addShared(_sharedname); // store the name for later cleanup

                close(shm_fd);
#endif // Linux/OSX

                return ptr;
        }

        void closeShared(T* ptr) throw(MemException)
        {
                if(!ptr)
                        return;

#ifdef WIN32
                if(!UnmapViewOfFile(ptr))
                        throw MemException("Failed to unmap file view");

                if(!CloseHandle(mapFile))
                        throw MemException("Failed to close handle");
#else
                int mures=munmap(ptr,memSize());    
                if(mures==-1)
                        throw MemException("Failed to unmap memory section");
#endif // WIN32
                _sharedname="";
        }

#ifdef WIN32
        void createFileMapping()
        {
#ifdef UNICODE
                const char* name=_sharedname.c_str();
                wchar_t namebuff[1024];

                ::MultiByteToWideChar(CP_ACP, NULL,name, -1, namebuff,int(_sharedname.size()+1));
#else
                const char* namebuff=_sharedname.c_str();
#endif // UNICODE

                mapFile = CreateFileMapping(INVALID_HANDLE_VALUE, NULL, PAGE_READWRITE, 0, DWORD(memSize()), namebuff);
        }
#endif // WIN32
};

typedef Matrix<real> RealMatrix;
typedef Matrix<vec3> Vec3Matrix;
typedef Matrix<indexval> IndexMatrix;
typedef Matrix<color> ColorMatrix;

class DataSet : public MetaType
{
protected:
        std::string name;
        std::vector<std::string> indexnames;
        std::vector<std::string> fieldnames;

public:
        DataSet(const char* name): name(name) { }

        const char* getName() const { return name.c_str(); }
        
        virtual DataSet* clone(const char* name,bool cloneNodes=false){ return NULL;}

        virtual Vec3Matrix* getNodes() const { return NULL ;}
        virtual void setNodes(Vec3Matrix *nodes) { }
        
        virtual std::vector<std::string> getIndexNames() const { return indexnames; }
        
        virtual void setIndexNames(std::vector<std::string> names) { indexnames=names; }

        virtual IndexMatrix* getIndexSet(const char* name) const { return NULL; }

        virtual bool hasIndexSet(const char* name) const { return getIndexSet(name)!=NULL; }
        
        virtual void setIndexSet(IndexMatrix *indices,const char *alias=NULL) {}
        
        virtual std::vector<std::string> getFieldNames() const { return fieldnames; }
        
        virtual void setFieldNames(std::vector<std::string> names) { fieldnames=names; }
        
        virtual RealMatrix* getDataField(const char* name) const { return NULL; }
        
        virtual bool hasDataField(const char* name) const { return getDataField(name)!=NULL; }
        
        virtual void setDataField(RealMatrix *field,const char *alias=NULL) {}
};

/*****************************************************************************************************************************/
/* Scene Objects */
/*****************************************************************************************************************************/

class Material;

class Texture
{
public:
        Texture(){}
        virtual ~Texture() {}

        virtual const char* getName() const {return "";}
        virtual const char* getFilename() const {return "";}
        virtual sval getWidth() const { return 0;}
        virtual sval getHeight() const { return 0;}
        virtual sval getDepth() const { return 0; }
        virtual bool hasAlpha() const {return false;}
        virtual TextureFormat getFormat() const { return TF_UNKNOWN; }
        virtual void fillBlack() {}
        virtual void fillColor(color col) {}
        virtual void fillColor(const ColorMatrix *mat,indexval depth) {}
        virtual void fillColor(const RealMatrix *mat,indexval depth,real minval=0.0,real maxval=1.0, const Material* colormat=NULL,const RealMatrix *alphamat=NULL,bool mulAlpha=true) {}
};

class GPUProgram
{
public:
        virtual ~GPUProgram(){}

        virtual std::string getName() const {return "";}
        virtual void setType(ProgramType pt) {}
        virtual ProgramType getType() const { return PT_VERTEX; }
        virtual std::string getLanguage() const { return ""; }
        
        virtual void setLanguage(const std::string& lang) {}
        
        virtual void setSourceCode(const std::string& code){}
        
        virtual bool hasError() const { return false; }
        
        virtual std::string getSourceCode() const { return ""; }
        
        virtual bool setParameter(const std::string& param, const std::string& val) { return false; }
        virtual std::string getParameter(const std::string& param) const { return ""; }
        virtual std::string getEntryPoint() const { return getParameter("entry_point"); }
        virtual std::string getProfiles() const { return getParameter("profiles"); }
        virtual std::vector<std::string> getParameterNames() const { return std::vector<std::string>(); }
        virtual void setEntryPoint(const std::string main) { setParameter("entry_point",main); }
        virtual void setProfiles(const std::string profiles) { setParameter("profiles",profiles); }
};

template<typename T>
class PositionQueue
{
protected:
        typedef std::pair<real,T> ListValue;
        Matrix<ListValue> vals;

public:
        PositionQueue() : vals("vals",0) {}
        virtual ~PositionQueue() {}

        void copyFrom(const PositionQueue<T>* queue)
        {
                clear();
                vals.append(queue->vals);
        }

        void add(real pos,T val)
        {
                vals.append(ListValue(pos,val));
                sort();
        }

        void fill(const RealMatrix* pos, const Matrix<T>* ctrls)
        {
                vals.setN(pos->n());
                for(indexval i=0;i<pos->n();i++)
                        vals[i]=ListValue(pos->at(i),ctrls->at(i));
        }

        sval size() const { return vals.n();}

        void clear() { vals.setN(0); }

        T get(indexval index) const throw(IndexException)
        {
                return vals.getAt(index).second;
        }

        void set(indexval index,real pos, T value) throw(IndexException)
        {
                vals.setAt(ListValue(pos,value),index);
                sort();
        }

        real pos(indexval index) const throw(IndexException)
        {
                return vals.getAt(index).first;
        }

        void remove(indexval index) throw(IndexException)
        {
                vals.removeRow(index);
                sort();
        }

        indexval find(real pos, const T& value) const
        {
                indexval i=0;
                sval numVals=vals.n();

                while(i<numVals && (!equalsEpsilon(vals[i].first,pos) || vals[i].second!=value))
                        i++;

                return i;
        }

        void sort()
        {
                if(vals.n()>1)
                        qsort(vals.dataPtr(),vals.n(),sizeof(ListValue),sortTupleFirstCB<ListValue>);
        }
};

template<typename T>
class ControlCurve
{
protected:
        Matrix<T> ctrls;
        Matrix<T> derivs;

public:
        ControlCurve() : ctrls("ctrls","",0,1), derivs("derivs","",0,2) {}
        virtual ~ControlCurve() {} 

        void copyFrom(const ControlCurve<T> *con)
        {
                ctrls.setN(0);
                ctrls.append(con->ctrls);
                calculateDerivs();
        }

        void clear()
        {
                ctrls.setN(0);
                derivs.setN(0);
        }

        virtual void addCtrlPoint(const T& t) { ctrls.append(t); calculateDerivs(); }
        virtual void setCtrlPoint(const T& t,indexval index) throw(IndexException) { ctrls.setAt(t,index); calculateDerivs(); }
        virtual void removeCtrlPoint(indexval index) throw(IndexException) { ctrls.removeRow(index); calculateDerivs(); }
        virtual sval numPoints() const { return ctrls.n(); }
        virtual T getCtrlPoint(indexval index) const throw(IndexException) { return ctrls.getAt(index); }

        virtual void setCtrlPoints(const Matrix<T>* pts)
        {
                ctrls.setN(0);
                for(sval i=0;i<pts->n();i++)
                        ctrls.append(pts->at(i));

                calculateDerivs();
        }

        virtual void calculateDerivs()
        {
                sval n=ctrls.n();
                derivs.setN(n);

                if(n==1){
                        derivs(0,0)=ctrls[0];
                        derivs(0,1)=ctrls[0];
                }
                else if(n==2){
                        derivs(0,0)=ctrls[0];
                        derivs(0,1)=ctrls[1];
                        derivs(1,0)=ctrls[1];
                        derivs(1,1)=ctrls[0];
                }
                else if(n>2){
                        RealMatrix mat("",n,3);

                        Matrix<T> localderivs("localderivs",n,1);

                        localderivs[0]=ctrls[0];
                        localderivs[n-1]=ctrls[n-1];

                        if(n==3)
                                localderivs[1]=ctrls[1]*6-ctrls[0]-ctrls[n-1];
                        else{
                                localderivs[1]=ctrls[1]*6-ctrls[0];
                                localderivs[n-2]=ctrls[n-2]*6-ctrls[n-1];
                        
                                for(indexval i=2;i<n-2;i++)
                                        localderivs[i]=ctrls[i]*6;
                        }

                        for(indexval i=0;i<n;i++){
                                mat(i,0)=4.0;
                                mat(i,1)=1.0;
                                mat(i,2)=1.0;
                        }

                        // Gaussian elimination
                        for(indexval i=2;i<n-1;i++){
                                mat(i,1)=mat(i,1)/mat(i-1,0);
                                mat(i,0)=mat(i,0)-(mat(i,1)*mat(i-1,2));
                                localderivs[i]=localderivs[i]-(localderivs[i-1]*mat(i,1));
                        }

                        localderivs[n-2]=localderivs[n-2]/mat(n-2,0);

                        for(indexval i=n-3;i>0;i--)
                                localderivs[i]=(localderivs[i]-(localderivs[i+1]*mat(i,2)))/mat(i,0);

                        for(indexval s=0;s<n;s++){
                                indexval e=clamp<indexval>(s+1,0,n-1);
                                derivs(s,0)=(localderivs[s] * (2.0 / 3.0)) + (localderivs[e] / 3.0);
                                derivs(s,1)=(localderivs[s] / 3.0) + (localderivs[e] * (2.0 / 3.0));
                        }
                }
                // else n==0 do nothing
        }

        virtual T at(real tt) const
        {
                indexval n=ctrls.n();
                real tn=tt*(n-1);
                
                indexval s=clamp<indexval>(indexval(tn),0,n-1);
                indexval e=clamp<indexval>(indexval(tn)+1,0,n-1);
                
                real t=tn-s;
                real t1=1.0-t;
                
                T d1=derivs(s,0); 
                T d2=derivs(s,1); 

                return (ctrls[s]*(t1*t1*t1))+(ctrls[e]*(t*t*t))+(d1*(3*t1*t1*t))+(d2*(3*t1*t*t)); // cubic bezier spline
        }
};


/*template<typename T>
class ControlCurve
{
protected:
        Matrix<T> ctrls;

public:
        ControlCurve() : ctrls("ctrls","",0,1) {}
        virtual ~ControlCurve() {} 

        virtual void addCtrlPoint(const T& t) { ctrls.append(t); }
        virtual void setCtrlPoint(const T& t,indexval index) throw(IndexException) { ctrls.setAt(t,index); }
        virtual void removeCtrlPoint(indexval index) throw(IndexException) { ctrls.removeRow(index); }
        virtual sval numPoints() const { return ctrls.n(); }
        virtual T getCtrlPoint(indexval index) const throw(IndexException) { return ctrls.getAt(index); }

        virtual void setCtrlPoints(const Matrix<T>* pts)
        {
                ctrls.setN(0);
                for(sval i=0;i<pts->n();i++)
                        ctrls.append(pts->at(i));
        }

        virtual void calculateDerivs(){}

        virtual T at(real tt) const
        {
                indexval n=ctrls.n();
                real tn=tt*(n-1);

                indexval s=clamp<indexval>(indexval(tn),0,n-1);
                indexval e=clamp<indexval>(indexval(tn)+1,0,n-1);
                indexval ds=clamp<indexval>(indexval(tn)-1,0,n-1);
                indexval de=clamp<indexval>(indexval(tn)+2,0,n-1);

                real t=tn-s;

                real t2=t*t;
                real t3=t2*t;
                real t3_05=t3*0.5;
                real t3_15=t3*1.5;
                real t_05=t*0.5;

                real cs=t3_15-2.5*t2+1;
                real ce=2*t2+t_05-t3_15;
                real dcs=t2-t_05-t3_05;
                real dce=t3_05-0.5*t2;

                T d1=ctrls[ds];
                T d2=ctrls[de];

                return ctrls[s]*cs+ctrls[e]*ce+d1*dcs+d2*dce;
        }
};*/

class Vec3Curve : public ControlCurve<vec3>
{
        bool isXFunc;
        bool isLinear;
public:
        Vec3Curve(bool isXFunc) : ControlCurve<vec3>(), isXFunc(isXFunc),isLinear(false) {}

        void setLinear(bool b) { isLinear=b; }
        bool isLinearFunc() const { return isLinear; }

        virtual void addCtrlPoint(const vec3& t)
        {
                vec3 tt=t;
                if(isXFunc){
                        real minx=ctrls.n()==0 ? 0.0 : ctrls[ctrls.n()-1].x();
                        tt=vec3(clamp<real>(t.x(),minx,1.0),clamp<real>(t.y(),0.0,1.0),0);
                }
                ControlCurve<vec3>::addCtrlPoint(tt);
        }

        virtual void setCtrlPoint(const vec3& t,indexval index) throw(IndexException)
        {
                vec3 tt=t;
                if(isXFunc){
                        real minx=index==0 ? 0.0 : ctrls[index-1].x();
                        tt=vec3(clamp<real>(t.x(),minx,1.0),clamp<real>(t.y(),0.0,1.0),0);
                }

                ctrls.setAt(tt,index); 

                if(isXFunc)
                        for(sval i=index+1;i<ctrls.n();i++)
                                ctrls[i].x(clamp<real>(ctrls[i].x(),ctrls[i-1].x(),1.0));

                calculateDerivs();
        }

        virtual void calculateDerivs()
        {
                if(isXFunc)
                        qsort(ctrls.dataPtr(),ctrls.n(),sizeof(vec3),vec3::compX);
                
                ControlCurve<vec3>::calculateDerivs();
        }

        real atX(real x, real threshold=0.0001) const
        {
                if(x<=ctrls[0].x())
                        return ctrls[0].y();

                if(x>=ctrls[ctrls.n()-1].x())
                        return ctrls[ctrls.n()-1].y();

                if(isLinear){
                        sval i=1;
                        while(ctrls[i].x()<x)
                                i++;

                        real xi=lerpXi<real>(x,ctrls[i-1].x(),ctrls[i].x());
                        return lerp<real>(xi,ctrls[i-1].y(),ctrls[i].y());
                }
                else{
                        real start=0.0, end=1.0, mid=0.5;

                        vec3 val=at(mid);
                        real diff=val.x()-x;

                        // bsearch for interpolated value whose X coord is close to `x'
                        while(fabs(diff)>threshold && (end-start)>threshold){
                                if(diff>0)
                                        end=mid;
                                else
                                        start=mid;

                                mid=start+(end-start)*0.5;
                                val=at(mid);
                                diff=val.x()-x;
                        }

                        return val.y();
                }
        }
};


class Spectrum
{
protected:
        PositionQueue<color> spec;

        Vec3Curve alphacurve;

        std::string name;

public:
        Spectrum(const std::string& name="") : alphacurve(true), name(name) {}

        virtual const char* getName() const { return name.c_str(); }

        virtual void clearSpectrum()
        {
                spec.clear();
                alphacurve.clear();
                updateSpectrum();
        }

        virtual void copySpectrumFrom(const Spectrum* s)
        {
                spec.copyFrom(&s->spec);
                alphacurve.copyFrom(&s->alphacurve);
                updateSpectrum();
        }
        
        virtual color getDefaultColor() const { return color(); }

        virtual real getAlpha() const { return 1.0; }

        virtual void updateSpectrum() {}

        virtual void addSpectrumValue(real pos,color value)
        {
                spec.add(pos,value);
                updateSpectrum();
        }

        virtual real getSpectrumPos(indexval index) const throw(IndexException)
        {
                return spec.pos(index);
        }

        virtual color getSpectrumValue(indexval index) const throw(IndexException)
        {
                return spec.get(index);
        }

        virtual indexval getSpectrumIndex(real pos,color value) const
        {
                return spec.find(pos,value);
        }

        virtual void setSpectrumValue(sval index, real pos,color value) throw(IndexException)
        {
                spec.set(index,pos,value);
                updateSpectrum();
        }

        virtual sval numSpectrumValues() const
        {
                return spec.size(); 
        }

        virtual void removeSpectrumValue(indexval index) throw(IndexException)
        {
                spec.remove(index);
                updateSpectrum();
        }

        virtual sval numAlphaCtrls() const 
        { 
                return alphacurve.numPoints(); 
        }

        virtual vec3 getAlphaCtrl(indexval index) const throw(IndexException) 
        { 
                return alphacurve.getCtrlPoint(index); 
        }

        virtual void addAlphaCtrl(vec3 v) 
        { 
                alphacurve.addCtrlPoint(v); 
                updateSpectrum();
        }

        virtual void removeAlphaCtrl(indexval index) throw(IndexException) 
        {
                alphacurve.removeCtrlPoint(index);
                updateSpectrum();
        }

        virtual void setAlphaCtrl(vec3 v, indexval index) throw(IndexException) 
        {
                alphacurve.setCtrlPoint(v,index);
                updateSpectrum();
        }

        virtual void setAlphaCurve(const Vec3Matrix* pts)
        {
                alphacurve.setCtrlPoints(pts);
                updateSpectrum();
        }

        virtual void setLinearAlpha(bool b)
        {
                alphacurve.setLinear(b);
                updateSpectrum();
        }

        virtual bool isLinearAlpha() const 
        {
                return alphacurve.isLinearFunc();
        }

        virtual color interpolateColor(real pos) const
        {
                color result;
                float alpha=getAlpha();
                indexval specsize=spec.size();

                if(specsize==0)
                        result=getDefaultColor();
                else if(pos<=spec.pos(0))
                        result= spec.get(0);
                else if(pos>=spec.pos(specsize-1))
                        result=spec.get(specsize-1);
                else{
                        sval index=0;
                        while(index<specsize-1 && spec.pos(index+1)<pos)
                                index++;

                        color cmin=spec.get(index);
                        color cmax=spec.get(index+1);

                        real interp=lerpXi(pos,spec.pos(index),spec.pos(index+1));

                        result= cmin.interpolate(interp,cmax);
                }

                if(alphacurve.numPoints()>1){
                        real a=clamp<real>(alphacurve.atX(pos),0.0,1.0);
                        if(alpha>=0.0 && specsize>0)
                                a*=alpha;

                        result.a(a);
                }
                else if(alpha>=0.0)
                        result.a(alpha);

                return result;
        }
        
        virtual void fillColorMatrix(ColorMatrix *col,const RealMatrix *mat,bool useValAsAlpha=false) throw(IndexException) 
        {
                bool hasMatAlpha=mat->m()>=col->m()*2;
                sval len=_min(col->n(),mat->n());
                sval width=_min(col->m(),mat->m());
                
                for(sval i=0;i<len;i++){
                        for(sval j=0;j<width;j++){
                                float val=mat->getAt(i,j);
                                color c=interpolateColor(val);
                                
                                if(hasMatAlpha)
                                        c.a(mat->getAt(i,col->m()+j)*c.a());
                                else if(useValAsAlpha)
                                        c.a(val*c.a());
                                
                                col->setAt(c,i,j);
                        }
                }
        }
};


class Material : public Spectrum
{
public:
        float alpha;
        bool useAlpha;

        Material() : alpha(1.0), useAlpha(true) {} 
        virtual ~Material(){}

        virtual Material* clone(const char* name) const { return NULL; }
        
        virtual void copyTo(Material* mat,bool copyTex=false,bool copySpec=false,bool copyProgs=false) const {}
        
        virtual real getAlpha() const { return useAlpha ? alpha : -1.0; }

        virtual color getDefaultColor() const { return getDiffuse(); }
        
        virtual void setAlpha(real alpha)
        {
                this->alpha=(float)alpha;
                setDiffuse(getDiffuse());
                setSpecular(getSpecular());
        }

        virtual bool usesInternalAlpha() const { return useAlpha; }
        
        virtual void useInternalAlpha(bool val)
        {
                useAlpha=val;
                setDiffuse(getDiffuse());
                setSpecular(getSpecular());
        }

        virtual color getAmbient() const { return color(); }
        virtual color getDiffuse() const { return color(); }
        virtual color getSpecular() const { return color(); }
        virtual color getEmissive() const { return color(); }

        virtual real getShininess() const { return 0.0f; }
        virtual real getPointSizeMin() const { return 0.0f; }
        virtual real getPointSizeMax() const { return 0.0f; }
        virtual real getPointSizeAbs() const { return 0.0f; }
        virtual bool usesPointAttenuation() const { return false; }
        virtual BlendMode getBlendMode() const { return BM_ALPHA; }

        virtual bool usesVertexColor() const { return false; }
        virtual bool usesLighting() const { return false; }
        virtual bool usesFlatShading() const { return false; }
        virtual bool usesDepthCheck() const { return false; }
        virtual bool usesDepthWrite() const { return false; }
        virtual bool usesTexFiltering() const { return false; }
        virtual bool isClampTexAddress() const { return false;}
        virtual bool isCullBackfaces() const { return false; }
        virtual bool usesPointSprites() const { return false; }
        virtual const char* getTexture() const { return ""; }
        virtual const char* getGPUProgram(ProgramType pt) const {  return ""; }

        virtual int getGPUParamInt(ProgramType pt, const std::string& name) { return 0; }
        virtual real getGPUParamReal(ProgramType pt, const std::string& name) { return 0; }
        virtual vec3 getGPUParamVec3(ProgramType pt, const std::string& name) { return vec3(); }
        virtual color getGPUParamColor(ProgramType pt, const std::string& name) { return color(); }

        virtual bool isTransparentColor() const { return useAlpha && alpha<1.0 && !usesVertexColor() && !strcmp(getTexture(),""); }

        virtual void setAmbient(const color & c) {}
        virtual void setDiffuse(const color & c){}
        virtual void setSpecular(const color & c){}
        virtual void setEmissive(const color & c) {}
        virtual void setShininess(real c){}
        
        virtual void setPointSize(real min,real max) {}
        virtual void setPointSizeAbs(real size) {}
        virtual void setPointAttenuation(bool enabled,real constant=0.0f,real linear=1.0f, real quad=0.0f) {}
        
        virtual void setBlendMode(BlendMode bm) {}
        virtual void useVertexColor(bool use) {}
        virtual void useLighting(bool use) {}
        virtual void useFlatShading(bool use) {}
        virtual void useDepthCheck(bool use) {}
        virtual void useDepthWrite(bool use) {}
        virtual void useTexFiltering(bool use) {}
        virtual void clampTexAddress(bool use) {}
        virtual void cullBackfaces(bool cull) {}
        virtual void usePointSprites(bool useSprites){}
        virtual void setTexture(const char* name){}
        virtual void setTexture(const Texture* tex) { setTexture(tex->getName()); }

        virtual void useSpectrumTexture(bool use) {}
        
        virtual void setGPUProgram(const std::string& name, ProgramType pt) {}
        virtual void setGPUProgram(const GPUProgram *prog) { setGPUProgram(prog->getName(),prog->getType()); }

        virtual bool setGPUParamInt(ProgramType pt,const std::string& name, int val) { return false; }
        virtual bool setGPUParamReal(ProgramType pt,const std::string& name, real val) { return false; }
        virtual bool setGPUParamVec3(ProgramType pt,const std::string& name, vec3 val) { return false; }
        virtual bool setGPUParamColor(ProgramType pt,const std::string& name, color val) { return false; }

        virtual void updateSpectrum() {}
};

class Light
{
public:
        Light(){}
        virtual ~Light() {}

        virtual void setPosition(vec3 &v){}
        
        virtual void setDirection(vec3 &v){}
        
        virtual void setDiffuse(const color & c){}
        
        virtual void setSpecular(const color & c){}

        virtual void setDirectional() {}
        
        virtual void setPoint() {}
        
        virtual void setSpotlight(real radsInner, real radsOuter, real falloff=1.0f) {}
        
        virtual void setAttenuation(real range, real constant=0.0f,real linear=1.0f, real quad=0.0f) {}

        virtual void setVisible(bool isVisible){}
        
        virtual bool isVisible() const { return false; }
};

class VertexBuffer
{
public:
        virtual ~VertexBuffer() {}

        virtual vec3 getVertex(int i) const { return vec3(0,0,0); }
        virtual vec3 getNormal(int i) const { return vec3(0,0,0); }
        virtual color getColor(int i) const { return color(0,0,0,0); }
        virtual vec3 getUVWCoord(int i) const { return vec3(0,0,0); }

        virtual sval numVertices() const { return 0; }
        virtual bool hasNormal() const { return false; }
        virtual bool hasColor() const { return false; }
        virtual bool hasUVWCoord() const { return false; }
};

class IndexBuffer
{
public:
        virtual ~IndexBuffer() {}
        
        virtual sval numIndices() const { return 0; }
        virtual sval indexWidth(int i) const { return 0; }
        virtual sval getIndex(int i,int w) const { return 0; }
};

template<typename Ctx>
class CallbackVertexBuffer : public VertexBuffer
{
public:
        typedef vec3 (*vecfunc)(Ctx,int);
        typedef color (*colfunc)(Ctx,int);

        vecfunc vertfunc;
        vecfunc normalfunc;
        colfunc colorfunc;
        vecfunc uvwfunc;
        sval numvertices;
        Ctx context;

        CallbackVertexBuffer(Ctx context, sval numvertices, vecfunc vertfunc, vecfunc normalfunc=NULL, colfunc colorfunc=NULL, vecfunc uvwfunc=NULL):
                vertfunc(vertfunc),normalfunc(normalfunc),colorfunc(colorfunc),uvwfunc(uvwfunc),numvertices(numvertices),context(context)
        { }

        virtual vec3 getVertex(int i) const { return vertfunc(context,i); }     
        virtual vec3 getNormal(int i) const { return normalfunc(context,i); }   
        virtual color getColor(int i) const { return colorfunc(context,i); }    
        virtual vec3 getUVWCoord(int i) const { return uvwfunc(context,i); }

        virtual sval numVertices() const { return numvertices; }
        virtual bool hasNormal() const { return normalfunc!=NULL; }
        virtual bool hasColor() const { return colorfunc!=NULL; }
        virtual bool hasUVWCoord() const { return uvwfunc!=NULL; }
};

template<typename Ctx>
class CallbackIndexBuffer : public IndexBuffer
{
public:
        typedef sval (*wfunc)(Ctx,int);
        typedef sval (*ifunc)(Ctx,int,int);

        wfunc widthfunc;
        ifunc indexfunc;

        sval numindices;
        Ctx context;

        CallbackIndexBuffer(Ctx context, sval numindices, wfunc widthfunc,ifunc indexfunc):
                widthfunc(widthfunc), indexfunc(indexfunc), numindices(numindices), context(context)
        {}

        virtual sval numIndices() const { return numindices; }
        virtual sval indexWidth(int i) const { return widthfunc(context,i); }
        virtual sval getIndex(int i,int w) const { return indexfunc(context,i,w); }
};

class MatrixVertexBuffer : public VertexBuffer
{
        Vec3Matrix* vecs;
        ColorMatrix* cols;
        IndexMatrix* extinds;
        sval numverts;
        bool deleteMatrices;

public:
        MatrixVertexBuffer(Vec3Matrix* vecs,ColorMatrix* cols=NULL,IndexMatrix* extinds=NULL) throw(RenderException) 
                : vecs(vecs), cols(cols),extinds(extinds), deleteMatrices(false)
        {
                if(vecs==NULL)
                        throw RenderException("Matrix 'vecs' must be provided");

                numverts=sval(extinds!=NULL ? extinds->n() : vecs->n());
        }
        
        MatrixVertexBuffer(const VertexBuffer* buf) throw(RenderException) : numverts(buf ? buf->numVertices() : 0), cols(NULL), extinds(NULL), deleteMatrices(true)
        {
                if(!numverts)
                        throw RenderException("VertexBuffer 'buf' must be provided");
                
                int columns=1;
                if(buf->hasUVWCoord())
                        columns=4;
                else if(buf->hasNormal())
                        columns=2;
                
                vecs=new Vec3Matrix("copyvecs",buf->numVertices(),columns);
                
                if(buf->hasColor())
                        cols=new ColorMatrix("copycols",buf->numVertices());
                
                for(sval i=0;i<numverts;i++){
                        vecs->at(i,0)=buf->getVertex(i);
                        if(buf->hasNormal())
                                vecs->at(i,1)=buf->getNormal(i);
                        if(buf->hasUVWCoord())
                                vecs->at(i,3)=buf->getUVWCoord(i);
                        
                        if(buf->hasColor())
                                cols->at(i)=buf->getColor(i);
                }
        }
        
        virtual ~MatrixVertexBuffer()
        {
                if(deleteMatrices){
                        SAFE_DELETE(vecs);
                        SAFE_DELETE(cols);
                        SAFE_DELETE(extinds);
                }
        }

        sval getIndex(sval i) const { return extinds!=NULL ? extinds->at(i) : i; }

        virtual vec3 getVertex(int i) const { return vecs->at(getIndex(i)); }
        virtual vec3 getNormal(int i) const { return vecs->at(getIndex(i),1); }
        virtual color getColor(int i) const { return cols->at(getIndex(i)); }
        virtual vec3 getUVWCoord(int i) const { return vecs->at(getIndex(i),3); }

        virtual sval numVertices() const { return numverts; }
        virtual bool hasNormal() const { return vecs->m()>1; }
        virtual bool hasColor() const { return cols!=NULL; }
        virtual bool hasUVWCoord() const { return vecs->m()>3; }
};

class MatrixIndexBuffer : public IndexBuffer
{
        IndexMatrix* indices;
        IndexMatrix* extinds;
        bool deleteMatrices;
        
public:
        MatrixIndexBuffer(IndexMatrix* indices,IndexMatrix* extinds=NULL) : indices(indices), extinds(extinds), deleteMatrices(false)
        {}
        
        MatrixIndexBuffer(const IndexBuffer* buf) throw(RenderException) : extinds(NULL), deleteMatrices(true)
        {
                if(!buf)
                        throw RenderException("IndexBuffer 'buf' must be provided");
                
                indices=new IndexMatrix("copyinds",buf->numIndices(),buf->indexWidth(0));
                
                for(sval i=0;i<buf->numIndices();i++)
                        for(sval j=0;j<_min(indices->m(),buf->indexWidth(i));j++)
                                indices->at(i,j)=buf->getIndex(i,j);
        }
        
        virtual ~MatrixIndexBuffer()
        {
                if(deleteMatrices){
                        SAFE_DELETE(indices);
                        SAFE_DELETE(extinds);
                }
        }

        virtual sval numIndices() const
        {
                if(indices==NULL)
                        return 0;
                else if(extinds!=NULL)
                        return sval(extinds->n());
                else
                        return sval(indices->n());
        }
        virtual sval indexWidth(int i) const { return sval(indices!=NULL ? indices->m() : 0); }
        virtual sval getIndex(int i,int j) const { return sval(indices->getAt(extinds!=NULL ? extinds->getAt(i) : i,j)); }
};

class Ray
{
        vec3 pos,dir,invdir;
        bool signx,signy,signz;
public:

        Ray()
        {}

        Ray(const vec3 &pos, const vec3 &dir): pos(pos)
        {
                setDirection(dir);
        }

        Ray(const Ray& r) : pos(r.pos)
        {
                setDirection(r.dir);
        }
        
        vec3 getPosition(real t=0) const { return pos+(dir*t); }

        vec3 getDirection() const { return dir; }

        void setPosition(const vec3 &v) { pos=v; }

        void setDirection(const vec3 &v)
        {
                if(v.isZero())
                        throw std::invalid_argument("Direction vector is zero length.");
                
                dir=v.norm();
                
                // store the direction inverse for bound box check efficiency
                invdir=dir.inv();
                
                // store component signs for bound box check efficiency
                signx=invdir.x()<0;
                signy=invdir.y()<0;
                signz=invdir.z()<0;
        }

        real distTo(const vec3 v) const { return dir.dot(v-pos); }
        
        bool onPlane(const vec3& planept, const vec3& planenorm) const 
        { 
            return pos.onPlane(planept,planenorm) && (pos+dir).onPlane(planept,planenorm); 
        }
        
        real intersectsPlane(const vec3 & planepos, const vec3 & planenorm) const
        {
                return planenorm.dot(planepos-pos)/planenorm.dot(dir);
        }
        
        realpair intersectsAABB(const vec3& minv, const vec3& maxv) const
        {
                real tmin, tmax, tymin, tymax, tzmin, tzmax;

                tmin  = ((signx ? maxv : minv).x() - pos.x()) * invdir.x();
                tmax  = ((signx ? minv : maxv).x() - pos.x()) * invdir.x();
                tymin = ((signy ? maxv : minv).y() - pos.y()) * invdir.y();
                tymax = ((signy ? minv : maxv).y() - pos.y()) * invdir.y();

                if ( (tmin > tymax) || (tymin > tmax) ) 
                        return realpair(-1,-1);

                if (tymin > tmin)
                        tmin = tymin;
                if (tymax < tmax)
                        tmax = tymax;

                tzmin = ((signz ? maxv : minv).z() - pos.z()) * invdir.z();
                tzmax = ((signz ? minv : maxv).z() - pos.z()) * invdir.z();

                if ( (tmin > tzmax) || (tzmin > tmax) ) 
                        return realpair(-1,-1);

                if (tzmin > tmin)
                        tmin = tzmin;
                if (tzmax < tmax)
                        tmax = tzmax;

                return realpair(tmin,tmax);//tmin < to && tmax > from;
        }
        
        realpair intersectsSphere(const vec3& center, real rad) const
        {
                real tca = distTo(center); // distance from `pos' to projection of `center' on the ray
                real thc=0; // offset from project of `center' on the ray where the sphere is intersected
        
                if(tca>0){ // if tca is positive, the sphere is in front of the ray, otherwise the result is where the sphere center projects on the ray
                        real r2=rad*rad;
                        vec3 L = center - pos;
                        real d2 = L.dot(L) - tca * tca;
                
                        if (d2 < r2) // ray intersects sphere at 2 points offset by thc, otherwise ray passes outside of sphere
                                thc = sqrt(r2 - d2);
                }               

                return realpair(tca - thc, tca + thc); // thc==0 if ray does not pass through sphere, so tca is projection distance either in front or behind
        }

        realpair intersectsRay(const Ray& ray) const
        {
                real t=0,s=0;
                vec3 p1=pos;
                vec3 p2=ray.getPosition();
                real t1=distTo(p2);
                real t2=ray.distTo(p1);
                vec3 pt1=getPosition(t1);
                vec3 pt2=ray.getPosition(t2);

                if(p2==pt1 || p1==pt2){ // check if either ray is pointing at the other
                        if(p2==pt1) // this was pointing at `ray'
                                t=t1;
                        if(p1==pt2) // `ray' was pointing at this
                                s=t2;
                }
                else{ 
                        vec3 norm=p1.planeNorm(p2,getPosition(1.0)); // calculate a normal based on the two origins at the point at t=1.0
                        if(ray.getPosition(1.0).onPlane(p1,norm)){ // if the point for `ray' at t=1.0 is on the plane then the rays intersect if they aren't parallel
                                vec3 rd=ray.getDirection();
                                real angle=dir.angleTo(rd);

                                if(angle>dEPSILON && angle<(dPI-dEPSILON)){ // rays are not parallel
                                        t=intersectsPlane(p2,rd.cross(norm)); // find where this intersects the plane defined by `ray' at right angle to the plane this and `ray' lie on
                                        s=ray.distTo(getPosition(t)); // find the intersect distance for `ray' by seeing how far away the interect point for this is
                                }
                        }
                }

                return realpair(t,s);
        }
        
        real intersectsLineSeg(const vec3& v1, const vec3& v2) const
        {
                real dist=intersectsPlane(v1,pos.planeNorm(v1,v2).cross(v2-v1));
                
                if(dist<0 || dist==realInf)
                    return -1;
                
                vec3 lpos=getPosition(dist);
                
                if(lpos==v1 || lpos==v2 || equalsEpsilon(lpos.lineDist(v1,v2),0))
                    return dist;
                
                //if(dist>=0 && dist<realInf && equalsEpsilon(lpos.lineDist(v1,v2),0))
                //              return dist;
                
                return -1;
        }

        realtriple intersectsTri(const vec3& v0, const vec3& v1, const vec3& v2) const 
        {
                vec3 e1=v1-v0; // triangle edge 1
                vec3 e2=v2-v0; // triangle edge 2
                vec3 p=dir.cross(e2); // direction perpendicular to ray line and edge 2
                real det=e1.dot(p); // determinant, will be 0 if edge 1 and p are perpendicular

                if(equalsEpsilon(det,0.0)) // ray is parallel with triangle's plane
                        return realtriple(-1,-1,-1);

                real invdet=1.0/det;

                vec3 t=pos-v0;
                real u=p.dot(t)*invdet;

                if(u < 0 || u > 1) // ray intersects triangle's plane outside the band parallel with e1
                        return realtriple(-1,-1,-1);

                vec3 q=t.cross(e1);
                real v=dir.dot(q)*invdet;

                if(v < 0 || u + v  > 1) // ray intersects triangle's plane outside the band parallel with e2 or outside the triangle area
                        return realtriple(-1,-1,-1);
 
                real len=e2.dot(q)*invdet;

                if(len>dEPSILON) // ray points away from triangle
                        return realtriple(len,u,v);
        
                return realtriple(-1,-1,-1);
        }
        
        std::vector<indextriple> intersectsTriMesh(const Vec3Matrix* const nodes, const IndexMatrix* const inds,
                const Vec3Matrix* const centers, const RealMatrix* const radii2, sval numResults=0,sval excludeInd=-1) const  throw(IndexException)
        {
                std::vector<indextriple> results;
                sval len=inds->n();
                
                if(centers && radii2)
                        len=_min(len,_min(centers->n(),radii2->n()));

                for(sval n=0;n<len && (numResults==0 || results.size()<numResults);n++){
                        if(n==excludeInd)
                                continue;

                        vec3 ncenter,npos,v0,v1,v2;
                        real nrad=0,cdist;
                        
                        if(centers && radii2){
                                ncenter=centers->at(n);
                                nrad=radii2->at(n);
                                //npos=getPosition(ncenter.distTo(pos));
        
                                //if(npos.distToSq(ncenter)>nrad)
                                //      continue;
                                
                                cdist=distTo(ncenter);
                                if((cdist*cdist)>nrad)
                                        continue;
                                
                                v0=nodes->getAt(inds->at(n,0));
                                v1=nodes->getAt(inds->at(n,1));
                                v2=nodes->getAt(inds->at(n,2));
                        }
                        else{
                                v0=nodes->getAt(inds->at(n,0));
                                v1=nodes->getAt(inds->at(n,1));
                                v2=nodes->getAt(inds->at(n,2));
                                
                                ncenter=(v0+v1+v2)/3.0;
                                nrad=_max(ncenter.distToSq(v0),_max(ncenter.distToSq(v1),ncenter.distToSq(v2)));
                                npos=getPosition(ncenter.distTo(pos));
        
                                if(npos.distToSq(ncenter)>nrad)
                                        continue;
                        }
                        
                        realtriple inter=intersectsTri(v0,v1,v2);
                        
                        if(inter.first>=0)
                                results.push_back(indextriple(indexval(n),inter));
                }
                
                return results;
        }
};

class transform
{
public:
        vec3 trans;
        vec3 scale;
        rotator rot;
        bool _isInverse;

        transform(const vec3& trans=vec3(),const vec3& scale=vec3(1.0),const rotator& rot=rotator(),bool isInv=false) : 
                        trans(trans),scale(scale),rot(rot),_isInverse(isInv)
        {}

        transform(const transform & t) : trans(t.trans),scale(t.scale),rot(t.rot),_isInverse(t.isInverse())
        {}

        transform(real x, real y, real z, real sx, real sy, real sz, real yaw, real pitch, real roll, bool isInv=false) :
                trans(x,y,z), scale(sx,sy,sz),rot(yaw,pitch,roll), _isInverse(isInv)
        {}

        vec3 getTranslation() const { return trans; }
        vec3 getScale() const { return scale; }
        rotator getRotation() const { return rot; }

        bool isInverse() const { return _isInverse; }

        void setTranslation(const vec3 &v){ trans=v; }
        void setScale(const vec3 &v) { scale=v; }
        void setRotation(const rotator &r) { rot=r; }

        vec3 operator * (const vec3 &v) const
        {
                if(_isInverse)
                        return scale*(rot*(v+trans));
                else
                        return trans+(rot*(v*scale));
        }

        friend vec3 operator * (const vec3 &v, const transform &t)
        {
                return t*v;
        }

        vec3 operator / (const vec3 &v) const
        {
                return inverse()*v;
        }

        friend vec3 operator / (const vec3 &v, const transform &t)
        {
                return t/v;
        }

        friend vec3& operator *= (vec3 &v, const transform &t)
        {
                vec3 vv=t*v;
                v.x(vv.x());
                v.y(vv.y());
                v.z(vv.z());
                return v;
        }
        
        Ray operator * (const Ray &r) const
        {
                transform t=*this;
                return Ray(t*r.getPosition(),t.directional()*r.getDirection());
        }
        
        friend Ray operator * (const Ray &r, const transform &t)
        {
                return t*r;
        }

        transform operator * (const transform& t) const
        {
                /*vec3 nscale=t.getScale()*scale;
                rotator nrot=rot*t.getRotation();
                vec3 ntrans=t.getTranslation()+(trans*t.directional());
                return transform(ntrans,nscale,nrot);*/
                transform _t=*this; 
                vec3 mincorner=_t*(t*vec3(0)); 
                vec3 maxcorner=_t*(t*vec3(1)); 
                vec3 xcorner=_t*(t*vec3(1,0,0)); 
                vec3 ycorner=_t*(t*vec3(0,1,0)); 
                rotator rot((xcorner-mincorner).norm(),(ycorner-mincorner).norm(),vec3(1,0,0),vec3(0,1,0));
                vec3 scale=rot/(maxcorner-mincorner);
                return transform(mincorner,scale,rot); 
        }

        bool operator == (const transform & t) const
        {
                return trans==t.trans && scale==t.scale && rot==t.rot && _isInverse==t._isInverse;
        }

        transform inverse() const
        {
                return transform(trans*-1,scale.inv(),rot.inverse(),!_isInverse);
        }

        transform directional() const
        {
                return transform(vec3(),scale,rot,_isInverse);
        }

        void toMatrix(real* mat) const
        {
                rot.toMatrix(mat);
                if(_isInverse){
                        mat[ 0]*=scale.x();
                        mat[ 1]*=scale.x();
                        mat[ 2]*=scale.x();
                        mat[ 4]*=scale.y();
                        mat[ 5]*=scale.y();
                        mat[ 6]*=scale.y();
                        mat[ 8]*=scale.z();
                        mat[ 9]*=scale.z();
                        mat[10]*=scale.z();
                        mat[ 3]=trans.dot(vec3(mat[0],mat[1],mat[2]));
                        mat[ 7]=trans.dot(vec3(mat[4],mat[5],mat[6]));
                        mat[11]=trans.dot(vec3(mat[8],mat[9],mat[10]));
                }
                else{
                        mat[ 0]*=scale.x();
                        mat[ 1]*=scale.y();
                        mat[ 2]*=scale.z();
                        mat[ 4]*=scale.x();
                        mat[ 5]*=scale.y();
                        mat[ 6]*=scale.z();
                        mat[ 8]*=scale.x();
                        mat[ 9]*=scale.y();
                        mat[10]*=scale.z();
                        mat[ 3]=trans.x();
                        mat[ 7]=trans.y();
                        mat[11]=trans.z();
                }
        }

        mat4 toMatrix() const
        {
                mat4 m;
                toMatrix((real*)m.m);
                return m;
        }

        friend std::ostream& operator << (std::ostream &out, const transform &t)
        {
                return out << "transform(" << t.getTranslation() << ", " << t.getScale() << ", " << t.getRotation() << ", " << (t.isInverse() ? "true" : "false") << ")";
        }
};

class Image
{
public:
        virtual ~Image() {} 

        virtual TextureFormat getFormat() const { return TF_UNKNOWN; }
        virtual sval getWidth() const { return 0; }
        virtual sval getHeight() const { return 0; }
        virtual sval getDepth() const { return 0; }
        virtual size_t getDataSize() const { return 0; }
        virtual u8* getData() { return 0; }
        virtual std::string encode(const std::string& format) { return ""; }
        virtual void fillRealMatrix(RealMatrix* mat) throw(IndexException) {}
        virtual void fillColorMatrix(ColorMatrix* mat) throw(IndexException) {} 
};

class Camera
{
public:
        Camera(){}
        virtual ~Camera() {}

        virtual const char* getName() const { return ""; }

        virtual real getAspectRatio() const { return 0; }
        
        virtual Ray* getProjectedRay(real x, real y, bool isAbsolute=true) const { return NULL; }

        virtual vec3 getPosition() const { return vec3();}
        virtual vec3 getLookAt() const { return vec3();}
        virtual rotator getRotation() const { return rotator(); }
        
        virtual mat4 getViewMatrix() const { return mat4(); }
        virtual mat4 getProjMatrix() const { return mat4(); }

        virtual vec3 getScreenPosition(vec3 pos) const { return vec3(); }

        virtual vec3 getWorldPosition(real x, real y, bool isAbsolute=true) const 
        {
                Ray *r=getProjectedRay(x,y,isAbsolute);
                vec3 pos=r->getPosition();
                delete r;
                return pos;
        }

        virtual void setPosition(const vec3 &v){}
        virtual void setLookAt(const vec3 &v) {}
        virtual void setUp(const vec3 & v){}
        virtual void setZUp(){}
        virtual void rotate(const rotator & r) {}
        virtual void setRotation(const rotator& r){}

        virtual void setNearClip(real dist) {}
        virtual void setFarClip(real dist) {}
        virtual void setVertFOV(real rads) {}
        virtual void setBGColor(const color & c) {}
        virtual void setAspectRatio(real rat) {}
        virtual void setViewport(real left=0.0f,real top=0.0f,real width=1.0f,real height=1.0f) {}
        virtual void setOrtho(bool isOrtho){}
        virtual void setWireframe(bool isWireframe){}
        virtual void setSecondaryCamera(bool selective){}

        virtual real getVertFOV() const { return 0; }
        
        virtual real getNearClip() const {return 0; }
        virtual real getFarClip() const {return 0; }

        virtual sval getWidth() const { return 0; }
        virtual sval getHeight() const { return 0; }

        virtual bool isPointInViewport(int x, int y) const { return false;}
        
        virtual bool isSecondaryCamera() { return false; }

        virtual void renderToFile(const std::string& filename,sval width,sval height, TextureFormat format=TF_RGB24,real stereoOffset=0.0) throw(RenderException) {}
        virtual void renderToStream(void* stream,sval width,sval height, TextureFormat format=TF_RGB24,real stereoOffset=0.0) throw(RenderException) {}
        virtual Image* renderToImage(sval width,sval height, TextureFormat format=TF_RGB24,real stereoOffset=0.0) throw(RenderException) { return 0; }
};

class Figure
{
public:
        virtual ~Figure(){}

        virtual const char* getName() {return "";}
        virtual void setPosition(const vec3 &v) {}
        virtual void setScale(const vec3 &v) {}
        virtual void setRotation(const rotator& r){}
        virtual void setTransform(const transform &t) { setTransform(t.trans,t.scale,t.rot); }
        virtual void setTransform(const vec3 &trans,const vec3 &scale, const rotator &rot)
        {
                setPosition(trans);
                setRotation(rot);
                setScale(scale);
        }

        virtual vec3 getPosition(bool isDerived=false) const { return vec3(); }
        virtual vec3 getScale(bool isDerived=false) const { return vec3(); }
        virtual rotator getRotation(bool isDerived=false) const { return rotator(); }
        virtual transform getTransform(bool isDerived=false) const { return transform(getPosition(isDerived),getScale(isDerived),getRotation(isDerived),false); }

        virtual void setMaterial(const char* mat) throw(RenderException) {}
        virtual void setMaterial(const Material *mat) throw(RenderException) { setMaterial(mat->getName()); }

        virtual const char* getMaterial() const { return ""; }

        virtual std::pair<vec3,vec3> getAABB() const { return std::pair<vec3,vec3>(vec3(),vec3()); }

        virtual void fillData(const VertexBuffer* vb, const IndexBuffer* ib,bool deferFill=false,bool doubleSided=false) throw(RenderException) {}
        
        virtual void setVisible(bool isVisible){}
        virtual bool isVisible() const { return false;}

        virtual bool isTransparent() const { return false; }
        virtual bool isOverlay() const { return false; }

        virtual void setTransparent(bool isTrans){}
        virtual void setOverlay(bool isOverlay){}

        virtual void setRenderQueue(sval queue){}
        virtual sval getRenderQueue() const { return 0; }

        virtual void setCameraVisibility(const Camera* cam, bool isVisible){}

        virtual void setParent(Figure *fig){}
};

class BBSetFigure : public Figure
{
public:
        virtual ~BBSetFigure(){}

        virtual void setDimension(real width, real height){}

        virtual real getWidth() const { return 0;}
        virtual real getHeight() const { return 0;}

        virtual void setUpVector(const vec3& v){}

        virtual int numBillboards() const {return 0;}

        virtual void fillData(const VertexBuffer* vb, const IndexBuffer* ib,bool deferFill=false,bool doubleSided=false) throw(RenderException) {}

        virtual void setBillboardPos(indexval index, const vec3& pos) throw(IndexException) {}
        virtual void setBillboardDir(indexval index, const vec3& dir) throw(IndexException) {}
        virtual void setBillboardColor(indexval index, const color& col) throw(IndexException) {}
};

class TextureVolumeFigure : public Figure
{
public:
        virtual ~TextureVolumeFigure() {}
        
        virtual void setNumPlanes(sval num){}
        virtual sval getNumPlanes() const { return 0; }
        
        virtual void setAlpha(real a) {}
        virtual real getAlpha() const { return 0; }

        virtual void setTexAABB(const vec3& minv, const vec3& maxv) {}

        virtual void setAABB(const vec3& minv, const vec3& maxv) {}
        //virtual void setTexCoordDir(bool invertX, bool invertY, bool invertZ) {}

        virtual vec3 getTexXiPos(vec3 pos) const { return vec3();}

        virtual vec3 getTexXiDir(vec3 pos) const { return vec3();}

        virtual sval getPlaneIntersects(vec3 planept, vec3 planenorm,vec3 buffer[6][2],bool transformPlane=false,bool isXiPoint=false)
        {
                return 0;
        }
};

class GlyphFigure : public Figure
{
public:
        virtual ~GlyphFigure() {}

        virtual void setGlyphScale(vec3 v) {}
        virtual vec3 getGlyphScale() const { return vec3(); }
        virtual void setGlyphName(const std::string& name) {}
        virtual std::string getGlyphName() const {return ""; }
        virtual void addGlyphMesh(const std::string& name,const Vec3Matrix* nodes, const Vec3Matrix* norms, const IndexMatrix* inds) {}
};

class RibbonFigure : public Figure
{
public:
        virtual ~RibbonFigure() {}
        virtual void setOrientation(const vec3& orient) {}
        virtual bool isCameraOriented() const { return true; }
        virtual vec3 getOrientation() const { return vec3(); }

        virtual void setNumRibbons(sval num) {}
        virtual sval numRibbons() const { return 0; }
        virtual sval numNodes(sval ribbon) const throw(IndexException) { return 0; }
        virtual void setMaxNodes(sval num) {}
        virtual sval getMaxNodes() const { return 0; }
        
        virtual void clearRibbons() {}
        virtual void removeRibbon(sval ribbon) throw(IndexException) {}
        virtual void removeNode(sval ribbon) throw(IndexException) {}
        virtual void addNode(sval ribbon,const vec3& pos, const color& col,real width, const rotator& rot=rotator(), real tex=0.0) throw(IndexException) {}
        virtual void setNode(sval ribbon,sval node,const vec3& pos, const color& col,real width, const rotator& rot=rotator(), real tex=0.0) throw(IndexException) {}
        virtual vec3 getNode(sval ribbon,sval node) throw(IndexException) { return vec3(); }
        virtual quadruple<color,real,rotator,real> getNodeProps(sval ribbon,sval node) throw(IndexException) { return quadruple<color,real,rotator,real>(); }
};

class TextFigure : public Figure
{
public:
        virtual ~TextFigure() {}
        virtual void setText(const std::string& text) {}
        virtual void setFont(const std::string& fontname) {}
        virtual void setColor(const color& col) {}
        
        virtual void setVAlign(VAlignType align){}
        virtual void setHAlign(HAlignType align){}
        virtual void setTextHeight(real height){}
        virtual void setSpaceWidth(real width) {}
        
        virtual void setCameraAlign(bool align) {}
        
        virtual std::string getText() const { return "";}
        virtual std::string getFont() const { return "";}
        virtual color getColor() const { return color(); }
        
        virtual VAlignType getVAlign() const { return V_CENTER; }
        virtual HAlignType getHAlign() const { return H_CENTER; }
        virtual real getTextHeight() const { return 0; }
        virtual real getSpaceWidth() const { return 0; }
        
        virtual bool isCameraAligned() const { return false; }
};

class Config
{
protected:
        typedef std::pair<std::string,std::string> strpair;
        typedef std::map<strpair,std::string> configmap;

        configmap map;

public:
        Config(){}

        void set(const char* group, const char* name, const char* value) { map[getPair(group,name)]=value; }
        void set(const char* name, const char* value) { map[getPair("",name)]=value; }
        bool hasValue(const char* group, const char* name) const { return map.find(getPair(group,name))!=map.end(); }
        bool hasValue(const char* name) const { return hasValue("",name); }
        const char* get(const char* group, const char* name) { return hasValue(group,name) ? map[getPair(group,name)].c_str() : ""; }
        const char* get(const char* name) { return get("",name); }

        std::string toString()
        {
                std::ostringstream out;
                for(configmap::const_iterator i=map.begin();i!=map.end();i++)
                        out << "(" <<(*i).first.first << ", " << (*i).first.second << ") = " << (*i).second << std::endl;

                return out.str();
        }

        friend std::ostream& operator << (std::ostream &out, const Config *c)
        {
                for(configmap::const_iterator i=c->map.begin();i!=c->map.end();i++)
                        out << "(" <<(*i).first.first << ", " << (*i).first.second << ") = " << (*i).second << std::endl;

                return out;
        }
        
private:
        strpair getPair(const char* group, const char* name) const
        {
                std::string sgroup=group,sname=name;
                std::transform(sgroup.begin(), sgroup.end(), sgroup.begin(), ::tolower);
                std::transform(sname.begin(), sname.end(), sname.begin(), ::tolower);
                return strpair(sgroup,sname);
        }
};

class RenderScene
{
        bool renderHighQuality;
        bool alwaysHighQuality;

public:
        RenderScene() : renderHighQuality(false), alwaysHighQuality(false) {}

        virtual ~RenderScene() {}

        virtual Camera* createCamera(const char* name,real left=0.0f,real top=0.0f,real width=1.0f,real height=1.0f) throw(RenderException) { return NULL; }
        
        virtual void setAmbientLight(const color & c) {}

        virtual void addResourceDir(const char* dir) {}
        
        virtual void initializeResources() {}

        virtual Material* createMaterial(const char* name) throw(RenderException) {return NULL;}

        virtual Figure* createFigure(const char* name, const char* mat,FigureType type) throw(RenderException) {return NULL;}

        virtual Light* createLight() throw(RenderException) { return NULL; }

        virtual Image* loadImageFile(const std::string &filename) throw(RenderException) { return NULL; }

        virtual Texture* createTexture(const char* name,sval width, sval height, sval depth, TextureFormat format) throw(RenderException) { return NULL; }

        virtual Texture* loadTextureFile(const char* name,const char* absFilename) throw(RenderException) { return NULL; }
        
        virtual GPUProgram* createGPUProgram(const char* name,ProgramType ptype,const char* language) throw(RenderException) { return NULL; }

        virtual void saveScreenshot(const char* filename,Camera* c=NULL,int width=0,int height=0,real stereoOffset=0.0,TextureFormat tf=TF_RGB24) throw(RenderException) {}

        virtual Config* getConfig() const {return NULL; }

        virtual void logMessage(const char* msg){}
        
        virtual void setBGObject(color col,bool enabled){}

        void setRenderHighQuality(bool val) { renderHighQuality=val; }
        
        void setAlwaysHighQuality(bool val) { alwaysHighQuality=val; }

        bool getRenderHighQuality() const { return renderHighQuality || alwaysHighQuality;}
        
        bool getAlwaysHighQuality() const  { return alwaysHighQuality; }
};

class RenderAdapter
{
public:
        virtual ~RenderAdapter(){}

        virtual u64 createWindow(int width, int height)  throw(RenderException) { return 0; }
        virtual void paint(){}
        virtual void resize(int x, int y,int width, int height){}
        virtual RenderScene* getRenderScene() { return NULL; }
};

RenderAdapter* getRenderAdapter(Config* config) throw(RenderException);

/*****************************************************************************************************************************/
/* Algorithms */
/*****************************************************************************************************************************/

template<typename T,typename V> void setMatrixMinMax(Matrix<T>* mat,const V& minv,const V& maxv)
{
        std::ostringstream os;
        os << V(minv);
        mat->meta("min",os.str().c_str());
        os.str("");
        os << V(maxv);
        mat->meta("max",os.str().c_str());
}

template<typename T> struct StrConvert { static T conv(const char* str) { return T(); } };
template<> struct StrConvert<real> { static real conv(const char* str) { return atof(str); } };
template<> struct StrConvert<indexval> { static indexval conv(const char* str) { return atol(str); } };

template<typename T> struct ParseLine 
{ 
        static void parse(const char* line,sval numvals,T* list) 
        {
                char *buf=new char[strlen(line)+1];
                strcpy(buf,line);
                        
                char* p=strtok(buf," \t\r\n");
        
                for(sval x=0;x<numvals && p;x++){
                        list[x]=StrConvert<T>::conv(p);
                        p=strtok(NULL," \t\r\n");
                }
                        
                delete buf;
        } 
};

template<> struct ParseLine<vec3>
{ 
        static void parse(const char* line,sval numvals,vec3* list) 
        {
                char *buf=new char[strlen(line)+1];
                strcpy(buf,line);
                        
                char* p=strtok(buf," \t\r\n");
        
                for(sval i=0;i<numvals && p;i++){
                        real x=StrConvert<real>::conv(p);
                        p=strtok(NULL," \t\r\n");
                        real y=(p ? StrConvert<real>::conv(p) : 0);
                        p=strtok(NULL," \t\r\n");
                        real z=(p ? StrConvert<real>::conv(p) : 0);
                        p=strtok(NULL," \t\r\n");
                        list[i]=vec3(x,y,z);
                }
                        
                delete buf;
        } 
};

template<typename T>
void readTextFileMatrix(const std::string & filename, sval numHeaders, Matrix<T>* mat)
{
        std::ifstream in(filename.c_str());

        if(!in)
                return;

        std::string line;
        sval numvals=sval(mat->m());
        T* entry=new T[numvals];

        if(numHeaders>0){
                sval* header=new sval[numHeaders];
                std::getline(in,line);  
                ParseLine<sval>::parse(line.c_str(),numHeaders,header);
                delete header;
        }

        while(std::getline(in,line)){           
                if(line.find_first_not_of(" \n\t\r")==std::string::npos)
                        continue;

                ParseLine<T>::parse(line.c_str(),numvals,entry);

                mat->append(entry[0]);
                sval pos=mat->n()-1;
                for(sval x=1;x<numvals;x++)
                        mat->at(pos,x)=entry[x];
        }

        delete entry;
}


//template<typename T> void convertStreamToRealMatrix(const T* stream, RealMatrix* mat);
template<typename T>
void convertStreamToRealMatrix(const T* stream, RealMatrix* mat)
{
        T minval=stream[0];
        T maxval=minval;
        sval mn=mat->n(),mm=mat->m();

        for(sval n=0;n<mn;n++)
                for(sval m=0;m<mm;m++){
                        T val=stream[n*mm+m]; // totally unsafe if stream is too short
                        minval=_min(val,minval);
                        maxval=_max(val,maxval);

                        mat->at(n,m)=val;
                }

        setMatrixMinMax<real,real>(mat,minval,maxval);
}

void convertUByteStreamToRealMatrix(const char* stream, RealMatrix* mat);
void convertUShortStreamToRealMatrix(const char* stream, RealMatrix* mat);
void convertByteStreamToRealMatrix(const char* stream, RealMatrix* mat);
void convertShortStreamToRealMatrix(const char* stream, RealMatrix* mat);
void convertUIntStreamToRealMatrix(const char* stream, RealMatrix* mat);
void convertIntStreamToRealMatrix(const char* stream, RealMatrix* mat);
void convertFloatStreamToRealMatrix(const char* stream, RealMatrix* mat);
void convertRealStreamToRealMatrix(const char* stream, RealMatrix* mat);
//void convertRGBA32StreamToRealMatrix(const u8* stream, RealMatrix* mat);

//RealMatrix* readImageFile(const std::string& filename);

std::pair<vec3,vec3> calculateBoundBox(const Vec3Matrix* mat);


realtriple calculateTriPlaneSlice(const vec3& planept, const vec3& planenorm, const vec3& a, const vec3& b, const vec3& c);
real calculateLinePlaneSlice(const vec3& planept, const vec3& planenorm, const vec3& a, const vec3& b);

void calculateTetValueIntersects(real val, real a, real b, real c, real d, real* results);

sval calculateHexValueIntersects(real val,const real* vals,intersect* results);

void basis_Tet1NL(real xi0, real xi1, real xi2, real* coeffs);

void basis_Hex1NL(real xi0, real xi1, real xi2, real* coeffs);

real basis_n_NURBS(sval ctrlpt,sval degree, real xi,const RealMatrix *knots);

void basis_NURBS_default(real u, real v, real w,sval ul, sval vl, sval wl, sval udegree, sval vdegree, sval wdegree,real* coeffs);

void catmullRomSpline(real t, real* coeffs);

real mat4Det(real a, real b, real c, real d, real e, real f, real g, real h, real i, real j, real k, real l, real m, real n, real o, real p);

bool pointInTet(vec3 pt, vec3 n1, vec3 n2, vec3 n3, vec3 n4);

bool pointInHex(vec3 pt, vec3 n1, vec3 n2, vec3 n3, vec3 n4, vec3 n5, vec3 n6, vec3 n7, vec3 n8);

inline float calculateTetVolume(vec3 a, vec3 b, vec3 c, vec3 d)
{
        return -mat4Det(a.x(),b.x(),c.x(),d.x(),a.y(),b.y(),c.y(),d.y(),a.z(),b.z(),c.z(),d.z(),1,1,1,1)/6.0;
}

vec3 pointSearchLinTet(vec3 pt, vec3 n1, vec3 n2, vec3 n3, vec3 n4);

vec3 pointSearchLinHex(vec3 pt, vec3 n1, vec3 n2, vec3 n3, vec3 n4, vec3 n5, vec3 n6, vec3 n7, vec3 n8);

template<typename T>
void cubicInterpMatrices(real t,const Matrix<T>* v1,const Matrix<T>* v2,const Matrix<T>* m1,const Matrix<T>* m2,const Matrix<T>* result)
{
        sval rows=_min(v1->n(),_min(v2->n(),_min(m1->n(),_min(m2->n(),result->n()))));
        sval cols=_min(v1->m(),_min(v2->m(),_min(m1->m(),_min(m2->m(),result->m()))));
        real coeffs[4];
        catmullRomSpline(t,coeffs);

        for(sval i=0;i<rows;i++)
                for(sval j=0;j<cols;j++){
                        T a=coeffs[0]*v1->at(i,j);
                        T b=coeffs[1]*v2->at(i,j);
                        T c=coeffs[2]*m1->at(i,j);
                        T d=coeffs[3]*m2->at(i,j);
                        result->at(i,j)=a+b+c+d;
                }
}

template<typename T>
quadruple<int,int,int,int> calculateBoundSquare(const Matrix<T>* const mat,const T& threshold)
{
        int minx=-1,maxx=-1,miny=-1,maxy=-1;

        sval rows=mat->n();
        sval cols=mat->m();

        for(sval i=0;i<rows && miny<0;i++)
                for(sval j=0;j<cols;j++)
                        if(mat->at(i,j)>threshold){
                                miny=int(i);
                                break;
                        }

        for(sval i=rows;i>0 && maxy<0;i--)
                for(sval j=0;j<cols;j++)
                        if(mat->at(i-1,j)>threshold){
                                maxy=int(i-1);
                                break;
                        }

        for(sval j=0;j<cols && minx<0;j++)
                for(sval i=0;i<rows;i++)
                        if(mat->at(i,j)>threshold){
                                minx=int(j);
                                break;
                        }

        for(sval j=cols;j>0 && maxx<0;j--)
                for(sval i=0;i<rows;i++)
                        if(mat->at(i,j-1)>threshold){
                                maxx=int(j-1);
                                break;
                        }

        return quadruple<int,int,int,int>(minx,miny,maxx,maxy);
}

template<typename T>
sval countValuesInRange(const Matrix<T> *mat, const T& minv, const T& maxv)
{
        sval count=0, rows=mat->n(), cols=mat->m();

        for(sval i=0;i<rows;i++)
                for(sval j=0;j<cols;j++){
                        sval val=mat->at(i,j);
                        if(val>=minv && val<=maxv)
                                count++;
                }

        return count;
}

template<typename T>
std::vector<vec3> findBoundaryPoints(const Matrix<T>* mat, const T& threshold)
{
        std::vector<vec3> result;
        
        sval rows=mat->n();
        sval cols=mat->m();
        
        for(sval i=0;i<rows;i++)
                for(sval j=0;j<cols;j++){
                        T val=mat->atc(i,j);
                        if(val<threshold)
                                continue;
                        
                        bool allInternal=true;
                        for(sval n=_max<sval>(0,i-1);allInternal && n<_min(rows,i+1);n++)
                                for(sval m=_max<sval>(0,j-1);allInternal && m<_min(cols,j+1);m++)
                                        if(n!=i || m!=j)
                                                allInternal=allInternal && mat->getAt(n,m)>=threshold;
                                
                        if(!allInternal)
                                result.push_back(vec3(i,j,0));
                }
                
        return result;
}

template<typename T>
T sumMatrix(const Matrix<T> *mat)
{
        T result=T();
        sval rows=mat->n();
        sval cols=mat->m();
        
        for(sval i=0;i<rows;i++)
                for(sval j=0;j<cols;j++)
                        result+=mat->atc(i,j);

        return result;
}

template<typename T>
std::pair<T,T> minmaxMatrix(const Matrix<T>* mat) throw(ValueException)
{
        sval rows=mat->n();
        sval cols=mat->m();
        
        if(rows==0)
                throw ValueException("mat","Matrix must be non-empty");

        std::pair<T,T> result(mat->atc(0,0),mat->atc(0,0));

        for(sval i=0;i<rows;i++)
                for(sval j=0;j<cols;j++){
                        T val=mat->atc(i,j);
                        if(val<result.first)
                                result.first=val;
                        else if(val>result.second)
                                result.second=val;
                }

        return result;
}

template<typename T>
T bilerpMatrix(const Matrix<T> *mat,real x, real y) throw(ValueException)
{
        if(mat->n()==0)
                throw ValueException("mat","Matrix must be non-empty");
        
        if(x<0.0 || x>1.0 || y<0.0 || y>1.0)
                return T();
                
        x*=mat->m()-1;
        y*=mat->n()-1;
                
        sval sx=(sval)floor(x);
        sval sy=(sval)floor(y);

        real dx=x-sx;
        real dy=y-sy;
        real dx1=1.0-dx;
        real dy1=1.0-dy;

        return dx*(dy*mat->atc(sy+1,sx+1)+dy1*mat->atc(sy,sx+1)) + dx1*(dy*mat->atc(sy+1,sx)+dy1*mat->atc(sy,sx));
}

template<typename T>
T trilerpMatrices(const Matrix<T>* mat1, const Matrix<T>* mat2, vec3 v1, vec3 v2) throw(ValueException)
{
        T val1=bilerpMatrix<T>(mat1,v1.x(),v1.y());
        T val2=bilerpMatrix<T>(mat2,v2.x(),v2.y());
        
        real absz=fabs(v1.z());
        real lerpval=lerpXi<real>(absz,0,absz+fabs(v2.z()));
        
        return lerp(lerpval,val1,val2);
}

inline vec3 getPlaneXi(const vec3& pos, const vec3& planepos, const rotator& orientinv, const vec3& dimvec)
{
        return (orientinv*(pos-planepos))/vec3(dimvec.x(),dimvec.y(),1);
}

void interpolateImageStack(const std::vector<RealMatrix*>& stack,const transform& stacktransinv,RealMatrix *out,const transform& outtrans);

real getImageStackValue(const std::vector<RealMatrix*>& stack,const vec3& pos);

void calculateImageHistogram(const RealMatrix* img, RealMatrix* hist, i32 minv); 

vec3* calculateTriNorms(vec3* nodes, sval numnodes, indexval* inds, sval numinds);

} // namespace RenderTypes
#endif // RENDERTYPES_H