gcg.h

Go to the documentation of this file.

#ifndef _GCG_H_

#ifdef __cplusplus
  extern "C++" {    // All names have C++ style
#else
  #error "GCC library must be compiled with C++ compilers"
#endif

#if !defined(_COUNTER_GCG_)
  #define _COUNTER_GCG_ 1
#endif

// Next definitions will be included once
#if _COUNTER_GCG_ == 1

#if defined(GCG_WINDOWS32) || defined(GCG_WINDOWS64)
  #ifndef VAPIWIN_API
    // GCGlib is deployed as a Dynamic Linked Library for Windows
    #define GCG_API_FUNCTION extern __declspec(dllexport)
    #define GCG_API_DATA     extern __declspec(dllexport)
    #define GCG_API_CLASS    __declspec(dllexport)
    #define GCG_API_TEMPLATE __declspec(dllexport)
  #else
    // VAPIWIN is a Dynamic Linked Library for use with non-Visual C++ compilers for video input
    #define GCG_API_FUNCTION extern
    #define GCG_API_DATA extern
    #define GCG_API_CLASS
    #define GCG_API_TEMPLATE
  #endif
#else
  #ifdef WIN32
    // GCGlib names need special modifiers in Windows
    #define GCG_API_FUNCTION extern __declspec(dllimport)
    #define GCG_API_DATA     extern __declspec(dllimport)
    #define GCG_API_CLASS
    #define GCG_API_TEMPLATE
  #else
    // GCGlib names defined as C and C++ standards
    #define GCG_API_FUNCTION extern
    #define GCG_API_DATA     extern
    #define GCG_API_CLASS
    #define GCG_API_TEMPLATE
  #endif
#endif

#include <stdio.h>
#include <math.h>
#include <new>

#if defined(__GNUC__) || defined(MINGW)
  #include <stdint.h>
#endif

// MSVC compatible 64 bits integer
#ifdef _MSC_VER
  typedef __int64   int64_t;
  typedef unsigned __int64  uint64_t;
#endif



GCG_API_DATA    const char *GCG_VERSION;
GCG_API_DATA    const char *GCG_BUILD_DATE;
GCG_API_DATA    const int GCG_VERSION_MAJOR;
GCG_API_DATA    const int GCG_VERSION_MINOR;
GCG_API_DATA    const int GCG_VERSION_BUILD;




/*******************************************************************
  GENERAL MACROS AND INLINE FUNCTIONS
********************************************************************/



#ifndef SQR
static inline int SQR(int x) { return x * x; }

static inline short SQR(short x) { return x * x; }

static inline long SQR(long x) { return x * x; }

static inline float SQR(float x) { return x * x; }

static inline double SQR(double x) { return x * x; }
#endif


#ifndef MIN
#define MIN(a, b) ((a > b) ? (b) : (a))
#endif


#ifndef MAX
#define MAX(a, b) ((a < b) ? (b) : (a))
#endif


#ifndef ABS
#define ABS(a) (((a) < 0) ? -(a) : (a))
#endif



#ifndef SWAP
#define SWAP(a, b) {(a) ^= (b); (b) ^= (a);  (a) ^= (b);}
#endif


#ifndef DEG2RAD
#define DEG2RAD(a)  (((a) * M_PI) / 180.0)
#endif


#ifndef RAD2DEG
#define RAD2DEG(a)  (((a) * 180) / M_PI)
#endif


#ifndef TRUE
#define TRUE  1
#endif



#ifndef FALSE
#define FALSE 0
#endif


#ifndef EPSILON
#define EPSILON 0.0000077
#endif


#ifndef INF
#define INF 1.0e+16
#endif


#ifndef FEQUAL
#define FEQUAL(a, b) (fabs((a)-(b)) < EPSILON)
#endif


#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif


#ifndef NULL
#define NULL 0
#endif


#ifndef ISPOWEROF2
#define ISPOWEROF2(n) (!(n & (n-1)))
#endif





GCG_API_FUNCTION    bool gcgReport(int code, const char *format = NULL, ...);

GCG_API_FUNCTION    int gcgGetReport(int sizestring = 0, char *extrastring = NULL);

GCG_API_FUNCTION    const char *gcgGetReportString(char **typestr = NULL, char **domainstr = NULL);

GCG_API_FUNCTION    bool gcgEnableReportLog(int code, bool enable);




GCG_API_FUNCTION    FILE *gcgGetLogStream();

GCG_API_FUNCTION    FILE *gcgSetLogStream(FILE *logstream);

GCG_API_FUNCTION    void gcgOutputLog(const char *format, ...);




#define GCG_REPORT_CODE(t, d, m) (((((int) d) & (int) 0xff) << 3) | ((((int) m) & (int) 0x3ff) << 11) | ((((int) t) & (int) 0x07) << 0))

#define GCG_REPORT_TYPE(c) ((((int) c) >> 0) & (int) 0x07)      // 3 bits

#define GCG_REPORT_DOMAIN(c) ((((int) c) >> 3) & (int) 0xff)    // 8 bits

#define GCG_REPORT_MESSAGE(c) ((((int) c) >> 11) & (int) 0x3ff)  // 10 bits




#define GCG_SUCCESS                 0 // Must be zero
#define GCG_ERROR                   1
#define GCG_WARNING                 2
#define GCG_INFORMATION             3 // Must be the greatest code


#define GCG_DOMAIN_MEMORY             8 // Must be the smallest code
#define GCG_DOMAIN_NETWORK            9
#define GCG_DOMAIN_FILE               10
#define GCG_DOMAIN_IMAGE              11
#define GCG_DOMAIN_CAMERA             12
#define GCG_DOMAIN_VIDEO              13
#define GCG_DOMAIN_THREAD             14
#define GCG_DOMAIN_GRAPHICS           15
#define GCG_DOMAIN_GEOMETRY           16
#define GCG_DOMAIN_ALGEBRA            17
#define GCG_DOMAIN_CONFIG             18
#define GCG_DOMAIN_DATASTRUCTURE      19
#define GCG_DOMAIN_SIGNAL             20

#define GCG_DOMAIN_USER               21   // Must be the greatest code

#define GCG_REPORT_ALL                -1




#define GCG_NONE                    0     // Must be the smallest code
#define GCG_ALLOCATIONERROR         1
#define GCG_INVALIDOBJECT           2
#define GCG_BADPARAMETERS           3
#define GCG_INITERROR               4
#define GCG_INTERNALERROR           5
#define GCG_TIMEOUT                 6
#define GCG_IOERROR                 7
#define GCG_FORMATMISMATCH          8
#define GCG_UNRECOGNIZEDFORMAT      9
#define GCG_UNSUPPORTEDFORMAT       10
#define GCG_UNSUPPORTEDOPERATION    11
#define GCG_INFORMATIONPROBLEM      12
#define GCG_READERROR               13
#define GCG_WRITEERROR              14
#define GCG_NOTCONNECTED            15
#define GCG_SOCKETERROR             16
#define GCG_SELECTERROR             17
#define GCG_CONNECTIONCLOSED        18
#define GCG_HOSTNOTFOUND            19
#define GCG_CONNECTERROR            20
#define GCG_SETSOCKETOPTERROR       21
#define GCG_BINDERROR               22
#define GCG_LISTENERROR             23
#define GCG_ACCEPTERROR             24
#define GCG_EXCEPTION               25
#define GCG_ACCEPTBADADDRESS        26
#define GCG_CONNECTIONREFUSED       27
#define GCG_CONNECTIONABORTED       28
#define GCG_WRITEBADADDRESS         29
#define GCG_READBADADDRESS          30
#define GCG_SENDSHUTDOWN            31
#define GCG_RECEIVESHUTDOWN         32
#define GCG_INVALIDSOCKET           33
#define GCG_STREAMCONTROLERROR      34
#define GCG_ALREADYRUNNING          35
#define GCG_CREATIONFAILED          36
#define GCG_MUTEXFAILURE            37
#define GCG_SEMAPHOREFAILURE        38
#define GCG_NOBUFFER                39
#define GCG_JOBNOTEXECUTED          40
#define GCG_INTERRUPTED             41
#define GCG_DISCONNECTED            42
#define GCG_STOPERROR               43
#define GCG_NOTRUNNING              44
#define GCG_EXITFAILURE             45
#define GCG_RELEASEERROR            46
#define GCG_LOCKERROR               47
#define GCG_OWNERSHIPERROR          48
#define GCG_REMOVALERROR            49
#define GCG_STOREFAILURE            50
#define GCG_OPENERROR               51
#define GCG_DECODEERROR             52
#define GCG_FUNCTIONCALLFAILED      53
#define GCG_GRAPHICSERROR           54
#define GCG_INCOMPATIBILITYERROR    55
#define GCG_CONFIGERROR             56
#define GCG_UNAVAILABLERESOURCE     57
#define GCG_STARTERROR              58
#define GCG_OUTOFBOUNDS             59
#define GCG_ENDOFSTREAM             60
#define GCG_OPERATIONFAILED         61
#define GCG_MEMORYLEAK              62
#define GCG_INSERTIONFAILED         63
#define GCG_ADAPTING                64

#define GCG_UNDEFINED               65 // Must be the greatest code




class    GCG_API_CLASS    gcgCLASS {
  public:
    void *operator new(size_t size);

    void *operator new(size_t size, const std::nothrow_t&) throw();

    void *operator new[](size_t size);

    void *operator new[](size_t size, const std::nothrow_t&) throw();

    void operator delete(void *p);

    void operator delete(void *p, const std::nothrow_t&) throw();

    void operator delete[](void *p);

    void operator delete[](void *p, const std::nothrow_t&) throw();
};


GCG_API_FUNCTION    uintptr_t getAllocatedMemoryCounter();




// GCG primitive types.

typedef float VECTOR2[2];
typedef float VECTOR3[3];
typedef float VECTOR4[4];

typedef float MATRIX2[4];
typedef float MATRIX3[9];
typedef float MATRIX34[12];
typedef float MATRIX4[16];

typedef double VECTOR2d[2];
typedef double VECTOR3d[3];
typedef double VECTOR4d[4];

typedef double MATRIX2d[4];
typedef double MATRIX3d[9];
typedef double MATRIX34d[12];
typedef double MATRIX4d[16];


// GCG algebra macros. More efficient than functions. But implies in a larger code.

// Expressions ///////////////////////////////////////////////
#define gcgDOTVECTOR2(v1, v2) ((v1)[0] * (v2)[0] + (v1)[1] * (v2)[1])
#define gcgDOTIMEDVECTOR2(u, a, b) ((u)[0] * (a) + (u)[1] * (b)}
#define gcgLENGTHVECTOR2(v) (sqrt(gcgDOTVECTOR2((v),(v))))

#define gcgDOTVECTOR3(v1, v2) ((v1)[0] * (v2)[0] + (v1)[1] * (v2)[1] + (v1)[2] * (v2)[2])
#define gcgDOTIMEDVECTOR3(u, a, b, c) ((u)[0] * (a) + (u)[1] * (b) + (u)[2] * (c))
#define gcgLENGTHVECTOR3(v) (sqrt(gcgDOTVECTOR3((v),(v))))

#define gcgDOTVECTOR4(v1, v2) ((v1)[0] * (v2)[0] + (v1)[1] * (v2)[1] + (v1)[2] * (v2)[2] + (v1)[3] * (v2)[3])
#define gcgDOTIMEDVECTOR4(u, a, b, c, d) ((u)[0] * (a) + (u)[1] * (b) + (u)[2] * (c) + (u)[3] * (d))
#define gcgLENGTHVECTOR4(v) (sqrt(gcgDOTVECTOR4((v),(v))))

#define gcgISEQUALVECTOR2(v1, v2) (fabs(v1[0] - v2[0]) > EPSILON && fabs((v1)[1] - (v2)[1]) > EPSILON)
#define gcgISEQUALVECTOR3(v1, v2) (fabs(v1[0] - v2[0]) > EPSILON && fabs((v1)[1] - (v2)[1]) > EPSILON && fabs((v1)[2] - (v2)[2]) > EPSILON)
#define gcgISEQUALVECTOR4(v1, v2) (fabs(v1[0] - v2[0]) > EPSILON && fabs((v1)[1] - (v2)[1]) > EPSILON && fabs((v1)[2] - (v2)[2]) > EPSILON && fabs((v1)[3] - (v2)[3]) > EPSILON)

// Computes the determinant of matrices
#define gcgDETERMINANT3x3_(m0, m1, m2, m3, m4, m5, m6, m7, m8)  ((m0)*((m8)*(m4)-(m7)*(m5)) - (m3)*((m8)*(m1)-(m7)*(m2)) + (m6)*((m5)*(m1)-(m4)*(m2)))
#define gcgDETERMINANT4x4_(m, det0, det4, det8, det12) \
( (m)[0]  * (det0 = gcgDETERMINANT3x3_((m)[5], (m)[6], (m)[7], (m)[9], (m)[10], (m)[11], (m)[13], (m)[14], (m)[15])) \
 -(m)[4]  * (det4 = gcgDETERMINANT3x3_((m)[1], (m)[2], (m)[3], (m)[9], (m)[10], (m)[11], (m)[13], (m)[14], (m)[15])) \
 +(m)[8]  * (det8 = gcgDETERMINANT3x3_((m)[1], (m)[2], (m)[3], (m)[5], (m)[6],  (m)[7],  (m)[13], (m)[14], (m)[15])) \
 -(m)[12] * (det12= gcgDETERMINANT3x3_((m)[1], (m)[2], (m)[3], (m)[5], (m)[6],  (m)[7],  (m)[9],  (m)[10], (m)[11])) )

#define gcgDETERMINANTMATRIX3(m) (gcgDETERMINANT3x3_((m)[0], (m)[1], (m)[2], (m)[3], (m)[4], (m)[5], (m)[6], (m)[7], (m)[8]))
#define gcgDETERMINANTMATRIX4(m) \
( (m)[0]  * gcgDETERMINANT3x3_((m)[5], (m)[6], (m)[7], (m)[9], (m)[10], (m)[11], (m)[13], (m)[14], (m)[15]) \
 -(m)[4]  * gcgDETERMINANT3x3_((m)[1], (m)[2], (m)[3], (m)[9], (m)[10], (m)[11], (m)[13], (m)[14], (m)[15]) \
 +(m)[8]  * gcgDETERMINANT3x3_((m)[1], (m)[2], (m)[3], (m)[5], (m)[6],  (m)[7],  (m)[13], (m)[14], (m)[15]) \
 -(m)[12] * gcgDETERMINANT3x3_((m)[1], (m)[2], (m)[3], (m)[5], (m)[6],  (m)[7],  (m)[9],  (m)[10], (m)[11]) )


// Commands ///////////////////////////////////////////////
#define gcgZEROVECTOR2(v) { (v)[0] = (v)[1] = 0.0; }
#define gcgZEROVECTOR3(v) { (v)[0] = (v)[1] = (v)[2] = 0.0; }
#define gcgZEROVECTOR4(v) { (v)[0] = (v)[1] = (v)[2] = (v)[3] = 0.0; }

#define gcgSETVECTOR2(dest, a, b) {(dest)[0] = a; (dest)[1] = b; }
#define gcgSETVECTOR3(dest, a, b, c) {(dest)[0] = a; (dest)[1] = b; (dest)[2] = c; }
#define gcgSETVECTOR4(dest, a, b, c, d) {(dest)[0] = a; (dest)[1] = b; (dest)[2] = c; (dest)[3]=d; }

static inline void gcgCOPYVECTOR2(VECTOR2 dest, VECTOR2 src) {dest[0] = src[0]; dest[1] = src[1]; }
static inline void gcgCOPYVECTOR2(VECTOR2d dest, VECTOR2d src) {dest[0] = src[0]; dest[1] = src[1]; }
static inline void gcgCOPYVECTOR2(VECTOR2 dest, VECTOR2d src) {dest[0] = (float) src[0]; dest[1] = (float) src[1]; }
static inline void gcgCOPYVECTOR2(VECTOR2d dest, VECTOR2 src) {dest[0] = src[0]; dest[1] = src[1]; }
static inline void gcgCOPYVECTOR3(VECTOR3 dest, VECTOR3 src) { dest[0] = src[0]; dest[1] = src[1]; dest[2] = src[2]; }
static inline void gcgCOPYVECTOR3(VECTOR3 dest, VECTOR3d src) { dest[0] = (float) src[0]; dest[1] = (float) src[1]; dest[2] = (float) src[2]; }
static inline void gcgCOPYVECTOR3(VECTOR3d dest, VECTOR3 src) { dest[0] = src[0]; dest[1] = src[1]; dest[2] = src[2]; }
static inline void gcgCOPYVECTOR3(VECTOR3d dest, VECTOR3d src) { dest[0] = src[0]; dest[1] = src[1]; dest[2] = src[2]; }
static inline void gcgCOPYVECTOR4(VECTOR4 dest, VECTOR4 src) { dest[0] = src[0]; dest[1] = src[1]; dest[2] = src[2]; dest[3] = src[3];}
static inline void gcgCOPYVECTOR4(VECTOR4 dest, VECTOR4d src) { dest[0] = (float) src[0]; dest[1] = (float) src[1]; dest[2] = (float) src[2]; dest[3] = (float) src[3];}
static inline void gcgCOPYVECTOR4(VECTOR4d dest, VECTOR4 src) { dest[0] = src[0]; dest[1] = src[1]; dest[2] = src[2]; dest[3] = src[3];}
static inline void gcgCOPYVECTOR4(VECTOR4d dest, VECTOR4d src) { dest[0] = src[0]; dest[1] = src[1]; dest[2] = src[2]; dest[3] = src[3];}

static inline void gcgSCALEVECTOR2(VECTOR2 dest, VECTOR2 v, float factor) { dest[0]=factor*v[0]; (dest)[1]=factor*v[1]; }
static inline void gcgSCALEVECTOR2(VECTOR2d dest, VECTOR2d v, double factor) { dest[0]=factor*v[0]; (dest)[1]=factor*v[1]; }
static inline void gcgSCALEVECTOR3(VECTOR3 dest, VECTOR3 v, float factor) {dest[0]=factor*v[0]; dest[1]=factor*v[1]; dest[2]=factor*v[2];}
static inline void gcgSCALEVECTOR3(VECTOR3d dest, VECTOR3d v, double factor) {dest[0]=factor*v[0]; dest[1]=factor*v[1]; dest[2]=factor*v[2];}
static inline void gcgSCALEVECTOR4(VECTOR4 dest, VECTOR4 v, float factor) {dest[0]=factor*v[0]; dest[1]=factor*v[1]; dest[2]=factor*v[2]; dest[3]=factor*v[3];}
static inline void gcgSCALEVECTOR4(VECTOR4d dest, VECTOR4d v, double factor) {dest[0]=factor*v[0]; dest[1]=factor*v[1]; dest[2]=factor*v[2]; dest[3]=factor*v[3];}

#define gcgADDVECTOR2(dest,v1,v2) {(dest)[0]=(v1)[0]+(v2)[0]; (dest)[1]=(v1)[1]+(v2)[1]; }
#define gcgADDVECTOR3(dest,v1,v2) {(dest)[0]=(v1)[0]+(v2)[0]; (dest)[1]=(v1)[1]+(v2)[1]; (dest)[2]=(v1)[2]+(v2)[2];}
#define gcgADDVECTOR4(dest,v1,v2) {(dest)[0]=(v1)[0]+(v2)[0]; (dest)[1]=(v1)[1]+(v2)[1]; (dest)[2]=(v1)[2]+(v2)[2]; (dest)[3]=(v1)[3]+(v2)[3];}

#define gcgSUBVECTOR2(dest, v1, v2) {(dest)[0]=(v1)[0]-(v2)[0]; (dest)[1]=(v1)[1]-(v2)[1]; }
#define gcgSUBVECTOR3(dest, v1, v2) {(dest)[0]=(v1)[0]-(v2)[0]; (dest)[1]=(v1)[1]-(v2)[1]; (dest)[2]=(v1)[2]-(v2)[2];}
#define gcgSUBVECTOR4(dest, v1, v2) {(dest)[0]=(v1)[0]-(v2)[0]; (dest)[1]=(v1)[1]-(v2)[1]; (dest)[2]=(v1)[2]-(v2)[2]; (dest)[3]=(v1)[3]-(v2)[3];}

// Normalize VECTOR2
static inline void gcgNORMALIZEVECTOR2(VECTOR2 dest, VECTOR2 v) {
  register float factor = (float) (1.0 / gcgLENGTHVECTOR2((v)));
  gcgSCALEVECTOR2((dest), (v), factor);
}

static inline void gcgNORMALIZEVECTOR2(VECTOR2d dest, VECTOR2d v) {
  register double factor = 1.0 / gcgLENGTHVECTOR2((v));
  gcgSCALEVECTOR2((dest), (v), factor);
}

// Normalize VECTOR3
static inline void gcgNORMALIZEVECTOR3(VECTOR3 dest, VECTOR3 v) {
  register float factor = (float) (1.0 / gcgLENGTHVECTOR3((v)));
  gcgSCALEVECTOR3((dest), (v), factor);
}

static inline void gcgNORMALIZEVECTOR3(VECTOR3d dest, VECTOR3d v) {
  register double factor = 1.0 / gcgLENGTHVECTOR3((v));
  gcgSCALEVECTOR3((dest), (v), factor);
}

// Normalize VECTOR4
static inline void gcgNORMALIZEVECTOR4(VECTOR4 dest, VECTOR4 v) {
  register float factor = (float) (1.0 / gcgLENGTHVECTOR4((v)));
  gcgSCALEVECTOR4((dest), (v), factor);
}

static inline void gcgNORMALIZEVECTOR4(VECTOR4d dest, VECTOR4d v) {
  register double factor = 1.0 / gcgLENGTHVECTOR4((v));
  gcgSCALEVECTOR4((dest), (v), factor);
}


#define gcgCROSSVECTOR3(dest, v1, v2) {                          \
              (dest)[0] = (v1)[1] * (v2)[2] - (v1)[2] * (v2)[1]; \
              (dest)[1] = (v1)[2] * (v2)[0] - (v1)[0] * (v2)[2]; \
              (dest)[2] = (v1)[0] * (v2)[1] - (v1)[1] * (v2)[0];}

// Creates a tensor from two vectors: T = v1 * transpose(v2)
#define gcgVECTOR3TOTENSOR(t, v1, v2) { \
  (t)[0] = (v1)[0]*(v2)[0]; (t)[1] = (v1)[0]*(v2)[1]; (t)[2] = (v1)[0]*(v2)[2]; \
  (t)[3] = (v1)[1]*(v2)[0]; (t)[4] = (v1)[1]*(v2)[1]; (t)[5] = (v1)[1]*(v2)[2]; \
  (t)[6] = (v1)[2]*(v2)[0]; (t)[7] = (v1)[2]*(v2)[1]; (t)[8] = (v1)[2]*(v2)[2]; }


// MATRIX COMMANDS MACROS

#define gcgSETMATRIX2(dest, a, b, c, d) {   \
  (dest)[ 0] = a; (dest)[ 1] = b; (dest)[ 2] = c; (dest)[ 3] = d; }

#define gcgSETMATRIX3(dest, a, b, c, d, e, f, g, h, i) {   \
  (dest)[ 0] = a; (dest)[ 1] = b; (dest)[ 2] = c; \
  (dest)[ 3] = d; (dest)[ 4] = e; (dest)[ 5] = f; \
  (dest)[ 6] = g; (dest)[ 7] = h; (dest)[ 8] = i; }

#define gcgSETMATRIX4(dest, a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p) {   \
  (dest)[ 0] = a; (dest)[ 1] = b; (dest)[ 2] = c; (dest)[ 3] = d; \
  (dest)[ 4] = e; (dest)[ 5] = f; (dest)[ 6] = g; (dest)[ 7] = h; \
  (dest)[ 8] = i; (dest)[ 9] = j; (dest)[10] = k; (dest)[11] = l; \
  (dest)[12] = m; (dest)[13] = n; (dest)[14] = o; (dest)[15] = p; }

