home *** CD-ROM | disk | FTP | other *** search
/ Photo CD Demo 1 / Demo.bin / gems / gemsii / noise3.c < prev    next >
C/C++ Source or Header  |  1991-09-22  |  5KB  |  229 lines
  1.  
  2. /* Copyright (c) 1988 Regents of the University of California */
  3.  
  4. #ifndef lint
  5. static char SCCSid[] = "@(#)noise3.c 1.6 8/2/91 LBL";
  6. #endif
  7.  
  8. /*
  9.  *  noise3.c - noise functions for random textures.
  10.  *
  11.  *     Credit for the smooth algorithm goes to Ken Perlin.
  12.  *     (ref. SIGGRAPH Vol 19, No 3, pp 287-96)
  13.  *
  14.  *     4/15/86
  15.  *     5/19/88    Added fractal noise function
  16.  */
  17.  
  18.  
  19. #define  A        0
  20. #define  B        1
  21. #define  C        2
  22. #define  D        3
  23.  
  24. #define  rand3a(x,y,z)    frand(67*(x)+59*(y)+71*(z))
  25. #define  rand3b(x,y,z)    frand(73*(x)+79*(y)+83*(z))
  26. #define  rand3c(x,y,z)    frand(89*(x)+97*(y)+101*(z))
  27. #define  rand3d(x,y,z)    frand(103*(x)+107*(y)+109*(z))
  28.  
  29. #define  hermite(p0,p1,r0,r1,t)  (    p0*((2.0*t-3.0)*t*t+1.0) + \
  30.                     p1*(-2.0*t+3.0)*t*t + \
  31.                     r0*((t-2.0)*t+1.0)*t + \
  32.                     r1*(t-1.0)*t*t )
  33.  
  34. static char  noise_name[4][8] = {"noise3a", "noise3b", "noise3c", "noise3"};
  35. static char  fnoise_name[] = "fnoise3";
  36. static char  hermite_name[] = "hermite";
  37.  
  38. double  *noise3(), fnoise3(), argument(), frand();
  39.  
  40. static long  xlim[3][2];
  41. static double  xarg[3];
  42.  
  43. #define  EPSILON    .0001        /* error allowed in fractal */
  44.  
  45. #define  frand3(x,y,z)    frand(17*(x)+23*(y)+29*(z))
  46.  
  47.  
  48. static double
  49. l_noise3(nam)            /* compute a noise function */
  50. register char  *nam;
  51. {
  52.     register int  i;
  53.     double  x[3];
  54.                     /* get point */
  55.     x[0] = argument(1);
  56.     x[1] = argument(2);
  57.     x[2] = argument(3);
  58.                     /* make appropriate call */
  59.     if (nam == fnoise_name)
  60.         return(fnoise3(x));
  61.     i = 4;
  62.     while (i--)
  63.         if (nam == noise_name[i])
  64.             return(noise3(x)[i]);
  65.     eputs(nam);
  66.     eputs(": called l_noise3!\n");
  67.     quit(1);
  68. }
  69.  
  70.  
  71. double
  72. l_hermite()            /* library call for hermite interpolation */
  73. {
  74.     double  t;
  75.     
  76.     t = argument(5);
  77.     return( hermite(argument(1), argument(2),
  78.             argument(3), argument(4), t) );
  79. }
  80.  
  81.  
  82. setnoisefuncs()            /* add noise functions to library */
  83. {
  84.     register int  i;
  85.  
  86.     funset(hermite_name, 5, ':', l_hermite);
  87.     funset(fnoise_name, 3, ':', l_noise3);
  88.     i = 4;
  89.     while (i--)
  90.         funset(noise_name[i], 3, ':', l_noise3);
  91. }
  92.  
  93.  
  94. double *
  95. noise3(xnew)            /* compute the noise function */
  96. register double  xnew[3];
  97. {
  98.     extern double  floor();
  99.     static double  x[3] = {-100000.0, -100000.0, -100000.0};
  100.     static double  f[4];
  101.  
  102.     if (x[0]==xnew[0] && x[1]==xnew[1] && x[2]==xnew[2])
  103.         return(f);
  104.     x[0] = xnew[0]; x[1] = xnew[1]; x[2] = xnew[2];
  105.     xlim[0][0] = floor(x[0]); xlim[0][1] = xlim[0][0] + 1;
  106.     xlim[1][0] = floor(x[1]); xlim[1][1] = xlim[1][0] + 1;
  107.     xlim[2][0] = floor(x[2]); xlim[2][1] = xlim[2][0] + 1;
  108.     xarg[0] = x[0] - xlim[0][0];
  109.     xarg[1] = x[1] - xlim[1][0];
  110.     xarg[2] = x[2] - xlim[2][0];
  111.     interpolate(f, 0, 3);
  112.     return(f);
  113. }
  114.  
  115.  
  116. static
  117. interpolate(f, i, n)
  118. double  f[4];
  119. register int  i, n;
  120. {
  121.     double  f0[4], f1[4];
  122.  
  123.     if (n == 0) {
  124.         f[A] = rand3a(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
  125.         f[B] = rand3b(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
  126.         f[C] = rand3c(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
  127.         f[D] = rand3d(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
  128.     } else {
  129.         n--;
  130.         interpolate(f0, i, n);
  131.         interpolate(f1, i | 1<<n, n);
  132.         f[A] = (1.0-xarg[n])*f0[A] + xarg[n]*f1[A];
  133.         f[B] = (1.0-xarg[n])*f0[B] + xarg[n]*f1[B];
  134.         f[C] = (1.0-xarg[n])*f0[C] + xarg[n]*f1[C];
  135.         f[D] = hermite(f0[D], f1[D], f0[n], f1[n], xarg[n]);
  136.     }
  137. }
  138.  
  139.  
  140. double
  141. frand(s)            /* get random number from seed */
  142. register long  s;
  143. {
  144.     s = s<<13 ^ s;
  145.     return(1.0-((s*(s*s*15731+789221)+1376312589)&0x7fffffff)/1073741824.0);
  146. }
  147.  
  148.  
  149. double
  150. fnoise3(p)            /* compute fractal noise function */
  151. double  p[3];
  152. {
  153.     double  floor();
  154.     long  t[3], v[3], beg[3];
  155.     double  fval[8], fc;
  156.     int  branch;
  157.     register long  s;
  158.     register int  i, j;
  159.                         /* get starting cube */
  160.     s = (long)(1.0/EPSILON);
  161.     for (i = 0; i < 3; i++) {
  162.         t[i] = s*p[i];
  163.         beg[i] = s*floor(p[i]);
  164.     }
  165.     for (j = 0; j < 8; j++) {
  166.         for (i = 0; i < 3; i++) {
  167.             v[i] = beg[i];
  168.             if (j & 1<<i)
  169.                 v[i] += s;
  170.         }
  171.         fval[j] = frand3(v[0],v[1],v[2]);
  172.     }
  173.                         /* compute fractal */
  174.     for ( ; ; ) {
  175.         fc = 0.0;
  176.         for (j = 0; j < 8; j++)
  177.             fc += fval[j];
  178.         fc *= 0.125;
  179.         if ((s >>= 1) == 0)
  180.             return(fc);        /* close enough */
  181.         branch = 0;
  182.         for (i = 0; i < 3; i++) {    /* do center */
  183.             v[i] = beg[i] + s;
  184.             if (t[i] > v[i]) {
  185.                 branch |= 1<<i;
  186.             }
  187.         }
  188.         fc += s*EPSILON*frand3(v[0],v[1],v[2]);
  189.         fval[~branch & 7] = fc;
  190.         for (i = 0; i < 3; i++) {    /* do faces */
  191.             if (branch & 1<<i)
  192.                 v[i] += s;
  193.             else
  194.                 v[i] -= s;
  195.             fc = 0.0;
  196.             for (j = 0; j < 8; j++)
  197.                 if (~(j^branch) & 1<<i)
  198.                     fc += fval[j];
  199.             fc = 0.25*fc + s*EPSILON*frand3(v[0],v[1],v[2]);
  200.             fval[~(branch^1<<i) & 7] = fc;
  201.             v[i] = beg[i] + s;
  202.         }
  203.         for (i = 0; i < 3; i++) {    /* do edges */
  204.             j = (i+1)%3;
  205.             if (branch & 1<<j)
  206.                 v[j] += s;
  207.             else
  208.                 v[j] -= s;
  209.             j = (i+2)%3;
  210.             if (branch & 1<<j)
  211.                 v[j] += s;
  212.             else
  213.                 v[j] -= s;
  214.             fc = fval[branch & ~(1<<i)];
  215.             fc += fval[branch | 1<<i];
  216.             fc = 0.5*fc + s*EPSILON*frand3(v[0],v[1],v[2]);
  217.             fval[branch^1<<i] = fc;
  218.             j = (i+1)%3;
  219.             v[j] = beg[j] + s;
  220.             j = (i+2)%3;
  221.             v[j] = beg[j] + s;
  222.         }
  223.         for (i = 0; i < 3; i++)        /* new cube */
  224.             if (branch & 1<<i)
  225.                 beg[i] += s;
  226.     }
  227. }
  228.  
  229.