TEXTURES.C

Go to the documentation of this file.

/* --
OpenFX Modelling, Animation and Rendering Package
-- */

// file textures.c

#define MODULE_TEXTUR2

/* surface shaders - most code scavanged from elsewhere and adapted     */
/* by converting to the approporate geometry                            */
/* Perlins SIGGRAPH 86 paper textures are very evident                  */
/* The CLOUD texture is called BOZO   see Perlin + DKB trace            */


#include "render.h"

#define MINX    -10000
#define MINY    MINX
#define MINZ    MINX

#define max_noise 10
#define NRIPPLES 10
#define MAXSIZE 267
#define RNDMASK 0x7FFF
#define RNDDIVISOR (double) RNDMASK
#define NUMBER_OF_WAVES 10
#define SINTABSIZE 1000
#define EPSILON   (double) 0.00001

#define FLOOR(x) ((x) >= 0.0 ? floor(x) : (0.0 - floor(0.0 - (x)) - 1.0))
#define FABS(x) ((x) < 0.0 ? (0.0 - x) : (x))
#define SCURVE(a) ((a)*(a)*(3.0-2.0*(a)))
#define REALSCALE ( 2.0 / 65535.0 )
#define Hash3d(a,b,c) hashTable[(int)(hashTable[(int)( \
     hashTable[(int)((a) & 0xfffL)] ^ ((b) & 0xfffL))] ^ ((c) & 0xfffL))]
#define INCRSUM(m,s,x,y,z)   ((s)*(RTable[m]*0.5+RTable[m+1]*(x) \
                                  + RTable[m+2]*(y) + RTable[m+3]*(z)))

vector Wave_Sources[NUMBER_OF_WAVES];
double Rip_Freq=10000.0;

static vector  ripple_o[NRIPPLES];
static double  ripple_l[NRIPPLES],ripple_h[NRIPPLES];
static unsigned short noise_table[max_noise+1][max_noise+1][max_noise+1];

static double sintab [SINTABSIZE];
static double frequency[NUMBER_OF_WAVES];
static double RTable[MAXSIZE];
static short *hashTable=NULL;
static unsigned short crctab[256] =
{
   0x0000, 0xc0c1, 0xc181, 0x0140, 0xc301, 0x03c0, 0x0280, 0xc241,
   0xc601, 0x06c0, 0x0780, 0xc741, 0x0500, 0xc5c1, 0xc481, 0x0440,
   0xcc01, 0x0cc0, 0x0d80, 0xcd41, 0x0f00, 0xcfc1, 0xce81, 0x0e40,
   0x0a00, 0xcac1, 0xcb81, 0x0b40, 0xc901, 0x09c0, 0x0880, 0xc841,
   0xd801, 0x18c0, 0x1980, 0xd941, 0x1b00, 0xdbc1, 0xda81, 0x1a40,
   0x1e00, 0xdec1, 0xdf81, 0x1f40, 0xdd01, 0x1dc0, 0x1c80, 0xdc41,
   0x1400, 0xd4c1, 0xd581, 0x1540, 0xd701, 0x17c0, 0x1680, 0xd641,
   0xd201, 0x12c0, 0x1380, 0xd341, 0x1100, 0xd1c1, 0xd081, 0x1040,
   0xf001, 0x30c0, 0x3180, 0xf141, 0x3300, 0xf3c1, 0xf281, 0x3240,
   0x3600, 0xf6c1, 0xf781, 0x3740, 0xf501, 0x35c0, 0x3480, 0xf441,
   0x3c00, 0xfcc1, 0xfd81, 0x3d40, 0xff01, 0x3fc0, 0x3e80, 0xfe41,
   0xfa01, 0x3ac0, 0x3b80, 0xfb41, 0x3900, 0xf9c1, 0xf881, 0x3840,
   0x2800, 0xe8c1, 0xe981, 0x2940, 0xeb01, 0x2bc0, 0x2a80, 0xea41,
   0xee01, 0x2ec0, 0x2f80, 0xef41, 0x2d00, 0xedc1, 0xec81, 0x2c40,
   0xe401, 0x24c0, 0x2580, 0xe541, 0x2700, 0xe7c1, 0xe681, 0x2640,
   0x2200, 0xe2c1, 0xe381, 0x2340, 0xe101, 0x21c0, 0x2080, 0xe041,
   0xa001, 0x60c0, 0x6180, 0xa141, 0x6300, 0xa3c1, 0xa281, 0x6240,
   0x6600, 0xa6c1, 0xa781, 0x6740, 0xa501, 0x65c0, 0x6480, 0xa441,
   0x6c00, 0xacc1, 0xad81, 0x6d40, 0xaf01, 0x6fc0, 0x6e80, 0xae41,
   0xaa01, 0x6ac0, 0x6b80, 0xab41, 0x6900, 0xa9c1, 0xa881, 0x6840,
   0x7800, 0xb8c1, 0xb981, 0x7940, 0xbb01, 0x7bc0, 0x7a80, 0xba41,
   0xbe01, 0x7ec0, 0x7f80, 0xbf41, 0x7d00, 0xbdc1, 0xbc81, 0x7c40,
   0xb401, 0x74c0, 0x7580, 0xb541, 0x7700, 0xb7c1, 0xb681, 0x7640,
   0x7200, 0xb2c1, 0xb381, 0x7340, 0xb101, 0x71c0, 0x7080, 0xb041,
   0x5000, 0x90c1, 0x9181, 0x5140, 0x9301, 0x53c0, 0x5280, 0x9241,
   0x9601, 0x56c0, 0x5780, 0x9741, 0x5500, 0x95c1, 0x9481, 0x5440,
   0x9c01, 0x5cc0, 0x5d80, 0x9d41, 0x5f00, 0x9fc1, 0x9e81, 0x5e40,
   0x5a00, 0x9ac1, 0x9b81, 0x5b40, 0x9901, 0x59c0, 0x5880, 0x9841,
   0x8801, 0x48c0, 0x4980, 0x8941, 0x4b00, 0x8bc1, 0x8a81, 0x4a40,
   0x4e00, 0x8ec1, 0x8f81, 0x4f40, 0x8d01, 0x4dc0, 0x4c80, 0x8c41,
   0x4400, 0x84c1, 0x8581, 0x4540, 0x8701, 0x47c0, 0x4680, 0x8641,
   0x8201, 0x42c0, 0x4380, 0x8341, 0x4100, 0x81c1, 0x8081, 0x4040
};