static inline void gcgCOPYMATRIX2(MATRIX2 dest, MATRIX2 m) {
  dest[ 0] = m[ 0]; dest[ 1] = m[ 1]; dest[ 2] = m[ 2]; dest[ 3] = m[ 3];
}

static inline void gcgCOPYMATRIX2(MATRIX2d dest, MATRIX2d m) {
  dest[ 0] = m[ 0]; dest[ 1] = m[ 1]; dest[ 2] = m[ 2]; dest[ 3] = m[ 3];
}

static inline void gcgCOPYMATRIX2(MATRIX2d dest, MATRIX2 m) {
  dest[ 0] = m[ 0]; dest[ 1] = m[ 1]; dest[ 2] = m[ 2]; dest[ 3] = m[ 3];
}

static inline void gcgCOPYMATRIX2(MATRIX2 dest, MATRIX2d m) {
  dest[ 0] = (float) m[ 0]; dest[ 1] = (float) m[ 1]; dest[ 2] = (float) m[ 2]; dest[ 3] = (float) m[ 3];
}

static inline void gcgCOPYMATRIX3(MATRIX3 dest, MATRIX3 m) {
  dest[ 0] = m[ 0]; dest[ 1] = m[ 1]; dest[ 2] = m[ 2];
  dest[ 3] = m[ 3]; dest[ 4] = m[ 4]; dest[ 5] = m[ 5];
  dest[ 6] = m[ 6]; dest[ 7] = m[ 7]; dest[ 8] = m[ 8];
}

static inline void gcgCOPYMATRIX3(MATRIX3d dest, MATRIX3d m) {
  dest[ 0] = m[ 0]; dest[ 1] = m[ 1]; dest[ 2] = m[ 2];
  dest[ 3] = m[ 3]; dest[ 4] = m[ 4]; dest[ 5] = m[ 5];
  dest[ 6] = m[ 6]; dest[ 7] = m[ 7]; dest[ 8] = m[ 8];
}

static inline void gcgCOPYMATRIX3(MATRIX3d dest, MATRIX3 m) {
  dest[ 0] = m[ 0]; dest[ 1] = m[ 1]; dest[ 2] = m[ 2];
  dest[ 3] = m[ 3]; dest[ 4] = m[ 4]; dest[ 5] = m[ 5];
  dest[ 6] = m[ 6]; dest[ 7] = m[ 7]; dest[ 8] = m[ 8];
}

static inline void gcgCOPYMATRIX3(MATRIX3 dest, MATRIX3d m) {
  dest[ 0] = (float) m[ 0]; dest[ 1] = (float) m[ 1]; dest[ 2] = (float) m[ 2];
  dest[ 3] = (float) m[ 3]; dest[ 4] = (float) m[ 4]; dest[ 5] = (float) m[ 5];
  dest[ 6] = (float) m[ 6]; dest[ 7] = (float) m[ 7]; dest[ 8] = (float) m[ 8];
}


static inline void gcgCOPYMATRIX4(MATRIX4 dest, MATRIX4 m) {
  dest[ 0] = m[ 0]; dest[ 1] = m[ 1]; dest[ 2] = m[ 2]; dest[ 3] = m[ 3];
  dest[ 4] = m[ 4]; dest[ 5] = m[ 5]; dest[ 6] = m[ 6]; dest[ 7] = m[ 7];
  dest[ 8] = m[ 8]; dest[ 9] = m[ 9]; dest[10] = m[10]; dest[11] = m[11];
  dest[12] = m[12]; dest[13] = m[13]; dest[14] = m[14]; dest[15] = m[15];
}

static inline void gcgCOPYMATRIX4(MATRIX4d dest, MATRIX4 m) {
  dest[ 0] = m[ 0]; dest[ 1] = m[ 1]; dest[ 2] = m[ 2]; dest[ 3] = m[ 3];
  dest[ 4] = m[ 4]; dest[ 5] = m[ 5]; dest[ 6] = m[ 6]; dest[ 7] = m[ 7];
  dest[ 8] = m[ 8]; dest[ 9] = m[ 9]; dest[10] = m[10]; dest[11] = m[11];
  dest[12] = m[12]; dest[13] = m[13]; dest[14] = m[14]; dest[15] = m[15];
}

static inline void gcgCOPYMATRIX4(MATRIX4d dest, MATRIX4d m) {
  dest[ 0] = m[ 0]; dest[ 1] = m[ 1]; dest[ 2] = m[ 2]; dest[ 3] = m[ 3];
  dest[ 4] = m[ 4]; dest[ 5] = m[ 5]; dest[ 6] = m[ 6]; dest[ 7] = m[ 7];
  dest[ 8] = m[ 8]; dest[ 9] = m[ 9]; dest[10] = m[10]; dest[11] = m[11];
  dest[12] = m[12]; dest[13] = m[13]; dest[14] = m[14]; dest[15] = m[15];
}

static inline void gcgCOPYMATRIX4(MATRIX4 dest, MATRIX4d m) {
  dest[ 0] = (float) m[ 0]; dest[ 1] = (float) m[ 1]; dest[ 2] = (float) m[ 2]; dest[ 3] = (float) m[ 3];
  dest[ 4] = (float) m[ 4]; dest[ 5] = (float) m[ 5]; dest[ 6] = (float) m[ 6]; dest[ 7] = (float) m[ 7];
  dest[ 8] = (float) m[ 8]; dest[ 9] = (float) m[ 9]; dest[10] = (float) m[10]; dest[11] = (float) m[11];
  dest[12] = (float) m[12]; dest[13] = (float) m[13]; dest[14] = (float) m[14]; dest[15] = (float) m[15];
}

#define gcgIDENTITYMATRIX2(matriz) {  \
        (matriz)[0]  = (matriz)[3] = 1.f; \
        (matriz)[1]  = (matriz)[2] = 0.0f; }

#define gcgIDENTITYMATRIX3(matriz) {  \
        (matriz)[0]  = (matriz)[4]  = (matriz)[8]  = 1.f; \
        (matriz)[1]  = (matriz)[2]  = (matriz)[3]  = (matriz)[5]  = (matriz)[6]  = (matriz)[7]  = 0.0f; }

#define gcgIDENTITYMATRIX4(matriz) {  \
        (matriz)[0]  = (matriz)[5]  = (matriz)[10] = (matriz)[15] = 1.f; \
        (matriz)[1]  = (matriz)[2]  = (matriz)[3]  = (matriz)[4]  = (matriz)[6]  = (matriz)[7]  = \
        (matriz)[8]  = (matriz)[9]  = (matriz)[11] = (matriz)[12] = (matriz)[13] = (matriz)[14] = 0.f; }

//Transpose matrix 3x3
#define gcgTRANSPOSEMATRIX3(tm, m) {   \
  (tm)[ 0] = (m)[0];   (tm)[1] = (m)[3];   (tm)[2] = (m)[6]; \
  (tm)[ 3] = (m)[1];   (tm)[4] = (m)[4];   (tm)[5] = (m)[7]; \
  (tm)[ 6] = (m)[2];   (tm)[7] = (m)[5];   (tm)[8] = (m)[8]; }

//Transpose matrix 4x4
#define gcgTRANSPOSEMATRIX4(tm, m) {   \
  (tm)[ 0] = (m)[ 0];   (tm)[ 1] = (m)[ 4];   (tm)[ 2] = (m)[ 8];   (tm)[ 3] = (m)[12]; \
  (tm)[ 4] = (m)[ 1];   (tm)[ 5] = (m)[ 5];   (tm)[ 6] = (m)[ 9];   (tm)[ 7] = (m)[13]; \
  (tm)[ 8] = (m)[ 2];   (tm)[ 9] = (m)[ 6];   (tm)[10] = (m)[10];   (tm)[11] = (m)[14]; \
  (tm)[12] = (m)[ 3];   (tm)[13] = (m)[ 7];   (tm)[14] = (m)[11];   (tm)[15] = (m)[15]; }

// Scale matrix 2x2
#define gcgSCALEMATRIX2(dest, m, s) \
  {(dest)[0] = (m)[0] * (s); (dest)[1] = (m)[1] * (s); (dest)[2] = (m)[2] * (s); (dest)[3] = (m)[3] * (s); }

// Scale matrix 3x3
#define gcgSCALEMATRIX3(dest, m, s) {  \
  (dest)[0] = (m)[0] * (s); (dest)[1] = (m)[1] * (s); (dest)[2] = (m)[2] * (s); \
  (dest)[3] = (m)[3] * (s); (dest)[4] = (m)[4] * (s); (dest)[5] = (m)[5] * (s); \
  (dest)[6] = (m)[6] * (s); (dest)[7] = (m)[7] * (s); (dest)[8] = (m)[8] * (s); }

// Scale matrix 4x4
#define gcgSCALEMATRIX4(dest, m, s) {   \
  (dest)[ 0] = (m)[ 0] * (s); (dest)[ 1] = (m)[ 1] * (s); (dest)[ 2] = (m)[ 2] * (s); (dest)[ 3] = (m)[ 3] * (s); \
  (dest)[ 4] = (m)[ 4] * (s); (dest)[ 5] = (m)[ 5] * (s); (dest)[ 6] = (m)[ 6] * (s); (dest)[ 7] = (m)[ 7] * (s); \
  (dest)[ 8] = (m)[ 8] * (s); (dest)[ 9] = (m)[ 9] * (s); (dest)[10] = (m)[10] * (s); (dest)[11] = (m)[11] * (s); \
  (dest)[12] = (m)[12] * (s); (dest)[13] = (m)[13] * (s); (dest)[14] = (m)[14] * (s); (dest)[15] = (m)[15] * (s); }

// Additions and subtractions
#define gcgADDMATRIX2(dest, m1, m2) {   \
  (dest)[ 0] = (m1)[ 0] + (m2)[ 0]; (dest)[ 1] = (m1)[ 1] + (m2)[ 1]; \
  (dest)[ 2] = (m1)[ 2] + (m2)[ 2]; (dest)[ 3] = (m1)[ 3] + (m2)[ 3]; }

#define gcgADDMATRIX3(dest, m1, m2) {   \
  (dest)[ 0] = (m1)[ 0] + (m2)[ 0]; (dest)[ 1] = (m1)[ 1] + (m2)[ 1]; (dest)[ 2] = (m1)[ 2] + (m2)[ 2]; \
  (dest)[ 3] = (m1)[ 3] + (m2)[ 3]; (dest)[ 4] = (m1)[ 4] + (m2)[ 4]; (dest)[ 5] = (m1)[ 5] + (m2)[ 5]; \
  (dest)[ 6] = (m1)[ 6] + (m2)[ 6]; (dest)[ 7] = (m1)[ 7] + (m2)[ 7]; (dest)[ 8] = (m1)[ 8] + (m2)[ 8]; }

#define gcgADDMATRIX4(dest, m1, m2) {   \
  (dest)[ 0] = (m1)[ 0] + (m2)[ 0]; (dest)[ 1] = (m1)[ 1] + (m2)[ 1]; (dest)[ 2] = (m1)[ 2] + (m2)[ 2]; (dest)[ 3] = (m1)[ 3] + (m2)[ 3]; \
  (dest)[ 4] = (m1)[ 4] + (m2)[ 4]; (dest)[ 5] = (m1)[ 5] + (m2)[ 5]; (dest)[ 6] = (m1)[ 6] + (m2)[ 6]; (dest)[ 7] = (m1)[ 7] + (m2)[ 7]; \
  (dest)[ 8] = (m1)[ 8] + (m2)[ 8]; (dest)[ 9] = (m1)[ 9] + (m2)[ 9]; (dest)[10] = (m1)[10] + (m2)[10]; (dest)[11] = (m1)[11] + (m2)[11]; \
  (dest)[12] = (m1)[12] + (m2)[12]; (dest)[13] = (m1)[13] + (m2)[13]; (dest)[14] = (m1)[14] + (m2)[14]; (dest)[15] = (m1)[15] + (m2)[15]; }

#define gcgSUBMATRIX2(dest, m1, m2) {   \
  (dest)[ 0] = (m1)[ 0] - (m2)[ 0]; (dest)[ 1] = (m1)[ 1] - (m2)[ 1]; \
  (dest)[ 2] = (m1)[ 2] - (m2)[ 2]; (dest)[ 3] = (m1)[ 3] - (m2)[ 3]; }

#define gcgSUBMATRIX3(dest, m1, m2) {   \
  (dest)[ 0] = (m1)[ 0] - (m2)[ 0]; (dest)[ 1] = (m1)[ 1] - (m2)[ 1]; (dest)[ 2] = (m1)[ 2] - (m2)[ 2]; \
  (dest)[ 3] = (m1)[ 3] - (m2)[ 3]; (dest)[ 4] = (m1)[ 4] - (m2)[ 4]; (dest)[ 5] = (m1)[ 5] - (m2)[ 5]; \
  (dest)[ 6] = (m1)[ 6] - (m2)[ 6]; (dest)[ 7] = (m1)[ 7] - (m2)[ 7]; (dest)[ 8] = (m1)[ 8] - (m2)[ 8]; }

#define gcgSUBMATRIX4(dest, m1, m2) {   \
  (dest)[ 0] = (m1)[ 0] - (m2)[ 0]; (dest)[ 1] = (m1)[ 1] - (m2)[ 1]; (dest)[ 2] = (m1)[ 2] - (m2)[ 2]; (dest)[ 3] = (m1)[ 3] - (m2)[ 3]; \
  (dest)[ 4] = (m1)[ 4] - (m2)[ 4]; (dest)[ 5] = (m1)[ 5] - (m2)[ 5]; (dest)[ 6] = (m1)[ 6] - (m2)[ 6]; (dest)[ 7] = (m1)[ 7] - (m2)[ 7]; \
  (dest)[ 8] = (m1)[ 8] - (m2)[ 8]; (dest)[ 9] = (m1)[ 9] - (m2)[ 9]; (dest)[10] = (m1)[10] - (m2)[10]; (dest)[11] = (m1)[11] - (m2)[11]; \
  (dest)[12] = (m1)[12] - (m2)[12]; (dest)[13] = (m1)[13] - (m2)[13]; (dest)[14] = (m1)[14] - (m2)[14]; (dest)[15] = (m1)[15] - (m2)[15]; }

//Multiply matrix 3x3
#define gcgMULTMATRIX3(result, m1, m2) {   \
        (result)[ 0] = (m1)[ 0] * (m2)[ 0] + (m1)[ 1] * (m2)[ 3] + (m1)[ 2] * (m2)[ 6]; \
        (result)[ 1] = (m1)[ 0] * (m2)[ 1] + (m1)[ 1] * (m2)[ 4] + (m1)[ 2] * (m2)[ 7]; \
        (result)[ 2] = (m1)[ 0] * (m2)[ 2] + (m1)[ 1] * (m2)[ 5] + (m1)[ 2] * (m2)[ 8]; \
        (result)[ 3] = (m1)[ 3] * (m2)[ 0] + (m1)[ 4] * (m2)[ 3] + (m1)[ 5] * (m2)[ 6]; \
        (result)[ 4] = (m1)[ 3] * (m2)[ 1] + (m1)[ 4] * (m2)[ 4] + (m1)[ 5] * (m2)[ 7]; \
        (result)[ 5] = (m1)[ 3] * (m2)[ 2] + (m1)[ 4] * (m2)[ 5] + (m1)[ 5] * (m2)[ 8]; \
        (result)[ 6] = (m1)[ 6] * (m2)[ 0] + (m1)[ 7] * (m2)[ 3] + (m1)[ 8] * (m2)[ 6]; \
        (result)[ 7] = (m1)[ 6] * (m2)[ 1] + (m1)[ 7] * (m2)[ 4] + (m1)[ 8] * (m2)[ 7]; \
        (result)[ 8] = (m1)[ 6] * (m2)[ 2] + (m1)[ 7] * (m2)[ 5] + (m1)[ 8] * (m2)[ 8]; }

// Multiplica as matrizes 4x4 m1 e m2 sem utilizar laco
#define gcgMULTMATRIX4(result, m1, m2) {   \
        (result)[ 0] = (m1)[ 0] * (m2)[ 0] + (m1)[ 1] * (m2)[ 4] + (m1)[ 2] * (m2)[ 8] + (m1)[ 3] * (m2)[12]; \
        (result)[ 1] = (m1)[ 0] * (m2)[ 1] + (m1)[ 1] * (m2)[ 5] + (m1)[ 2] * (m2)[ 9] + (m1)[ 3] * (m2)[13]; \
        (result)[ 2] = (m1)[ 0] * (m2)[ 2] + (m1)[ 1] * (m2)[ 6] + (m1)[ 2] * (m2)[10] + (m1)[ 3] * (m2)[14]; \
        (result)[ 3] = (m1)[ 0] * (m2)[ 3] + (m1)[ 1] * (m2)[ 7] + (m1)[ 2] * (m2)[11] + (m1)[ 3] * (m2)[15]; \
        (result)[ 4] = (m1)[ 4] * (m2)[ 0] + (m1)[ 5] * (m2)[ 4] + (m1)[ 6] * (m2)[ 8] + (m1)[ 7] * (m2)[12]; \
        (result)[ 5] = (m1)[ 4] * (m2)[ 1] + (m1)[ 5] * (m2)[ 5] + (m1)[ 6] * (m2)[ 9] + (m1)[ 7] * (m2)[13]; \
        (result)[ 6] = (m1)[ 4] * (m2)[ 2] + (m1)[ 5] * (m2)[ 6] + (m1)[ 6] * (m2)[10] + (m1)[ 7] * (m2)[14]; \
        (result)[ 7] = (m1)[ 4] * (m2)[ 3] + (m1)[ 5] * (m2)[ 7] + (m1)[ 6] * (m2)[11] + (m1)[ 7] * (m2)[15]; \
        (result)[ 8] = (m1)[ 8] * (m2)[ 0] + (m1)[ 9] * (m2)[ 4] + (m1)[10] * (m2)[ 8] + (m1)[11] * (m2)[12]; \
        (result)[ 9] = (m1)[ 8] * (m2)[ 1] + (m1)[ 9] * (m2)[ 5] + (m1)[10] * (m2)[ 9] + (m1)[11] * (m2)[13]; \
        (result)[10] = (m1)[ 8] * (m2)[ 2] + (m1)[ 9] * (m2)[ 6] + (m1)[10] * (m2)[10] + (m1)[11] * (m2)[14]; \
        (result)[11] = (m1)[ 8] * (m2)[ 3] + (m1)[ 9] * (m2)[ 7] + (m1)[10] * (m2)[11] + (m1)[11] * (m2)[15]; \
        (result)[12] = (m1)[12] * (m2)[ 0] + (m1)[13] * (m2)[ 4] + (m1)[14] * (m2)[ 8] + (m1)[15] * (m2)[12]; \
        (result)[13] = (m1)[12] * (m2)[ 1] + (m1)[13] * (m2)[ 5] + (m1)[14] * (m2)[ 9] + (m1)[15] * (m2)[13]; \
        (result)[14] = (m1)[12] * (m2)[ 2] + (m1)[13] * (m2)[ 6] + (m1)[14] * (m2)[10] + (m1)[15] * (m2)[14]; \
        (result)[15] = (m1)[12] * (m2)[ 3] + (m1)[13] * (m2)[ 7] + (m1)[14] * (m2)[11] + (m1)[15] * (m2)[15]; }

// Multiply the 3x4 affine matrices m1 e m2.
#define gcgMULTMATRIX34(result, m1, m2) {   \
        (result)[0]  = (m1)[0] * (m2)[0] + (m1)[1] * (m2)[4] + (m1)[2] * (m2)[8];            \
        (result)[1]  = (m1)[0] * (m2)[1] + (m1)[1] * (m2)[5] + (m1)[2] * (m2)[9];            \
        (result)[2]  = (m1)[0] * (m2)[2] + (m1)[1] * (m2)[6] + (m1)[2] * (m2)[10];           \
        (result)[3]  = (m1)[0] * (m2)[3] + (m1)[1] * (m2)[7] + (m1)[2] * (m2)[11] + (m1)[3]; \
        (result)[4]  = (m1)[4] * (m2)[0] + (m1)[5] * (m2)[4] + (m1)[6] * (m2)[8];            \
        (result)[5]  = (m1)[4] * (m2)[1] + (m1)[5] * (m2)[5] + (m1)[6] * (m2)[9];            \
        (result)[6]  = (m1)[4] * (m2)[2] + (m1)[5] * (m2)[6] + (m1)[6] * (m2)[10];           \
        (result)[7]  = (m1)[4] * (m2)[3] + (m1)[5] * (m2)[7] + (m1)[6] * (m2)[11] + (m1)[7]; \
        (result)[8]  = (m1)[8] * (m2)[0] + (m1)[9] * (m2)[4] + (m1)[10] * (m2)[8];           \
        (result)[9]  = (m1)[8] * (m2)[1] + (m1)[9] * (m2)[5] + (m1)[10] * (m2)[9];           \
        (result)[10] = (m1)[8] * (m2)[2] + (m1)[9] * (m2)[6] + (m1)[10] * (m2)[10];          \
        (result)[11] = (m1)[8] * (m2)[3] + (m1)[9] * (m2)[7] + (m1)[10] * (m2)[11] + (m1)[11]; }


//Multiply a 3x1 vector by a 1x3 vector resulting in a 3x3 matrix
#define gcgMATRIX3VECTOR3VECTOR3(result, v1, v2) {   \
  (result)[0] = (v1)[0] * (v2)[0]; (result)[1] = (v1)[0] *  (v2)[1]; (result)[2] =  (v1)[0] *  (v2)[2]; \
  (result)[3] = (v1)[1] * (v2)[0]; (result)[4] = (v1)[1] *  (v2)[1]; (result)[5] =  (v1)[1] *  (v2)[2]; \
  (result)[6] = (v1)[2] * (v2)[0]; (result)[7] = (v1)[2] * (v2)[1]; (result)[8] =  (v1)[2] *  (v2)[2]; }


#define gcgAPPLYMATRIX3VECTOR3(r, t, vetor) { \
        (r)[0] =  (vetor)[0] * (t)[0]  +  (vetor)[1] * (t)[1]  +  (vetor)[2] * (t)[2]; \
        (r)[1] =  (vetor)[0] * (t)[3]  +  (vetor)[1] * (t)[4]  +  (vetor)[2] * (t)[5]; \
        (r)[2] =  (vetor)[0] * (t)[6]  +  (vetor)[1] * (t)[7]  +  (vetor)[2] * (t)[8]; }


#define gcgAPPLYMATRIX4VECTOR4(r, t, vetor) { \
        (r)[0] =  (vetor)[0] * (t)[0]  +  (vetor)[1] * (t)[1]  +  (vetor)[2] * (t)[2]  +  (vetor)[3] * (t)[3];  \
        (r)[1] =  (vetor)[0] * (t)[4]  +  (vetor)[1] * (t)[5]  +  (vetor)[2] * (t)[6]  +  (vetor)[3] * (t)[7];  \
        (r)[2] =  (vetor)[0] * (t)[8]  +  (vetor)[1] * (t)[9]  +  (vetor)[2] * (t)[10] +  (vetor)[3] * (t)[11]; \
        (r)[3] =  (vetor)[0] * (t)[12] +  (vetor)[1] * (t)[13] +  (vetor)[2] * (t)[14] +  (vetor)[3] * (t)[15]; }

#define gcgAPPLYMATRIX34VECTOR3(r, t, vetor) { \
        (r)[0] =  (vetor)[0] * (t)[0]  +  (vetor)[1] * (t)[1]  +  (vetor)[2] * (t)[2]  + (t)[3];  \
        (r)[1] =  (vetor)[0] * (t)[4]  +  (vetor)[1] * (t)[5]  +  (vetor)[2] * (t)[6]  + (t)[7];  \
        (r)[2] =  (vetor)[0] * (t)[8]  +  (vetor)[1] * (t)[9]  +  (vetor)[2] * (t)[10] + (t)[11]; }

