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