static short InitRTable(void);
static short InitTextureTable(void);
static int RTableNoise(vector v);
static int Crc16(char *buf, int count);
static void InitialiseSimpleNoise(void);
static double SimpleNoise(double x, double y, double z);

void FreeTextureHashTable(void){
  if(hashTable != NULL)X__Free(hashTable);
  hashTable=NULL;
  return;
}

short InitialiseNoise(void){
   register int i = 0;
   InitialiseSimpleNoise();
   if(InitRTable() == FAIL)return FAIL;
   for (i = 0 ; i < SINTABSIZE ; i++)
      sintab[i] = sin(i/(double)SINTABSIZE * (3.14159265359 * 2.0));
   for (i = 0 ; i < NUMBER_OF_WAVES ; i++){
     Dnoise2(Wave_Sources[i], (double) i, 0.0, 0.0);
     normalize(Wave_Sources[i]);
     frequency[i] = (rand() & RNDMASK) / RNDDIVISOR + 0.01;
   }
   return OK;
}

static short InitRTable(void){
   int i;
   vector rp;
   if(InitTextureTable() == FAIL)return FAIL;
   for (i = 0; i < MAXSIZE; i++){
     rp[0] = rp[1] = rp[2] = (double)i;
     RTable[i] = (unsigned int) RTableNoise(rp) * REALSCALE - 1.0;
   }
   return OK;
}

static short InitTextureTable(void){
   int i, j, temp;
   srand(0);
   if((hashTable = (short int *) X__Malloc(4096*sizeof(short int))) == NULL)
     return FAIL;
   for (i = 0; i < 4096; i++)
      hashTable[i] = i;
   for (i = 4095; i >= 0; i--)
      {
      j = rand() % 4096;
      temp = hashTable[i];
      hashTable[i] = hashTable[j];
      hashTable[j] = temp;
      }
  return OK;
}

static int RTableNoise(vector v){
   v[0] *= .12345;
   v[1] *= .12345;
   v[2] *= .12345;
   return (Crc16((char *) v, sizeof(vector)));
}

static int Crc16(char *buf, int count){
   register unsigned short crc = 0;
   while (count--)
      crc = (crc >> 8) ^ crctab[ (unsigned char) (crc ^ *buf++) ];
   return ((int) crc);
}