#define gcgLINEARMATRIX4VECTOR3(r, t, vetor) { \
        (r)[0] =  (vetor)[0] * (t)[0]  +  (vetor)[1] * (t)[1]  +  (vetor)[2] * (t)[2];  \
        (r)[1] =  (vetor)[0] * (t)[4]  +  (vetor)[1] * (t)[5]  +  (vetor)[2] * (t)[6];  \
        (r)[2] =  (vetor)[0] * (t)[8]  +  (vetor)[1] * (t)[9]  +  (vetor)[2] * (t)[10]; }

// QUATERNION ALGEBRA ///////////////////////////////////////////////

// Expressions //////////////////////////////////////////////////////
// Commands /////////////////////////////////////////////////////////

// Makes an unit quaternion
#define gcgNORMALIZEQUATERNION(qf, q) { \
    register double mag = 1.0/(q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3]); \
    gcgSCALEVECTOR3(qf, q, mag); }

// Quaternion conjugate
#define gcgCONJUGATEQUATERNION(qi, q) {   \
  (qi)[0] = -(q)[0];                      \
  (qi)[1] = -(q)[1];                      \
  (qi)[2] = -(q)[2];                      \
  (qi)[3] =  (q)[3]; }


// Multiply two quaternions (Graham's product)
static inline void gcgMULTQUATERNION(VECTOR4 dest, VECTOR4 q2, VECTOR4 q1) {
  VECTOR4 t1;
  t1[3] = (q1)[3]*(q2)[3] - (q1)[0]*(q2)[0] - (q1)[1]*(q2)[1] - (q1)[2]*(q2)[2];
        t1[0] = (q1)[3]*(q2)[0] + (q1)[0]*(q2)[3] + (q1)[1]*(q2)[2] - (q1)[2]*(q2)[1];
        t1[1] = (q1)[3]*(q2)[1] + (q1)[1]*(q2)[3] + (q1)[2]*(q2)[0] - (q1)[0]*(q2)[2];
        t1[2] = (q1)[3]*(q2)[2] + (q1)[2]*(q2)[3] + (q1)[0]*(q2)[1] - (q1)[1]*(q2)[0];
  gcgCOPYVECTOR4(dest, t1);
}

static inline void gcgMULTQUATERNION(VECTOR4d dest, VECTOR4d q2, VECTOR4d q1) {
  VECTOR4d t1;
  t1[3] = (q1)[3]*(q2)[3] - (q1)[0]*(q2)[0] - (q1)[1]*(q2)[1] - (q1)[2]*(q2)[2];
        t1[0] = (q1)[3]*(q2)[0] + (q1)[0]*(q2)[3] + (q1)[1]*(q2)[2] - (q1)[2]*(q2)[1];
        t1[1] = (q1)[3]*(q2)[1] + (q1)[1]*(q2)[3] + (q1)[2]*(q2)[0] - (q1)[0]*(q2)[2];
        t1[2] = (q1)[3]*(q2)[2] + (q1)[2]*(q2)[3] + (q1)[0]*(q2)[1] - (q1)[1]*(q2)[0];
  gcgCOPYVECTOR4(dest, t1);
}

// Build a 4x4 rotation matrix from a quaternion
static inline void gcgQUATERNIONTOMATRIX4(MATRIX4 m, VECTOR4 q) {
    register double xx = (q)[0] * (q)[0], yy = (q)[1] * (q)[1], zz = (q)[2] * (q)[2];
    register double xy = (q)[0] * (q)[1], xz = (q)[0] * (q)[2], yz = (q)[1] * (q)[2];
    register double xw = (q)[0] * (q)[3], yw = (q)[1] * (q)[3], zw = (q)[2] * (q)[3];
    m[0] =  (float) (1.0 - (yy + yy + zz + zz));
    m[1] =  (float) (xy + xy - zw - zw);
    m[2] =  (float) (xz + xz + yw + yw);
    m[3] = 0.0f;
    m[4] =  (float) (xy + xy + zw + zw);
    m[5] =  (float) (1.0 - (zz + zz + xx + xx));
    m[6] =  (float) (yz + yz - xw - xw);
    m[7] = 0.0f;
    m[8] =  (float) (xz + xz - yw - yw);
    m[9] =  (float) (yz + yz + xw + xw);
    m[10] =  (float) (1.0 - (yy + yy + xx + xx));
    m[11] =  m[12] = m[13] = m[14] = 0.0f;
    m[15] = 1.0f;
}

// Build a 4x4 rotation matrix from a quaternion
static inline void gcgQUATERNIONTOMATRIX4(MATRIX4d m, VECTOR4d q) {
    register double xx = (q)[0] * (q)[0], yy = (q)[1] * (q)[1], zz = (q)[2] * (q)[2];
    register double xy = (q)[0] * (q)[1], xz = (q)[0] * (q)[2], yz = (q)[1] * (q)[2];
    register double xw = (q)[0] * (q)[3], yw = (q)[1] * (q)[3], zw = (q)[2] * (q)[3];
    m[0] =  (1.0 - (yy + yy + zz + zz));
    m[1] =  (xy + xy - zw - zw);
    m[2] =  (xz + xz + yw + yw);
    m[3] = 0.0;
    m[4] =  (xy + xy + zw + zw);
    m[5] =  (1.0 - (zz + zz + xx + xx));
    m[6] =  (yz + yz - xw - xw);
    m[7] = 0.0;
    m[8] =  (xz + xz - yw - yw);
    m[9] =  (yz + yz + xw + xw);
    m[10] =  (1.0 - (yy + yy + xx + xx));
    m[11] =  m[12] = m[13] = m[14] = 0.0f;
    m[15] = 1.0f;
}

// Build a quaternion that forms a rotation around an axis
static inline void gcgAXISTOQUATERNION(VECTOR4 q, float phi, VECTOR3 a) {
    gcgSCALEVECTOR3(q, a, sinf(phi * 0.5f));
    q[3] = cosf(phi * 0.5f);
}

static inline void gcgAXISTOQUATERNION(VECTOR4d q, double phi, VECTOR3d a) {
    register double factor =  sin(phi * 0.5);
    gcgSCALEVECTOR3(q, a, factor);
    q[3] =  cos(phi * 0.5);
}


// Obtains the angle and rotation axis from a quaternion
#define gcgQUATERNIONTOAXIS(angle, a, q) {    \
    register double cos_a = q[3];                       \
    angle = acos(cos_a) * 2;                  \
    register double sin_a = sqrt(1.0 - cos_a * cos_a ); \
    if(fabs(sin_a) < 0.0005) sin_a = 1.0;     \
    gcgSETVECTOR3(a, q[0] / sin_a, q[1] / sin_a, q[2] / sin_a); }



// Build a 3x3 rotation matrix from a quaternion
#define gcgQUATERNIONTOMATRIX3(m, q) {          \
    register double xx = (q)[0] * (q)[0], yy = (q)[1] * (q)[1], zz = (q)[2] * (q)[2];\
    register double xy = (q)[0] * (q)[1], xz = (q)[0] * (q)[2], yz = (q)[1] * (q)[2];\
    register double xw = (q)[0] * (q)[3], yw = (q)[1] * (q)[3], zw = (q)[2] * (q)[3];\
    m[0] =  (1.0 - 2.0 * ( yy + zz));    \
    m[1] =  (2.0 * (xy - zw));           \
    m[2] =  (2.0 * (xz + yw));           \
    m[3] =  (2.0 * (xy + zw));           \
    m[4] =  (1.0 - 2.0 * (zz + xx));     \
    m[5] =  (2.0 * (yz - xw));           \
    m[6] =  (2.0 * (xz - yw));           \
    m[7] =  (2.0 * (yz + xw));           \
    m[8] =  (1.0 - 2.0 * (yy + xx)); }

// Build a 3x4 rotation matrix from a quaternion
#define gcgQUATERNIONTOMATRIX34(m, q) {          \
    register double xx = (q)[0] * (q)[0], yy = (q)[1] * (q)[1], zz = (q)[2] * (q)[2];\
    register double xy = (q)[0] * (q)[1], xz = (q)[0] * (q)[2], yz = (q)[1] * (q)[2];\
    register double xw = (q)[0] * (q)[3], yw = (q)[1] * (q)[3], zw = (q)[2] * (q)[3];\
    m[0] =  (1.0 - 2.0 * ( yy + zz));    \
    m[1] =  (2.0 * (xy - zw));           \
    m[2] =  (2.0 * (xz + yw));           \
    m[3] = 0.0;                                 \
    m[4] =  (2.0 * (xy + zw));           \
    m[5] =  (1.0 - 2.0 * (zz + xx));     \
    m[6] =  (2.0 * (yz - xw));           \
    m[7] = 0.0;                                 \
    m[8] =  (2.0 * (xz - yw));           \
    m[9] =  (2.0 * (yz + xw));           \
    m[10] =  (1.0 - 2.0 * (yy + xx));    \
    m[11] = 0.0;    }


// Build a quaternion from a 4x4 rotation matrix
#define gcgMATRIX4TOQUATERNION(q, mat) {                                  \
  register double trace = 1.0 + (mat)[0] + (mat)[5] + (mat)[10];          \
  if(trace > EPSILON) {                                                   \
    double S = sqrt(trace) * 2.0;                                          \
    (q)[0] = ((mat)[9] - (mat)[6]) / S;                                   \
    (q)[1] = ((mat)[2] - (mat)[8]) / S;                                   \
    (q)[2] = ((mat)[4] - (mat)[1]) / S;                                   \
    (q)[3] = 0.25 * S;                                                    \
  } else                                                                  \
     if((mat)[0] > (mat)[5] && (mat)[0] > (mat)[10]) {                    \
        double S  = sqrt(1.0 + (mat)[0] - (mat)[5] - (mat)[10]) * 2.0;     \
        (q)[0] = 0.25 * S;                                                \
        (q)[1] = ((mat)[4] + (mat)[1] ) / S;                              \
        (q)[2] = ((mat)[2] + (mat)[8] ) / S;                              \
        (q)[3] = ((mat)[9] - (mat)[6] ) / S;                              \
    } else if((mat)[5] > (mat)[10]) {                                                   \
        double S  = sqrt(1.0 + (mat)[5] - (mat)[0] - (mat)[10] ) * 2.0;    \
        (q)[0] = ((mat)[4] + (mat)[1] ) / S;                              \
        (q)[1] = 0.25 * S;                                                \
        (q)[2] = ((mat)[9] + (mat)[6] ) / S;                              \
        (q)[3] = ((mat)[2] - (mat)[8] ) / S;                              \
    } else {                                                                \
        double S  = sqrt(1.0 + (mat)[10] - (mat)[0] - (mat)[5] ) * 2.0;    \
        (q)[0] = ((mat)[2] + (mat)[8] ) / S;                              \
        (q)[1] = ((mat)[9] + (mat)[6] ) / S;                              \
        (q)[2] = 0.25 * S;                                                \
        (q)[3] = ((mat)[4] - (mat)[1] ) / S;                              \
    }                                                                     \
    gcgNORMALIZEQUATERNION(q, q); }


//Diadic tensor product of two 3D vectors
#define gcgTENSORPRODUCT3(result, v1, v2) {   \
        (result)[ 0] = (v1)[ 0] * (v2)[ 0]; \
        (result)[ 1] = (v1)[ 0] * (v2)[ 1]; \
        (result)[ 2] = (v1)[ 0] * (v2)[ 2]; \
        (result)[ 3] = (v1)[ 1] * (v2)[ 0]; \
        (result)[ 4] = (v1)[ 1] * (v2)[ 1]; \
        (result)[ 5] = (v1)[ 1] * (v2)[ 2]; \
        (result)[ 6] = (v1)[ 2] * (v2)[ 0]; \
        (result)[ 7] = (v1)[ 2] * (v2)[ 1]; \
        (result)[ 8] = (v1)[ 2] * (v2)[ 2]; }



// Build a quaternion from Euler angles
static inline void gcgEULERTOQUATERNION(VECTOR4 q, float pitch, float yaw, float roll) {
  VECTOR4 qx = {sinf(pitch * -0.5f), 0.0f, 0.0f,  cosf(pitch  * -0.5f)};
  VECTOR4 qy = {0.0f,  sinf(yaw * -0.5f), 0.0f,  cosf(yaw * -0.5f)};
  VECTOR4 qz = {0.0f, 0.0f,  sinf(roll * -0.5f),  cosf(roll * -0.5f)};
  gcgMULTQUATERNION(qy, qy, qx);
  gcgMULTQUATERNION(q,  qz, qy);
}

static inline void gcgEULERTOQUATERNION(VECTOR4d q, double pitch, double yaw, double roll) {
  VECTOR4d qx = {sin(pitch * -0.5), 0, 0, cos(pitch  * -0.5)};
  VECTOR4d qy = {0,  sin(yaw * -0.5), 0,  cos(yaw * -0.5)};
  VECTOR4d qz = {0, 0,  sin(roll * -0.5), cos(roll * -0.5)};
  gcgMULTQUATERNION(qy, qy, qx);
  gcgMULTQUATERNION(q,  qz, qy);
}


// Vector rotation by quaternion application:  vi = qr * v * qr-1
static inline void gcgAPPLYQUATERNIONVECTOR3(VECTOR4 vi, VECTOR4 qr, VECTOR3 va) {
  VECTOR4 qi, kr;
  gcgCONJUGATEQUATERNION(qi, qr);
  gcgCOPYVECTOR3((float*) kr, va);
  (kr)[3] = 0.0;
  gcgMULTQUATERNION(kr, qr, kr);
  gcgMULTQUATERNION(kr, kr, qi);
  gcgCOPYVECTOR3(vi, kr);
}

static inline void gcgAPPLYQUATERNIONVECTOR3(VECTOR4d vi, VECTOR4d qr, VECTOR3d va) {
  VECTOR4d qi, kr;
  gcgCONJUGATEQUATERNION(qi, qr);
  gcgCOPYVECTOR3((double*) kr, va);
  (kr)[3] = 0.0;
  gcgMULTQUATERNION(kr, qr, kr);
  gcgMULTQUATERNION(kr, kr, qi);
  gcgCOPYVECTOR3(vi, kr);
}

// FUNCTION PROTOTYPES





GCG_API_FUNCTION    bool gcgInverseMatrix2(MATRIX2 inverse, MATRIX2 matrix);

GCG_API_FUNCTION    bool gcgInverseMatrix3(MATRIX3 inverse, MATRIX3 matrix);

GCG_API_FUNCTION    bool gcgInverseMatrix4(MATRIX4 inverse, MATRIX4 matrix);

GCG_API_FUNCTION    bool gcgInverseMatrix2(MATRIX2d inverse, MATRIX2d matrix);

GCG_API_FUNCTION    bool gcgInverseMatrix3(MATRIX3d inverse, MATRIX3d matrix);

GCG_API_FUNCTION    bool gcgInverseMatrix4(MATRIX4d inverse, MATRIX4d matrix);

GCG_API_FUNCTION    bool gcgEigenSymmetric(int norder, float matrix[], float eigenvectors[], float eigenvalues[]);

GCG_API_FUNCTION    bool gcgEigenSymmetric(int norder, double matrix[], double eigenvectors[], double eigenvalues[]);


GCG_API_FUNCTION    bool gcgEigenSymmetricMatrix3(MATRIX3 matrix, MATRIX3 eigenvectors, VECTOR3 eigenvalues);

GCG_API_FUNCTION    bool gcgEigenSymmetricMatrix3(MATRIX3d matrix, MATRIX3d eigenvectors, VECTOR3d eigenvalues);



template<class NUMTYPE> class    GCG_API_CLASS    gcgMATRIX;






GCG_API_FUNCTION    bool gcgDrawLogo(float scale = 0.67);

GCG_API_FUNCTION    bool gcgDrawAABox(VECTOR3 bbmin, VECTOR3 bbmax);

GCG_API_FUNCTION    bool gcgDrawAABox(VECTOR3d bbmin, VECTOR3d bbmax);

GCG_API_FUNCTION    bool gcgDrawOBox(VECTOR3 obbangles, VECTOR3 obbcenter, VECTOR3 obbsize);

GCG_API_FUNCTION    bool gcgDrawOBox(VECTOR3d obbangles, VECTOR3d obbcenter, VECTOR3d obbsize);

GCG_API_FUNCTION    bool gcgDrawVectorPyramid(float x, float y, float z, VECTOR3 vector, float scale);

GCG_API_FUNCTION    bool gcgDrawVectorPyramid(float x, float y, float z, VECTOR3d vector, float scale);

GCG_API_FUNCTION    bool gcgDrawVectorPyramidClosed(float x, float y, float z, VECTOR3 vector, float scale);

GCG_API_FUNCTION    bool gcgDrawVectorPyramidClosed(float x, float y, float z, VECTOR3d vector, float scale);

GCG_API_FUNCTION    bool gcgDrawVectorThetrahedron(float x, float y, float z, VECTOR3 vector, float scale);

GCG_API_FUNCTION    bool gcgDrawVectorThetrahedron(float x, float y, float z, VECTOR3d vector, float scale);

GCG_API_FUNCTION    bool gcgDrawVectorThetrahedronClosed(float x, float y, float z, VECTOR3 vector, float scale);

GCG_API_FUNCTION    bool gcgDrawVectorThetrahedronClosed(float x, float y, float z, VECTOR3d vector, float scale);




GCG_API_FUNCTION    float gcgBlinkingAlpha();

GCG_API_FUNCTION    bool gcgHeatColor(float normheat, VECTOR3 color);

GCG_API_FUNCTION    bool gcgHeatColor(double normheat, VECTOR3d color);



// GCGFONT Versions
#define GCG_FONT_VERSION 0xf010

// Compression modes
#define GCG_FONT_UNCOMPRESSED    0
#define GCG_FONT_COMPRESSED_RLE  1

#pragma pack(push, 1)  // Must have 1 byte alignment

// Structure for the GCG font file format .gff version 1.0.
typedef struct _GCGFONTDATA {
    // Header
    unsigned short  version;        // First two bytes indicate version: 1.0 = 0xf010
    unsigned int    fontdatasize;   // Font data size: header + information + maps
    unsigned char   compression;    // Font image compression mode

    // Font information dependent on version number
    unsigned int    width;          // Font image width
    unsigned int    height;         // Font image height
    unsigned int    cellwidth;      // Character cell width
    unsigned int    cellheight;     // Character cell height
    unsigned char   firstASCII;     // ASCII value of the first character in the font map

    // Information about the chars inside the cell
    unsigned int    xoffset;        // Offset in x for all characters in pixels
    unsigned int    yoffset;        // Offset in y for all characters in pixels
    unsigned int    charheight;     // Height of all characters in font

    // Font data maps as a sequence of chars.
    // nchars = ((unsigned int) width / cellwidth) * ((unsigned int) height / cellheight);
    // unsigned char charwidths[nchars];       // Width of each character texture in pixels
    // unsigned char basewidths[nchars];       // Character bases in pixels
    // unsigned char fontmap[remainingbytes];  // Font image
} GCGFONTDATA;

// Class for internal font handling
class     GCG_API_CLASS    gcgFONT  :  public gcgCLASS {
  public:
    unsigned int  textureID;      // The font texture ID for OpenGL
    unsigned int  nchars;         // Number of characters in this font
    unsigned char firstASCII;     // ASCII code of the first character in this font
    int           font_base;      // Character OpenGL list base.

    unsigned int  charheight;     // Character height in pixels
    unsigned char *basewidths;    // Character widths in pixels.

  public:
    gcgFONT();
    virtual ~gcgFONT();

    // Object initiation
    bool loadGFF(char *fontname);  // Loads a font from a .gff file
    bool createFromData(GCGFONTDATA *font);
    void releaseFont(); // Free all font resources

    // Functions to retrieve string information for rendering
    void textLength(unsigned int *strwidth, unsigned int *strheight, char *str, int strlength); // Computes the box spanned by this string
    int  charsThatFit(float width, char *str, int strlength); // Computes how many chars fit to the width.
};


// Utilities for font creation
GCG_API_FUNCTION    bool gcgSaveGFF(char *name, GCGFONTDATA *gcgfont);
GCG_API_FUNCTION    GCGFONTDATA *gcgImportBFF(char *name); // Imports a .bff font file

#pragma pack(pop)


// Number of fonts internally implemented
#define GCG_SYSTEM_FONTS   12

// Registered system fonts
#define GCG_FONT_SANSSERIF_11_NORMAL        0
#define GCG_FONT_TAHOMA_12_NORMAL           1
#define GCG_FONT_TAHOMA_12_BOLD             2
#define GCG_FONT_TAHOMA_13_NORMAL           3
#define GCG_FONT_COURIERNEW_15_NORMAL       4
#define GCG_FONT_COURIERNEW_15_BOLD         5
#define GCG_FONT_COURIERNEW_15_ITALIC       6
#define GCG_FONT_COURIERNEW_15_BOLDITALIC   7
#define GCG_FONT_TAHOMA_17_NORMAL           8
#define GCG_FONT_TAHOMA_17_BOLD             9
#define GCG_FONT_TAHOMA_17_ITALIC           10
#define GCG_FONT_TAHOMA_17_BOLDITALIC       11

// Class definition
class    GCG_API_CLASS    gcgTEXT   :  public gcgCLASS {
  private:
    float   x, y, z;          // Next print position
    gcgFONT *currentfont;     // Selected font object
    float   scalex, scaley;   // Font scale.
    float   textangle;        // Text orientation
    int     borderpixels[2];  // Border pixels separation

    int     viewport[4];      // OpenGL viewport
    float   textbox[4];       // Output text box

    int     wrap;             // Wrap text to box?
    bool    noTextBox;        // Text box enabled?

  public:
    // Construction
    gcgTEXT();
    virtual ~gcgTEXT();

    // Font parameters
    bool setFont(gcgFONT *font);
    bool setSystemFont(int fontindex);
    void fontScale(float scalex, float scaley);
    float setBestFitSystemFont(char *str, float width); // Selects the best normal system font and scaling that fits to the width passed.
    gcgFONT* getCurrentFont();

    // Text box functions. Use screen-space coordinates (pixels).
    void enableTextBox(float x, float y, float width, float height);  // Restricts the output area to this box.
    void enableTextBoxAt3DPos(float xpos, float ypos, float zpos, float width, float height);
    void adjustTextBox(float xoffset, float yoffset, float widthoffset, float heightoffset);
    void setBorder(int npixelsX, int npixelsY); // Number of pixels for text separation
    void disableTextBox(); // The output area is the entire viewport.
    void drawTextBox(float boxR, float boxG, float boxB, float boxA, float frameR, float frameG, float frameB, float frameA, float linewidth);

    // Text output functions. Use screen-space coordinates (pixels).
    void wrapText(bool wraptobox);  // Text is wrapped to the box if wraptobox is TRUE.
    void textPosition(float x, float y); // Position inside the text box.
    void textOrientation(float angle);   // Position inside the text box and orientation.
    void gcgprintf(const char *format, ...);

    // Information

  private:
    void initTextOutput();   // Initialization for al text output functions
    void finishTextOutput(); // Finishing the text output for all functions
    void outputText(int length, char *str); // Prints text to the output area. Implements the wrap functionality.
};


// Generates C code for internal registration as system font
bool gcgGenerateSystemFontCode(char *name, GCGFONTDATA *gcgfont);



// GCG geometry macros. More efficient than functions. But results in a larger code.

// Add a point to a Axis Aligned Box min/max vectors
#define gcgADDPPOINTAABox(bbmin, bbmax, p) { \
  if((bbmin)[0] > (p)[0]) (bbmin)[0] = (p)[0]; \
  if((bbmin)[1] > (p)[1]) (bbmin)[1] = (p)[1]; \
  if((bbmin)[2] > (p)[2]) (bbmin)[2] = (p)[2]; \
  if((bbmax)[0] < (p)[0]) (bbmax)[0] = (p)[0]; \
  if((bbmax)[1] < (p)[1]) (bbmax)[1] = (p)[1]; \
  if((bbmax)[2] < (p)[2]) (bbmax)[2] = (p)[2]; }


