#include "voxeldata.h"
#include <qdatastream.h>
#include <QPixmap>
#include <math.h>


bool cosTableInited = false;
int cosTable[720];
int *VoxelData::cosT = &(cosTable[0]);
int *VoxelData::sinT = &(cosTable[90]); // sin(x) is basically same as cos(x + 90) (in degrees)
double sqrt2;               // square root of 2
const int scale = 65536/2;  // fixed point maths constants
const int scaleBits = 15;


// Windows math.h header lacks M_PI definition, it's available almost everywhere else as a #define
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif


VoxelData::VoxelData( int w, int h, int d )
{
    vw = w;
    vh = h;
    vd = d;
    init();
    if ( !cosTableInited ) { 
        cosTableInited = true;
        for ( int i = 0; i < 720; i++ ) 
            cosT[i] = int(scale*cos(i*M_PI/180.0));
        sqrt2 = sqrt(2);
    }
}


VoxelData::~VoxelData()
{
    for ( int z = 0; z < vd; z++ ) {
        for ( int x = 0; x < vw; x++ ) {
            delete [] (layers[z])[x];
        }
        delete [] layers[z];
    }
}


QRgb **VoxelData::layer( int z )
{
    return layers[z];
}


void VoxelData::init()
{
    layers.resize(vd);
    for ( int z = 0; z < vd; z++ ) {
        layers[z] = new QRgb * [vh];
        for ( int x = 0; x < vw; x++ ) {
            (layers[z])[x] = new QRgb [vh];
            memset( layer(z)[x], 0, sizeof(QRgb)*vh );
        }
    }
}


// clamp a value to between 0 and 255
// probably better if clamp doesn't need to be used as branching code doesn't perform
// well on some processors
#define clamp(x) \
                        ((x) > 255) ? 255 : ((x) < 0) ? 0 : (x)

// Idea for a clamp macro that doesn't use branching code (requires some temporary variables
// which may use more registers on some CPUs)
#define fastClampIdea(x) \
                        b = 0 - ((x & ~0xFF) == 0); \
                        y = (x & b) | (~b & 0xFF);


QPixmap VoxelData::render( int size, int rx, int ry, int rz )
{
    // Use a fixed width and height near to max possible for different rotations
    int imgSize = int((vw + vh + vd) / sqrt2);
    int dimension = imgSize << size; // size used as a power of 2 scale factor

    QImage img(dimension, dimension, 32);
    uint **imgData = (uint**)img.jumpTable();

    // Use DDA algorithm to trace a ray through the voxel data
    int sx = sinT[rx], cx = cosT[rx];
    int sy = sinT[ry], cy = cosT[ry];
    int sz = sinT[rz], cz = cosT[rz];

    // s and c -> sin and cos
    // d -> delta
    // u, v and w are voxel data coordinates
    int sysx = (sy*sx)>>scaleBits;
    int cxsz = (cx*sz)>>scaleBits;
    int cxcz = (cx*cz)>>scaleBits;

    int du0 = cxcz-((sysx*sz)>>scaleBits);
    int du1 = (cy*sz)>>scaleBits;
    int du2 = (-cz*sx-sy*cxsz)>>scaleBits;

    int dv0 = -cxsz-((sysx*cz)>>scaleBits);
    int dv1 = (cz*cy)>>scaleBits;
    int dv2 = (sz*sx-sy*cxcz)>>scaleBits;

    int dw0 = (sx*cy)>>scaleBits;
    int dw1 = sy;
    int dw2 = (cx*cy)>>scaleBits;

    // voxel coordinate deltas for each change in y of image coordinate system
    du0 >>= size;
    dv0 >>= size;
    dw0 >>= size;

    // voxel coordinate deltas for each change in x of image coordinate system
    du1 >>= size;
    dv1 >>= size;
    dw1 >>= size;

    // voxel coordinate deltas for each change in z of image coordinate system
    du2 >>= size;
    dv2 >>= size;
    dw2 >>= size;

    // voxel coordinates when x, y and z are zero
    int u0 = (vw*scale - dimension*(du0+du1+du2))/2;
    int v0 = (vh*scale - dimension*(dv0+dv1+dv2))/2;
    int w0 = (vd*scale - dimension*(dw0+dw1+dw2))/2;

    // x, y and z are in the coordinate system of the output image
    for ( int y = 0; y < dimension; y++ ) {
        uint *imgRow = imgData[y];
        int u1 = u0;
        int v1 = v0;
        int w1 = w0;
        for ( int x = 0; x < dimension; x++ ) {
            int u2 = u1;
            int v2 = v1;
            int w2 = w1;
            bool foundVoxel = false;
            bool haveProbedInsideModel = false; // true when u, v and w have been in range 

            for ( int z = 0; z < dimension; z++ ) {
                int u = u2 >> scaleBits;
                int v = v2 >> scaleBits;
                int w = w2 >> scaleBits;
                // ensure our u, v and w array indexes are within range
                if ( u >= 0 && u < vw && v >= 0 && v < vh && w >= 0 && w < vd ) {
                    haveProbedInsideModel = true; 
                    QRgb pix = layer(w)[u][v];
                    if ( pix ) {
                        // cheap lighting based on distance from one of the corners
                        double light = (x + (dimension - y) + z) / (double)dimension - 1.0;
                        int r = int(light * qRed(pix));
                        int g = int(light * qGreen(pix));
                        int b = int(light * qBlue(pix));
                        // keep the colour components in the range of 0-255 (after applying lighting)
                        //*imgRow = qRgb(clamp(r), clamp(g), clamp(b));
                        *imgRow = qRgb(r, g, b);
                        foundVoxel = true;
                        break; // leave the z for-loop
                    }
                } else if ( haveProbedInsideModel ) {
                    // u, v and w were in range and now are not. Means that we have now stepped
                    // through to the other side of the model as projected on to the image
                    // coordinates and as we keep stepping through we won't find anything else.
                    break; // leave the z for-loop
                } else {
                    // Skip ahead with faster loops to more quickly get u, v and w within range
                    // probably with some algerbra might be able to work it out
                    do {
                        u2 += du2;
                        v2 += dv2;
                        w2 += dw2;
                        u = u2 >> scaleBits;
                        z++;
                    } while ( ( u < 0 || u > vw ) && z < dimension );
                    do {
                        u2 += du2;
                        v2 += dv2;
                        w2 += dw2;
                        v = v2 >> scaleBits;
                        z++;
                    } while ( ( v < 0 || v > vh ) && z < dimension );
                    do {
                        u2 += du2;
                        v2 += dv2;
                        w2 += dw2;
                        w = w2 >> scaleBits;
                        z++;
                    } while ( ( w < 0 || w > vd ) && z < dimension );
                }
                u2 += du2;
                v2 += dv2;
                w2 += dw2;
            }
            // if we don't find a voxel, paint it black
            if (!foundVoxel) 
                *imgRow = qRgb( 0, 0, 0 );
            imgRow++;
            u1 += du1;
            v1 += dv1;
            w1 += dw1;
        }
        u0 += du0;
        v0 += dv0;
        w0 += dw0;
    }

    QPixmap pix;
    pix.convertFromImage( img );    
    return pix;
}