void Snoise2(double x, double y, double z, double *r){
   register long ix, iy, iz, jx, jy, jz;
   double sx, sy, sz, tx, ty, tz;
   double sum;
   short m;
   /* ensures the values are positive. */
   x -= MINX;
   y -= MINY;
   z -= MINZ;
   /* its equivalent integer lattice point. */
   ix = (long)x; iy = (long)y; iz = (long)z;
   jx = ix + 1; jy = iy + 1; jz = iz + 1;
   sx = SCURVE(x - ix); sy = SCURVE(y - iy); sz = SCURVE(z - iz);
   /* the complement values of sx,sy,sz */
   tx = 1.0 - sx; ty = 1.0 - sy; tz = 1.0 - sz;
   /*
    *  interpolate!
    */
   m = Hash3d( ix, iy, iz ) & 0xFF;
   sum = INCRSUM(m,(tx*ty*tz),(x-ix),(y-iy),(z-iz));
   m = Hash3d( jx, iy, iz ) & 0xFF;
   sum += INCRSUM(m,(sx*ty*tz),(x-jx),(y-iy),(z-iz));
   m = Hash3d( ix, jy, iz ) & 0xFF;
   sum += INCRSUM(m,(tx*sy*tz),(x-ix),(y-jy),(z-iz));
   m = Hash3d( jx, jy, iz ) & 0xFF;
   sum += INCRSUM(m,(sx*sy*tz),(x-jx),(y-jy),(z-iz));
   m = Hash3d( ix, iy, jz ) & 0xFF;
   sum += INCRSUM(m,(tx*ty*sz),(x-ix),(y-iy),(z-jz));
   m = Hash3d( jx, iy, jz ) & 0xFF;
   sum += INCRSUM(m,(sx*ty*sz),(x-jx),(y-iy),(z-jz));
   m = Hash3d( ix, jy, jz ) & 0xFF;
   sum += INCRSUM(m,(tx*sy*sz),(x-ix),(y-jy),(z-jz));
   m = Hash3d( jx, jy, jz ) & 0xFF;
   sum += INCRSUM(m,(sx*sy*sz),(x-jx),(y-jy),(z-jz));
   sum = sum + 0.5;          /* range at this point -0.5 - 0.5... */
   if (sum < 0.0)
      sum = 0.0;
   if (sum > 1.0)
      sum = 1.0;
   *r = sum;
   return;
}

void Dnoise2(vector result, double x, double y, double z){
   register long ix, iy, iz, jx, jy, jz;
   double px, py, pz, s;
   double sx, sy, sz, tx, ty, tz;
   short m;
   /* ensures the values are positive. */
   x -= MINX;
   y -= MINY;
   z -= MINZ;
   /* its equivalent integer lattice point. */
   ix = (long)x; iy = (long)y; iz = (long)z;
   jx = ix+1; jy = iy + 1; jz = iz + 1;
   sx = SCURVE(x - ix); sy = SCURVE(y - iy); sz = SCURVE(z - iz);
   /* the complement values of sx,sy,sz */
   tx = 1.0 - sx; ty = 1.0 - sy; tz = 1.0 - sz;
   m = Hash3d( ix, iy, iz ) & 0xFF;
   px = x-ix;  py = y-iy;  pz = z-iz;
   s = tx*ty*tz;
   result[0] = INCRSUM(m,s,px,py,pz);
   result[1] = INCRSUM(m+4,s,px,py,pz);
   result[2] = INCRSUM(m+8,s,px,py,pz);
   m = Hash3d( jx, iy, iz ) & 0xFF;
   px = x-jx;
   s = sx*ty*tz;
   result[0] += INCRSUM(m,s,px,py,pz);
   result[1] += INCRSUM(m+4,s,px,py,pz);
   result[2] += INCRSUM(m+8,s,px,py,pz);
   m = Hash3d( jx, jy, iz ) & 0xFF;
   py = y-jy;
   s = sx*sy*tz;
   result[0] += INCRSUM(m,s,px,py,pz);
   result[1] += INCRSUM(m+4,s,px,py,pz);
   result[2] += INCRSUM(m+8,s,px,py,pz);
   m = Hash3d( ix, jy, iz ) & 0xFF;
   px = x-ix;
   s = tx*sy*tz;
   result[0] += INCRSUM(m,s,px,py,pz);
   result[1] += INCRSUM(m+4,s,px,py,pz);
   result[2] += INCRSUM(m+8,s,px,py,pz);
   m = Hash3d( ix, jy, jz ) & 0xFF;
   pz = z-jz;
   s = tx*sy*sz;
   result[0] += INCRSUM(m,s,px,py,pz);
   result[1] += INCRSUM(m+4,s,px,py,pz);
   result[2] += INCRSUM(m+8,s,px,py,pz);
   m = Hash3d( jx, jy, jz ) & 0xFF;
   px = x-jx;
   s = sx*sy*sz;
   result[0] += INCRSUM(m,s,px,py,pz);
   result[1] += INCRSUM(m+4,s,px,py,pz);
   result[2] += INCRSUM(m+8,s,px,py,pz);
   m = Hash3d( jx, iy, jz ) & 0xFF;
   py = y-iy;
   s = sx*ty*sz;
   result[0] += INCRSUM(m,s,px,py,pz);
   result[1] += INCRSUM(m+4,s,px,py,pz);
   result[2] += INCRSUM(m+8,s,px,py,pz);
   m = Hash3d( ix, iy, jz ) & 0xFF;
   px = x-ix;
   s = tx*ty*sz;
   result[0] += INCRSUM(m,s,px,py,pz);
   result[1] += INCRSUM(m+4,s,px,py,pz);
   result[2] += INCRSUM(m+8,s,px,py,pz);
}