// GCG geometry prototypes.

// Functions for quaternions
GCG_API_FUNCTION    void gcgQuaternionSlerp(VECTOR4 dest, VECTOR4 p, VECTOR4 q, float t);

// Functions for alignment tranformations
// Computes the rotation matrix that aligns the z-axis with the direction specified by vector dir.
GCG_API_FUNCTION    bool gcgComputeAlignMatrix(MATRIX4 matrix, VECTOR3 dir);

// Trackball simulation
GCG_API_FUNCTION    void gcgTrackball(VECTOR4 q, float p1x, float p1y, float p2x, float p2y);

// Functions for bounding boxes
GCG_API_FUNCTION    void gcgAABoxFromOBox(VECTOR3 aabbmin, VECTOR3 aabbmax, VECTOR3 obbangles, VECTOR3 obbposition, VECTOR3 obbmin, VECTOR3 obbmax);

// Functions for computing intersections
GCG_API_FUNCTION    int  gcgIntersectTriangleTriangle(VECTOR3 V0, VECTOR3 V1, VECTOR3 V2, VECTOR3 U0, VECTOR3 U1, VECTOR3 U2);
GCG_API_FUNCTION    int  gcgIntersectTriangleTriangleLine(VECTOR3 V0, VECTOR3 V1, VECTOR3 V2, VECTOR3 U0, VECTOR3 U1, VECTOR3 U2, int *coplanar, VECTOR3 isectpt1, VECTOR3 isectpt2);
GCG_API_FUNCTION    int  gcgIntersectLineLine(VECTOR3 p1, VECTOR3 p2, VECTOR3 p3, VECTOR3 p4, VECTOR3 pa, VECTOR3 pb,  double *mua, double *mub);
GCG_API_FUNCTION    int  gcgIntersectPointLine(VECTOR3 point, VECTOR3 edgeA, VECTOR3 edgeB);
GCG_API_FUNCTION    int  gcgIntersectPointTriangle(VECTOR3 p, VECTOR3 A, VECTOR3 B, VECTOR3 C);
GCG_API_FUNCTION    bool gcgIntersectLineBox2D(VECTOR2 p1, VECTOR2 p2, VECTOR2 box2dMin, VECTOR2 box2dMax);

// Functions for triangles

GCG_API_FUNCTION    float gcgTriangleArea(VECTOR3 v0, VECTOR3 v1, VECTOR3 v2);

GCG_API_FUNCTION    void gcgTriangleGradient(VECTOR3 grad_u, VECTOR3 v0, VECTOR3 v1, VECTOR3 v2, float u0, float u1, float u2);

GCG_API_FUNCTION    float gcgSphericalDistance(float latitude1, float longitude1, float latitude2, float longitude2);


//************* Double precision versions *************
GCG_API_FUNCTION    void gcgQuaternionSlerp(VECTOR4d dest, VECTOR4d p, VECTOR4d q, double t);
GCG_API_FUNCTION    bool gcgComputeAlignMatrix(MATRIX4d matrix, VECTOR3d dir);
GCG_API_FUNCTION    void gcgTrackball(VECTOR4d q, double p1x, double p1y, double p2x, double p2y);
GCG_API_FUNCTION    void gcgAABoxFromOBox(VECTOR3d aabbmin, VECTOR3d aabbmax, VECTOR3d obbangles, VECTOR3d obbposition, VECTOR3d obbmin, VECTOR3d obbmax);
GCG_API_FUNCTION    int  gcgIntersectTriangleTriangle(VECTOR3d V0, VECTOR3d V1, VECTOR3d V2, VECTOR3d U0, VECTOR3d U1, VECTOR3d U2);
GCG_API_FUNCTION    int  gcgIntersectTriangleTriangleLine(VECTOR3d V0, VECTOR3d V1, VECTOR3d V2, VECTOR3d U0, VECTOR3d U1, VECTOR3d U2, int *coplanar, VECTOR3d isectpt1, VECTOR3d isectpt2);
GCG_API_FUNCTION    int  gcgIntersectLineLine(VECTOR3d p1, VECTOR3d p2, VECTOR3d p3, VECTOR3d p4, VECTOR3d pa, VECTOR3d pb,  double *mua, double *mub);
GCG_API_FUNCTION    int  gcgIntersectPointLine(VECTOR3d point, VECTOR3d edgeA, VECTOR3d edgeB);
GCG_API_FUNCTION    int  gcgIntersectPointTriangle(VECTOR3d p, VECTOR3d A, VECTOR3d B, VECTOR3d C);
GCG_API_FUNCTION    bool gcgIntersectLineBox2D(VECTOR2d p1, VECTOR2d p2, VECTOR2d box2dMin, VECTOR2d box2dMax);
GCG_API_FUNCTION    double gcgTriangleArea(VECTOR3d v0, VECTOR3d v1, VECTOR3d v2);
GCG_API_FUNCTION    void gcgTriangleGradient(VECTOR3d grad_u, VECTOR3d v0, VECTOR3d v1, VECTOR3d v2, double u0, double u1, double u2);
GCG_API_FUNCTION    double gcgSphericalDistance(double latitude1, double longitude1, double latitude2, double longitude2);


// Internal structures for mesh representation.

// Class gcgPOLYGONVERTEX: main abstract class for vertex handling.
class gcgPOLYGONVERTEX   :  public gcgCLASS {
  public:
    // Internal use in gcgPOLYGON class. Do not modify.
    unsigned int        nhalf_edges;    // Number of half-edges containing this vertex
    unsigned int        maxhalf_edges;  // Capacity of the half_edges array
    struct _GCGHALFEDGE **half_edges;   // Array of half-edges containing this vertex
    gcgPOLYGONVERTEX    *prevVertex;    // Previous vertex of a mesh
    gcgPOLYGONVERTEX    *nextVertex;    // Next vertex of a mesh

  public:
    virtual ~gcgPOLYGONVERTEX() {}

    // Returns true if this vextex is equal to vertex2
    virtual bool isEqual(gcgPOLYGONVERTEX *vertex2);

    // Copies the contents of vertex2 to vertex
    virtual void copy(gcgPOLYGONVERTEX *vertex2);
};

// Class gcgVERTEX3: main abstract template class for
//                   3D vertex handling.
// - use:
//      gcgVERTEX3<float> for VECTOR3
//      gcgVERTEX3<double> for VECTOR3d
template <typename T>
class gcgVERTEX3 : public gcgPOLYGONVERTEX {
  public:
    T position[3];

  public:
    virtual bool isEqual(T vertex2[3]) {
      return FEQUAL(this->position[0], vertex2[0]) &&
             FEQUAL(this->position[1], vertex2[1]) &&
             FEQUAL(this->position[2], vertex2[2]);
    };

    virtual void copy(T vertex2[3]) {
      this->position[0] = vertex2[0];
      this->position[1] = vertex2[1];
      this->position[2] = vertex2[2];
    };
};


// Internal polygon information
class    GCG_API_CLASS    gcgPOLYGON   :  public gcgCLASS {
  public:
   unsigned int nvertices;         // Number of vertices of this polygon
   struct _GCGHALFEDGE *half_edge; // Half-edges forming this polygon. One per vertex.
   gcgPOLYGON   *prevPolygon;      // Previous polygon of a mesh
   gcgPOLYGON   *nextPolygon;      // Next polygon of a mesh
};


// Internal Half-Edge information
typedef struct _GCGHALFEDGE {
   gcgPOLYGONVERTEX    *srcVertex;           // Vertex to which this half-edge belongs
   gcgPOLYGON          *face;                // Face to which this half-edge belongs
   struct _GCGHALFEDGE *nextHalfEdge;        // Next half-edge of this face
   struct _GCGHALFEDGE *prevHalfEdge;        // Previous half-edge of this face
   struct _GCGHALFEDGE *nextBrotherHalfEdge; // Next face half-edge
   struct _GCGHALFEDGE *prevBrotherHalfEdge; // Previous face half-edge
} GCGHALFEDGE;


// Class gcgPOLYGONMESH: container for polygons and vertices.
class    GCG_API_CLASS    gcgPOLYGONMESH    :  public gcgCLASS {
  public:
    // Polygon information
    unsigned int     npolygons;   // Number of polygons of this object
    gcgPOLYGON       *polygons;   // Linked list of polygons

    // Vertex information
    unsigned int     nvertices;   // Number of vertices of this object
    gcgPOLYGONVERTEX *vertices;   // Linked list of vertices

    // Acess information
    GCGHALFEDGE*     last_half_edge;
    GCGHALFEDGE*     first_half_edge_brother;
    GCGHALFEDGE*     last_half_edge_brother;

  public:
    // Constructors and destructors
    gcgPOLYGONMESH();
    virtual ~gcgPOLYGONMESH();

    // Resources handling.
    void destroyPolygon();

    // Object data access and initialization
    bool addPolygon(gcgPOLYGON *polygon);
    bool addVertex(gcgPOLYGONVERTEX *vertex);
    gcgPOLYGONVERTEX* findSimilarVertex(gcgPOLYGONVERTEX *vertex);

    // Functions to attach vertices to polygons and edges
    bool  addPolygonVertex(gcgPOLYGON *ipolygon, gcgPOLYGONVERTEX *vertex, unsigned int newvertexIndex = 0);

    // Functions to retrieve information
    gcgPOLYGONVERTEX* getPolygonVertex(gcgPOLYGON *polygon, unsigned int vertexIndex);
    gcgPOLYGONVERTEX* getPolygonNextVertex(); // Gets the next edge in the last polygon acessed by getPolygonVertex.
    gcgPOLYGON*       getPolygonFromEdge(gcgPOLYGONVERTEX *vertexID1, gcgPOLYGONVERTEX *vertexID2);
    gcgPOLYGON*       getNextPolygonFromEdge();

    // Methods to remove objects. Keep remaining polygons and edges.
    bool removeVertex(gcgPOLYGONVERTEX *vertex);
    bool removeEdge(gcgPOLYGONVERTEX *vertex1, gcgPOLYGONVERTEX *vertex2);
    bool removePolygon(gcgPOLYGON *ipolygon);
    bool removeHalfEdge(GCGHALFEDGE* ihalf_edge);
    unsigned int removeIsolatedVertices();

    // Stellar Operations
    bool edgeCollapse(gcgPOLYGONVERTEX *vertexID1, gcgPOLYGONVERTEX *vertex2);
    bool halfEdgeCollapse(gcgPOLYGONVERTEX *vertexID1, gcgPOLYGONVERTEX *vertex2, gcgPOLYGONVERTEX *newvertex);
    bool edgeFlip(gcgPOLYGONVERTEX *vertex1, gcgPOLYGONVERTEX *vertex2);
    bool edgeWeld(gcgPOLYGONVERTEX *vertex);
    bool edgeSplit(gcgPOLYGONVERTEX *vertex1, gcgPOLYGONVERTEX *vertex2, gcgPOLYGONVERTEX *newvertex);
    bool polygonSplit(gcgPOLYGON *ipolygon, unsigned int vertexIndex, gcgPOLYGONVERTEX *newvertex);

    // Test link condition for edge collapse operation
    bool testLinkCondition(gcgPOLYGONVERTEX *vertex1, gcgPOLYGONVERTEX *vertex2);

    // Return the valence of ivertex
    unsigned int getValence(gcgPOLYGONVERTEX *vertexID);
    // Get the star of a vertex
    // starVertices - Returns a list of vertices
    // starPolygons - Returns a list of the faces
    bool getVertexStar(gcgPOLYGONVERTEX *vertex, gcgPOLYGONVERTEX **starVertexVertices, unsigned int* numVertices, gcgPOLYGON **starVertexPolygons, unsigned int* numPolygons);
    unsigned int getVerticesVertexStar(gcgPOLYGONVERTEX *vertexIndex, gcgPOLYGONVERTEX **starVertexVertices, unsigned int* numVertices);
    unsigned int getPolygonsVertexStar(gcgPOLYGONVERTEX *vertexID, gcgPOLYGON **starVertexPolygons, unsigned int* numPolygons);
    // Get the star of an edge
    unsigned int getEdgeStar(gcgPOLYGONVERTEX *vertex1, gcgPOLYGONVERTEX *vertex2, gcgPOLYGONVERTEX **starEdgeVertices, unsigned int* numVertices, gcgPOLYGON **starEdgePolygons, unsigned int* numPolygons);
    // Return the intersection between star1 and star2
    bool getStarIntersection(gcgPOLYGONVERTEX** starVertexVertices1, gcgPOLYGON** starVertexPolygons1, gcgPOLYGONVERTEX** starVertexVertices2, gcgPOLYGON** starVertexPolygons2, gcgPOLYGONVERTEX** starIntersectionVertices, gcgPOLYGON** starIntersectionPolygons, unsigned int numVertices1, unsigned int numPolygons1, unsigned int numVertices2, unsigned int numPolygons2, unsigned int* numVerticesIntersection, unsigned int* numPolygonsIntersection);
    // Compare two stars
    bool starComparison(gcgPOLYGONVERTEX** starVertices1, gcgPOLYGON** starPolygons1, gcgPOLYGONVERTEX** starVertices2, gcgPOLYGON** starPolygons2, unsigned int numVertices1, unsigned int numPolygons1, unsigned int numVertices2, unsigned int numPolygons2);

  private:// Check functions
    void testVertex();
    void testVertex(unsigned int vertexIndex);
    void testHalfEdgeBrother(GCGHALFEDGE* ihalf_edge);
    void testPolygons();
    void testStar(unsigned int *starVertex[2], unsigned int numVertices, unsigned int numPolygons);
    void testStarVertices(gcgPOLYGONVERTEX** starVertexVertices, unsigned int numVertices);
    void testStarVerticesRemove(gcgPOLYGONVERTEX** starVertexVertices, unsigned int numVertices);
    void testStarPolygons(gcgPOLYGON** starVertexPolygons, unsigned int numPolygons);
    void testStarPolygonsRemove(gcgPOLYGON** starVertexPolygons, unsigned int numPolygons);
    void testHalfEdgeBrother2();
    void testHalfEdgeBrother3();
};



// Vertex indices
#define GCG_VAPEX   0
#define GCG_VESQ    1
#define GCG_VDIR    2

// Sides for connections on level 2:
//    Defining the base of a side (NORTH,WEST) as the base of the opposite side,
//  (SOUTH, EAST) we force the subdivisions of the first be propagated to the
//   other.
//    As such, we can change the topology of the generated mesh as follows:
//
//        Plane:    no forced connection
//        Cylinder: WEST is connected to EAST
//        Torus:    WEST is connected to EAST and NORTH is connected to SOUTH
//
//            NORTH
//          p3------p2
//           | \  / |
//    WEST   |  \/  |   EAST
//           |  /\  |
//           | /  \ |
//          p0------p1
//            SOUTH
//
#define GCG_TOPOLOGY_PLANE      0     // Default topology
#define GCG_TOPOLOGY_CYLINDER   1
#define GCG_TOPOLOGY_TORUS      2

// Flags for composing mask actions for vertex indexes of an adaptive object.
//            - GCG_REFINE: this triangle is too coarse and should be refined.
//            - GCG_SIMPLIFY: this triangle is too detailed and should be simplified.
//            - GCG_MAXIMUMLOD: this triangle reached its maximum level of detail.
//                              Avoid refining the sons of this node. Note that these
//                              may be refined to keep the mesh consistent.
//            - GCG_MINIMUMLOD: this triangle reached its minimum level of detail.
//                              Avoid simplifying the parents of this node. Note
//                              that these may be simplified to keep the mesh consistent.
//            - GCG_INVISIBLE: informs that this triangle is completely invisible
//            - GCG_TOTALLYVISIBLE: informs that this triangle is completely visible
//

#define GCG_REFINE           ((unsigned char) 0x01)
#define GCG_SIMPLIFY         ((unsigned char) 0x02)
#define GCG_MAXIMUMLOD       ((unsigned char) 0x04)
#define GCG_MINIMUMLOD       ((unsigned char) 0x08)
#define GCG_INVISIBLE        ((unsigned char) 0x10)
#define GCG_TOTALLYVISIBLE   ((unsigned char) 0x20)
#define GCG_SELECT0          ((unsigned char) 0x40)
#define GCG_SELECT1          ((unsigned char) 0x80)



//  gcgADAPTIVE: abstract interface to support refinable meshes
// - informs the geometry of a refinable object
// - must allocate all vertices.
// - all methods are abstract
class    GCG_API_CLASS    gcgADAPTIVE    :  public gcgCLASS {
  public:

    // beginTesselation() is called in the beginning of a tesselation to
    // obtain the domain for subdivision in the form of four vertices. Its call
    // indicates that the previous tesselation was discarded and all vertices
    // may be freed.
    //  - For the semi-regular A48 mesh, the 4 vertices form a losangle
    //    as follows:
    //                 d3-----d2
                //                  |   /  |
                //                  |  /   |
                //                  | /    |
                //                 d0-----d1
                //
    //  - For every tesselation, this method must allocate the vertices as
    //    described above.
    virtual void beginTesselation() = 0;

    // action() is called to query if the triangle formed by counterclockwise
    // indexes iVertices[VAPEX], iVertices[VESQ], iVertices[DIR], should be
    // REFINED or SIMPLIFIED. It must return a byte mask with bits composed by
    // any combination of the flags:
    //      - GCG_REFINE: this triangle is too coarse and should be refined.
    //      - GCG_SIMPLIFY: this triangle is too detailed and should be simplified.
    //      - GCG_MAXIMUMLOD: this triangle reached its maximum level of detail.
    //                        Avoid refining the sons of this node. Note that these
    //                        may be refined to keep the mesh consistent.
    //      - GCG_MINIMUMLOD: this triangle reached its minimum level of detail.
    //                        Avoid simplifying the parents of this node. Note
    //                        that these may be simplified to keep the mesh consistent.
    //      - GCG_INVISIBLE: informs that this triangle is completely invisible
    //      - GCG_TOTALLYVISIBLE: informs that this triangle is completely visible
    //
    // It must return 0 if no action should be applied.
    virtual unsigned char action(unsigned char nivel, int *iVertices, unsigned char oldmask) = 0;

    // sample() must compute a NEW POINT for triangle subdivision
    // (splitting), and return its integer index. If the parameter apex2 is
    // positive, then it's the index of the apex of the triangle that forms
    // a losangle (diamond) with the triangle vertices pointed by iVertices.
    //
                // iVertices[VAPEX] +------+ iVertices[VDIR]
                //                  |    / |
                //                  |   /  |                         - VAPEX = 0
                //                  |  /   |                         - VESQ  = 1
                //                  | /    |                         - VDIR  = 2
                //  iVertices[VESQ] +------+ apex2 (if positive)
                //
                // - apex2 may be negative. In this case only one triangle is used in the
    //   subdivision.
    //
                // - it must return the interpolated point using the four points, case apex2 >= 0,
    //   or using the triangle vertices, otherwise. This index will be used for
    //   further subdivision queries or sampling.
                //
    virtual int sample(int *iVertices, int apex2) = 0;

    // endTesselation() is called when a tesselation has finished. Can be used
    // to adequate the object status.
    //
    virtual void endTesselation() = 0;
};


// Class gcgSEMIREGULAR48: restricted unbalanced binary trees
//    - forms a semi-regular A48 mesh.
//    - controls the top-down subdivision of the mesh
//    - define the mesh topology
class    GCG_API_CLASS    gcgSEMIREGULAR48   :  public gcgCLASS {
  private:
    // Para formar a estrutura básica da árvore
    typedef struct __nodo {
      // Most used fields first
      int           iInfo;            // Index for front info, if this node
                                      // is at tree front then iInfo >= 0

      unsigned char nivel;            // Nivel (detalhe) desse nodo na árvore
      //unsigned int  indice;         // Índice único (estrutural) que
                                      //  identifica esse triangulo

      // Controle da estrutura de subdivisao
      int           filhoesq;         // Filho esquerdo do nodo
                                      // Filho direito do nodo = filhoesq + 1 sempre

      int           vizesq,
                    vizdir, vizbase;  // Indices para os vizinhos do nodo
    } NODO;

    // Information for front nodes
    typedef struct __frontinfo {
      unsigned char lastmask;       // Last flags returned by an action call
    } FRONTINFO;

   // Connection control for topology
    unsigned char topology;

    gcgADAPTIVE *geometria;   // Objeto corrente que fornece a geometria

    int           nnodos;       // Quantidade de nodos atual
    NODO          *nodos;       // Buffer dinamico para controle dos nodos
    int           maxnodos;     // Capacidade do buffer de nodos

    int           ntriangulos;  // Quantidade de triangulos atual
    int           maxtriang;    // Capacidade do buffer de triangulos
    FRONTINFO     *front;       // Front information: one for each triangle
    int           *vertices;    // Buffer dinamico para controle dos vertices dos triangulos

    int           *stack;       // Stack of integers for depth search
    int           top;          // Top of the stack = next node to be expanded
    int           maxstack;     // Stack capacity

    // For triangle selection
    unsigned char hidemask;     // Front triangle selection: if a logical AND of
                                // this mask with the lastmask is non-zero,
                                // then it's not returned in a getIndexArray call.
    unsigned char showmask;     // A logical AND with this mask and the lastmask,
                                // determines those triangles that must be shown.
                                // Only used if hidemask is non-zero.
    int           nhidden;      // Number of hidden elements in a getIndexArray call.

  // Metodos privados /////////////////////////////////////////////////////////
  private:
          void divide_node(int inodo);
        void triangle_bissection(int inodo, int novovertice, int ifilhoesq);
        int  initTesselation(gcgADAPTIVE *geo);


        // Metodos públicos ////////////////////////////////////////////////////////
  public:
    // Creates a semiregular 4-8 mesh which domain and geometry is
    // defined by a gcgADAPTIVE object
    gcgSEMIREGULAR48();
    virtual ~gcgSEMIREGULAR48();

    void setTopology(int topologycode);

    // Generates a mesh WITH a restriction on the number of triangles using
    // a breath-first search.
    //    - uses the gcgADAPTIVE interface to build the mesh.
    //    - Returns the number of generated triangles.
    //
          int tesselate(int maxtriang, gcgADAPTIVE *geometry);

    // Generates a mesh WITHOUT a restriction on the number of triangles using
    // a depth-first search. It is faster than the other version, but the stop
    // criteria are exclusively implemented by the refinable object.
    //
    //    - uses the gcgADAPTIVE interface to build the mesh.
    //    - Returns the number of generated triangles.
    //
          int tesselate(gcgADAPTIVE *geometry);

          // Returns the vertex array of last tesselation. Stores the number of
    // triangles (vextex triplets) in the 'ntriang' parameter if it is non-null.
          int* getIndexArray(int *ntriang);

          // Selects the triangles that will be delivered by getIndexArray.
          int hideTriangles(unsigned char hidemask, unsigned char mustshowmask);
};




template<class NUMTYPE> class    GCG_API_CLASS    gcgDISCRETE1D;
template<class NUMTYPE> class    GCG_API_CLASS    gcgDISCRETE2D;



#define GCG_BORDER_EXTENSION_ZERO               1

#define GCG_BORDER_EXTENSION_CLAMP              2

#define GCG_BORDER_EXTENSION_PERIODIC           3

#define GCG_BORDER_EXTENSION_SYMMETRIC_NOREPEAT 4

#define GCG_BORDER_EXTENSION_SYMMETRIC_REPEAT   5




template<class NUMTYPE>    GCG_API_FUNCTION    bool gcgLinearSolver(gcgDISCRETE2D<NUMTYPE> *A, gcgDISCRETE1D<NUMTYPE> *b, gcgDISCRETE1D<NUMTYPE> *x);


template <class NUMTYPE = float> class    GCG_API_CLASS    gcgSEQUENCEDERIVATIVE2D;
template <class NUMTYPE = float> class    GCG_API_CLASS    gcgFILTERMASK1D;

// Signal comparison: error computation.

