/*----------------------------------------------------------------------
 * texture3d.c GL example of the mesa 3d-texture extention to simulate 
 *             procedural texturing, it uses a lattice gradient noise and
 *             perlin turbulence functions.
 *
 * Author:   Daniel Barrero
 *           barrero@irit.fr
 *
 * cc texture3d.c -o texture3d -lglut -lMesaGLU -lMesaGL -lX11 -lXext -lm
 *
 *----------------------------------------------------------------------- */

#include <GL/gl.h>
#include <GL/glu.h>
#include <GL/glut.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
 
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#define B 256    
#define SIZE 16

float   noise[SIZE+1][SIZE+1][SIZE+1];

float rotX=0,rotY=0,rotZ=30;
int cube=1;

/*****
    system indepedant well behaved noise function, based on code in:
   "Random Number Generators: Good Ones are Hard to Find" 
    by Stephen K. Park and Keith W. Miller in Communications of the ACM,
   31, 10 (Oct. 1988) pp. 1192-1201.
*******/
#define K_A 16807
#define K_M 2147483647                  /* Mersenne prime 2^31 -1 */
#define K_Q 127773                      /* K_M div K_A */
#define K_R 2836                        /* K_M mod K_A */
double drand(long s)
{
 static long seed;
 long hi, lo;
 if(s!=0) {
    seed=s;
    drand(0);
    drand(0);
    drand(0);
 }
 hi = seed / K_Q;
 lo = seed % K_Q;
 if ((seed = K_A * lo - K_R * hi) <= 0)  seed += K_M;
 return ((float) seed / K_M);
}

