Uniform Spatial Subdivision

#include <math.h>
#include <stdio.h>
#include <assert.h>
#include <sys/stat.h>
#include <stdlib.h>
#include <stdio.h>
#include <vector>
#include <GL/freeglut.h>

const int IMAGESIZE=128;
const int Rows=IMAGESIZE;
const int Cols=IMAGESIZE;

struct Float3
{
    float X;
    float Y;
    float Z;
    Float3(float x=0, float y=0, float z=0)
        :
    X(x), Y(y), Z(z){}

    Float3(const Float3& rhs)
        :
    X(rhs.X), Y(rhs.Y), Z(rhs.Z){}

   inline Float3& operator+=(const Float3& rhs)
   {    X += rhs.X; Y += rhs.Y; Z += rhs.Z; return *this; }

   inline Float3& operator-=(const Float3& rhs)
    {X -= rhs.X; Y -= rhs.Y; Z -= rhs.Z; return *this;}

   inline Float3& operator*=(const float& rhs)
    {X *= rhs; Y *= rhs; Z *= rhs; return *this;}

   inline Float3& operator/=(const float& rhs)
    {X /= rhs; Y /= rhs; Z /= rhs; return *this;}

   inline bool operator!=(const Float3& rhs)
    {return X!=rhs.X || Y!=rhs.Y || Z!=rhs.Z; }

   inline Float3 operator*(float rhs) const
    {return Float3(X*rhs, Y*rhs, Z*rhs);   }

   inline Float3 operator+(Float3& rhs) const
    {return Float3(X+rhs.X, Y+rhs.Y, Z+rhs.Z);}

   inline Float3 operator-(Float3& rhs) const
    {return Float3(X-rhs.X, Y-rhs.Y, Z-rhs.Z);   }

   inline Float3 operator/(Float3& rhs) const
    {return Float3(X/rhs.X, Y/rhs.Y, Z/rhs.Z);   }

    inline float length()
    {return sqrt(X*X +Y*Y + Z*Z); }

   inline float lengthsq()
    {return X*X +Y*Y + Z*Z; }
};

typedef struct
{
    Float3 p1;
    Float3 p2;
    Float3 p3;
    Float3 n1;
    Float3 n2;
    Float3 n3;
}triangle;

int nTriangles;
triangle *T;
Float3 *normalArray;

struct Bound
{
    Float3 minBound;
    Float3 maxBound;
    Float3 center;
    float dia;
}B;

struct uchar4 {
    unsigned char r,g,b,a;
};

uchar4*raster;
uchar4 *colorArray;

struct Camera {
    Float3 From,At,Up;
};

Camera Cam = {Float3(0.0f,1.0f/3.0f,2.0f),
              Float3(0.0f,1.0f/3.0f,0.0f),
              Float3(0.0f,1.0f,0.f)};

struct  Ray{
    Float3 ray_d;
    Float3 ray_o;
    float tMin, tMax;
};

typedef struct {
    int nx,ny,nz;
    Bound *box;
    int Ghit(Ray ray, int hitrec, Bound *box);
    triangle tri;
}Grid;
Grid G;
std::vector<Grid> cell;


struct Co_ordinateFrame
{
    Float3 x_axis;
    Float3 y_axis;
    Float3 z_axis;
    Float3 origin;
};
Co_ordinateFrame CF;

float txmin, tymin, tzmin;
float txmax, tymax, tzmax;