void STurbulence(double x, double y, double z, double *r){
   register double pixelSize = 0.1;
   register double t = 0.0;
   register double scale;
   double   value;
   for (scale = 1.0 ; scale > pixelSize ; scale *= 0.5) {
      Snoise2(x/scale, y/scale, z/scale, &value);
      t += FABS (value) * scale;
      }
   *r=t;
   return;
}

void DTurbulence(vector result, double x, double y, double z){
   register double pixelSize = 0.01;
   register double scale;
   vector value;
   result[0] = 0.0;
   result[1] = 0.0;
   result[2] = 0.0;
   value[0] = value[1] = value[2] = 0.0;
   for (scale = 1.0 ; scale > pixelSize ; scale *= 0.5) {
      Dnoise2(value, x/scale, y/scale, z/scale);
      result[0] += value[0] * scale;
      result[1] += value[1] * scale;
      result[2] += value[2] * scale;
      }
}

void Cycloidal(double value, double *r){
 if (value >= 0.0)
    *r = (sintab [(int)((value - floor (value)) * SINTABSIZE)]);
 else
    *r = (0.0 - sintab [(int)((0.0 - (value + floor (0.0 - value)))
                                  * SINTABSIZE)]);
 return;
}

static void InitialiseSimpleNoise(void){
 long x,y,z,xx,yy,zz;
 for(x=0;x<=max_noise;x++)
 for(y=0;y<=max_noise;y++)
 for(z=0;z<=max_noise;z++){
   noise_table[x][y][z]=(unsigned short)rand();
   if(x == max_noise)xx=0;
   else              xx=x;
   if(y == max_noise)yy=0;
   else              yy=y;
   if(z == max_noise)zz=0;
   else              zz=z;
   noise_table[x][y][z]=noise_table[xx][yy][zz];
 }
}

static double SimpleNoise(double x, double y, double z){
 long ix,iy,iz,n;
 double ox,oy,oz,n00,n01,n10,n11,n0,n1;
 ix=((long)x+30000)%max_noise;
 iy=((long)y+30000)%max_noise;
 iz=((long)z+30000)%max_noise;
 ox=x-floor(x); oy=y-floor(y); oz=z-floor(z);
 n=noise_table[ix][iy][iz];
 n00=n+ox*(noise_table[ix+1][iy][iz]-n);
 n=noise_table[ix][iy][iz+1];
 n01=n+ox*(noise_table[ix+1][iy][iz+1]-n);
 n=noise_table[ix][iy+1][iz];
 n10=n+ox*(noise_table[ix+1][iy+1][iz]-n);
 n=noise_table[ix][iy+1][iz+1];
 n11=n+ox*(noise_table[ix+1][iy+1][iz+1]-n);
 n0=n00+oy*(n10-n00);
 n1=n01+oy*(n11-n01);
 return (n0+oz*(n1-n0))*3.052e-5;
}


double Triangle_Wave(double value){
  register double offset;
  if (value >= 0.0) offset = value - floor(value);
  else offset = value - (-1.0 - floor(FABS(value)));
  if (offset >= 0.5) return (2.0 * (1.0 - offset));
  else return (2.0 * offset);
}

static void bozo_grain(double alpha, double beta, double gamma, double *color, 
               unsigned char *txcol, double turb, double zscale){
   double noise;
   vector BozoTurbulence;
   double x,y,z,fraction;
   x=alpha; y=beta; z=gamma;
   z = z*zscale;
   x=alpha; y=beta; z=gamma;
   if (turb != 0.0)
      {
      DTurbulence (BozoTurbulence, x, y, z);
      x += BozoTurbulence[0] * turb;
      y += BozoTurbulence[1] * turb;
      z += BozoTurbulence[2] * turb;
      }
   Snoise2 (x, y, z, &noise);
   if (noise < 0.5)return;      /* basic colour */
   else if (noise < 0.6) {      /* blend to white */
      fraction=(noise - 0.5)/(0.6 - 0.5);
      color[0] = color[0] + fraction*(255.0 - color[0]);
      color[1] = color[1] + fraction*(255.0 - color[1]);
      color[2] = color[2] + fraction*(255.0 - color[2]);
   }
   else {                       /* blend to grey */
      fraction=1.0 - 0.5*(noise - 0.6)/(1.0 - 0.6);
      color[0] = fraction*255.0;
      color[1] = fraction*255.0;
      color[2] = fraction*255.0;
   }
   return;
}