template<class NUMTYPE>    GCG_API_FUNCTION    double computeMSEwith(gcgDISCRETE1D<NUMTYPE> *src1, gcgDISCRETE1D<NUMTYPE> *src2); // Computes the Mean Squared Error.
template<class NUMTYPE>    GCG_API_FUNCTION    double computePSNRwith(gcgDISCRETE1D<NUMTYPE> *src1, gcgDISCRETE1D<NUMTYPE> *src2);// Computes the Peak Signal-to-Noise Ratio in dB.
template<class NUMTYPE>    GCG_API_FUNCTION    double computeMAEwith(gcgDISCRETE1D<NUMTYPE> *src1, gcgDISCRETE1D<NUMTYPE> *src2); // Computes the Mean Absolute Error.

template<class NUMTYPE>    GCG_API_FUNCTION    double computeMSEwith(gcgDISCRETE2D<NUMTYPE> *src1, gcgDISCRETE2D<NUMTYPE> *src2); // Computes the Mean Squared Error.
template<class NUMTYPE>    GCG_API_FUNCTION    double computePSNRwith(gcgDISCRETE2D<NUMTYPE> *src1, gcgDISCRETE2D<NUMTYPE> *src2);// Computes the Peak Signal-to-Noise Ratio in dB.
template<class NUMTYPE>    GCG_API_FUNCTION    double computeMAEwith(gcgDISCRETE2D<NUMTYPE> *src1, gcgDISCRETE2D<NUMTYPE> *src2); // Computes the Mean Absolute Error.


// Signal analysis.

// Computes the gradient of 1D signal src. The derivative is computed using well-known Sobel operator and
// the result is stored in dx.
template<class NUMTYPE>    GCG_API_FUNCTION    bool gcgFirstDerivativeSobel1D(gcgDISCRETE1D<NUMTYPE> *src, gcgDISCRETE1D<NUMTYPE> *dx);

// Computes the gradient of 2D signal src. The partial derivatives are computed using well-known Sobel
// operator and the components are stored in dx and dy. The pointer for dx or dy may be NULL.
template<class NUMTYPE>    GCG_API_FUNCTION    bool gcgGradientSobel2D(gcgDISCRETE2D<NUMTYPE> *src, gcgDISCRETE2D<NUMTYPE> *dx, gcgDISCRETE2D<NUMTYPE> *dy);

// Computes the (optical) flow given the partial derivatives in X, Y and time of up to three channels.
// Uses the Augereau method (Bertrand Augereau, Benoît Tremblais, Christine Fernandez-Maloigne,
// "Vectorial computation of the optical flow in colorimage sequences", CIC05).
template<class NUMTYPE>    GCG_API_FUNCTION    bool gcgOpticalFlowAugereau2D(gcgDISCRETE2D<NUMTYPE> *dx1, gcgDISCRETE2D<NUMTYPE> *dy1, gcgDISCRETE2D<NUMTYPE> *dt1,
                                                         gcgDISCRETE2D<NUMTYPE> *dx2, gcgDISCRETE2D<NUMTYPE> *dy2, gcgDISCRETE2D<NUMTYPE> *dt2,
                                                         gcgDISCRETE2D<NUMTYPE> *dx3, gcgDISCRETE2D<NUMTYPE> *dy3, gcgDISCRETE2D<NUMTYPE> *dt3,
                                                         gcgDISCRETE2D<float> *outflowX, gcgDISCRETE2D<float> *outflowY);
template<class NUMTYPE>    GCG_API_FUNCTION    bool gcgOpticalFlowAugereau2D(gcgDISCRETE2D<NUMTYPE> *dx1, gcgDISCRETE2D<NUMTYPE> *dy1, gcgDISCRETE2D<NUMTYPE> *dt1,
                                                         gcgDISCRETE2D<NUMTYPE> *dx2, gcgDISCRETE2D<NUMTYPE> *dy2, gcgDISCRETE2D<NUMTYPE> *dt2,
                                                         gcgDISCRETE2D<NUMTYPE> *dx3, gcgDISCRETE2D<NUMTYPE> *dy3, gcgDISCRETE2D<NUMTYPE> *dt3,
                                                         gcgDISCRETE2D<double> *outflowX, gcgDISCRETE2D<double> *outflowY);




template <class NUMTYPE = float>
class    GCG_API_CLASS    gcgBASIS1D   :  public gcgCLASS {
  virtual unsigned int getNumberOfCoefficients() = 0;

  virtual bool setNumberOfSamples(unsigned int nsamples) = 0;

  virtual bool projectSignal(int atX, gcgDISCRETE1D<float> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef) = 0;
  virtual bool projectSignal(int atX, gcgDISCRETE1D<double> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef) = 0;
  virtual bool projectSignal(int atX, gcgDISCRETE1D<short> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef) = 0;
  virtual bool projectSignal(int atX, gcgDISCRETE1D<int> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef) = 0;
  virtual bool projectSignal(int atX, gcgDISCRETE1D<long> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef) = 0;

  virtual bool reconstructSignal(int atX, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE1D<float> *outputvector) = 0;
  virtual bool reconstructSignal(int atX, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE1D<double> *outputvector) = 0;
  virtual bool reconstructSignal(int atX, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE1D<short> *outputvector) = 0;
  virtual bool reconstructSignal(int atX, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE1D<int> *outputvector) = 0;
  virtual bool reconstructSignal(int atX, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE1D<long> *outputvector) = 0;
};

template <class NUMTYPE = float>
class    GCG_API_CLASS    gcgBASIS2D   :  public gcgCLASS {
  virtual unsigned int getNumberOfCoefficients() = 0;

  virtual bool projectSignal(int atX, int atY, gcgDISCRETE2D<float> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef) = 0;
  virtual bool projectSignal(int atX, int atY, gcgDISCRETE2D<double> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef) = 0;
  virtual bool projectSignal(int atX, int atY, gcgDISCRETE2D<short> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef) = 0;
  virtual bool projectSignal(int atX, int atY, gcgDISCRETE2D<int> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef) = 0;
  virtual bool projectSignal(int atX, int atY, gcgDISCRETE2D<long> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef) = 0;

  virtual bool reconstructSignal(int atX, int atY, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE2D<float> *outputvector) = 0;
  virtual bool reconstructSignal(int atX, int atY, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE2D<double> *outputvector) = 0;
  virtual bool reconstructSignal(int atX, int atY, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE2D<short> *outputvector) = 0;
  virtual bool reconstructSignal(int atX, int atY, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE2D<int> *outputvector) = 0;
  virtual bool reconstructSignal(int atX, int atY, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE2D<long> *outputvector) = 0;
};

template<class NUMTYPE> class    GCG_API_CLASS    gcgLEGENDREBASIS1D;
template<class NUMTYPE> class    GCG_API_CLASS    gcgLEGENDREBASIS2D;
template<class NUMTYPE> class    GCG_API_CLASS    gcgDWTBASIS1D;







class    GCG_API_CLASS    gcgIMAGE   :  public gcgCLASS {
  public:
    unsigned char  *data;         
    unsigned int   width;         
    unsigned int   height;        
    unsigned int   bpp;           
    unsigned int   rowsize;       

    unsigned int   palettecolors; 
    unsigned char  *palette;      

    unsigned int    colormasks[4];
    unsigned int    max[4];       

  private:
    unsigned int    rightshift[4];
    unsigned int    nbits[4];     

  public:
    gcgIMAGE();

    virtual ~gcgIMAGE();

    bool createImage(unsigned int width, unsigned int height, unsigned int bpp, bool usealpha = false);


    bool createSimilar(gcgIMAGE *source);

    bool destroyImage();

    bool isCompatibleWith(gcgIMAGE *source);

    bool isValid();

    bool duplicateImage(gcgIMAGE *source);

    bool duplicateSubimage(gcgIMAGE *source, int left, int top, unsigned int width, unsigned int height);

    bool copyImage(gcgIMAGE *source);

        bool loadImage(const char *filename);

    bool loadJPG(const char *filename);

    bool loadBMP(const char *filename);

    bool loadPCX(const char *filename);

    bool loadTGA(const char *filename);

    bool saveBMP(const char *filename);

    bool saveJPG(const char *filename , int quality);

    bool getPixelColor(unsigned int i, unsigned int j, VECTOR4 color);

    bool setPixelColor(unsigned int i, unsigned int j, VECTOR4 color);

    int getPixelIndex(unsigned int i, unsigned int j);

    bool setPixelIndex(unsigned int i, unsigned int j, int newindex);

    bool getPaletteColor(int index, VECTOR4 color);

    bool setPaletteColor(int index, VECTOR4 color);

    bool convolutionX(gcgIMAGE *source, gcgDISCRETE1D<float> *mask, float addthis = 0);

    bool convolutionY(gcgIMAGE *source, gcgDISCRETE1D<float> *mask, float addthis = 0);

    bool convolutionXY(gcgIMAGE *source, gcgDISCRETE2D<float> *mask, float addthis = 0);

    bool templateMatching(unsigned int imgleft, unsigned int imgbottom, unsigned int imgwidth, unsigned int imgheight, gcgDISCRETE2D<float> *mask, VECTOR2 position);

    bool transformColorSpace(gcgIMAGE *source, MATRIX4 matrix);

    bool convertBits(gcgIMAGE *source, unsigned int newbpp);

    bool convertGrayScale(gcgIMAGE *source);

    bool convertGrayScale8bits(gcgIMAGE *source);

    bool isGrayScale();

    bool forceLinearPalette();

    bool scale(float weight, gcgIMAGE *source, float addthis = 0.0);

    bool combineAdd(gcgIMAGE *source1, gcgIMAGE *source2, float weight1 = 1.0, float weight2 = 1.0, float addthis = 0.0);

    bool combineMult(gcgIMAGE *source1, gcgIMAGE *source2, float add1 = 0.0, float add2 = 0.0, float addthis = 0.0);

    bool splitChannels(gcgIMAGE *red, gcgIMAGE *green, gcgIMAGE *blue, gcgIMAGE *alpha = NULL);

    bool mergeChannels(gcgIMAGE *red, gcgIMAGE *green, gcgIMAGE *blue, gcgIMAGE *alpha = NULL);

    bool exportChannels(gcgDISCRETE2D<float> *red, gcgDISCRETE2D<float> *green, gcgDISCRETE2D<float> *blue, gcgDISCRETE2D<float> *alpha);

    bool exportChannels(gcgDISCRETE2D<double> *red, gcgDISCRETE2D<double> *green, gcgDISCRETE2D<double> *blue, gcgDISCRETE2D<double> *alpha);

    bool exportChannels(gcgDISCRETE2D<short> *red, gcgDISCRETE2D<short> *green, gcgDISCRETE2D<short> *blue, gcgDISCRETE2D<short> *alpha);

    bool exportChannels(gcgDISCRETE2D<int> *red, gcgDISCRETE2D<int> *green, gcgDISCRETE2D<int> *blue, gcgDISCRETE2D<int> *alpha);

    bool exportChannels(gcgDISCRETE2D<long> *red, gcgDISCRETE2D<long> *green, gcgDISCRETE2D<long> *blue, gcgDISCRETE2D<long> *alpha);

    bool importChannels(gcgDISCRETE2D<float> *red, gcgDISCRETE2D<float> *green, gcgDISCRETE2D<float> *blue, gcgDISCRETE2D<float> *alpha);

    bool importChannels(gcgDISCRETE2D<double> *red, gcgDISCRETE2D<double> *green, gcgDISCRETE2D<double> *blue, gcgDISCRETE2D<double> *alpha);

    bool importChannels(gcgDISCRETE2D<short> *red, gcgDISCRETE2D<short> *green, gcgDISCRETE2D<short> *blue, gcgDISCRETE2D<short> *alpha);

    bool importChannels(gcgDISCRETE2D<int> *red, gcgDISCRETE2D<int> *green, gcgDISCRETE2D<int> *blue, gcgDISCRETE2D<int> *alpha);

    bool importChannels(gcgDISCRETE2D<long> *red, gcgDISCRETE2D<long> *green, gcgDISCRETE2D<long> *blue, gcgDISCRETE2D<long> *alpha);

    // Histogram processing.
    // The srcimage can be the destination (this).
    bool histogramRGBA(unsigned int hvectorR[256], unsigned int hvectorG[256], unsigned int hvectorB[256], unsigned int hvectorA[256]); // Computes the histogram for each 8bit channel.
    bool histogramGray(unsigned int hvector[256]); // Computes the grayscale histogram (but does not convert the image to grayscale).
    int  histogramIndex(unsigned int vectorsize, unsigned int *hvector); // Paletted images only. Returns the number of bins or zero if image is invalid.
    bool equalizeHistogram(gcgIMAGE *source);


    // Binarization.
    // The srcimage can be the destination (this).
    bool binarizeRGBA(gcgIMAGE *srcimage, int threR, int threG, int threB, int threA); // Binarize RGBA channels independently. Negative values indicate that the channel should not be binarized.
    bool binarizeGray(gcgIMAGE *srcimage, int threshold);  // Binarize this image using the RGB average.
    bool binarizeIndex(gcgIMAGE *srcimage, int threshold); // Binarize the index values. For paletted images only.

    // Image bidimensional support processing functions.
    // The srcimage can be the destination (this).
    bool verticalFlip();  // Mirrors the image in the vertical direction.

    // Utilities
    bool unpackBMP(void *bmpHeader, unsigned char *bmpdata, unsigned int npalcolors, unsigned char *pal);

        private:
    // Load specific formats
    bool loadUncompressed8BitTGA(FILE *file);
    bool loadUncompressedTrueColorTGA(FILE *file);
    bool loadCompressedTrueColorTGA(FILE *file);

    // Unpacking BITMAP
    bool unpackPalettedBMP(void *bmpHeader, unsigned char *bmpdata, unsigned int npalcolors, unsigned char *pal);
    bool unpack16BitBMP(void *bmpHeader, unsigned char *bmpdata, unsigned int npalcolors, unsigned char *pal);
    bool unpack24BitBMP(void *bmpHeader, unsigned char *bmpdata);
    bool unpack32BitBMP(void *bmpHeader, unsigned char *bmpdata, unsigned int npalcolors, unsigned char *pal);

    // Specific converters
    bool convertPalettedto8Bits(gcgIMAGE *srcimage);
    bool convertPalettedto24Bits(gcgIMAGE *srcimage);
    bool convert16Bitsto24Bits(gcgIMAGE *srcimage);
    bool convert32Bitsto24Bits(gcgIMAGE *srcimage);
    bool convert16Bitsto32Bits(gcgIMAGE *srcimage);
};





GCG_API_FUNCTION    unsigned int gcgPackRLE8(unsigned int fullsize, unsigned char *srcdata, unsigned char *RLEdata);

GCG_API_FUNCTION    unsigned int gcgUnpackRLE8(unsigned char *RLEdata, unsigned char *dstdata);

GCG_API_FUNCTION    unsigned int gcgPackRLE32(unsigned int fullsize, unsigned char *srcdata, unsigned char *RLEdata);

GCG_API_FUNCTION    unsigned int gcgUnpackRLE32(unsigned char *RLEdata, unsigned char *dstdata);


GCG_API_FUNCTION    unsigned int gcgCompressImageRLE8(unsigned int width, unsigned int height, unsigned char *srcdata, unsigned char *RLEdata);

GCG_API_FUNCTION    unsigned int gcgDecompressImageRLE8(unsigned int width, unsigned int height, unsigned char *RLEdata, unsigned char *dstdata);

GCG_API_FUNCTION    unsigned int gcgDecompressImageRLE4(unsigned int width, unsigned int height, unsigned char *RLEdata, unsigned char *dstdata);



GCG_API_FUNCTION    bool gcgSaveBytesAsText(unsigned int size, unsigned char *data, char *outputname);

GCG_API_FUNCTION    bool gcgSaveFloatsAsText(unsigned int size, float *data, char *outputname);

GCG_API_FUNCTION    bool gcgSaveBMPcode(char *inputname, char *outputname);




// Macros for plotting graphs
#define   GCG_PLOT_MAX_LINES      16  // In any implementation, supremum is 255 lines
#define   GCG_PLOT_MAX_AXIS        3  // Prepared for 3D plots

// Macros for range definition
#define   GCG_PLOT_X      0
#define   GCG_PLOT_Y      1
#define   GCG_PLOT_Z      2

#define   GCG_PLOT_RANGE_INI   0
#define   GCG_PLOT_RANGE_END   1

class    GCG_API_CLASS    gcgPLOT   :  public gcgCLASS {
  private:
    // Basic struct for independent lines
    typedef struct __line {
      char            *linename;    // Name to identify this line
      float           r, g, b, a;   // Color for plotting this line
      float           linewidth;    // Line width
      unsigned short  bitpattern;   // Pattern for lines

      VECTOR3         *data;        // Pointer to point data buffer
      unsigned int    maxdata;      // Data buffer capacity
      unsigned int    ndata;        // Number of points in buffer
    } LINE;

    // Basic struct for axis
    typedef struct __axis {
      char            *title;       // Title to identify this axis
      float           r, g, b, a;   // Color for plotting this title
      float           linewidth;    // Axis width
      unsigned short  bitpattern;   // Pattern for lines
      float           ini;          // Range start
      float           end;          // Range end
    } AXIS;

    // Lines buffer
    LINE lines[GCG_PLOT_MAX_LINES];
    unsigned int nlines;     // Current number of lines

    // Plot axis
    AXIS  axis[GCG_PLOT_MAX_AXIS];  // For each axis, information and a real interval.

    // Plot info
    float   framebox[4];                    // Output frame box
    float   linewidth;                      // Frame line width
    float   frameR, frameG, frameB, frameA; // Frame color
    float   plotR, plotG, plotB, plotA;     // Plot rectangle color

    unsigned int gridpixelsX;               // Interval in pixels between two consecutive X values
    unsigned int gridpixelsY;               // Interval in pixels between two consecutive Y values
    float   gridR, gridG, gridB, gridA;     // Grid color

    // Text output
    gcgTEXT textoutput;                     // Text output object
    unsigned int textWidthX, textHeightX;   // Text extent in X
    unsigned int textWidthY, textHeightY;   // Text extent in Y
    float        factorX;                   // Text factor for X extents

  public:
    // Constructor
    gcgPLOT();
    virtual ~gcgPLOT();

    // INPUT DATA SPECIFICATION
    int newLine();  // Adds a new line to this plot
    bool addPoint2D(unsigned int lineindex, float x, float y);
    void setAxisRange(unsigned int axis, float ini, float end);
    void clipOutOfRange(unsigned int lineindex, unsigned int iaxis);
    void rangeFromData(unsigned int iaxis);
    void clearData(unsigned int lineindex);
    void clearAllData();

    // DRAW PARAMETERS SPECIFICATION
    void setFrame(float x, float y, float width, float height, float linewidth);
    void adjustFrame(float xoffset, float yoffset, float widthoffset, float heightoffset);
    void setFrameColor(float r, float g, float b, float a);

    void setPlotColor(float r, float g, float b, float a);
    void setGrid(unsigned int pixelsx, unsigned int pixelsy);
    void setGridColor(float r, float g, float b, float a);

    void setAxis(unsigned int axis, char *title, float r, float g, float b, float a, float linewidth, unsigned short bitpattern);
    void setLine(unsigned char lineindex, char *linename, float r, float g, float b, float a, float linewidth, unsigned short bitpattern);

    // DRAW FUNCTIONS
    void draw();
};




// Generate random numbers uniformly distributed with period (2^19937) - 1
//
//  M. Matsumoto & T. Nishimura, "Mersenne Twister: A 623-Dimensionally
//   Equidistributed Uniform Pseudo-Random Number Generator", ACM Transactions
//   on Modeling and Computer Simulation, vol. 8, no. 1, 1998, pp. 3-30.
class    GCG_API_CLASS    gcgRANDOM   :  public gcgCLASS {
  #define GCG_MERSENNE_N          624
  #define GCG_MERSENNE_MATRIX_A   0x9908b0df  // Constant vector a

  private:
    unsigned long mt[GCG_MERSENNE_N]; // The array for the state vector
    unsigned long mag01[2];           // mag01[x] = x * MATRIX_A for x=0,1
    int mti;                          // Index for Mersenne Twisted state vector
    unsigned int lastInterval;        // Last interval length for intRandom
    unsigned int RLimit;              // Rejection limit used by intRandom

  public:
    // Constructor
    gcgRANDOM();                      // Use C conventional rand() to define seed
    gcgRANDOM(unsigned int seed);
    virtual ~gcgRANDOM();

    // SERIES SPECIFICATION
    void  setSeed(unsigned int seed);

    // Random number generation
    unsigned int bitRandom();         // Generate 32 random bits
    double random();                  // Output random float in U[0, 1)

    // Random number generation in an interval
    int intRandom(int min, int max);  // Output random integer in U[min, max]
    float floatRandom(float min, float max);  // Output random float in U[min, max]
    double doubleRandom(double min, double max);  // Output random double in U[min, max]
};



// Generate random numbers with Gaussian distribution with 0 mean and 1.0 as
// standard-deviation N(0, 1).
class    GCG_API_CLASS    gcgRANDOMGAUSSIAN   :  public gcgCLASS {
  private:
    gcgRANDOM ngen1;    // Independent random generators
    gcgRANDOM ngen2;
    bool hasValue;      // Flags the existence of a second value.
    double secondValue; // Keeps the second value

  public:
    // Constructor
    gcgRANDOMGAUSSIAN();  // Use C conventional rand() to define seed
    gcgRANDOMGAUSSIAN(unsigned int seed1, unsigned int seed2); // seed1 must be different than seed2.
    virtual ~gcgRANDOMGAUSSIAN();

    // SERIES SPECIFICATION: seed1 must be different than seed2.
    void  setSeed(unsigned int seed1, unsigned int seed2);

    // RANDOM NUMBER GENERATION
    double random();                     // Output random float in N(0, 1)
};

// Class to generate vectors with isotropic distribution.
class    GCG_API_CLASS    gcgRANDOMVECTOR   :  public gcgCLASS {
  private:
    gcgRANDOM ngen1;  // Independent random generators
    gcgRANDOM ngen2;

  public:
    // Constructor
    gcgRANDOMVECTOR(); // Use C conventional rand() to define seed
    gcgRANDOMVECTOR(unsigned int seed1, unsigned int seed2); // seed1 must be different than seed2.
    virtual ~gcgRANDOMVECTOR();

    // SERIES SPECIFICATION: seed1 must be different than seed2.
    void  setSeed(unsigned int seed1, unsigned int seed2);

    // RANDOM VECTOR GENERATION
    void random(VECTOR3 v);          // Output isotropic float vector in S^2.
    void random(VECTOR3d v);         // Output isotropic double vector in S^2.
};


class    GCG_API_CLASS    gcgEVENTSPERSECOND   :  public gcgCLASS {
  private:
    void  *handle; 

  public:
    gcgEVENTSPERSECOND();

    virtual ~gcgEVENTSPERSECOND();

    void  startingEvents();

    float finishedEvents(int nevents = 1);


    float getEventsPerSecond();            // Retorna o ultimo EPS
};


// FUNCTION PROTOTYPES

// Compute the Probability Distribution Function from a histogram.
GCG_API_FUNCTION    int gcgPDFfromHistogram(float *pdf, unsigned int *histogram, unsigned int Nbins);


// Find the optimal threshold using Otsu's method.
// Nobuyuki Otsu, "A threshold selection method from gray level histograms",
// IEEE Transactions on Systems, Man and Cybernetics, vol.9, n.1, pp.62-66, jan/1979.
GCG_API_FUNCTION    unsigned int gcgOptimalThresholdOtsu(float *pdf, unsigned int left, unsigned int right);


//************* Double precision versions *************
GCG_API_FUNCTION    int gcgPDFfromHistogram(double *pdf, unsigned int *histogram, unsigned int Nbins);
GCG_API_FUNCTION    unsigned int gcgOptimalThresholdOtsu(double *pdf, unsigned int left, unsigned int right);




#define GCG_FRUSTUM_PARTIAL             -1
#define GCG_FRUSTUM_OUTSIDE        0
#define GCG_FRUSTUM_INSIDE         1