void paintBox ( GLdouble x0, GLdouble x1, GLdouble y0, GLdouble y1,
        GLdouble z0, GLdouble z1, GLenum type)
{
    static GLdouble n[6][3] = {
        {-1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {1.0, 0.0, 0.0},
        {0.0, -1.0, 0.0}, {0.0, 0.0, 1.0}, {0.0, 0.0, -1.0}
    };
    static GLint faces[6][4] = {
        { 0, 1, 2, 3 }, { 3, 2, 6, 7 }, { 7, 6, 5, 4 },
        { 4, 5, 1, 0 }, { 5, 6, 2, 1 }, { 7, 4, 0, 3 }
    };
    GLdouble v[8][3], tmp;
    GLdouble vt[8][3];
    GLint i;

    if (x0 > x1) {
        tmp = x0; x0 = x1; x1 = tmp;
    }
    if (y0 > y1) {
        tmp = y0; y0 = y1; y1 = tmp;
    }
 if (z0 > z1) {
        tmp = z0; z0 = z1; z1 = tmp;
    }
    v[0][0] = v[1][0] = v[2][0] = v[3][0] = x0;
    v[4][0] = v[5][0] = v[6][0] = v[7][0] = x1;
    v[0][1] = v[1][1] = v[4][1] = v[5][1] = y0;
    v[2][1] = v[3][1] = v[6][1] = v[7][1] = y1;
    v[0][2] = v[3][2] = v[4][2] = v[7][2] = z0;
    v[1][2] = v[2][2] = v[5][2] = v[6][2] = z1;

    vt[0][0] = vt[1][0] = vt[2][0] = vt[3][0] = 0.0;
    vt[4][0] = vt[5][0] = vt[6][0] = vt[7][0] = 1.0;
    vt[0][1] = vt[1][1] = vt[4][1] = vt[5][1] = 0.0;
    vt[2][1] = vt[3][1] = vt[6][1] = vt[7][1] = 1.0;
    vt[0][2] = vt[3][2] = vt[4][2] = vt[7][2] = 0.0;
    vt[1][2] = vt[2][2] = vt[5][2] = vt[6][2] = 1.0;

    for (i = 0; i < 6; i++) {
      glBegin(type);
        glNormal3dv(&n[i][0]);
        glVertex3dv(&v[faces[i][0]][0]);
        glTexCoord3dv(&vt[faces[i][0]][0]);

        glNormal3dv(&n[i][0]);
        glVertex3dv(&v[faces[i][1]][0]);
        glTexCoord3dv(&vt[faces[i][1]][0]);
 
        glNormal3dv(&n[i][0]);
        glVertex3dv(&v[faces[i][2]][0]);
        glTexCoord3dv(&vt[faces[i][2]][0]);
 
        glNormal3dv(&n[i][0]);
        glVertex3dv(&v[faces[i][3]][0]);
        glTexCoord3dv(&vt[faces[i][3]][0]);
      glEnd();
    }
}

void drawCube(float size)
{
 paintBox(-size/2., size/2., -size/2., size/2.,-size/2., size/2.,GL_QUADS);
}                                         

float noise31(float vec[3])
{
 long tmp;
 tmp= 103*vec[0]+107*vec[1]+109*vec[2];
 return (float)(((rand()*tmp) % ( B + B)) - B) / B;
}

void initNoise()
{
  float   tmp,u,v;
  long    i,j, k, ii,jj,kk;
  srand(1);
  drand(1107);
  for (i=0; i<SIZE; i++)
    for (j=0; j<SIZE; j++)
      for (k=0; k<SIZE; k++)       {
        noise[i][j][k] = (float)drand(0);
      } 
   for (i=0; i<SIZE+1; i++)
     for (j=0; j<SIZE+1; j++)
       for (k=0; k<SIZE+1; k++)       { 
          ii = (i == SIZE)? 0:  i;
          jj = (j == SIZE)? 0:  j;   
          kk = (k == SIZE)? 0:  k;
          noise[i][j][k] = noise[ii][jj][kk];
       } 
}
float noise3(float pnt[3])
{
 float  t1;  
 float  p_l,p_l2,    /* value lerped down left side of face1 & face 2 */ 
        p_r,p_r2,    /* value lerped down left side of face1 & face 2 */
        p_face1,     /* value lerped across face 1 (x-y plane ceil of z) */
       p_face2,     /* value lerped across face 2 (x-y plane floor of z) */
        p_final;     /* value lerped through cube (in z)                  */
 float   tnoise;  
 register int      x, y, z,px,py,pz;    
 px = (int)pnt[0];  
 py = (int)pnt[1];
 pz = (int)pnt[2];
 x = px &(SIZE-1); /* make sure the values are in the table           */
 y = py &(SIZE-1); /* Effectively, replicates the table thought space */
 z = pz &(SIZE-1);
 t1 = pnt[1] - py; 
 p_l  = noise[x][y][z+1]+t1*(noise[x][y+1][z+1]-noise[x][y][z+1]);
 p_r  =noise[x+1][y][z+1]+t1*(noise[x+1][y+1][z+1]-noise[x+1][y][z+1]);
 p_l2 = noise[x][y][z]+ t1*( noise[x][y+1][z] - noise[x][y][z]);
 p_r2 = noise[x+1][y][z]+ t1*(noise[x+1][y+1][z] - noise[x+1][y][z]);
 t1 = pnt[0] - px;  
 p_face1 = p_l + t1 * (p_r - p_l);
 p_face2 = p_l2 + t1 * (p_r2 -p_l2);
 t1 = pnt[2] - pz;
 p_final =  p_face2 + t1*(p_face1 -p_face2); 
 return(p_final);
}

float turbulence(float point[3], float lofreq, float hifreq)
{
 float freq, t, p[3];

 p[0] = point[0] + 123.456;
 p[1] = point[1];
 p[2] = point[2];
 t = 0;
 for (freq = lofreq ; freq < hifreq ; freq *= 2.) {
     t += (fabs(noise3(p)) / freq ) ;
     p[0] *= 2.;
     p[1] *= 2.;
     p[2] *= 2.;
 }
 return t - 0.3; /* readjust to make mean value - 0.0 */
}


#define clamp(x,a,b)  ((x<a)?a:((x>b)?b:x))
void calcMarble(unsigned char *voxels,int tex_width,int tex_height,int tex_depth)
{
 unsigned char *vp;
 int i,j,k,l,m;
 float vec[3],vecs[3],vect[3];
 float tmp,tmp1,marble,f;
 initNoise();
 vp=voxels;
 vecs[2]=1.0/tex_depth;
 vecs[1]=1.0/tex_height;
 vecs[0]=1.0/tex_width;
 for (i=0,vec[0]=0;i<tex_width;i++){
    vec[0]+=vecs[0];
    for (j=0,vec[1]=0;j<tex_height;j++) {
       vec[1]+=vecs[1];
       for (k=0,vec[2]=0;k<tex_depth;k++) {
           vec[2]+=vecs[2];
           /* basic marble */
           /*tmp=(vec[2]*vec[1]*vec[0]*2.0*turbulence(vec,0.01,1));*/
           tmp=(vec[1]+3.0*turbulence(vec,0.01,1));
           marble= 0; f = 1;
           for (l = 0; l < 4; l += 1) {
                for(m=0; m<3;m++) { vect[m]=(vec[m]*3.7)*f; }
                marble += noise3(vect) * 1/f;
                f *= 2;
           }
           marble = clamp(4*marble - 3, 0, 1);
           /*tmp=marble;*/

           tmp=sin(tmp*M_PI)+1;
           tmp=sqrt(tmp)*.7071;
           tmp1=sqrt(tmp);
/*
           *vp++=0.3*255+0.6*tmp*255;
           *vp++=0.3*255+0.8*tmp1*255;
           *vp++=0.6*255+0.4*tmp*255;
           *vp++=255;
*/
           *vp++=(64+tmp1*92)*marble+(1-marble)*164;
           *vp++=(64+tmp*64)*marble+(1-marble)*164;
           *vp++=(102+tmp1*192)*marble+(1-marble)*255;
           *vp++=255;
       }
    }
 }
}


void createTexture3D(int tex_width,int tex_height,int tex_depth)
{
 unsigned char *voxels;                       /* texture data ptr */

 /* creates a 3D texture image of marble */
 printf("creating 3d textures...\n");
 voxels = (unsigned char  *) malloc((4*tex_width*tex_height*tex_depth));
 calcMarble(voxels,tex_width,tex_height,tex_depth);

 printf("setting up 3d texture...\n");
 glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
 glTexParameteri(GL_TEXTURE_3D_EXT, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
 glTexParameteri(GL_TEXTURE_3D_EXT, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
 glTexParameteri(GL_TEXTURE_3D_EXT, GL_TEXTURE_WRAP_S,     GL_REPEAT);
 glTexParameteri(GL_TEXTURE_3D_EXT, GL_TEXTURE_WRAP_T,     GL_REPEAT);
 glTexParameteri(GL_TEXTURE_3D_EXT, GL_TEXTURE_WRAP_R_EXT, GL_REPEAT);
 glTexImage3DEXT(GL_TEXTURE_3D_EXT, 0, GL_RGBA,
                    tex_width, tex_height, tex_depth,
                    0, GL_RGBA, GL_UNSIGNED_BYTE, voxels);   
 glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
 printf("finished setting up 3d texture image...\n");       
}

void init(void)
{
 /* Init lights */
 GLfloat mat_specular[]  = { 1.0, 1.0, 1.0, 1.0 };
 GLfloat mat_shininess[] = { 64.0 };  
 GLfloat gray[]          = { 0.6, 0.6, 0.6, 0.0 };
 GLfloat white[]         = { 1.0, 1.0, 1.0, 0.0 };
 GLfloat light_position[]= { 0.0, 0.5, 3.0, 0.0 };
 /* Texture reference planes */
 GLfloat sPlaneEqn[] = { 1.0, 0.0, 0.0, 0.5 };
 GLfloat tPlaneEqn[] = { 0.0, 1.0, 0.0, 0.5 };
 GLfloat rPlaneEqn[] = { 0.0, 0.0, 1.0, 0.5 };
 
 glLightfv(GL_LIGHT1, GL_POSITION, light_position);
 glLightfv(GL_LIGHT1, GL_AMBIENT,  white);
 glLightfv(GL_LIGHT1, GL_DIFFUSE,  white);
 glLightfv(GL_LIGHT1, GL_SPECULAR, white);
 glEnable(GL_LIGHTING);
 glEnable(GL_LIGHT1);

 /* Object underlying material properties */
 glMaterialfv(GL_FRONT, GL_SPECULAR,  mat_specular);
 glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess);

 /* See if the 3D texture extension is supported */
 if(!glutExtensionSupported("GL_EXT_texture3D")) {
   printf("Sorry this GL implementation (%s) does not support 3D texture extensions\n", (char*)(glGetString(GL_RENDERER)));
   exit(0);
 } 
 /* If the extension is supported generate the texture */
 createTexture3D(32,32,32);

/* Enable 3D texturing and atomatic texture coords generation */
 glEnable(GL_TEXTURE_3D_EXT);
 glTexGeni(GL_S, GL_TEXTURE_GEN_MODE, GL_OBJECT_LINEAR); 
 glTexGeni(GL_T, GL_TEXTURE_GEN_MODE, GL_OBJECT_LINEAR); 
 glTexGeni(GL_R, GL_TEXTURE_GEN_MODE, GL_OBJECT_LINEAR); 
 glTexGenfv(GL_S, GL_OBJECT_PLANE, sPlaneEqn);
 glTexGenfv(GL_T, GL_OBJECT_PLANE, tPlaneEqn);
 glTexGenfv(GL_R, GL_OBJECT_PLANE, rPlaneEqn);
 glEnable(GL_TEXTURE_GEN_S);
 glEnable(GL_TEXTURE_GEN_T);
 glEnable(GL_TEXTURE_GEN_R);

 /* enable blending functions so the texture can be lighted */ 
 glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
 glEnable(GL_BLEND); 
 glEnable(GL_AUTO_NORMAL);
 glEnable(GL_NORMALIZE);

 glClearColor(0.0, 0.0, 0.0, 0.0);
 glClearIndex( 0.0 );

 glShadeModel(GL_SMOOTH);  
 glDepthFunc(GL_LEQUAL);
 glEnable(GL_DEPTH_TEST);
 glEnable(GL_CULL_FACE);
/*
 glFrontFace(GL_CW);
*/
 glCullFace(GL_BACK);
}

void resize ( int w, int h)
{
 glViewport(0,0,w,h);
 glMatrixMode(GL_PROJECTION);
 glLoadIdentity();
 glFrustum(-1.0, 1.0, -1.0, 1.0, 1.5, 40.0);
 gluLookAt(0.0, 2.0, 5.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
 glMatrixMode(GL_MODELVIEW); 
 glLoadIdentity();
}

void drawStuff(void)
{
 glPushMatrix();
 glRotatef(rotZ,1.0,0.0,0.0);
 glRotatef(rotY,0.0,1.0,0.0);
 glRotatef(rotX,0.0,0.0,1.0);
 if(!cube) {
  glFrontFace(GL_CW);
  glutSolidTeapot(1.2);
 } else { 
  glFrontFace(GL_CCW);
  drawCube(1.0);
 }
 glPopMatrix();
}

void display(void)
{
 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
 drawStuff();
 glutSwapBuffers();
}

void printhelp(void) 
{
 printf("Texture 3d extension demo keyboard  options:\n");
 printf("  ESC       Exit\n");
 printf("   h        Show help\n"); 
 printf("   s        Smooth shading\n");
 printf("   f        Flat shading\n");
 printf("  y/Y       Rotate around Y\n");
 printf("  x/X       Rotate around X\n");
 printf("  z/Z       Rotate around Z\n"); 
}

void keyboard(unsigned char c, int x, int y) 
{
 switch(c) {
   case 27: exit(0);
            break;
   case 's': 
            glShadeModel(GL_SMOOTH);
            glutPostRedisplay();
            break;
   case 'f':
            glShadeModel(GL_FLAT);
            glutPostRedisplay();
            break;
   case 'h':
            printhelp();
            break;
   case 'x': rotX-=15;
            glutPostRedisplay();
            break;
   case 'X': rotX+=15;
            glutPostRedisplay();
            break;
   case 'y': rotY-=15;
            glutPostRedisplay();
            break;
   case 'Y': rotY+=15;
            glutPostRedisplay();
            break;
   case 'z': rotZ-=15;
            glutPostRedisplay();
            break;
   case 'Z': rotZ+=15;
            glutPostRedisplay();
            break;
   case 'e': exit(0);
            break;
  default: break;
 }
}

void menu( int choice ) 
{
 switch(choice) {
   case 0:   
            glShadeModel(GL_SMOOTH);
            glutPostRedisplay();
        glutChangeToMenuEntry(1,"Smooth --> Flat",1);
        break;
   case 1:   
        glShadeModel(GL_FLAT);
        glutPostRedisplay();
        glutChangeToMenuEntry(1,"Flat ---> Smooth",0);
        break;
   case 2: glDisable(GL_TEXTURE_3D_EXT); 
        glutPostRedisplay();
        glutChangeToMenuEntry(2,"Flat ---> Texture3D",3);
        break;
   case 3: glEnable(GL_TEXTURE_3D_EXT); 
        glutPostRedisplay();
        glutChangeToMenuEntry(2,"Texture 3D --> Flat",2);
        break;
   case 4:
        glDisable(GL_TEXTURE_GEN_S);
        glDisable(GL_TEXTURE_GEN_T);
        glDisable(GL_TEXTURE_GEN_R);
        glutPostRedisplay();
        glutChangeToMenuEntry(3,"Texture coord Fixed --> Gen",5);
        break;
   case 5:
        glEnable(GL_TEXTURE_GEN_S);
        glEnable(GL_TEXTURE_GEN_T);
        glEnable(GL_TEXTURE_GEN_R);
        glutPostRedisplay();
        glutChangeToMenuEntry(3,"Texture coord Gen --> fixed",4);
        break;
   case 6: cube=0;
        glutPostRedisplay();
        glutChangeToMenuEntry(4,"Teapot --> Cube",7);
        break;
   case 7: cube=1;
        glutPostRedisplay();
        glutChangeToMenuEntry(4,"Cube --> Teapot",6);
        break;
   case 20: exit(0);
        break;
   default:
        break;
 }
}

void main( int argc, char **argv)
{
 glutInit(&argc,argv);
 glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);
 glutCreateWindow("3d TextureEXT Test");
 init();
 glutDisplayFunc(display);
 glutReshapeFunc(resize);
 glutKeyboardFunc(keyboard);
 glutCreateMenu(menu);
 glutAddMenuEntry("Smooth --> Flat ",1);
 glutAddMenuEntry("Texture3D --> Flat ",2);
 glutAddMenuEntry("Texture Gen --> Fixed",4);
 glutAddMenuEntry("Cube --> Teapot",6);
 glutAddMenuEntry("Exit",20);
 glutAttachMenu(GLUT_RIGHT_BUTTON);
 printhelp();
 glutMainLoop();
}