void VoxelData::write( QDataStream &ds )
{
    ds << vw;
    ds << vh;
    ds << vd;
    for ( int z = 0; z < vd; z++ ) {
        QRgb **l = layer(z);
        for ( int y = 0; y < vh; y++ ) {
            for ( int x = 0; x < vw; x+=8 ) {

                // Send a byte mask of the next 8 colors
                unsigned short v = 0;
                for ( int i = 0; i < 8 && x+i<vw; i++, v <<= 1 ) 
                    if ( l[x+i][y] ) 
                        v |= 0x1;
                v >>= 1;
                unsigned char v2 = (unsigned char)v;
                ds << v2;
                for ( int i = 0; i < 8 && x+i<vw; i++ ) 
                    if ( l[x+i][y] ) 
                        ds << l[x+i][y];

            }
        }
    }
}


void VoxelData::read( QDataStream &ds )
{
    ds >> vw;
    ds >> vh;
    ds >> vd;
    init();

    for ( int z = 0; z < vd; z++ ) {
        QRgb **l = layer(z);
        for ( int y = 0; y < vh; y++ ) {
            for ( int x = 0; x < vw; x+=8 ) {

                // Read a byte mask of the next 8 colors
                unsigned char v = 0xFF;
                ds >> v;
                for ( int i = 0; i < 8 && x+i<vw; i++, v <<= 1 ) 
                    if ( v & 0x80 )
                        ds >> l[x+i][y];
                    else
                        l[x+i][y] = 0;

            }
        }
    }
}


QRgb VoxelData::colorOf( int x, int y, int z, int rx, int ry, int rz )
{
    QRgb rgb = 0;

    if ( layer(z)[x][y] ) {
        double cosz = cos(-rz*M_PI/180.0);
        double sinz = sin(-rz*M_PI/180.0);
        double cosy = cos(-ry*M_PI/180.0);
        double siny = sin(-ry*M_PI/180.0);
        double cosx = cos(-rx*M_PI/180.0);
        double sinx = sin(-rx*M_PI/180.0);

        int tx = x-vw/2;
        int ty = y-vh/2;
        int tz = z;

        // y
        int vx = int(tx*cosy + z*siny);
        int vz = int(tx*siny + z*cosy);
        tx = vx;

        // z
        vx = int(tx*cosz - ty*sinz);
        int vy = int(tx*sinz + ty*cosz);
        ty = vy;
        tz = vz;

        // x
        vy = int(ty*cosx - tz*sinx);
        vz = int(ty*sinx + tz*cosx);

        vx += vw/2;
        vy += vh/2;

        rgb = layer(z)[x][y];
        int r = qRed(rgb);
        int g = qGreen(rgb);
        int b = qBlue(rgb);
        r -= r*vy/vh/2;
        g -= g*vy/vh/2;
        b -= b*vy/vh/2;

        rgb = qRgb( QMAX(r,0), QMAX(g,0), QMAX(b,0) );
    }

    return rgb;
}