class    GCG_API_CLASS    gcgFRUSTUM   :  public gcgCLASS {
      public: // But read-only
        // CAMERA POSITION AND ORIENTATION: in world coordinates
        float   x, y, z;    // Observer's position

        VECTOR3 view;     // Viewing direction: given by the quaternion
        VECTOR3 up;       // Head up direction: given by the quaternion
        VECTOR3 right;    // Right side direction: given by the quaternion

      private:
        // CAMERA PARAMETERS (always updated by changes in frustum)
        VECTOR4d orientation; // QUATERNION reflecting the viewers orientation.
        MATRIX4d view_matrix;

        // PARAMETROS DA PROJECAO
        MATRIX4d proj_matrix;

        // Equacoes dos planos do FRUSTUM corrente (redundante pois é calculado toda hora)
        // Coeficientes do plano superior do frustum
        float topA, topB, topC, topD;

        // Coeficientes do plano inferior do frustum
        float bottomA, bottomB, bottomC, bottomD;

        // Coeficientes do plano posterior do frustum
        float farA, farB, farC, farD;

        // Coeficientes do plano anterior do frustum
        float nearA, nearB, nearC, nearD;

        // Coeficientes do plano direito do frustum
        float rightA, rightB, rightC, rightD;

        // Coeficientes do plano esquerdo do frustum
        float leftA, leftB, leftC, leftD;

      public:
        //metodos construtores
        gcgFRUSTUM();
        virtual ~gcgFRUSTUM();

        // FREE VIEWER

        // Posicao do observador
        void setPosition(float x, float y, float z); // Posicao do observador
        void setPosition(double x, double y, double z); // Posicao do observador
        void setPosition(VECTOR3 pos); // Posicao do observador
        void setPosition(VECTOR3d pos); // Posicao do observador

        void movePosition(float dx, float dy, float dz);    // Move o observador, dado um vetor
        void movePosition(double dx, double dy, double dz); // Move o observador, dado um vetor
        void movePosition(VECTOR3 pos); // Posicao do observador
        void movePosition(VECTOR3d pos); // Posicao do observador

        void advancePosition(float forward, float upward, float rightward); // Avanca usando as direcoes de visualizacao
        void advancePosition(double forward, double upward, double rightward); // Avanca usando as direcoes de visualizacao

        // Direcoes de observacao
        void resetOrientation();    // Loads the canonical reference system

        void setViewVector(float dirX, float dirY, float dirZ); // Configura vetor de direcao
        void setViewVector(double dirX, double dirY, double dirZ); // Configura vetor de direcao
        void setViewVector(VECTOR3 dir); // Configura vetor de direcao
        void setViewVector(VECTOR3d dir); // Configura vetor de direcao

        void setViewVectorSpherical(float latit, float longit); // latit em [-90, 90]  longit em [-180, 180]
        void setViewVectorSpherical(double latit, double longit); // latit em [-90, 90]  longit em [-180, 180]

        void setUpVector(float cimaX, float cimaY, float cimaZ); // Configura vetor acima
        void setUpVector(double cimaX, double cimaY, double cimaZ); // Configura vetor acima
        void setUpVector(VECTOR3 cima); // Configura vetor acima
        void setUpVector(VECTOR3d cima); // Configura vetor acima

        void setTarget(float alvoX, float alvoY, float alvoZ); // Ajusta as direcoes para visualizar o alvo.
        void setTarget(double alvoX, double alvoY, double alvoZ); // Ajusta as direcoes para visualizar o alvo.
        void setTarget(VECTOR3 alvo); // Ajusta as direcoes para visualizar o alvo.
        void setTarget(VECTOR3d alvo); // Ajusta as direcoes para visualizar o alvo.

        // OBSERVER IN THE CENTER OF AN UNITARY SPHERE (planetarium)

        // Rotates the unit ball where the observer is centered (planetarium)
        //  - all functions use the local viewing system to apply rotations
        void rotateOrientationTrackball(double x1, double y1, double x2, double y2); // Trackball: p1,p2 em [-1,1]x[-1,1] plano no qual a bola é projetada. rotacao livre.
        void rotateOrientationSpherical(double dlat, double dlong);     // Free rotation of the planetarium using spherical angles
        void rotateOrientationEuler(double pitch, double yaw, double row); // Rotation of the planetarium using Euler angles but keeping main axis aligned with Y e X
        void rotateOrientationAxis(double ang, double axisX, double axisY, double axisZ); // Rotation of the planetarium around a given axis.

        void alignAxis(); // Aligns up vector with Y, and right vector with X

        // OBSERVER IS FREE ROTATING AROUND AN ARBITRARY SPHERE (orbital)

        // Rotation of the observer's system around a sphere
        //  - all functions align the local viewing system to the canonical system,
        //    apply rotations and then align the canonical system back to the local system
        void rotateOrbitTrackball(double centroX, double centroY, double centroZ, double x1, double y1, double x2, double y2); // Trackball: p1,p2 em [-1,1]x[-1,1] plano no qual a esfera é projetada. rotacao livre.
        void rotateOrbitSpherical(double centroX, double centroY, double centroZ, double dlat, double dlong); // Spherical coordinates
        void rotateOrbitEuler(double centroX, double centroY, double centroZ, double pitch, double yaw, double roll); // Euler angles

        // OBSERVER IS FREE ROTATING AROUND ARBITRARY AXIS

        // Rotation of the observer's system around an axis
        void rotateAxis(double ang, double ox, double oy, double oz, double dirX, double dirY, double dirZ);


        // VISIBILITY TESTS

        int isPointVisible(float x, float y, float z);
        int isSphereVisible(float x, float y, float z, float raio);
        int isBoxVisible(float xmin, float ymin, float zmin, float xmax, float ymax, float zmax);

        //  GEOMETRIC PROJECTIONS
        void setPerspective(float fovh, float aspect, float prox, float dist);
        void setOrthographic(float  left, float right, float  top, float bottom, float near, float far);

        //  GP Interfaces

        // Loads the projection and modelview matrix to OpenGL
        void exportOpenGL();

        //  CAMERA INFORMATION RETRIEVAL

        void getViewMatrix(MATRIX4 matrot);
        void getViewMatrix(MATRIX4d matrot);

        void getProjectionMatrix(MATRIX4 matproj);
        void getProjectionMatrix(MATRIX4d matproj);

      private:
        // Calcula as direcoes do observador
        void computeParameters();
};




class gcgINTERNAL_VIDEO; // Internal class depends on platform

class    GCG_API_CLASS    gcgVIDEO   :  public gcgCLASS {
  public:
    unsigned int  width;       //< Frame width in pixels. Read-only.
    unsigned int  height;      //< Frame height in pixels. Read-only.
    unsigned char bpp;         //< Bits per pixel. Read-only.
    float         fps;         //< Frames per second. Read-only.
    unsigned int  bitrate;         //< Bits per second. Read-only.

    gcgINTERNAL_VIDEO *handle; //< Handle for internal video object or platform specific API. Read-only.

  public:

    gcgVIDEO();

    virtual ~gcgVIDEO();
};


class    GCG_API_CLASS    gcgVIDEOCAPTURE : public gcgVIDEO {
        unsigned int  id;   //< Current camera device identifier. Read-only.

  public:
                gcgVIDEOCAPTURE();

                virtual ~gcgVIDEOCAPTURE();

    virtual bool destroyVideo();

                virtual int getNumberOfCameras();

                virtual bool getCameraName(unsigned int id, char* pCameraName, unsigned int maxName);

                virtual bool openCamera(unsigned int id, unsigned int width, unsigned int height, unsigned char bpp, float fps);

                virtual bool start();

                virtual bool stop();

                virtual bool resume();

                virtual bool pause();

                virtual bool isCapturing();

                virtual bool setCallBackFunction(void (*callback) (gcgVIDEOCAPTURE*));

                virtual double copyFrameTo(gcgIMAGE *dstimg);
};


class    GCG_API_CLASS    gcgVIDEOFILE : public gcgVIDEO {
        public:
    gcgVIDEOFILE();

    virtual ~gcgVIDEOFILE();

    virtual bool destroyVideo();

    virtual bool openVideoFile(const char *filename);

    virtual long getFrameCount();

    virtual long getCurrentPosition();

    virtual bool setCurrentPosition(long newposition);

    virtual bool copyFrame(gcgIMAGE *dstimg);
};





class    GCG_API_CLASS    gcgDATA  :  public gcgCLASS {
  public:
    virtual ~gcgDATA(){};
};



class    GCG_API_CLASS    gcgITERATOR   :  public gcgCLASS {
  public:
    virtual ~gcgITERATOR();

    virtual gcgDATA *next() = 0;
};


class    GCG_API_CLASS    gcgDATASTRUCTURE   :  public gcgCLASS {
    virtual uintptr_t getCounter() = 0;

    virtual bool deleteAll() = 0;

    virtual gcgITERATOR *getIterator(int traversemode = 0) = 0;

    virtual gcgITERATOR *detach(int traversemode = 0) = 0;
};

class    GCG_API_CLASS    gcgLINK  :  public gcgDATA {
  public:
    union {
      gcgLINK *next;      
      gcgLINK *right;     
    };
};

class    GCG_API_CLASS    gcgDOUBLELINK  :  public gcgLINK {
  public:
    union {
      gcgLINK *previous;  
      gcgLINK *left;      
    };
};


class    GCG_API_CLASS    gcgLISTSTRUCTURE   :  public gcgDATASTRUCTURE {
  public:
    virtual bool insertFirst(gcgLINK *newnode) = 0;

    virtual bool insertLast(gcgLINK *newnode) = 0;

    virtual bool insertAfter(gcgLINK *newnode, gcgLINK *node) = 0;

    virtual bool insertBefore(gcgLINK *newnode, gcgLINK *node) = 0;

    virtual bool moveToFirst(gcgLINK *node) = 0;

    virtual bool moveToLast(gcgLINK *node) = 0;

    virtual bool switchNodes(gcgLINK *node1, gcgLINK *node2) = 0;

    virtual bool remove(gcgLINK *node) = 0;
};






class    GCG_API_CLASS    gcgLINKEDLIST   :  public gcgLISTSTRUCTURE {
  public:
    gcgLINK       *first;   
    gcgLINK       *last;    
    uintptr_t     counter;  

  public:
    gcgLINKEDLIST();

    virtual ~gcgLINKEDLIST();

    // gcgLINKEDLIST specific methods

    gcgLINK *dequeueFirst();


    gcgLINK *dequeueLast();

    // gcgDATASTRUCTURE interface

    uintptr_t getCounter();

    bool deleteAll();

    gcgITERATOR *getIterator(int traversemode = 0);

    gcgITERATOR *detach(int traversemode = 0);

    // gcgLISTSTRUCTURE interface

    bool insertFirst(gcgLINK *node);

    bool insertLast(gcgLINK *node);

    bool insertAfter(gcgLINK *newnode, gcgLINK *node);

    bool insertBefore(gcgLINK *newnode, gcgLINK *node);

    bool moveToFirst(gcgLINK *node);

    bool moveToLast(gcgLINK *node);

    bool switchNodes(gcgLINK *node1, gcgLINK *node2);

    bool remove(gcgLINK *node);
};




class    GCG_API_CLASS    gcgDOUBLELINKEDLIST   :  public gcgLISTSTRUCTURE {
  public:
    gcgLINK   *first;   
    gcgLINK   *last;    
    uintptr_t counter;  

  public:
    gcgDOUBLELINKEDLIST();

    virtual ~gcgDOUBLELINKEDLIST();

    // gcgDOUBLELINKEDLIST specific methods

    gcgLINK *dequeueFirst();


    gcgLINK *dequeueLast();

    // gcgDATASTRUCTURE interface

    uintptr_t getCounter();

    bool deleteAll();

    gcgITERATOR *getIterator(int traversemode = 0);

    gcgITERATOR *detach(int traversemode = 0);


    // gcgLISTSTRUCTURE interface

    bool insertFirst(gcgLINK *node);

    bool insertLast(gcgLINK *node);

    bool insertAfter(gcgLINK *newnode, gcgLINK *node);

    bool insertBefore(gcgLINK *newnode, gcgLINK *node);

    bool moveToFirst(gcgLINK *node);

    bool moveToLast(gcgLINK *node);

    bool switchNodes(gcgLINK *node1, gcgLINK *node2);

    bool remove(gcgLINK *node);
};










class    GCG_API_CLASS    gcgORDEREDNODE  :  public gcgDOUBLELINK {
  public:
    union {
      signed char status;   
            signed char balance;  
    };
};


class    GCG_API_CLASS    gcgORDEREDSTRUCTURE   :  public gcgDATASTRUCTURE {
  public:
    virtual gcgORDEREDNODE *insert(gcgORDEREDNODE *newnode) = 0;

    virtual gcgORDEREDNODE *search(gcgORDEREDNODE *key) = 0;

    virtual gcgORDEREDNODE *remove(gcgORDEREDNODE *key) = 0;

    virtual int compare(gcgORDEREDNODE *refnode1, gcgORDEREDNODE *refnode2) = 0;
};



class    GCG_API_CLASS    gcgAVLTREE   :  public gcgDATASTRUCTURE {
  public:
          gcgORDEREDNODE  *root;    
    uintptr_t       counter;  

  private:
          struct _AVL_LEVEL_PATH  *path;
          int                     levels;
          int                     maxlevels;

  public:
    gcgAVLTREE();

    virtual ~gcgAVLTREE();

    // gcgORDEREDSTRUCTURE interface

    gcgORDEREDNODE *insert(gcgORDEREDNODE *newnode);

    gcgORDEREDNODE *search(gcgORDEREDNODE *key);

    gcgORDEREDNODE *remove(gcgORDEREDNODE *key);

    // gcgDATASTRUCTURE interface

    uintptr_t getCounter();

    virtual bool deleteAll();

    gcgITERATOR *getIterator(int traversemode = 0);

    gcgITERATOR *detach(int traversemode = 0);

    virtual int compare(gcgORDEREDNODE *refnode1, gcgORDEREDNODE *refnode2) = 0;
};


class    GCG_API_CLASS    gcgHASHTABLE   :  public gcgORDEREDSTRUCTURE {
  public:
    uintptr_t     counter;      
    uintptr_t     capacity;     
    float         minimumload;  
    float         maximumload;  
    gcgCLASS      **table;      

  public:
    gcgHASHTABLE(uintptr_t capacity = 0, float minload = 0.8, float maxload = 16.0);

    virtual ~gcgHASHTABLE();

    // gcgHASHTABLE specific methods

    bool setCapacity(uintptr_t capacity);

    bool setLoadLimits(float minload, float maxload);

    // gcgDATASTRUCTURE interface

    uintptr_t getCounter();

    bool deleteAll();

    gcgITERATOR *getIterator(int traversemode = 0);

    gcgITERATOR *detach(int traversemode = 0);

    // gcgORDEREDSTRUCTURE interface

    gcgORDEREDNODE *insert(gcgORDEREDNODE *newnode);

    gcgORDEREDNODE *search(gcgORDEREDNODE *key);

    gcgORDEREDNODE *remove(gcgORDEREDNODE *key);

    virtual int compare(gcgORDEREDNODE *refnode1, gcgORDEREDNODE *refnode2) = 0;

    virtual long hash(gcgORDEREDNODE *node) = 0;
};




typedef union _GCGDIGITALKEY {
  void        *pkey;    
  intptr_t    intkey;   
  long        lkey;     
  int         ikey;     
  short       skey;     
  char        ckey;     

  // Constructors for implicit initialization
  _GCGDIGITALKEY()       : pkey(NULL) {}
  _GCGDIGITALKEY(int64_t l) : intkey((intptr_t) l) {}
  _GCGDIGITALKEY(int l) : intkey((intptr_t) l) {}
  _GCGDIGITALKEY(short l) : intkey((intptr_t) l) {}
  _GCGDIGITALKEY(char l) : intkey((intptr_t) l) {}
  _GCGDIGITALKEY(uint64_t l) : intkey((intptr_t) l) {}
  _GCGDIGITALKEY(unsigned int l)  : intkey((intptr_t) l) {}
  _GCGDIGITALKEY(unsigned short l)  : intkey((intptr_t) l) {}
  _GCGDIGITALKEY(unsigned char l)  : intkey((intptr_t) l) {}
  _GCGDIGITALKEY(void *p) : pkey(p) {}
  _GCGDIGITALKEY(char *p) : pkey((void*) p) {}
  _GCGDIGITALKEY(const char *p) : pkey((void*) p) {}
} GCGDIGITALKEY;



class    GCG_API_CLASS    gcgDATAMAP   :  public gcgDATASTRUCTURE {
  public:

    virtual gcgDATA *insert(GCGDIGITALKEY key, gcgDATA *value) = 0;

    virtual gcgDATA *search(GCGDIGITALKEY key) = 0;

    virtual gcgDATA *remove(GCGDIGITALKEY key) = 0;
};







class    GCG_API_CLASS    gcgPATRICIAMAP   :  public gcgDATAMAP {
  public:
    uintptr_t     counter;    
    unsigned int  keysize;    

  private:
    struct _GCG_PATRICIA_NODE *root;  
    gcgDATA *insertByValue(GCGDIGITALKEY key, gcgDATA *value); // Optimization: faster insertion for keys with size fixed up to sizeof(intptr_t)
    gcgDATA *removeByValue(GCGDIGITALKEY key);                 // Optimization: faster removals for keys with size fixed up to sizeof(intptr_t)

  public:
    gcgPATRICIAMAP(unsigned int keysize);

    virtual ~gcgPATRICIAMAP();

    bool setKeySize(unsigned int keysize);

    gcgITERATOR *suffixIterator(GCGDIGITALKEY prefixkey, unsigned int prefixsize = 0, int traversemode = 0);


    // gcgDATASTRUCTURE interface

    uintptr_t getCounter();

    bool deleteAll();

    gcgITERATOR *getIterator(int traversemode = 0);

    gcgITERATOR *detach(int traversemode = 0);

    // gcgORDEREDSTRUCTURE interface

    virtual gcgDATA *insert(GCGDIGITALKEY key, gcgDATA *value);

    virtual gcgDATA *search(GCGDIGITALKEY key);

    virtual gcgDATA *remove(GCGDIGITALKEY key);
};





GCG_API_FUNCTION    unsigned int gcgGetNumberOfProcessors();

class    GCG_API_CLASS    gcgTHREAD   :  public gcgCLASS {
  public:
    void    *handle;    

  public:
    gcgTHREAD();

    virtual ~gcgTHREAD();

    bool startThread();

    bool abortThread();

    bool waitThread();

    bool isRunning();

  public:
    virtual void run() = 0;
};

class    GCG_API_CLASS    gcgLOCK   :  public gcgCLASS {
  public:
    void *handle;    

    gcgLOCK();

    virtual ~gcgLOCK();

    bool lock(long timeoutUsec = -1);

    bool unlock();

    bool wait(long timeoutUsec = -1);

    bool wakeUp(int numberOfThreads = -1);
};

class    GCG_API_CLASS    gcgQUEUE   :  public gcgCLASS {
  private:
    void *handle; 

  public:
    gcgQUEUE();

    virtual ~gcgQUEUE();

    bool enqueueTail(gcgLINK *node);

    bool enqueueHead(gcgLINK *node);

    gcgLINK *dequeue();

    gcgLINK *waitDequeue(long timeoutUsec = -1);

    bool isEmpty();

    uintptr_t getCounter();

    bool deleteAll();

    gcgLOCK *wakeUpWhenEnqueued(gcgLOCK *notifyme, int wakeupNumThreads);


    gcgLOCK *wakeUpWhenDequeued(gcgLOCK *notifyme, int wakeupNumThreads);
};


class    GCG_API_CLASS    gcgCACHE   :  public gcgCLASS {
  public:
    unsigned long totalcost; 
    unsigned long capacity;  
    unsigned int  keysize;   

  private:
    void *handle;       

  public:
    gcgCACHE(unsigned int keysizebytes, unsigned long cachecapacity);

    virtual ~gcgCACHE();

    // Application interface: retrieving data from the cache

    void *get(void *key);

    bool flush(void *key);

    bool setCapacity(unsigned long maximumtotalcost);

    bool setKeySize(unsigned int keysizebytes);

    // Data interface: managing data creation, removal and persistence.

    virtual void *retrieveData(void *key, unsigned long *datacost) = 0;

    virtual void discardData(void *key, void *data, unsigned long datacost) = 0;
};



class    GCG_API_CLASS    gcgJOB : public gcgORDEREDNODE {
  public:
    gcgJOB() {}

    virtual ~gcgJOB() {}

    virtual void run() = 0;
};


class    GCG_API_CLASS    gcgTHREADPOOL   :  public gcgCLASS {
  public:
    void *handle;   

    gcgTHREADPOOL(unsigned int numThreads = gcgGetNumberOfProcessors());

    virtual ~gcgTHREADPOOL();

    bool setNumberOfThreads(unsigned int numThreads = gcgGetNumberOfProcessors());

    unsigned int getNumberOfThreads();

    bool assignJob(gcgJOB *job);


    uintptr_t getNumberOfPendingJobs();

    bool wait();

    bool waitAndDestroy();

    bool destroyPool();

    bool setOutputQueue(gcgQUEUE *executed, gcgQUEUE *discarded = NULL);

    bool getOutputQueue(gcgQUEUE **pexecuted, gcgQUEUE **pdiscarded);

    bool isRunning();
};



class    GCG_API_CLASS    gcgPRODUCERCONSUMER   :  public gcgCLASS {
    private:
      void *handle; 

    public:
      gcgPRODUCERCONSUMER(unsigned int bufferlength = 4);

      virtual ~gcgPRODUCERCONSUMER();

      bool setBufferSize(unsigned int bufferlength);

      bool destroyProducerConsumer();

      bool put(void *val);

      void *get();

      uintptr_t getCounter();

      void setTimeout(long timeoutprod, long timeoutcons);

      void enable(bool enabled);
};



class    GCG_API_CLASS    gcgTEXTURE2D   :  public gcgCLASS {
  public:
    unsigned int    width;        
    unsigned int    height;       

    unsigned int    actualwidth;  
    unsigned int    actualheight; 

    unsigned int    idOpengl;     

    int             internalformat; 
    int             format; 
    int             type;   

  public:
    gcgTEXTURE2D();

    virtual ~gcgTEXTURE2D();

    bool destroyTexture();

    bool createTexture(unsigned int width, unsigned int height, unsigned char bpp, bool usealpha, bool mipmapped = false);

    bool uploadImage(gcgIMAGE *srcimage, int destposX = 0, int destposY = 0, unsigned int level = 0, bool mipmapped = false);

    bool isCompatibleWith(gcgIMAGE *srcimage);

    bool bind();

    bool unbind();
};






// Event codes for waitEvent() of gcgINETPEER

#define GCG_EVENT_NOTREADYORDISCONNECTED  0x0

#define GCG_EVENT_RECEIVEREADY            0x1

#define GCG_EVENT_SENDREADY               0x2

#define GCG_EVENT_RECEIVESENDREADY        (GCG_EVENT_RECEIVEREADY | GCG_EVENT_SENDREADY) // 3

// Extra event codes for eventHandler() of gcgINETSERVICE

#define GCG_EVENT_INSERTED                0x4

#define GCG_EVENT_REMOVED                 0x8

#define GCG_EVENT_VISIT                   0x10

#define GCG_EVENT_NEWCONNECTION           0x20


class    GCG_API_CLASS    gcgINETPEER   :  public gcgCLASS {
  public:
    int   local_ip[4];  
    int   local_port;   
    bool  isConnected;  

    int   remote_ip[4]; 
    int   remote_port;  

    void  *handle;      

  public:
    gcgINETPEER();

    virtual ~gcgINETPEER(); // Safely disconnects the peer

    bool connectToHost(const char *host, int port);

    bool send(void *buffer, unsigned int nbytes);

    int receive(void *buffer, unsigned int nbytes);