void readModel(char *filename, Bound *B, Camera *cam)//, Grid *G
{
    B->minBound.X = B->minBound.Y = B->minBound.Z = 1e5;
    B->maxBound.X = B->maxBound.Y = B->maxBound.Z = -1e5;

    txmin = tymin = tzmin = 1e5;
    txmax = tymax = tzmax = -1e5;

    int nVertices;
    FILE *fl=fopen(filename,"r");
    if (fl==NULL)
    {
        printf("Can not read file.\n");
        exit(1);
    }
    fscanf(fl,"%d", &nVertices);
    float *V = (float *) malloc(nVertices*3*sizeof(float));
    float *N = (float *) malloc(nVertices*3*sizeof(float));
    float *tex = (float *) malloc(nVertices*2*sizeof(float));
    for (int i=0; i<nVertices; i++){
        fscanf(fl,"%f%f%f%f%f%f%f%f",V+i*3,V+i*3+1,V+i*3+2,N+i*3,N+i*3+1,N+i*3+2,tex+i*2,tex+i*2+1);
        if (V[i*3] < B->minBound.X) B->minBound.X=V[i*3];
        if (V[i*3] > B->maxBound.X) B->maxBound.X=V[i*3];
        if (V[i*3+1] < B->minBound.Y) B->minBound.Y=V[i*3+1];
        if (V[i*3+1] > B->maxBound.Y) B->maxBound.Y=V[i*3+1];
        if (V[i*3+2] < B->minBound.Z) B->minBound.Z=V[i*3+2];
        if (V[i*3+2] > B->maxBound.Z) B->maxBound.Z=V[i*3+2];
    }
    B->center.X = (B->minBound.X + B->maxBound.X)/2.0f;
    B->center.Y = (B->minBound.Y + B->maxBound.Y)/2.0f;
    B->center.Z = (B->minBound.Z + B->maxBound.Z)/2.0f;

    float dX = (B->maxBound.X - B->minBound.X);
    float dY = (B->maxBound.Y - B->minBound.Y);
    float dZ = (B->maxBound.Z - B->minBound.Z);

    B->dia = sqrt(dX*dX+dY*dY+dZ*dZ);
    cam->At = B->center;
    if ((dX > dY) && (dX > dZ)){
    cam->Up.X = 1.0f; cam->Up.Y = cam->Up.Z = 0.0f;
    if (dY > dZ){
    cam->From.Z = B->center.Z + B->dia;
    cam->From.X = B->center.X; cam->From.Y = B->center.Y;
    }
    else{
    cam->From.Y = B->center.Y + B->dia;
    cam->From.X = B->center.X; cam->From.Z = B->center.Z;
    }
    }
    else if((dY > dX) && (dY > dZ)){
    cam->Up.Y = 1.0f; cam->Up.X = cam->Up.Z = 0.0f;
    if (dX > dZ){
    cam->From.Z = B->center.Z + B->dia;
    cam->From.Y = B->center.Y; cam->From.X = B->center.X;
    }
    else {
    cam->From.X = B->center.X + B->dia;
    cam->From.Y = B->center.Y; cam->From.Z = B->center.Z;
    }
    }
    else {
    cam->Up.Z = 1.0f; cam->Up.X = cam->Up.Y = 0.0f;
    if (dX > dY){
    cam->From.Y = B->center.Y + B->dia;
    cam->From.X = B->center.X; cam->From.Z = B->center.Z;
    }
    else{
    cam->From.X = B->center.X - B->dia;
    cam->From.Y = B->center.Y; cam->From.Z = B->center.Z;
    }
    }

    G.nx = B->maxBound.X/(IMAGESIZE-1);
    G.ny = B->maxBound.Y/(IMAGESIZE-1);
    G.nz = B->maxBound.Z/(IMAGESIZE-1);

    fscanf(fl,"%d", &nTriangles);
    char textureFilename[256];
    fscanf(fl,"%s",textureFilename);
    printf("[Vertices: %d], [Triangles: %d], [Texture File: %s]\n",
    nVertices, nTriangles, textureFilename);
    //E = (int *) malloc(nTriangles*3*sizeof(int));
    T = (triangle *)malloc(nTriangles*sizeof(triangle));
    normalArray = (Float3 *)malloc(nTriangles*sizeof(Float3));
    int E1, E2, E3;
    for (int i=0; i<nTriangles; i++){
        fscanf(fl,"%d%d%d",&E1,&E2,&E3);
        T[i].p1.X = V[E1*3];T[i].p1.Y = V[E1*3+1];T[i].p1.Z = V[E1*3+2];
        T[i].p2.X = V[E2*3];T[i].p2.Y = V[E2*3+1];T[i].p2.Z = V[E2*3+2];
        T[i].p3.X = V[E3*3];T[i].p3.Y = V[E3*3+1];T[i].p3.Z = V[E3*3+2];
        //-- Set n1, n2, n3
        T[i].n1.X = N[E1*3];T[i].n1.Y = N[E1*3+1];T[i].n1.Z = N[E1*3+2];
        T[i].n2.X = N[E1*3];T[i].n2.Y = N[E1*3+1];T[i].n2.Z = N[E1*3+2];
        T[i].n3.X = N[E1*3];T[i].n3.Y = N[E1*3+1];T[i].n3.Z = N[E1*3+2];
    }

    fclose(fl);
    return;
}