static void wood_grain2(double alpha, double beta, double gamma,
                double *color,
                unsigned char *txcol, double turb){
 double x,y,z,noise,length;
 vector WoodTurbulence,point;
 x=alpha,y=beta,z=gamma;
 x *= 10.0; y *= 10.0; z *= 10.0;
 DTurbulence (WoodTurbulence,x,y,z);
 Cycloidal((x + WoodTurbulence[0])*turb, &point[0]);
 Cycloidal((y + WoodTurbulence[1])*turb, &point[1]);
 point[2] = 0.0;
 point[0] += x;
 point[1] += y;
 point[2] += z;
 length=R_length(point);
 noise = Triangle_Wave(length);
 if(noise > 0.6){
   color[0] = (double)txcol[0];
   color[1] = (double)txcol[1];
   color[2] = (double)txcol[2];
 }
 return;
}

static void waves2(double alpha, double beta, double gamma,
   vector Vector, double Frequency, double Phase, double Bump_Amount){
   double x,y,z;
   register int i;
   vector point;
   register double length, scalar, index;
   double sinValue ;
   x=alpha; y=beta; z=gamma;
   x *= .030; y *= .030; z=0;
   for (i = 0 ; i < NUMBER_OF_WAVES ; i++) {
      point[0] = x;
      point[1] = y;
      point[2] = z;
      vecsub(point, Wave_Sources[i], point);
      length=dot(point, point);
      if (length == 0.0)length = 1.0;
      length = sqrt(length);
      index = (length * Frequency * frequency[i]) + Phase;
      Cycloidal (index, &sinValue);
      scalar =  sinValue * Bump_Amount / frequency[i];
      vecscale(scalar/length/(double)NUMBER_OF_WAVES,point,point);
      vecsum(Vector, point, Vector);
   }
   normalize(Vector);
}

static void spotted_grain(double alpha, double beta, double gamma,
                   double *color, unsigned char *txcol){
   double x,y,z;
   double noise;
   x=alpha,y=beta; z=gamma;
   Snoise2(x,y,z,&noise);
   color[0] += (noise*255);             /* "white (1.0) * noise" */
   color[1] += (noise*255);
   color[2] += (noise*255);
   return;
}

static void agate_grain(double alpha, double beta, double gamma,
                   double *color, unsigned char *txcol){
   double x,y,z,r;
   double noise, hue;
   x=alpha; y=beta; z=gamma; 
   STurbulence(x,y,z,&r);
   Cycloidal(1.3 * r + 1.1 * z, &noise);
   noise += 1;
   noise *= 0.5;
   noise = pow(noise, 0.77);
   hue = 1.0 - noise;
   if (noise < 0.5)
      {
      color[0] += (1.0 - (noise / 10))*255;
      color[1] += (1.0 - (noise / 5))*255;
      color[2] += hue*255;
      }
   else if (noise < 0.6)
      {
      color[0] += 0.9*255;
      color[1] += 0.7*255;
      color[2] += hue*255;
      }
   else
      {
      color[0] += (0.6 + hue)*255;
      color[1] += (0.3 + hue)*255;
      color[2] += hue*255;
      }
   return;
   }

static double new_bumpy(double alpha, double beta, double gamma,double Bump_Amount, vector normal){
   double x,y,z;
   vector result;
   x=alpha; y=beta; z=gamma; 
   Dnoise2(result,x*10.0,y*10.0,z*10.0);
   Bump_Amount *= 2.0;
   VECSCALE(Bump_Amount,result,result)
   vecsum(normal,result,normal);             /* displace "normal" */
   normalize(normal);                        /* normalize normal! */
   return 0.0;
}

static double dents(double alpha, double beta, double gamma, double Bump_Amount, vector normal){
   double x,y,z,noise;
   vector result;
   x=alpha; y=beta; z=gamma;
   x *= 10.0; y *= 10.0; z *= 10.0;
   Snoise2(x,y,z,&noise);
   noise=noise*noise*noise*Bump_Amount*2.0;
   Dnoise2(result,x,y,z);
   vecscale(noise,result,result);
   vecsum(normal,result,normal);             /* displace "normal" */
   normalize(normal);                        /* normalize normal! */
   return 0.0;
}