    bool getConnectionFrom(gcgINETPEER *oldpeerobject);

    bool shutdownSend();

    bool shutdownReceive();

    bool disconnect();

    int waitEvent(bool checkReceive = true, bool checkSend = false, int timeoutUsec = -1);
};


class    GCG_API_CLASS    gcgINETSERVICE   :  public gcgCLASS {
  public:
    bool  isListening; 
    void *handle;      

  public:
    gcgINETSERVICE();

    virtual ~gcgINETSERVICE();

    bool listenPort(int port);

    bool stopListening();

    bool insertPeer(gcgINETPEER *peer);

    int getPeerCount();

    bool flagEvents(gcgINETPEER *peer, bool checkReceive, bool checkSend);

    bool postRemove(gcgINETPEER *peer);

    bool postRemoveAll();

    bool postVisit(gcgINETPEER *peer);

    bool postVisitAll();

    bool postExit();

    bool processEvents(unsigned int numThreads = gcgGetNumberOfProcessors());

    bool setTimeout(int timeoutUsec);

    virtual void eventHandler(int events, gcgINETPEER *peer) = 0;
};


// End of definitions included once
#endif // #if _COUNTER_GCG_ == 1



// Next definitions are inserted more than once

// Define the type
#undef NUMTYPE
#undef COERSIONFLOAT
#undef COERSIONDOUBLE
#undef COERSIONSHORT
#undef COERSIONINT
#undef COERSIONLONG
#if   _COUNTER_GCG_ == 1
  #define NUMTYPE float
  #define COERSIONFLOAT float
  #define COERSIONDOUBLE double
  #define COERSIONSHORT float
  #define COERSIONINT float
  #define COERSIONLONG float
#elif _COUNTER_GCG_ == 2
  #define NUMTYPE double
  #define COERSIONFLOAT double
  #define COERSIONDOUBLE double
  #define COERSIONSHORT double
  #define COERSIONINT double
  #define COERSIONLONG double
#elif _COUNTER_GCG_ == 3
  #define NUMTYPE short
  #define COERSIONFLOAT float
  #define COERSIONDOUBLE double
  #define COERSIONSHORT short
  #define COERSIONINT int
  #define COERSIONLONG long
#elif _COUNTER_GCG_ == 4
  #define NUMTYPE int
  #define COERSIONFLOAT float
  #define COERSIONDOUBLE double
  #define COERSIONSHORT int
  #define COERSIONINT int
  #define COERSIONLONG long
#elif _COUNTER_GCG_ == 5
  #define NUMTYPE long
  #define COERSIONFLOAT float
  #define COERSIONDOUBLE double
  #define COERSIONSHORT long
  #define COERSIONINT long
  #define COERSIONLONG long
#endif



template <>
class     GCG_API_CLASS    gcgDISCRETE1D<NUMTYPE>   :  public gcgCLASS {
  public:
    NUMTYPE      *data;           
    unsigned int length;          
    bool         isstaticdata;    

    int          origin;          
    int          extension;       

  public:
    gcgDISCRETE1D();

    gcgDISCRETE1D(unsigned int length, int origin, NUMTYPE *sampledata = NULL, bool isstaticdata = false, int borderextension = GCG_BORDER_EXTENSION_CLAMP);

    virtual ~gcgDISCRETE1D();

    bool createSignal(unsigned int length, int origin, NUMTYPE *sampledata = NULL, bool isstaticdata = false, int borderextension = GCG_BORDER_EXTENSION_CLAMP);

    bool destroySignal();

    bool createSimilar(gcgDISCRETE1D<float> *source);

    bool createSimilar(gcgDISCRETE1D<double> *source);

    bool createSimilar(gcgDISCRETE1D<short> *source);

    bool createSimilar(gcgDISCRETE1D<int> *source);

    bool createSimilar(gcgDISCRETE1D<long> *source);

    bool isCompatibleWith(gcgDISCRETE1D<float> *source); // Checks if both signals are compatible: same dimensions, same origin.
    bool isCompatibleWith(gcgDISCRETE1D<double> *source);
    bool isCompatibleWith(gcgDISCRETE1D<short> *source);
    bool isCompatibleWith(gcgDISCRETE1D<int> *source);
    bool isCompatibleWith(gcgDISCRETE1D<long> *source);
    bool duplicateSignal(gcgDISCRETE1D<float> *source); // Makes an independent copy of the source signal.
    bool duplicateSignal(gcgDISCRETE1D<double> *source);
    bool duplicateSignal(gcgDISCRETE1D<short> *source);
    bool duplicateSignal(gcgDISCRETE1D<int> *source);
    bool duplicateSignal(gcgDISCRETE1D<long> *source);
    bool duplicateSubsignal(gcgDISCRETE1D<float> *source, int left, unsigned int width); // Makes an independent copy of the source subsignal.
    bool duplicateSubsignal(gcgDISCRETE1D<double> *source, int left, unsigned int width);
    bool duplicateSubsignal(gcgDISCRETE1D<short> *source, int left, unsigned int width);
    bool duplicateSubsignal(gcgDISCRETE1D<int> *source, int left, unsigned int width);
    bool duplicateSubsignal(gcgDISCRETE1D<long> *source, int left, unsigned int width);

    // Data access and modification. Not recommended for high
    // performance processes. Uses the origin and border extension.
    NUMTYPE getDataSample(int i);
    bool setDataSample(int i, NUMTYPE value);

    // Signal processing functions
    // The src can be the destination (this).
    bool normalize(gcgDISCRETE1D<float> *src, float floor = 0, float ceil = 1);
    bool normalize(gcgDISCRETE1D<double> *src, double floor = 0, double ceil = 1);
    bool normalize(gcgDISCRETE1D<short> *src, short floor = 0, short ceil = 1);
    bool normalize(gcgDISCRETE1D<int> *src, int floor = 0, int ceil = 1);
    bool normalize(gcgDISCRETE1D<long> *src, long floor = 0, long ceil = 1);
    bool convolution(gcgDISCRETE1D<NUMTYPE> *src, gcgDISCRETE1D<float> *mask);
    bool convolution(gcgDISCRETE1D<NUMTYPE> *src, gcgDISCRETE1D<double> *mask);
    bool convolution(gcgDISCRETE1D<NUMTYPE> *src, gcgDISCRETE1D<short> *mask);
    bool convolution(gcgDISCRETE1D<NUMTYPE> *src, gcgDISCRETE1D<int> *mask);
    bool convolution(gcgDISCRETE1D<NUMTYPE> *src, gcgDISCRETE1D<long> *mask);
    bool templateMatching(unsigned int imgleft, unsigned int imgwidth, gcgDISCRETE1D<float> *mask, int *position);
    bool templateMatching(unsigned int imgleft, unsigned int imgwidth, gcgDISCRETE1D<double> *mask, int *position);
    bool templateMatching(unsigned int imgleft, unsigned int imgwidth, gcgDISCRETE1D<short> *mask, int *position);
    bool templateMatching(unsigned int imgleft, unsigned int imgwidth, gcgDISCRETE1D<int> *mask, int *position);
    bool templateMatching(unsigned int imgleft, unsigned int imgwidth, gcgDISCRETE1D<long> *mask, int *position);

    // Computes the scalar product of two signals. At least one of the signals must have extension
    // GCG_BORDER_EXTENSION_ZERO. Aligns the origin of srcsignal2 with the sample of this object
    // with coordinate atX.
    COERSIONFLOAT  scalarProduct(int atX, gcgDISCRETE1D<float> *srcsignal2);
    COERSIONDOUBLE scalarProduct(int atX, gcgDISCRETE1D<double> *srcsignal2);
    COERSIONSHORT  scalarProduct(int atX, gcgDISCRETE1D<short> *srcsignal2);
    COERSIONINT    scalarProduct(int atX, gcgDISCRETE1D<int> *srcsignal2);
    COERSIONLONG   scalarProduct(int atX, gcgDISCRETE1D<long> *srcsignal2);

    // Produces a binary signal. The samples are mapped to: lessvalue, if sample is below the threshold; greaterequalvalue,
    // if the sample is greater or equal the threshold. The src can be the destination (this).
    bool binarize(gcgDISCRETE1D<float> *src, float threshold, NUMTYPE lessvalue = 0, NUMTYPE greaterequalvalue = 1);
    bool binarize(gcgDISCRETE1D<double> *src, double threshold, NUMTYPE lessvalue = 0, NUMTYPE greaterequalvalue = 1);
    bool binarize(gcgDISCRETE1D<short> *src, short threshold, NUMTYPE lessvalue = 0, NUMTYPE greaterequalvalue = 1);
    bool binarize(gcgDISCRETE1D<int> *src, int threshold, NUMTYPE lessvalue = 0, NUMTYPE greaterequalvalue = 1);
    bool binarize(gcgDISCRETE1D<long> *src, long threshold, NUMTYPE lessvalue = 0, NUMTYPE greaterequalvalue = 1);

    // Signal composition. Source signals must be compatible.
    // The destination can be one of the sources.
    bool scale(gcgDISCRETE1D<float> *srcsignal1, COERSIONFLOAT weight, COERSIONFLOAT addthis = 0); // Scales a signal: I = w * I1 + at
    bool scale(gcgDISCRETE1D<double> *srcsignal1, COERSIONDOUBLE weight, COERSIONDOUBLE addthis = 0);
    bool scale(gcgDISCRETE1D<short> *srcsignal1, COERSIONSHORT weight, COERSIONSHORT addthis = 0);
    bool scale(gcgDISCRETE1D<int> *srcsignal1, COERSIONINT weight, COERSIONINT addthis = 0);
    bool scale(gcgDISCRETE1D<long> *srcsignal1, COERSIONLONG weight, COERSIONLONG addthis = 0);
    bool power(gcgDISCRETE1D<float> *srcsignal1, COERSIONFLOAT power); // Powers a signal: I = I1^power
    bool power(gcgDISCRETE1D<double> *srcsignal1, COERSIONDOUBLE power);
    bool power(gcgDISCRETE1D<short> *srcsignal1, COERSIONSHORT power);
    bool power(gcgDISCRETE1D<int> *srcsignal1, COERSIONINT power);
    bool power(gcgDISCRETE1D<long> *srcsignal1, COERSIONLONG power);
    bool combineAdd(gcgDISCRETE1D<NUMTYPE> *srcsignal1, gcgDISCRETE1D<float> *srcsignal2, COERSIONFLOAT weight1 = 1, COERSIONFLOAT weight2 = 1);  // Combine two signals by addition: I = w1 * I1 + w2 * I2
    bool combineAdd(gcgDISCRETE1D<NUMTYPE> *srcsignal1, gcgDISCRETE1D<double> *srcsignal2, COERSIONDOUBLE weight1 = 1, COERSIONDOUBLE weight2 = 1);
    bool combineAdd(gcgDISCRETE1D<NUMTYPE> *srcsignal1, gcgDISCRETE1D<short> *srcsignal2, COERSIONSHORT weight1 = 1, COERSIONSHORT weight2 = 1);
    bool combineAdd(gcgDISCRETE1D<NUMTYPE> *srcsignal1, gcgDISCRETE1D<int> *srcsignal2, COERSIONINT weight1 = 1, COERSIONINT weight2 = 1);
    bool combineAdd(gcgDISCRETE1D<NUMTYPE> *srcsignal1, gcgDISCRETE1D<long> *srcsignal2, COERSIONLONG weight1 = 1, COERSIONLONG weight2 = 1);
    bool combineMult(gcgDISCRETE1D<NUMTYPE> *srcsignal1, gcgDISCRETE1D<float> *srcsignal2, COERSIONFLOAT add1 = 0, COERSIONFLOAT add2 = 0);  // Combine two signals by multiplication: I = (a1 + I1) * (a2 + I2)
    bool combineMult(gcgDISCRETE1D<NUMTYPE> *srcsignal1, gcgDISCRETE1D<double> *srcsignal2, COERSIONDOUBLE add1 = 0, COERSIONDOUBLE add2 = 0);
    bool combineMult(gcgDISCRETE1D<NUMTYPE> *srcsignal1, gcgDISCRETE1D<short> *srcsignal2, COERSIONSHORT add1 = 0, COERSIONSHORT add2 = 0);
    bool combineMult(gcgDISCRETE1D<NUMTYPE> *srcsignal1, gcgDISCRETE1D<int> *srcsignal2, COERSIONINT add1 = 0, COERSIONINT add2 = 0);
    bool combineMult(gcgDISCRETE1D<NUMTYPE> *srcsignal1, gcgDISCRETE1D<long> *srcsignal2, COERSIONLONG add1 = 0, COERSIONLONG add2 = 0);
    bool combineNorm(gcgDISCRETE1D<NUMTYPE> *srcsignal1, gcgDISCRETE1D<float> *srcsignal2); // Combine two signals by the L2 norm: I = sqrt(I1^2 + I2^2)
    bool combineNorm(gcgDISCRETE1D<NUMTYPE> *srcsignal1, gcgDISCRETE1D<double> *srcsignal2);
    bool combineNorm(gcgDISCRETE1D<NUMTYPE> *srcsignal1, gcgDISCRETE1D<short> *srcsignal2);
    bool combineNorm(gcgDISCRETE1D<NUMTYPE> *srcsignal1, gcgDISCRETE1D<int> *srcsignal2);
    bool combineNorm(gcgDISCRETE1D<NUMTYPE> *srcsignal1, gcgDISCRETE1D<long> *srcsignal2);

    // Interfaces for gcgIMAGE
    bool exportGrayScaleImage(gcgIMAGE *outimage, bool normalize);
};



template <>
class     GCG_API_CLASS    gcgDISCRETE2D<NUMTYPE>   :  public gcgCLASS {
  public:
    // Signal info and data
    NUMTYPE       *data;        
    unsigned int  width;        
    unsigned int  height;       
    bool          isstaticdata; 

    // Origin coordinate: indicates the position of the origin (x = 0, y = 0) of the signal.
    int    originX;     
    int    originY;     

    int    extensionX;  
    int    extensionY;  

    int    errorcode;   

  public:
    // Discrete signal construction
    gcgDISCRETE2D();
    gcgDISCRETE2D(unsigned int width, unsigned int length, int originX, int originY, NUMTYPE *sampledata = NULL,
                  bool isstaticdata = false, int borderextX = GCG_BORDER_EXTENSION_CLAMP, int borderextY = GCG_BORDER_EXTENSION_CLAMP);
    virtual ~gcgDISCRETE2D();

    // Empty discrete signal creation and destruction
    bool createSignal(unsigned int width, unsigned int height, int originX, int originY, NUMTYPE *sampledata = NULL,
                      bool isstaticdata = false, int borderextX = GCG_BORDER_EXTENSION_CLAMP, int borderextY = GCG_BORDER_EXTENSION_CLAMP);
    void destroySignal();  // Release all resources used by this discrete signal.

    // Discrete signal copying, duplication and conversion.
    bool createSimilar(gcgDISCRETE2D<float> *source); // Creates a signal similar as the source. The data is not copied.
    bool createSimilar(gcgDISCRETE2D<double> *source);
    bool createSimilar(gcgDISCRETE2D<short> *source);
    bool createSimilar(gcgDISCRETE2D<int> *source);
    bool createSimilar(gcgDISCRETE2D<long> *source);
    bool isCompatibleWith(gcgDISCRETE2D<float> *source); // Checks if both signals are compatible: same dimensions, same origin.
    bool isCompatibleWith(gcgDISCRETE2D<double> *source);
    bool isCompatibleWith(gcgDISCRETE2D<short> *source);
    bool isCompatibleWith(gcgDISCRETE2D<int> *source);
    bool isCompatibleWith(gcgDISCRETE2D<long> *source);
    bool duplicateSignal(gcgDISCRETE2D<float> *source); // Makes an independent copy of the source signal.
    bool duplicateSignal(gcgDISCRETE2D<double> *source);
    bool duplicateSignal(gcgDISCRETE2D<short> *source);
    bool duplicateSignal(gcgDISCRETE2D<int> *source);
    bool duplicateSignal(gcgDISCRETE2D<long> *source);
    bool duplicateSubsignal(gcgDISCRETE2D<float> *source, int left, int top, unsigned int width, unsigned int height); // Makes an independent copy of the source subsignal.
    bool duplicateSubsignal(gcgDISCRETE2D<double> *source, int left, int top, unsigned int width, unsigned int height);
    bool duplicateSubsignal(gcgDISCRETE2D<short> *source, int left, int top, unsigned int width, unsigned int height);
    bool duplicateSubsignal(gcgDISCRETE2D<int> *source, int left, int top, unsigned int width, unsigned int height);
    bool duplicateSubsignal(gcgDISCRETE2D<long> *source, int left, int top, unsigned int width, unsigned int height);

    // Data access and modification. Not recommended for high
    // performance processes. Uses the origin and border extension.
    NUMTYPE getDataSample(int i, int j);
    bool setDataSample(int i, int j, NUMTYPE value);

    // Signal processing functions
    // The src can be the destination (this).
    bool normalize(gcgDISCRETE2D<float> *src, float floor = 0, float ceil = 1);
    bool normalize(gcgDISCRETE2D<double> *src, double floor = 0, double ceil = 1);
    bool normalize(gcgDISCRETE2D<short> *src, short floor = 0, short ceil = 1);
    bool normalize(gcgDISCRETE2D<int> *src, int floor = 0, int ceil = 1);
    bool normalize(gcgDISCRETE2D<long> *src, long floor = 0, long ceil = 1);
    bool convolutionX(gcgDISCRETE2D<NUMTYPE> *src, gcgDISCRETE1D<float> *mask);
    bool convolutionX(gcgDISCRETE2D<NUMTYPE> *src, gcgDISCRETE1D<double> *mask);
    bool convolutionX(gcgDISCRETE2D<NUMTYPE> *src, gcgDISCRETE1D<short> *mask);
    bool convolutionX(gcgDISCRETE2D<NUMTYPE> *src, gcgDISCRETE1D<int> *mask);
    bool convolutionX(gcgDISCRETE2D<NUMTYPE> *src, gcgDISCRETE1D<long> *mask);
    bool convolutionY(gcgDISCRETE2D<NUMTYPE> *src, gcgDISCRETE1D<float> *mask);
    bool convolutionY(gcgDISCRETE2D<NUMTYPE> *src, gcgDISCRETE1D<double> *mask);
    bool convolutionY(gcgDISCRETE2D<NUMTYPE> *src, gcgDISCRETE1D<short> *mask);
    bool convolutionY(gcgDISCRETE2D<NUMTYPE> *src, gcgDISCRETE1D<int> *mask);
    bool convolutionY(gcgDISCRETE2D<NUMTYPE> *src, gcgDISCRETE1D<long> *mask);
    bool convolutionXY(gcgDISCRETE2D<NUMTYPE> *src, gcgDISCRETE2D<float> *mask);
    bool convolutionXY(gcgDISCRETE2D<NUMTYPE> *src, gcgDISCRETE2D<double> *mask);
    bool convolutionXY(gcgDISCRETE2D<NUMTYPE> *src, gcgDISCRETE2D<short> *mask);
    bool convolutionXY(gcgDISCRETE2D<NUMTYPE> *src, gcgDISCRETE2D<int> *mask);
    bool convolutionXY(gcgDISCRETE2D<NUMTYPE> *src, gcgDISCRETE2D<long> *mask);
    bool templateMatching(unsigned int imgleft, unsigned int imgbottom, unsigned int imgwidth, unsigned int imgheight, gcgDISCRETE2D<float> *mask, int *positionX, int *positionY);
    bool templateMatching(unsigned int imgleft, unsigned int imgbottom, unsigned int imgwidth, unsigned int imgheight, gcgDISCRETE2D<double> *mask, int *positionX, int *positionY);
    bool templateMatching(unsigned int imgleft, unsigned int imgbottom, unsigned int imgwidth, unsigned int imgheight, gcgDISCRETE2D<short> *mask, int *positionX, int *positionY);
    bool templateMatching(unsigned int imgleft, unsigned int imgbottom, unsigned int imgwidth, unsigned int imgheight, gcgDISCRETE2D<int> *mask, int *positionX, int *positionY);
    bool templateMatching(unsigned int imgleft, unsigned int imgbottom, unsigned int imgwidth, unsigned int imgheight, gcgDISCRETE2D<long> *mask, int *positionX, int *positionY);

    // Computes the scalar product of two signals. At least one of the signals must have extension
    // GCG_BORDER_EXTENSION_ZERO for X, and at least one of the signals must have extension
    // GCG_BORDER_EXTENSION_ZERO for Y. Aligns the origin of srcsignal2 with the sample of this object
    // with relative coordinates (atX, atY).
    COERSIONFLOAT scalarProduct(int atX, int atY, gcgDISCRETE2D<float> *srcsignal2);
    COERSIONDOUBLE scalarProduct(int atX, int atY, gcgDISCRETE2D<double> *srcsignal2);
    COERSIONSHORT scalarProduct(int atX, int atY, gcgDISCRETE2D<short> *srcsignal2);
    COERSIONINT scalarProduct(int atX, int atY, gcgDISCRETE2D<int> *srcsignal2);
    COERSIONLONG scalarProduct(int atX, int atY, gcgDISCRETE2D<long> *srcsignal2);

    // Produces a binary signal. The samples are mapped to: lessvalue, if sample is below the threshold; greaterequalvalue,
    // if the sample is greater or equal the threshold. The src can be the destination (this).
    bool binarize(gcgDISCRETE2D<float> *src, float threshold, NUMTYPE lessvalue = 0, NUMTYPE greaterequalvalue = 1);
    bool binarize(gcgDISCRETE2D<double> *src, double threshold, NUMTYPE lessvalue = 0, NUMTYPE greaterequalvalue = 1);
    bool binarize(gcgDISCRETE2D<short> *src, short threshold, NUMTYPE lessvalue = 0, NUMTYPE greaterequalvalue = 1);
    bool binarize(gcgDISCRETE2D<int> *src, int threshold, NUMTYPE lessvalue = 0, NUMTYPE greaterequalvalue = 1);
    bool binarize(gcgDISCRETE2D<long> *src, long threshold, NUMTYPE lessvalue = 0, NUMTYPE greaterequalvalue = 1);