Float3 Scale(Float3 vec1,float s)
{
    Float3 res;
    res.X = vec1.X*s;
    res.Y = vec1.Y*s;
    res.Z = vec1.Z*s;

    return res;
}
  float Dot(Float3 vec1,Float3 vec2)
 {
    float dot = vec1.X*vec2.X+vec1.Y*vec2.Y +vec1.Z*vec2.Z;

    return dot;
 }

 Float3 Cross(Float3 vec1,Float3 vec2)
 {
     Float3 cross;

     cross.X = (vec1.Y*vec2.Z) - (vec1.Z*vec2.Y);
     cross.Y = (vec1.Z*vec2.X) - (vec1.X*vec2.Z);
     cross.Z = (vec1.X*vec2.Y) - (vec1.Y*vec2.X);

     return cross;
 }
 Float3 normalize(Float3 vec)
 {

     float mag = sqrt(Dot(vec,vec));
     Float3 result;
     result.X = vec.X /mag;
     result.Y = vec.Y/mag;
     result.Z = vec.Z/mag;

     return result;
 }

 void  ComputeCameraframe()
 {
    CF.origin = Cam.From;
    CF.z_axis = normalize(Cam.From - Cam.At);
    CF.x_axis = normalize(Cross(Cam.Up,CF.z_axis));
    CF.y_axis = normalize(Cross(CF.z_axis,CF.x_axis));
 }

 Ray get_Ray(int row, int col)
 {
    Ray ray;
    ray.ray_o = CF.origin;//eye point

    Float3  scale_resx = Scale(CF.x_axis,1.0-(2.0*col)/Cols);
    Float3  scale_resy = Scale(CF.y_axis,1.0-(2.0*row)/Rows);
    ray.ray_d = scale_resy - CF.z_axis - scale_resx;
    ray.tMin = 0.0f;
    ray.tMax = 10000000000.0f;
    return ray;
}

bool ray_box_intersection(Ray ray, Bound box)
{

       if(ray.ray_d.X >= 0.0f)
       {
         txmin = (box.minBound.X - ray.ray_o.X)/ray.ray_d.X;
         txmax = (box.maxBound.X - ray.ray_o.X)/ray.ray_d.X;
       }
       else
       {
           txmin = (box.maxBound.X - ray.ray_o.X)/ray.ray_d.X;
           txmax = (box.minBound.X - ray.ray_o.X)/ray.ray_d.X;
       }

       if(ray.ray_d.Y >= 0.0f)
       {
            tymin = (box.minBound.Y - ray.ray_o.Y)/ray.ray_d.Y;
            tymax = (box.maxBound.Y - ray.ray_o.Y)/ray.ray_d.Y;
       }
       else
       {
           tymin = (box.maxBound.Y - ray.ray_o.Y)/ray.ray_d.Y;
           tymax = (box.minBound.Y - ray.ray_o.Y)/ray.ray_d.Y;
       }
       if(txmin > tymax || tymin > txmax)
           return false;
       if(tymin > txmin)
           txmin = tymin;
       if(tymax < txmax)
           txmax = tymax;

       if(ray.ray_d.Z >= 0.0f)
       {
             tzmin = (box.minBound.Z - ray.ray_o.Z)/ray.ray_d.Z;
             tzmax = (box.maxBound.Z - ray.ray_o.Z)/ray.ray_d.Z;
       }
       else
       {
           tzmin = (box.maxBound.Z - ray.ray_o.Z)/ray.ray_d.Z;
           tzmax = (box.minBound.Z - ray.ray_o.Z)/ray.ray_d.Z;
       }
       if(txmin > tzmax || tzmin > txmax)
           return false;
       if(tzmin > txmin)
           txmin = tzmin;
       if(tzmax < txmax)
           txmax = tzmax;
    if (txmin < ray.tMax && txmax > ray.tMin){
        if (txmin > ray.tMin) ray.tMin = txmin;
        if (txmax < ray.tMax) ray.tMax = txmax;
    }
    else return false;
}