static double wrinkles(double alpha, double beta, double gamma, double Bump_Amount, vector normal){
   double x,y,z;
   register int i;
   register double scale = 1.0;
   vector result, value;
   x=alpha; y=beta; z=gamma;
   result[0] = 0.0;
   result[1] = 0.0;
   result[2] = 0.0;
   x *= 10.0; y *= 10.0; z *= 10.0;
   for (i = 0; i < 10 ; scale *= 2.0, i++){
     Dnoise2(value, x * scale, y * scale, z * scale);   /* * scale,*/
     result[0] += FABS (value[0] / scale);
     result[1] += FABS (value[1] / scale);
     result[2] += FABS (value[2] / scale);
   }
   vecscale(Bump_Amount*2.0,result,result);
   vecsum(normal,result,normal);             /* displace "normal" */
   normalize(normal);                        /* normalize normal! */
   return 0.0;
}

/*
   Granite - union of the "spotted" and the "dented" textures,
   using a 1/f fractal noise function for color values.
*/
static void granite_grain(double alpha, double beta, double gamma,
                   double *color, unsigned char *txcol){
   double x,y,z,temp;
   register int i;
   register double noise = 0.0, freq = 1.0;
   x=alpha; y=beta; z=gamma;
   for (i = 0; i < 6 ; freq *= 2.0, i++)   {
      Snoise2 (x * 4 * freq, y * 4 * freq, z * 4 * freq, &temp);
      temp = 0.5 - temp;
      temp = FABS(temp);
      noise += temp / freq;
   }
   color[0] *= noise;                  /* "white (1.0) * noise" */
   color[1] *= noise;
   color[2] *= noise;
   return;
}

static void marble_grain(double alpha, double beta, double gamma,
                  double *color, unsigned char *txcol){
 vector noisevec;
 double i,sinval;
 alpha *= 20;
 beta  *= 20;
 gamma *= 10;
 Dnoise2(noisevec,alpha,beta,gamma);
 i = fabs(noisevec[0]*0.5  - 0.25) +
     fabs(noisevec[1]*0.5  - 0.125) +
     fabs(noisevec[2]*0.5  - 0.0626);
 if(i < 0.0)i=0.0;
 if(i > 0.999)i=0.999;
 i=sin(gamma*1.0+i*2.0)*0.5 + 0.5;
#if 0
 double d,dd,i;
 alpha *= 20;
 beta  *= 20;
 gamma *= 10;
 d=30000+beta*5+17*SimpleNoise(alpha,beta,gamma);
 dd=(double)(((long)floor(d))%17);
/* if(dd < 4)i=0.7+0.2*SimpleNoise(alpha,beta,gamma);    original */
 if(dd < 2)i=0.7+0.2*SimpleNoise(alpha,beta,gamma);
/* else if(dd < 9 || dd >= 12){                    original */
 else if(dd < 5 || dd >= 10){
   d=fabs(d-floor(d/17)*17-10.5)*0.1538462;
   i=0.4+0.3*d+0.2*SimpleNoise(alpha,beta,gamma);
 }
 else i=0.2+0.2*SimpleNoise(alpha,beta,gamma);
#endif
 color[0]=(double)txcol[0]*i+color[0]*(1.0-i);
 color[1]=(double)txcol[1]*i+color[1]*(1.0-i);
 color[2]=(double)txcol[2]*i+color[2]*(1.0-i);
}

static void fire_grain(double alpha, double beta, double gamma,
                  double *color, unsigned char *txcol){
 vector noisevec;
 double i;
 alpha *= 10;
 beta  *= 10;
 gamma *= 10;
 Dnoise2(noisevec,alpha,beta,gamma);
 i = fabs(noisevec[0]*0.5  - 0.25) +
     fabs(noisevec[1]*0.5  - 0.125) +
     fabs(noisevec[2]*0.5  - 0.0626);
 //i=fabs(noisevec[0]*0.5 - 0.25);
 if(i < 0.0)i=0.0;
 if(i > 0.999)i=0.999;
 color[0]=(double)txcol[0]*i+color[0]*(1.0-i);
 color[1]=(double)txcol[1]*i+color[1]*(1.0-i);
 color[2]=(double)txcol[2]*i+color[2]*(1.0-i);
}