    // Signal composition. Source signals must be compatible.
    // The destination can be one of the sources.
    bool scale(gcgDISCRETE2D<float> *srcsignal1, COERSIONFLOAT weight, COERSIONFLOAT addthis = 0); // Scales a signal: I = w * I1 + at
    bool scale(gcgDISCRETE2D<double> *srcsignal1, COERSIONDOUBLE weight, COERSIONDOUBLE addthis = 0);
    bool scale(gcgDISCRETE2D<short> *srcsignal1, COERSIONSHORT weight, COERSIONSHORT addthis = 0);
    bool scale(gcgDISCRETE2D<int> *srcsignal1, COERSIONINT weight, COERSIONINT addthis = 0);
    bool scale(gcgDISCRETE2D<long> *srcsignal1, COERSIONLONG weight, COERSIONLONG addthis = 0);
    bool power(gcgDISCRETE2D<float> *srcsignal1, COERSIONFLOAT power); // Powers a signal: I = I1^power
    bool power(gcgDISCRETE2D<double> *srcsignal1, COERSIONDOUBLE power);
    bool power(gcgDISCRETE2D<short> *srcsignal1, COERSIONSHORT power);
    bool power(gcgDISCRETE2D<int> *srcsignal1, COERSIONINT power);
    bool power(gcgDISCRETE2D<long> *srcsignal1, COERSIONLONG power);
    bool combineAdd(gcgDISCRETE2D<NUMTYPE> *srcsignal1, gcgDISCRETE2D<float> *srcsignal2, COERSIONFLOAT weight1 = 1, COERSIONFLOAT weight2 = 1);  // Combine two signals by addition: I = w1 * I1 + w2 * I2
    bool combineAdd(gcgDISCRETE2D<NUMTYPE> *srcsignal1, gcgDISCRETE2D<double> *srcsignal2, COERSIONDOUBLE weight1 = 1, COERSIONDOUBLE weight2 = 1);
    bool combineAdd(gcgDISCRETE2D<NUMTYPE> *srcsignal1, gcgDISCRETE2D<short> *srcsignal2, COERSIONSHORT weight1 = 1, COERSIONSHORT weight2 = 1);
    bool combineAdd(gcgDISCRETE2D<NUMTYPE> *srcsignal1, gcgDISCRETE2D<int> *srcsignal2, COERSIONINT weight1 = 1, COERSIONINT weight2 = 1);
    bool combineAdd(gcgDISCRETE2D<NUMTYPE> *srcsignal1, gcgDISCRETE2D<long> *srcsignal2, COERSIONLONG weight1 = 1, COERSIONLONG weight2 = 1);
    bool combineMult(gcgDISCRETE2D<NUMTYPE> *srcsignal1, gcgDISCRETE2D<float> *srcsignal2, COERSIONFLOAT add1 = 0, COERSIONFLOAT add2 = 0);  // Combine two signals by multiplication: I = (a1 + I1) * (a2 + I2)
    bool combineMult(gcgDISCRETE2D<NUMTYPE> *srcsignal1, gcgDISCRETE2D<double> *srcsignal2, COERSIONDOUBLE add1 = 0, COERSIONDOUBLE add2 = 0);
    bool combineMult(gcgDISCRETE2D<NUMTYPE> *srcsignal1, gcgDISCRETE2D<short> *srcsignal2, COERSIONSHORT add1 = 0, COERSIONSHORT add2 = 0);
    bool combineMult(gcgDISCRETE2D<NUMTYPE> *srcsignal1, gcgDISCRETE2D<int> *srcsignal2, COERSIONINT add1 = 0, COERSIONINT add2 = 0);
    bool combineMult(gcgDISCRETE2D<NUMTYPE> *srcsignal1, gcgDISCRETE2D<long> *srcsignal2, COERSIONLONG add1 = 0, COERSIONLONG add2 = 0);
    bool combineNorm(gcgDISCRETE2D<NUMTYPE> *srcsignal1, gcgDISCRETE2D<float> *srcsignal2); // Combine two signals by the L2 norm: I = sqrt(I1^2 + I2^2)
    bool combineNorm(gcgDISCRETE2D<NUMTYPE> *srcsignal1, gcgDISCRETE2D<double> *srcsignal2);
    bool combineNorm(gcgDISCRETE2D<NUMTYPE> *srcsignal1, gcgDISCRETE2D<short> *srcsignal2);
    bool combineNorm(gcgDISCRETE2D<NUMTYPE> *srcsignal1, gcgDISCRETE2D<int> *srcsignal2);
    bool combineNorm(gcgDISCRETE2D<NUMTYPE> *srcsignal1, gcgDISCRETE2D<long> *srcsignal2);

    // Domain processing (2D support)
    bool verticalFlip();  // Mirrors the signal in the vertical direction.

    // Interfaces for gcgIMAGE
    bool exportGrayScaleImage(gcgIMAGE *outimage, bool normalize);
};



// Next section are included only for floating point types
#if _COUNTER_GCG_ < 3



template <>
class     GCG_API_CLASS    gcgMATRIX<NUMTYPE>   :  public gcgCLASS {
  public:
    // Matrix info and data
    NUMTYPE       *data;        
    unsigned int  rows;         
    unsigned int  columns;      
    bool          isstaticdata; 

  public:
    // Matrix construction
    gcgMATRIX();
    gcgMATRIX(unsigned int rows, unsigned int columns, NUMTYPE *sampledata = NULL, bool isstaticdata = false);
    virtual ~gcgMATRIX();

    // Empty matrix creation and destruction
    bool createMatrix(unsigned int rows, unsigned int columns, NUMTYPE *sampledata = NULL, bool isstaticdata = false);
    void destroyMatrix();  // Release all resources used by this matrix

    // Matrix copying, duplication and conversion.
    bool createSimilar(gcgMATRIX<float> *source); // Creates a matrix similar as the source. The data is not copied.
    bool createSimilar(gcgMATRIX<double> *source);
    bool isCompatibleWith(gcgMATRIX<float> *source); // Checks if both matrices are compatible: same dimensions
    bool isCompatibleWith(gcgMATRIX<double> *source);
    bool duplicateMatrix(gcgMATRIX<float> *source); // Makes an independent copy of the source matrix.
    bool duplicateMatrix(gcgMATRIX<double> *source);
    bool duplicateSubMatrix(gcgMATRIX<float> *source, int firstrow, int firstcolumn, unsigned int rows, unsigned int columns); // Makes an independent copy of the source submatrix.
    bool duplicateSubMatrix(gcgMATRIX<double> *source, int firstrow, int firstcolumn, unsigned int rows, unsigned int columns);

    // Matrix access and modification. Not recommended for high performance operations.
    NUMTYPE getDataSample(int row, int column);
    bool setDataSample(int row, int column, NUMTYPE value);

    // Matrix composition. Source matrices must be compatible.
    // The destination can be one of the sources.
    bool scale(gcgMATRIX<float> *srcmatrix1, COERSIONFLOAT weight, COERSIONFLOAT addthis = 0); // Scales a matrix: M = w * M1 + at
    bool scale(gcgMATRIX<double> *srcmatrix1, COERSIONDOUBLE weight, COERSIONDOUBLE addthis = 0);
    bool power(gcgMATRIX<float> *srcmatrix1, COERSIONFLOAT power); // Powers the elements of the matrix: M = M1^power
    bool power(gcgMATRIX<double> *srcmatrix1, COERSIONDOUBLE power);
    bool combineAdd(gcgMATRIX<NUMTYPE> *srcmatrix1,  gcgMATRIX<float> *srcmatrix2, COERSIONFLOAT weight1 = 1, COERSIONFLOAT weight2 = 1);  // Combine two signals by addition: M = w1 * M1 + w2 * M2
    bool combineAdd(gcgMATRIX<NUMTYPE> *srcmatrix1,  gcgMATRIX<double> *srcmatrix2, COERSIONDOUBLE weight1 = 1, COERSIONDOUBLE weight2 = 1);
    bool combineMult(gcgMATRIX<NUMTYPE> *srcmatrix1, gcgMATRIX<float> *srcmatrix2, COERSIONFLOAT add1 = 0, COERSIONFLOAT add2 = 0);  // Combine two signals by multiplication: M = (a1 + M1) * (a2 + M2)
    bool combineMult(gcgMATRIX<NUMTYPE> *srcmatrix1, gcgMATRIX<double> *srcmatrix2, COERSIONDOUBLE add1 = 0, COERSIONDOUBLE add2 = 0);
    bool combineNorm(gcgMATRIX<NUMTYPE> *srcmatrix1, gcgMATRIX<float> *srcmatrix2); // Combine two signals by the L2 norm: M = sqrt(M1^2 + M2^2)
    bool combineNorm(gcgMATRIX<NUMTYPE> *srcmatrix1, gcgMATRIX<double> *srcmatrix2);
};





template <>
class     GCG_API_CLASS    gcgLEGENDREBASIS1D<NUMTYPE> : public gcgBASIS1D<NUMTYPE> {
  public:
    unsigned int degree;    
    unsigned int nsamples;  

  private:
    double      *coefficients;  // Coefficients of polynomials. Double precision reduces numerical problems.
    unsigned int npolynomials;  // number of polynomials = maximum degree of the polynomial basis + 1

    NUMTYPE *discretepolynomials; // Discrete version of polynomials

  public:
    gcgLEGENDREBASIS1D();
    virtual ~gcgLEGENDREBASIS1D();

    bool setBasisDegree(unsigned int degree);

    bool setNumberOfSamples(unsigned int nsamples);

    unsigned int getNumberOfCoefficients();

    NUMTYPE getPointValue(unsigned int idegree, NUMTYPE x);

    NUMTYPE getIntegralValue(unsigned int idegree, NUMTYPE a, NUMTYPE b);

    NUMTYPE getDiscreteValue(unsigned int idegree, unsigned int i);

    bool projectSignal(int atX, gcgDISCRETE1D<float> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef);
    bool projectSignal(int atX, gcgDISCRETE1D<double> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef);
    bool projectSignal(int atX, gcgDISCRETE1D<short> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef);
    bool projectSignal(int atX, gcgDISCRETE1D<int> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef);
    bool projectSignal(int atX, gcgDISCRETE1D<long> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef);

    bool reconstructSignal(int atX, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE1D<float> *outputvector);
    bool reconstructSignal(int atX, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE1D<double> *outputvector);
    bool reconstructSignal(int atX, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE1D<short> *outputvector);
    bool reconstructSignal(int atX, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE1D<int> *outputvector);
    bool reconstructSignal(int atX, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE1D<long> *outputvector);

    bool basisInformation(const char *filename);
};





template <>
class     GCG_API_CLASS    gcgLEGENDREBASIS2D<NUMTYPE> : public gcgBASIS2D<NUMTYPE> {
  public:
    unsigned int degree;       
    unsigned int nsamplesX;    
    unsigned int nsamplesY;    

  private:
    double       *coefficients; // Coefficients of 2D polynomials. Double precision reduces numerical problems.
    unsigned int npolynomials;  // number of polynomials for x and y = ((degree-2)*(degree-1))/2

    NUMTYPE *discretepolynomials; // Discrete version of polynomials

  public:
    gcgLEGENDREBASIS2D();
    virtual ~gcgLEGENDREBASIS2D();

    bool setBasisDegree(unsigned int degree);

    bool setNumberOfSamples(unsigned int nsamplesX, unsigned int nsamplesY);

    unsigned int getNumberOfCoefficients();

    bool synchBasisWithFile(const char *filename, unsigned int nsamplesX, unsigned int nsamplesY);

    NUMTYPE getPointValue(unsigned int i, unsigned int j, NUMTYPE x, NUMTYPE y);

    NUMTYPE getIntegralValue(unsigned int i, unsigned int j, NUMTYPE xmin, NUMTYPE xmax, NUMTYPE ymin, NUMTYPE ymax);

    NUMTYPE getDiscreteValue(unsigned int i, unsigned int j, unsigned int ix, unsigned int iy);

    bool projectSignal(int atX, int atY, gcgDISCRETE2D<float> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef);
    bool projectSignal(int atX, int atY, gcgDISCRETE2D<double> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef);
    bool projectSignal(int atX, int atY, gcgDISCRETE2D<short> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef);
    bool projectSignal(int atX, int atY, gcgDISCRETE2D<int> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef);
    bool projectSignal(int atX, int atY, gcgDISCRETE2D<long> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef);


    bool reconstructSignal(int atX, int atY, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE2D<float> *outputvector);
    bool reconstructSignal(int atX, int atY, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE2D<double> *outputvector);
    bool reconstructSignal(int atX, int atY, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE2D<short> *outputvector);
    bool reconstructSignal(int atX, int atY, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE2D<int> *outputvector);
    bool reconstructSignal(int atX, int atY, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE2D<long> *outputvector);

    bool basisInformation(const char *filename, bool checkorthogonality, bool savefunctionsasimages);
};




template <>
class     GCG_API_CLASS    gcgDWTBASIS1D<NUMTYPE> : public gcgBASIS1D<NUMTYPE> {
  public:
    gcgDISCRETE1D<NUMTYPE> H;       
    gcgDISCRETE1D<NUMTYPE> G;       
    gcgDISCRETE1D<NUMTYPE> rH;      
    gcgDISCRETE1D<NUMTYPE> rG;      

    unsigned int nsamples;          
    unsigned int ncoefficients;     

  private:
    void        *buffer;            

  public:
    gcgDWTBASIS1D();
    virtual ~gcgDWTBASIS1D();

    bool setWavelets(gcgDISCRETE1D<NUMTYPE> *H, gcgDISCRETE1D<NUMTYPE> *G, gcgDISCRETE1D<NUMTYPE> *rH, gcgDISCRETE1D<NUMTYPE> *rG);

    bool setNumberOfSamples(unsigned int nsamples);

    unsigned int getNumberOfCoefficients();

    bool projectSignal(int atX, gcgDISCRETE1D<float> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef);
    bool projectSignal(int atX, gcgDISCRETE1D<double> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef);
    bool projectSignal(int atX, gcgDISCRETE1D<short> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef);
    bool projectSignal(int atX, gcgDISCRETE1D<int> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef);
    bool projectSignal(int atX, gcgDISCRETE1D<long> *vector, gcgDISCRETE1D<NUMTYPE> *outputcoef);

    bool reconstructSignal(int atX, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE1D<float> *outputvector);
    bool reconstructSignal(int atX, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE1D<double> *outputvector);
    bool reconstructSignal(int atX, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE1D<short> *outputvector);
    bool reconstructSignal(int atX, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE1D<int> *outputvector);
    bool reconstructSignal(int atX, gcgDISCRETE1D<NUMTYPE> *inputcoef, gcgDISCRETE1D<long> *outputvector);

    bool basisInformation(const char *filename);
};


#endif



// Class for computing the partial derivative of a 2D signal Icurrent (given
// the sequence of signals {Iprevious, Icurrent, Inext} consecutively defined
// in a axis t, perpendicular to the signals) in respect to axis t. The axis
// t is frequently defined as time, over which the 2D signal vary.
//
// In each call of sequenceDerivative(), the parameters are compared with
// those from previous invocation, and already computed information is reused.
//

template <>
class     GCG_API_CLASS    gcgSEQUENCEDERIVATIVE2D<NUMTYPE>   :  public gcgCLASS {
  private:
    gcgDISCRETE2D<NUMTYPE> *current;
    gcgDISCRETE2D<NUMTYPE> *next;
    gcgDISCRETE2D<NUMTYPE> *lowcurrent;
    gcgDISCRETE2D<NUMTYPE> *lownext;
    gcgDISCRETE2D<NUMTYPE> temp;

  public:
    // Constructors and destructors
    gcgSEQUENCEDERIVATIVE2D<NUMTYPE>();
    virtual ~gcgSEQUENCEDERIVATIVE2D<NUMTYPE>();

    // Computes a new partial derivative.
    virtual bool sequenceDerivative(gcgDISCRETE2D<NUMTYPE> *current, gcgDISCRETE2D<NUMTYPE> *next, gcgDISCRETE2D<NUMTYPE> *outputDt);
};



// Class for generating filter masks.
//
// Use gcgFILTERMASK<float> or gcgFILTERMASK<> for single precision.
// Use gcgFILTERMASK<double> for double precision. Methods for projection
// and reconstruction are avalaible. The performance of the numeric methods can be
// analysed by using the basisInformation() method.
//

template<>
class     GCG_API_CLASS    gcgFILTERMASK1D<NUMTYPE> : public gcgDISCRETE1D<NUMTYPE> {
  public:
    gcgFILTERMASK1D();
    virtual ~gcgFILTERMASK1D();

    // Box filter construction
    bool createBoxFilter(unsigned int nsamples);

    // Triangle filter construction
    bool createTriangleFilter(unsigned int nsamples);
};



// Signal comparison: error computation.

template<>    GCG_API_TEMPLATE    double computeMSEwith<NUMTYPE>(gcgDISCRETE1D<NUMTYPE> *src1, gcgDISCRETE1D<NUMTYPE> *src2); // Computes the Mean Squared Error.
template<>    GCG_API_TEMPLATE    double computePSNRwith<NUMTYPE>(gcgDISCRETE1D<NUMTYPE> *src1, gcgDISCRETE1D<NUMTYPE> *src2);// Computes the Peak Signal-to-Noise Ratio in dB.
template<>    GCG_API_TEMPLATE    double computeMAEwith<NUMTYPE>(gcgDISCRETE1D<NUMTYPE> *src1, gcgDISCRETE1D<NUMTYPE> *src2); // Computes the Mean Absolute Error.

template<>    GCG_API_TEMPLATE    double computeMSEwith<NUMTYPE>(gcgDISCRETE2D<NUMTYPE> *src1, gcgDISCRETE2D<NUMTYPE> *src2); // Computes the Mean Squared Error.
template<>    GCG_API_TEMPLATE    double computePSNRwith<NUMTYPE>(gcgDISCRETE2D<NUMTYPE> *src1, gcgDISCRETE2D<NUMTYPE> *src2);// Computes the Peak Signal-to-Noise Ratio in dB.
template<>    GCG_API_TEMPLATE    double computeMAEwith<NUMTYPE>(gcgDISCRETE2D<NUMTYPE> *src1, gcgDISCRETE2D<NUMTYPE> *src2); // Computes the Mean Absolute Error.


// Signal analysis.

// Computes the gradient of 1D signal src. The derivative is computed using well-known Sobel operator and
// the result is stored in dx.
template<>    GCG_API_TEMPLATE    bool gcgFirstDerivativeSobel1D<NUMTYPE>(gcgDISCRETE1D<NUMTYPE> *src, gcgDISCRETE1D<NUMTYPE> *dx);

// Computes the gradient of 2D signal src. The partial derivatives are computed using well-known Sobel
// operator and the components are stored in dx and dy. The pointer for dx or dy may be NULL.
template<>    GCG_API_TEMPLATE    bool gcgGradientSobel2D<NUMTYPE>(gcgDISCRETE2D<NUMTYPE> *src, gcgDISCRETE2D<NUMTYPE> *dx, gcgDISCRETE2D<NUMTYPE> *dy);

// Computes the (optical) flow given the partial derivatives in X, Y and time of up to three channels.
// Uses the Augereau method (Bertrand Augereau, Benoît Tremblais, Christine Fernandez-Maloigne,
// "Vectorial computation of the optical flow in colorimage sequences", CIC05).
template<>    GCG_API_TEMPLATE    bool gcgOpticalFlowAugereau2D<NUMTYPE>(gcgDISCRETE2D<NUMTYPE> *dx1, gcgDISCRETE2D<NUMTYPE> *dy1, gcgDISCRETE2D<NUMTYPE> *dt1,
                                                         gcgDISCRETE2D<NUMTYPE> *dx2, gcgDISCRETE2D<NUMTYPE> *dy2, gcgDISCRETE2D<NUMTYPE> *dt2,
                                                         gcgDISCRETE2D<NUMTYPE> *dx3, gcgDISCRETE2D<NUMTYPE> *dy3, gcgDISCRETE2D<NUMTYPE> *dt3,
                                                         gcgDISCRETE2D<float> *outflowX, gcgDISCRETE2D<float> *outflowY);

template<>    GCG_API_TEMPLATE    bool gcgOpticalFlowAugereau2D<NUMTYPE>(gcgDISCRETE2D<NUMTYPE> *dx1, gcgDISCRETE2D<NUMTYPE> *dy1, gcgDISCRETE2D<NUMTYPE> *dt1,
                                                         gcgDISCRETE2D<NUMTYPE> *dx2, gcgDISCRETE2D<NUMTYPE> *dy2, gcgDISCRETE2D<NUMTYPE> *dt2,
                                                         gcgDISCRETE2D<NUMTYPE> *dx3, gcgDISCRETE2D<NUMTYPE> *dy3, gcgDISCRETE2D<NUMTYPE> *dt3,
                                                         gcgDISCRETE2D<double> *outflowX, gcgDISCRETE2D<double> *outflowY);



template<>     GCG_API_TEMPLATE    bool gcgLinearSolver(gcgDISCRETE2D<NUMTYPE> *A, gcgDISCRETE1D<NUMTYPE> *b, gcgDISCRETE1D<NUMTYPE> *x);

/*
// Tsai method

#define NONE 0
#define THREE_PARAM 1
#define FULL 2

class gcgCalibrator {
  public:
    double r1, r2, r3, r4, r5, r6, r7, r8, r9; //elementos da matriz de rotação
    double Tx, Ty, Tz; //elementos do vetor de translação
    double f; //distância focal
    double k1; //coeficiente de distorção
};

class gcgCameraCalibration: public gcgCalibrator{
    public:

        int n; //número de pontos para a calibração
        double Cx, Cy; //coordenadas da projeção do centro óptico
        double Rx, Ry, Rz; //ângulos da rotação

        VECTOR2d *worldPoints; //vetor de pontos do mundo
        VECTOR2d *imagePoints; //vetor de pontos da imagem
        VECTOR2d *distortedPoints; //vetor de pontos com distorção

        gcgCameraCalibration(int n, VECTOR2d *worldPoints, VECTOR2d *imagePoints, double Cx, double Cy);

        void rotationAngles();

        void rotationMatrix();

        int TsaiMethodCoplanar(int otimizationType);

        void otimization();

        void otimizationFull();
};

class gcgProjectorCalibration: public gcgCameraCalibration{
    public:

    double Rc[9], Tc[3], fc; //parâmetros da câmera calibrada

    gcgProjectorCalibration(int n, VECTOR2d *worldPoints, VECTOR2d *imagePoints, double Cx, double Cy,
                             double rc1, double rc2, double rc3, double rc4, double rc5,  double rc6, double rc7, double rc8, double rc9,
                             double Tcx, double Tcy, double Tcz, double fc);

    gcgProjectorCalibration(int n, VECTOR2d *worldPoints, VECTOR2d *imagePoints, double Cx, double Cy, double Rc[9], double Tc[3], double fc);

    gcgProjectorCalibration(int n, VECTOR2d *worldPoints, VECTOR2d *imagePoints, double Cx, double Cy, gcgCameraCalibration *camera);

    int projectorCalibration(int otimizationType);
};

class gcgCameraProjectorCalibration: public gcgCameraCalibration{
    public:

    gcgProjectorCalibration *projector;

    gcgCameraProjectorCalibration(int nc, VECTOR2d *worldPointsCamera, VECTOR2d *imagePointsCamera, double CxCamera, double CyCamera,
                                  int np, VECTOR2d *worldPointsProjetor, VECTOR2d *imagePointsProjetor, double CxProjetor, double CyProjetor);

    int cameraProjectorCalibration(int otimizationType);
};


// Levenberg-Marquadt method
class gcgLevMarq{
    public:

        int numPt; //número de pontos da função
        int numPar; //número de parâmetros a serem otimizados
        gcgDISCRETE1D<double> *p; //vetor que armazena os parâmetros
        gcgDISCRETE1D<double> *x; //vetor de tamanho >= numPt para armazenar a saída da função

        double epsilon1;
        double epsilon2;
        double epsilon3;
        double tau;
        int kmax;
        double jacobianInterval;

        gcgLevMarq();

        virtual void f(int numPt, int numPar, gcgDISCRETE1D<double> *p, gcgDISCRETE1D<double> *residuos) = 0;

        void multTransposeXMatrix(gcgDISCRETE2D<double> *matriz, gcgDISCRETE2D<double> *resultado);

        void multTransposeXMatrix(gcgDISCRETE2D<double> *matriz, gcgDISCRETE1D<double> *vetor, gcgDISCRETE1D<double> *resultado);

        double infinityNorm(gcgDISCRETE1D<double> *vetor);

        double euclidianNorm(gcgDISCRETE1D<double> *vetor);

        double multVectorTXVector(gcgDISCRETE1D<double> *vetor1, gcgDISCRETE1D<double> *vetor2);

        void jacobianCalculate(gcgDISCRETE1D<double> *p, double intervalo, gcgDISCRETE2D<double> *jacobiana);

        void otimization();
};

*/


// Update counter
#if _COUNTER_GCG_ == 1
    #undef _COUNTER_GCG_
    #define _COUNTER_GCG_ 2
#elif _COUNTER_GCG_ == 2
    #undef _COUNTER_GCG_
    #define _COUNTER_GCG_ 3
#elif _COUNTER_GCG_ == 3
    #undef _COUNTER_GCG_
    #define _COUNTER_GCG_ 4
#elif _COUNTER_GCG_ == 4
    #undef _COUNTER_GCG_
    #define _COUNTER_GCG_ 5
#else
    #define _GCG_H_ // Stop loading
#endif

// Loads several times to automatically include classes with distinct types
#if !defined(_GCG_H_)
  #include "gcg.h"
#endif

#ifdef __cplusplus
} // extern "C++"
#endif

#endif // #ifndef _GCG_H_