int ray_trianglelist_intersection(Ray ray,  triangle *triangles, float *ray_t, float *t_beta, float *t_gamma)
{
    float t,beta,gamma;
    float a,b,c,d,e, f, g, h, i,j,k,l;

    float nearest_t=ray.tMax;
    float nearest_beta, nearest_gamma;
    int hitId = -1;

    for(int nt = 0;nt<nTriangles;nt++)
    {
        Float3 col1 = triangles[nt].p1 - triangles[nt].p2;
        Float3 col2 = triangles[nt].p1 - triangles[nt].p3;

        a = col1.X;   d = col2.X;  g = ray.ray_d.X;
        b = col1.Y;   e = col2.Y;  h = ray.ray_d.Y;
        c = col1.Z;   f = col2.Z;  i = ray.ray_d.Z;

        Float3 col3 = triangles[nt].p1 - ray.ray_o;

        j = col3.X;
        k = col3.Y;
        l = col3.Z;

        float M = a*(e*i-h*f)+b*(g*f-d*i)+c*(d*h-e*g);
        //float M = a*(e*i-h*f)+d*(h*c-b*i)+g*(b*f-e*c);
        //compute t

        t = -(f*(a*k-j*b)+e*(j*c-a*l)+d*(b*l-k*c))/M;
        //t = (a*(e*l-f*k) + d*(k*c-b*l) + j*(b*f-c*e))/M;

        if(t<ray.tMin || t>nearest_t)
            continue;
        //compute gamma
        gamma = (i*(a*k-j*b)+h*(j*c-a*l)+g*(b*l-k*c))/M;
        //gamma = (a*(k*i-l*h)-j*(b*i-c*h)+g*(b*l-c*k))/M;
        if(gamma<0.0f || gamma>1.0f)
            continue;
        //compute beta
        beta = (j*(e*i-h*f)+k*(g*f-d*i)+l*(d*h-e*g))/M;
        //beta = (j*(e*i-h*f)-d*(k*i-l*h)+g*(k*f-l*e))/M;

        if(beta<0.0f || beta > 1.0f-gamma)
            continue;
        nearest_t = t; nearest_beta = beta; nearest_gamma = gamma;
        hitId = nt;
    }
    if (nearest_t < ray.tMax){
        *ray_t = nearest_t; *t_beta = nearest_beta; *t_gamma = nearest_gamma;
        return hitId;
    }
    else
        return -1;
}


