home *** CD-ROM | disk | FTP | other *** search
/ rtsi.com / 2014.01.www.rtsi.com.tar / www.rtsi.com / OS9 / OSK / GRAPHICS / rayshade.lzh / noise.c < prev    next >
Text File  |  1990-09-25  |  10KB  |  387 lines
  1. /*
  2.  * noise.c
  3.  *
  4.  * Copyright (C) 1989, Robert Skinner, Craig E. Kolb, F. Kenton Musgrave
  5.  *
  6.  * This software may be freely copied, modified, and redistributed,
  7.  * provided that this copyright notice is preserved on all copies.
  8.  *
  9.  * There is no warranty or other guarantee of fitness for this software,
  10.  * it is provided solely "as is".  Bug reports or fixes may be sent
  11.  * to the author, who may or may not act on them as he desires.
  12.  *
  13.  * You may not include this software in a program or other software product
  14.  * without supplying the source, or without informing the end-user that the
  15.  * source is available for no extra charge.
  16.  *
  17.  * If you modify this software, you should include a notice giving the
  18.  * name of the person performing the modification, the date of modification,
  19.  * and the reason for such modification.
  20.  *
  21.  * $Id: noise.c,v 3.0.1.1 90/04/04 18:59:22 craig Exp $
  22.  *
  23.  * $Log:    noise.c,v $
  24.  * Revision 3.0.1.1  90/04/04  18:59:22  craig
  25.  * patch5: Lint removal.
  26.  * 
  27.  * Revision 3.0  89/10/27  02:05:57  craig
  28.  * Baseline for first official release.
  29.  * 
  30.  */
  31. #include <stdio.h>
  32. #include <math.h>
  33. #include "constants.h"
  34. #include "typedefs.h"
  35. #include "funcdefs.h"
  36.  
  37. #define MINX        -10000
  38. #define MINY        MINX
  39. #define MINZ        MINX
  40.  
  41. #define SCURVE(a) ((a)*(a)*(3.0-2.0*(a)))
  42. #define REALSCALE ( 2.0 / 65536.0 )
  43. #define NREALSCALE ( 2.0 / 4096.0 )
  44. #define Hash3d(a,b,c) hashTable[hashTable[hashTable[(a) & 0xfff] ^ ((b) & 0xfff)] ^ ((c) & 0xfff)]
  45. #define Hash(a,b,c) (xtab[(xtab[(xtab[(a) & 0xff] ^ (b)) & 0xff] ^ (c)) & 0xff] & 0xff)
  46.  
  47. #define INCRSUM(m,s,x,y,z)    ((s)*(RTable[m]*0.5        \
  48.                     + RTable[m+1]*(x)    \
  49.                     + RTable[m+2]*(y)    \
  50.                     + RTable[m+3]*(z)))    \
  51.  
  52.  
  53. #define MAXSIZE 267
  54.  
  55. double        RTable[MAXSIZE];
  56. static short    *hashTable;
  57.  
  58. static unsigned short xtab[256] =
  59. {
  60.    0x0000, 0xc0c1, 0xc181, 0x0140, 0xc301, 0x03c0, 0x0280, 0xc241,
  61.    0xc601, 0x06c0, 0x0780, 0xc741, 0x0500, 0xc5c1, 0xc481, 0x0440,
  62.    0xcc01, 0x0cc0, 0x0d80, 0xcd41, 0x0f00, 0xcfc1, 0xce81, 0x0e40,
  63.    0x0a00, 0xcac1, 0xcb81, 0x0b40, 0xc901, 0x09c0, 0x0880, 0xc841,
  64.    0xd801, 0x18c0, 0x1980, 0xd941, 0x1b00, 0xdbc1, 0xda81, 0x1a40,
  65.    0x1e00, 0xdec1, 0xdf81, 0x1f40, 0xdd01, 0x1dc0, 0x1c80, 0xdc41,
  66.    0x1400, 0xd4c1, 0xd581, 0x1540, 0xd701, 0x17c0, 0x1680, 0xd641,
  67.    0xd201, 0x12c0, 0x1380, 0xd341, 0x1100, 0xd1c1, 0xd081, 0x1040,
  68.    0xf001, 0x30c0, 0x3180, 0xf141, 0x3300, 0xf3c1, 0xf281, 0x3240,
  69.    0x3600, 0xf6c1, 0xf781, 0x3740, 0xf501, 0x35c0, 0x3480, 0xf441,
  70.    0x3c00, 0xfcc1, 0xfd81, 0x3d40, 0xff01, 0x3fc0, 0x3e80, 0xfe41,
  71.    0xfa01, 0x3ac0, 0x3b80, 0xfb41, 0x3900, 0xf9c1, 0xf881, 0x3840,
  72.    0x2800, 0xe8c1, 0xe981, 0x2940, 0xeb01, 0x2bc0, 0x2a80, 0xea41,
  73.    0xee01, 0x2ec0, 0x2f80, 0xef41, 0x2d00, 0xedc1, 0xec81, 0x2c40,
  74.    0xe401, 0x24c0, 0x2580, 0xe541, 0x2700, 0xe7c1, 0xe681, 0x2640,
  75.    0x2200, 0xe2c1, 0xe381, 0x2340, 0xe101, 0x21c0, 0x2080, 0xe041,
  76.    0xa001, 0x60c0, 0x6180, 0xa141, 0x6300, 0xa3c1, 0xa281, 0x6240,
  77.    0x6600, 0xa6c1, 0xa781, 0x6740, 0xa501, 0x65c0, 0x6480, 0xa441,
  78.    0x6c00, 0xacc1, 0xad81, 0x6d40, 0xaf01, 0x6fc0, 0x6e80, 0xae41,
  79.    0xaa01, 0x6ac0, 0x6b80, 0xab41, 0x6900, 0xa9c1, 0xa881, 0x6840,
  80.    0x7800, 0xb8c1, 0xb981, 0x7940, 0xbb01, 0x7bc0, 0x7a80, 0xba41,
  81.    0xbe01, 0x7ec0, 0x7f80, 0xbf41, 0x7d00, 0xbdc1, 0xbc81, 0x7c40,
  82.    0xb401, 0x74c0, 0x7580, 0xb541, 0x7700, 0xb7c1, 0xb681, 0x7640,
  83.    0x7200, 0xb2c1, 0xb381, 0x7340, 0xb101, 0x71c0, 0x7080, 0xb041,
  84.    0x5000, 0x90c1, 0x9181, 0x5140, 0x9301, 0x53c0, 0x5280, 0x9241,
  85.    0x9601, 0x56c0, 0x5780, 0x9741, 0x5500, 0x95c1, 0x9481, 0x5440,
  86.    0x9c01, 0x5cc0, 0x5d80, 0x9d41, 0x5f00, 0x9fc1, 0x9e81, 0x5e40,
  87.    0x5a00, 0x9ac1, 0x9b81, 0x5b40, 0x9901, 0x59c0, 0x5880, 0x9841,
  88.    0x8801, 0x48c0, 0x4980, 0x8941, 0x4b00, 0x8bc1, 0x8a81, 0x4a40,
  89.    0x4e00, 0x8ec1, 0x8f81, 0x4f40, 0x8d01, 0x4dc0, 0x4c80, 0x8c41,
  90.    0x4400, 0x84c1, 0x8581, 0x4540, 0x8701, 0x47c0, 0x4680, 0x8641,
  91.    0x8201, 0x42c0, 0x4380, 0x8341, 0x4100, 0x81c1, 0x8081, 0x4040
  92. };
  93.  
  94. double Chaos(), Marble();
  95.  
  96. InitTextureTable()
  97. {
  98.     int i, j, temp;
  99.  
  100. #ifdef SYSV
  101.     (void)srand48(0);
  102. #else
  103.     (void)srandom(0);
  104. #endif
  105.  
  106.     hashTable = (short *) Malloc(4096*sizeof(short int));
  107.     for (i = 0; i < 4096; i++)
  108.         hashTable[i] = i;
  109.     for (i = 4095; i > 0; i--) {
  110.         j = (int)(nrand() * 4096);
  111.         temp = hashTable[i];
  112.         hashTable[i] = hashTable[j];
  113.         hashTable[j] = temp;
  114.     }
  115. }
  116.  
  117.  
  118. InitRTable()
  119. {
  120.     int i;
  121.     Vector rp;
  122.  
  123.     InitTextureTable();
  124.  
  125.     for (i = 0; i < MAXSIZE; i++) {
  126.            rp.x = rp.y = rp.z = (double)i;
  127.            RTable[i] = R(&rp)*REALSCALE - 1.0;
  128.      }
  129. }
  130.  
  131.  
  132. R(v)
  133. Vector *v;
  134. {
  135.     v->x *= .12345;
  136.     v->y *= .12345;
  137.     v->z *= .12345;
  138.  
  139.     return Crc16(v, sizeof(Vector));
  140. }
  141.  
  142. /*
  143.  * Note that passing a double to Crc16 and interpreting it as
  144.  * an array of chars means that machines with different floating-point
  145.  * representation schemes will evaluate Noise(point) differently.
  146.  */
  147. int
  148. Crc16(buf, count)
  149. register unsigned char *buf;
  150. register int  count;
  151. {
  152.     register unsigned int crc = 0;
  153.  
  154.     while (count--)
  155.         crc = (crc >> 8) ^ xtab[ (crc ^ *buf++) & 255];
  156.  
  157.     return crc;
  158. }
  159.  
  160.  
  161. /*
  162.  * Robert's Skinner's Perlin-style "Noise" function
  163.  */
  164. double
  165. Noise(point)
  166. Vector *point;
  167. {
  168.     register int    ix, iy, iz, jx, jy, jz;
  169.     double        x, y, z;
  170.     double    sx, sy, sz, tx, ty, tz;
  171.     double    sum;
  172.     short    m;
  173.  
  174.  
  175.     /* ensures the values are positive. */
  176.     x = point->x - MINX; y = point->y - MINY; z = point->z - MINZ;
  177.  
  178.     /* its equivalent integer lattice point. */
  179.     ix = (int)x; iy = (int)y; iz = (int)z;
  180.     jx = ix+1; jy = iy + 1; jz = iz + 1;
  181.  
  182.     sx = SCURVE(x - ix); sy = SCURVE(y - iy); sz = SCURVE(z - iz);
  183.  
  184.     /* the complement values of sx,sy,sz */
  185.     tx = 1.0 - sx; ty = 1.0 - sy; tz = 1.0 - sz;
  186.  
  187.     /*
  188.      *  interpolate!
  189.      */
  190.     m = Hash3d( ix, iy, iz ) & 0xFF;
  191.     sum = INCRSUM(m,(tx*ty*tz),(x-ix),(y-iy),(z-iz));
  192.  
  193.     m = Hash3d( jx, iy, iz ) & 0xFF;
  194.     sum += INCRSUM(m,(sx*ty*tz),(x-jx),(y-iy),(z-iz));
  195.  
  196.     m = Hash3d( ix, jy, iz ) & 0xFF;
  197.     sum += INCRSUM(m,(tx*sy*tz),(x-ix),(y-jy),(z-iz));
  198.  
  199.     m = Hash3d( jx, jy, iz ) & 0xFF;
  200.     sum += INCRSUM(m,(sx*sy*tz),(x-jx),(y-jy),(z-iz));
  201.  
  202.     m = Hash3d( ix, iy, jz ) & 0xFF;
  203.     sum += INCRSUM(m,(tx*ty*sz),(x-ix),(y-iy),(z-jz));
  204.  
  205.     m = Hash3d( jx, iy, jz ) & 0xFF;
  206.     sum += INCRSUM(m,(sx*ty*sz),(x-jx),(y-iy),(z-jz));
  207.  
  208.     m = Hash3d( ix, jy, jz ) & 0xFF;
  209.     sum += INCRSUM(m,(tx*sy*sz),(x-ix),(y-jy),(z-jz));
  210.  
  211.     m = Hash3d( jx, jy, jz ) & 0xFF;
  212.     sum += INCRSUM(m,(sx*sy*sz),(x-jx),(y-jy),(z-jz));
  213.  
  214.     return sum;
  215.  
  216. } /* Noise() */
  217.  
  218. /*
  219.  * Vector-valued "Noise"
  220.  */
  221. DNoise(point, result)
  222. Vector *point, *result;
  223. {
  224.     register int    ix, iy, iz, jx, jy, jz;
  225.     double        x, y, z;
  226.     double px, py, pz, s;
  227.     double    sx, sy, sz, tx, ty, tz;
  228.     short    m;
  229.  
  230.     /* ensures the values are positive. */
  231.     x = point->x - MINX; y = point->y - MINY; z = point->z - MINZ;
  232.  
  233.     /* its equivalent integer lattice point. */
  234.     ix = (int)x; iy = (int)y; iz = (int)z;
  235.     jx = ix+1; jy = iy + 1; jz = iz + 1;
  236.  
  237.     sx = SCURVE(x - ix); sy = SCURVE(y - iy); sz = SCURVE(z - iz);
  238.  
  239.     /* the complement values of sx,sy,sz */
  240.     tx = 1.0 - sx; ty = 1.0 - sy; tz = 1.0 - sz;
  241.  
  242.     /*
  243.      *  interpolate!
  244.      */
  245.     m = Hash3d( ix, iy, iz ) & 0xFF;
  246.     px = x-ix;  py = y-iy;  pz = z-iz;
  247.     s = tx*ty*tz;
  248.     result->x = INCRSUM(m,s,px,py,pz);
  249.     result->y = INCRSUM(m+4,s,px,py,pz);
  250.     result->z = INCRSUM(m+8,s,px,py,pz);
  251.  
  252.     m = Hash3d( jx, iy, iz ) & 0xFF;
  253.     px = x-jx;
  254.     s = sx*ty*tz;
  255.     result->x += INCRSUM(m,s,px,py,pz);
  256.     result->y += INCRSUM(m+4,s,px,py,pz);
  257.     result->z += INCRSUM(m+8,s,px,py,pz);
  258.  
  259.     m = Hash3d( jx, jy, iz ) & 0xFF;
  260.     py = y-jy;
  261.     s = sx*sy*tz;
  262.     result->x += INCRSUM(m,s,px,py,pz);
  263.     result->y += INCRSUM(m+4,s,px,py,pz);
  264.     result->z += INCRSUM(m+8,s,px,py,pz);
  265.  
  266.     m = Hash3d( ix, jy, iz ) & 0xFF;
  267.     px = x-ix;
  268.     s = tx*sy*tz;
  269.     result->x += INCRSUM(m,s,px,py,pz);
  270.     result->y += INCRSUM(m+4,s,px,py,pz);
  271.     result->z += INCRSUM(m+8,s,px,py,pz);
  272.  
  273.     m = Hash3d( ix, jy, jz ) & 0xFF;
  274.     pz = z-jz;
  275.     s = tx*sy*sz;
  276.     result->x += INCRSUM(m,s,px,py,pz);
  277.     result->y += INCRSUM(m+4,s,px,py,pz);
  278.     result->z += INCRSUM(m+8,s,px,py,pz);
  279.  
  280.     m = Hash3d( jx, jy, jz ) & 0xFF;
  281.     px = x-jx;
  282.     s = sx*sy*sz;
  283.     result->x += INCRSUM(m,s,px,py,pz);
  284.     result->y += INCRSUM(m+4,s,px,py,pz);
  285.     result->z += INCRSUM(m+8,s,px,py,pz);
  286.  
  287.     m = Hash3d( jx, iy, jz ) & 0xFF;
  288.     py = y-iy;
  289.     s = sx*ty*sz;
  290.     result->x += INCRSUM(m,s,px,py,pz);
  291.     result->y += INCRSUM(m+4,s,px,py,pz);
  292.     result->z += INCRSUM(m+8,s,px,py,pz);
  293.  
  294.     m = Hash3d( ix, iy, jz ) & 0xFF;
  295.     px = x-ix;
  296.     s = tx*ty*sz;
  297.     result->x += INCRSUM(m,s,px,py,pz);
  298.     result->y += INCRSUM(m+4,s,px,py,pz);
  299.     result->z += INCRSUM(m+8,s,px,py,pz);
  300. }
  301.  
  302. double
  303. Marble(vec)
  304. Vector *vec;
  305. {
  306.     double i;
  307.  
  308.     i = sin(8. * Chaos(vec, 6) + 7. * vec->z) + 1;
  309.     i *= 0.5;
  310.     i = pow(i, 0.77);
  311.     return i;
  312. }
  313.  
  314. double
  315. Chaos(vec, octaves)
  316. Vector *vec;
  317. int octaves;
  318. {
  319.     double f, s, t;
  320.     int n;
  321.     Vector tp;
  322.  
  323.     s = f = 1.0;
  324.     t = 0.;
  325.  
  326.     for (n = 0; n < octaves; n++) {
  327.         tp.x = f * vec->x;
  328.         tp.y = f * vec->y;
  329.         tp.z = f * vec->z;
  330.         t += Noise(&tp) * s;
  331.         f *= 2.0;
  332.         s *= 0.5;
  333.     }
  334.  
  335.     return t;
  336. }
  337.  
  338. VfBm(vec, omega, lambda, octaves, ans)
  339. Vector *vec, *ans;
  340. double omega, lambda;
  341. int octaves;
  342. {
  343.     register int i;
  344.     double l, o;
  345.     Vector tp, n;
  346.  
  347.     ans->x = ans->y = ans->z = 0.;
  348.  
  349.     l = o = 1.;
  350.     for (i = 0; i < octaves; i++) {
  351.         tp.x = l * vec->x;
  352.         tp.y = l * vec->y;
  353.         tp.z = l * vec->z;
  354.         DNoise(&tp, &n);
  355.         ans->x += o * n.x;
  356.         ans->y += o * n.y;
  357.         ans->z += o * n.z;
  358.         l *= lambda;
  359.         o *= omega;
  360.         if (o < EPSILON)
  361.             break;
  362.     }
  363. }
  364.  
  365. double
  366. fBm(vec, omega, lambda, octaves)
  367. register Vector *vec;
  368. double omega, lambda;
  369. int octaves;
  370. {
  371.     register int i;
  372.     double l, n, a, o;
  373.     Vector tp;
  374.  
  375.     a = 0; l = o = 1.;
  376.     for (i = 0; i < octaves; i++) {
  377.         tp.x = l * vec->x;
  378.         tp.y = l * vec->y;
  379.         tp.z = l * vec->z;
  380.         n = o * Noise(&tp);
  381.         a += n;
  382.         l *= lambda;
  383.         o *= omega;
  384.     }
  385.     return a;
  386. }
  387.