static double aztec_bumpy(double a1, double b1, double g1, double bumpyness, vector n){
 double alpha,beta,gamma,i1,i2,i3,c;
 short j1,j2,j3;
 vector v1,v2;
 vector unitx={1.0,0.0,0.0};
 vector unity={0.0,1.0,0.0};
 gamma =g1;
 c = 2.0*bumpyness;
 alpha  = a1*10;
 beta   = b1*10;
 gamma *= 30;
 j1=SimpleNoise(alpha,beta,gamma)*10;
 j2=SimpleNoise(alpha+0.15,beta,gamma)*10;
 j3=SimpleNoise(alpha,beta+0.15,gamma)*10;
 if(j1%2)i1=0; else i1=0.3;
 if(j2%2)i2=0; else i2=0.3;
 if(j3%2)i3=0; else i3=0.3;
 vecscale((i2-i1)*c,unitx,v1);
 vecscale((i3-i1)*c,unity,v2);
 vecsum(n,v1,n);
 vecsum(n,v2,n);
 normalize(n);
 return 0.0;
}

static void blend_grain(double alpha, double beta, double gamma, 
                 double *color,unsigned char *txcol){
 gamma=(gamma-floor(gamma))*2;
 if(gamma <= 1.0){
   color[0]=color[0] + ((double)txcol[0] - color[0])*gamma;
   color[1]=color[1] + ((double)txcol[1] - color[1])*gamma;
   color[2]=color[2] + ((double)txcol[2] - color[2])*gamma;
 }
 else{
   gamma -= 1.0;
   color[0]=(double)txcol[0] + (color[0] - (double)txcol[0])*gamma;
   color[1]=(double)txcol[1] + (color[1] - (double)txcol[1])*gamma;
   color[2]=(double)txcol[2] + (color[2] - (double)txcol[2])*gamma;
 }
}

void init_ripples(void){
  vector v1,v2,vy;
  double rn,rn1,rn2;
  short j;
  Rip_Freq=10000.0;
  veccopy(Ground.dy,vy);
  normalize(vy);
  vecscale(Ground.lengthx,vy,vy); /* x vector length in y direction */
  vecsum(Ground.p,Ground.dx,ripple_o[0]);  ripple_l[0]=Ground.lengthy;
  ripple_h[0]=Ground.height*0.4;  /* 0.4  is a good height */
  vecsum(Ground.p,vy,ripple_o[1]);  ripple_l[1]=Ground.lengthy;
  ripple_h[1]=Ground.height*0.4;
  if(NRIPPLES > 2)for(j=2;j<NRIPPLES;j++){
    rn1 = (double)(rand())/32767.0 * 10;
    rn2 = (double)(rand())/32767.0 * 10;
    vecscale(rn1,Ground.dx,v1);
    vecscale(rn2,vy,v2);
    vecsum(Ground.p,v1,v1);
    vecsum(v1,v2,ripple_o[j]);
    rn = 0.01 + (double)(rand())/32767.0;
    ripple_l[j]=Ground.lengthy*rn;
    ripple_h[j]=Ground.height*0.4;
  }
}

void calcripple(vector point, vector wave_centre,
                double wave_propagate,double wave_wavelength,
                double wave_drag, double wave_amplitude,
                vector ripple){
    double riprad,temp1,temp2,temp3;
    long iriprad;
    vecsub(point,wave_centre,ripple);
    normalize(ripple); /* unit vector pointing away from wave centre */
    temp1 = wave_centre[0] - point[0];
    temp2 = wave_centre[1] - point[1];
    temp3 = wave_centre[2] - point[2];

    temp1 = temp1 * temp1;  /* square it */
    temp2 = temp2 * temp2;  /* square it */
    temp3 = temp3 * temp3;  /* square it */

    riprad = (double)sqrt(temp1 + temp2 + temp3);  /* distance from centre to point */
/* distance damping */
/*    if (wave_drag < 1.01) damper = 1.0;
**    else damper = min(1.0,wave_drag/riprad);
*/
    riprad += wave_propagate * wave_wavelength;  /* move it */
    riprad /= wave_wavelength;  /* scale half-wave crest to 0..1 */

    iriprad = (long) riprad;
    if (iriprad & 1)  {/* negate it */
      ripple[0] = -ripple[0];
      ripple[1] = -ripple[1];
      ripple[2] = -ripple[2];
      }
    riprad -= (double) iriprad;  /* just get fraction 0..1 */
    riprad -= 0.5;  /* scale to -0.5 .. +0.5 */
    if (riprad < 0.0) riprad = -riprad;  /* absolute value */
    riprad  = 0.5 - riprad;  /* invert */
/*    riprad *= damper; */
    riprad *= wave_amplitude;
    vecscale(riprad,ripple,ripple);  /* scale bend */
}