int Grid::Ghit(Ray ray, int hitrec, Bound *box)//,
{
    int I_inc, I_stop, J_inc, J_stop, K_inc, K_stop;
    I_inc = J_inc = K_inc = 1;
    I_stop = J_stop =  K_stop = 1;
    float dtX,dtY,dtZ;
    dtX = dtY = dtZ = 0.0f;
    float tXnext, tYnext, tZnext , tlast;

    Float3 pointOnRay = ray.ray_o + Scale(ray.ray_d, ray.tMax);

    if(ray.ray_d.X >= 0.0f )
    {
        txmin = (box->minBound.X - ray.ray_o.X)/ray.ray_d.X;
        I_inc += 1;
        I_stop = nx;
    }
    else
    {
        txmin = (box->maxBound.X - ray.ray_o.X)/ray.ray_d.X;
        txmax = (box->minBound.X - ray.ray_o.X)/ray.ray_d.X;

        I_inc -= 1;
        I_stop = -1;
    }

    if(ray.ray_d.Y >= 0.0f)
        {

          tymin = (box->minBound.Y - ray.ray_o.Y)/ray.ray_d.Y;
          tymax = (box->maxBound.Y - ray.ray_o.Y)/ray.ray_d.Y;
          J_inc += 1;
          J_stop = ny;
        }
       else
        {
           tymin = (box->maxBound.Y - ray.ray_o.Y)/ray.ray_d.Y;
           tymax = (box->minBound.Y - ray.ray_o.Y)/ray.ray_d.Y;
           J_inc -= 1;
           J_stop = -1;
        }

    if(txmin > tymax || tymin > txmax) return -1;

    if(ray.ray_d.Z >= 0.0f)
    {
        tzmin = (box->minBound.Z - ray.ray_o.Z)/ray.ray_d.Z;
        tzmax = (box->maxBound.Z - ray.ray_o.Z)/ray.ray_d.Z;
        K_inc += 1;
        K_stop = nz;
    }
    else
    {
        tzmin = (box->maxBound.Z - ray.ray_o.Z)/ray.ray_d.Z;
        tzmax = (box->minBound.Z - ray.ray_o.Z)/ray.ray_d.Z;
        K_inc -= 1;
        K_stop = -1;
    }

    if(tzmin > txmax || txmin > tzmax) return -1;
    if(tzmin > tymax || tymin > tzmax) return -1;

    dtX = (txmax - txmin)/nx;
    dtY = (tymax - tymin)/ny;
    dtZ = (tzmax - tzmin)/nz;
    int  I, J, K;

    // if a point on the ray is inside the grid

    if(ray.tMin >= box->minBound.X)
    {
        I = floor(nx*(pointOnRay.X - box->minBound.X)/(box->maxBound.X - box->minBound.X));
        if(I < 0) I = 0;
        if(I > nx - 1) I = nx - 1;

        J = floor(ny*(pointOnRay.Y- box->minBound.Y)/(box->maxBound.Y - box->minBound.Y));
        if(J < 0) J = 0;
        if(J > ny-1) J = ny - 1;

        K = floor(nz*(pointOnRay.Z - box->minBound.Z)/(box->maxBound.Z - box->minBound.Z));
        if(K < 0 ) K = 0;
        if(K > nz-1) K = nz - 1;
        tXnext = txmin + (I + (I_inc + 1)/2)*dtX;
        tYnext = tymin + (J + (J_inc + 1)/2)*dtY;
        tZnext = tzmin + (K + (I_inc + 1)/2)*dtZ;
        tlast = ray.tMin;
    }
    else if (txmin > tymin && txmin > tzmin)
    {
        I = I_stop - I_inc;
        pointOnRay.Y = ray.ray_o.Y + ray.tMin;//txminyd;
        J = floor(ny*(pointOnRay.Y - box->minBound.Y)/(box->maxBound.Y- box->minBound.Y));
        if(J < 0) J = 0;
        if(J > ny-1) J = ny-1;

        pointOnRay.Z = ray.ray_o.Z + ray.tMin; //tyminZd;
        K = floor(nz*(pointOnRay.Z - box->minBound.Z)/(box->maxBound.Z- box->minBound.Z));
        if(K < 0 ) K = 0;
        if(K > nz-1) K = nz-1;

        tXnext = txmin + dtX;
        tYnext = tymin + (J + (J_inc + 1)/2)*dtY;
        tZnext = tzmin + (K + (K_inc + 1)/2)*dtZ;
        tlast = txmin;
    }
    else if(tymin > tzmin)
    {
        J = J_stop - J_inc;
        pointOnRay.X = ray.ray_o.X + ray.tMin; //tyminxd
        I = floor(nx*(pointOnRay.X - box->minBound.X)/(box->maxBound.X - box->minBound.X));
        if(I < 0) I = 0;
        if(I > nx-1) I = nx-1;

        pointOnRay.Z = ray.ray_o.Z + ray.tMin; //tyminzd
        K = floor(nz*(pointOnRay.Z - box->minBound.Z)/(box->maxBound.Z- box->minBound.Z));

        if(K < 0 ) K = 0;
        if(K > nz -1) K = nz -1;

        tYnext = tymin + dtY;
        tXnext = txmin + (I + (I_inc + 1)/2)*dtX;
        tZnext = tzmin + (K + (K_inc + 1)/2)*dtZ;
        tlast = tymin;
    }
    else K = K_stop - K_inc;

        pointOnRay.X = ray.ray_o.X + ray.tMin;//tzminxd
        I = floor(nx*(pointOnRay.X - box->minBound.X)/(box->maxBound.X- box->minBound.X));

        if (I < 0 ) I = 0;
        if (I > nx - 1) I = nx - 1;

        pointOnRay.Y = ray.ray_o.Y + ray.tMin; //tzminyd
        J = floor(nx*(pointOnRay.Y - box->minBound.Y)/(box->maxBound.Y- box->minBound.Y));

        if(J < 0) J=0;
        if(J > ny-1) J = ny-1;

        tZnext = tzmin + dtZ;
        tXnext = txmin + (I + (I_inc + 1)/2)*dtX;
        tYnext = tymin + (J + (J_inc + 1)/2)*dtZ;
        tlast = tzmin;

    if(tlast > ray.tMax) return -1;

    while(true)
    {
        if(tXnext < tYnext && tXnext < tZnext)
        {
            if(cell[I,J].Ghit(ray, hitrec, box) && hitrec < ray.tMax) //box //ray,tlast,tXnext,rec

                return hitrec;
                tlast = tXnext;
                tXnext += dtX;
                I += I_inc;
                if(I==I_stop) return -1;
        }
        else if(tYnext < tZnext)
        {
            if(cell[I,J].Ghit(ray, hitrec, box) && hitrec < ray.tMax) //box //ray,tlast,tXnext,rec
                    return hitrec;
                tlast = tYnext;
                tYnext += dtY;
                J += J_inc;
                if(J==J_stop) return -1;
        }
        else
            {
                if(cell[I,J].Ghit(ray,hitrec,box) && hitrec < ray.tMax) //box //ray,tlast,tXnext,rec
                    return hitrec;
                tlast = tZnext;
                tZnext += dtZ;
                K += K_inc;
                if(K==K_stop) return -1;
            }
    }
}