void Ocean(vector p0, vector Groundn, vector gn){
 int j;
 vector ripple;
 veccopy(Groundn,gn);
 for(j=0;j<NRIPPLES;j++){
   calcripple(p0,ripple_o[j],wave_phase,ripple_l[j],
             1.0,ripple_h[j],ripple);
   vecsum(ripple,gn,gn);
 }
 if(normalize(gn))veccopy(Groundn,gn);
}


BOOL GetInternalTexture(matl *Mat,double alpha, double beta, double gamma,
                        long textID, vector p, vector n, long obj, face *f,
                        double *colors, double *colort,
                        long *flag_trans, double tr, double *tr1,
                        long trace_relection_depth,long trace_refraction_depth,
                        long refractive_flag){
 double c;
 vector d;
 long t;
 unsigned char tExco[3]={128,128,0}; 
 unsigned char fParam=128;
 BOOL fullsat=FALSE;
 VECCOPY(Mat->texco,tExco)
 fParam=Mat->params[0];
 if(textID == 1){/* wood finish */
   wood_grain2(alpha,beta,gamma,colors,tExco,0.078125);
 }
 else if(textID ==  7)new_bumpy  (alpha,beta,gamma,(double)fParam*SHSCALE,n);
 else if(textID ==  8)aztec_bumpy(alpha,beta,gamma,(double)fParam*SHSCALE,n);
 else if(textID == 13)dents      (alpha,beta,gamma,(double)fParam*SHSCALE,n);
 else if(textID ==  9)wrinkles   (alpha,beta,gamma,(double)fParam*SHSCALE,n);
 else if(textID == 2){/* marble finish */
   marble_grain(alpha,beta,gamma,colors,tExco);
 }
 else if(textID == 12){ /* cloudy  */
   bozo_grain(alpha,beta,gamma,colors,tExco,
             (double)(fParam)*7.8125e-3,
//             (double)(f->matsp)*0.375);
             (double)(128)*0.375);
   fullsat=TRUE;
 }
 else if(textID == 3){ /* granite  */
   granite_grain(alpha,beta,gamma,colors,tExco);
 }
 if(textID == 4){   /* BRICK */  
   alpha=(alpha-floor(alpha));
   gamma=(gamma-floor(gamma));
   beta =(beta-floor(beta));
   if(beta > 0.95){
     colors[0]=(double)tExco[0];
     colors[1]=(double)tExco[1];
     colors[2]=(double)tExco[2];
   }
   else if((gamma > 0.95 || (gamma > 0.45 && gamma < 0.5)) ||
     (gamma < 0.5 && ( (alpha > 0.45 && alpha < 0.5) || alpha > 0.95) ) ||
     (gamma > 0.5 && ((alpha > 0.25 && alpha < 0.3) ||
     (alpha > 0.7 && alpha < 0.75)))
   ){
     colors[0]=(double)tExco[0];
     colors[1]=(double)tExco[1];
     colors[2]=(double)tExco[2];
   }
 }
 else if(textID == 5){ /* luminous */
   fullsat=TRUE;
 }
 else if(textID == 6){ /* FIRE */ 
   fire_grain(alpha,beta,gamma,colors,tExco);
   fullsat=TRUE;
 }
 else if(textID == 10){ /* spotted  */
   spotted_grain(alpha,beta,gamma,colors,tExco);
 }
 else if(textID == 11){ /* agate  */
   agate_grain(alpha,beta,gamma,colors,tExco);
 }
 else if(textID == 14){/* blend finish */
   blend_grain(alpha,beta,gamma,colors,tExco);
 }
#if 0
 else if(trace_refractions && textID == 15){/* refractive glass */  // THESE last two have been
   c=(double)((p[0]*p[0]) + (p[1]*p[1]) + (p[2]*p[2]));             // temporarily removed
   c=sqrt(c);
   VECSCALE((1.0/c),p,d)
   trace_REFR_ray(p,d,n,obj,f,(double)BIGP,colort,
                        trace_relection_depth,trace_refraction_depth,
                        refractive_flag);
   *flag_trans=1;// tr=1.0;
   *tr1=1.0-tr;
 }
 else if(trace_reflections && textID == 16){
  waves2(alpha,beta,gamma,n,Rip_Freq,wave_phase,         /* 0.2 is max   */
                (double)fParam * 0.000390625);/* 0.02 is good */
   c=(double)((p[0]*p[0]) + (p[1]*p[1]) + (p[2]*p[2]));
   c=sqrt(c);
   VECSCALE((1.0/c),p,d)
   trace_REFR_ray(p,d,n,obj,f,(double)BIGP,colort,
                        trace_relection_depth,trace_refraction_depth,
                        refractive_flag);

   *flag_trans=1;// tr=1.0;
   *tr1=1.0-tr;
 }
#endif 
 return fullsat;
}