int ray_uss_intersection(Ray ray, Grid *G, float *ray_t, float *t_beta, float *t_gamma )
{
   int USS_hitID = 0;
   int  USS_hit = G->Ghit(ray, USS_hitID, &B);
   int M, N;
   M = G->nx;
   N = G->ny;
   if(USS_hit < 0)
   {
       while (&cell!= NULL)
       {
           int Tri_hitID = ray_trianglelist_intersection(ray, &cell[M,N].tri, ray_t, t_beta, t_gamma);
            if (Tri_hitID >= 0 )//&& the point is inside the cell
            {
                //grid = cell[USS_hit];
                return Tri_hitID;
            }
       }
   }
   return -1;
}

int ray_scene_intersection(Ray ray, Grid *G, Bound *B, float *ray_t, float *t_beta, float *t_gamma)
{
    // Step I: Do bounding box test
    if (ray_box_intersection(ray, *B))
        //return ray_trianglelist_intersection(ray, T, ray_t, t_beta, t_gamma);
        return ray_uss_intersection(ray, G, ray_t, t_beta, t_gamma);
    else return -1;
    //
    // Step II:  Do bounding box test
    // if (ray_box_intersection(ray, *B))
    //    //return ray_uss_intersection(ray, G);
    ////
    //return ray_trianglelist_intersection(ray, T, ray_t, t_beta, t_gamma);
}
uchar4 get_Color(Ray r)
 {
     uchar4 background_color = {0,0,0,255};
     uchar4 foreground_color; //= {100,200,150,255};
     float t,b,g;
     int hitId = ray_scene_intersection(r,&G,&B,&t,&b,&g);

    if (hitId < 0)
         return background_color;
     else{

         //--compute foreground color
        // find the point of intersection
          Float3 hitvert1 = T[hitId].n1;
          Float3 hitvert2 = T[hitId].n2;
          Float3 hitvert3 = T[hitId].n3;
          Float3 interNorm1 = Scale(hitvert1,(1.0f-b-g))+Scale(hitvert2,b)+Scale(hitvert3,g);

         // Float3 normal = normalArray[hitId];

         // Assume light at camera from point.
          Float3 light_Pos = Cam.From;
          Float3 pI = r.ray_o + Scale(r.ray_d,t);
        // Compute and normalize the light direction: normalize(minus(from,point_of_intersection))
          Float3 light_dir = normalize(light_Pos -pI);
          float ndotl = Dot(interNorm1,light_dir);
          if (ndotl<0.0f)ndotl = 0.0f;
           foreground_color.r = ndotl*255;
           foreground_color.g = ndotl*255;
           foreground_color.b = ndotl*255;
           foreground_color.a = 255;
        return foreground_color;
     }
}


unsigned char *updateRasterTexture()
 {
    for (int row=0;row<Rows;row++)
        for(int col=0;col<Cols;col++)
        {
            Ray ray_r = get_Ray(row, col);//get_Ray(Rows/2, Cols/2);//
            int index = row*Cols+col;
            raster[index]= get_Color(ray_r);
        }
    return (unsigned char *)raster;

 }

GLuint textureId;
 static const int nVertices = 4;

 static GLfloat V[nVertices*3] = {
    -1.0f, -1.0f, 0.0f,
    1.0f, -1.0f, 0.0f,
    1.0f,1.0f,0.0f,
    -1.0f,1.0f,0.0f,
 };

 static GLfloat Te[nVertices*2] = {
    0.0f,0.0f,
    1.0f,0.0f,
    1.0f,1.0f,
    0.0f,1.0f
 };

 void displayHandler()
 {
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    //glActiveTexture(GL_TEXTURE0);
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    glEnable(GL_TEXTURE_2D);
    glBindTexture(GL_TEXTURE_2D,textureId);
    unsigned char *buffer=NULL;
    buffer = updateRasterTexture();
    glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, Cols, Rows,
                  GL_RGBA, GL_UNSIGNED_BYTE, buffer);
    glBegin(GL_QUADS);
        glTexCoord2f(Te[0],1.0-Te[1]);
        glVertex3f(V[0],V[1],V[2]);
        glTexCoord2f(Te[2],1.0-Te[3]);
        glVertex3f(V[3],V[4],V[5]);
        glTexCoord2f(Te[4],1.0-Te[5]);
        glVertex3f(V[6],V[7],V[8]);
        glTexCoord2f(Te[6],1.0-Te[7]);
        glVertex3f(V[9],V[10],V[11]);
    glEnd();
    glDisable(GL_TEXTURE_2D);
    glBindTexture(GL_TEXTURE_2D,0);
    glutSwapBuffers();
 }
 void kbdHandler( unsigned char key, int x, int y )
 {
    switch( key )
    {
       case 27:
          exit(0);
          break;
      default: return;
    }
    glutPostRedisplay();
 }

 void reshapeHandler(int w, int h)
 {
     glViewport(0, 0, w, h);
     glMatrixMode(GL_PROJECTION);
     glLoadIdentity();

    glutPostRedisplay();
 }

 void resize(int windowWidth, int windowHeight) {
    glViewport(0, 0, windowWidth, windowHeight);
 }
 void inittexture()
 {
    //glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
    // Initialize texture
    glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
    glGenTextures(1, &textureId);                       // Generate OpenGL texture IDs
    glBindTexture(GL_TEXTURE_2D, textureId);                    // Bind Our Texture
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
    glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_DECAL);
    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, Cols, Rows, 0,
        GL_RGBA, GL_UNSIGNED_BYTE, NULL);
 //build our texture mipmaps
    static int textureSize;
    //gluBuild2DMipmaps( GL_TEXTURE_2D, 4, textureSize, textureSize,
    //         GL_RGBA, GL_UNSIGNED_BYTE, NULL );
    //glGenerateMipmap(GL_TEXTURE_2D);
 }
 int initGraphics(int argc, char *argv[], int imageSize)
 {
    // Initalize OpenGL Context
    glutInit(&argc, argv);
    glutInitWindowPosition(100,150);
    glutInitWindowSize(imageSize,imageSize);
    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH);
    glutCreateWindow (argv[0]);
    glClearColor(0,0,0,0);
    glShadeModel(GL_FLAT);
    glEnable(GL_DEPTH_TEST);
    inittexture();
    // Set Callback Handler
    glutDisplayFunc(displayHandler);
    //glutIdleFunc(Idle);
    glutReshapeFunc(reshapeHandler);
    glutKeyboardFunc(kbdHandler);
   return 0;
 }

 int main(int argc,char** argv)
{
    raster = (uchar4 *)malloc(Rows*Cols*sizeof(uchar4));
    readModel("legoModel.txt", &B, &Cam );
    ComputeCameraframe();
    if (initGraphics(argc,argv,512))
        {
            printf("Can not initialize graphics.\n");
            return 0;
        }
    glutMainLoop();
    return 0;
}