home *** CD-ROM | disk | FTP | other *** search
/ The Fred Fish Collection 1.5 / ffcollection-1-5-1992-11.iso / ff_progs / prog-asm / talincod.lha / talincode.lha / perlin.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-08-06  |  6.8 KB  |  284 lines
  1.  
  2.  
  3. #include <math.h>
  4.  
  5.     /* Perlin Noise Function -- from Graphics Gems II */
  6.  
  7. #if 0
  8. #define hermite(p0,p1,r0,r1,t) (    p0*((2.*t-3.)*t*t+1.) + \
  9.                                     p1*(-2.*t+3.)*t*t + \
  10.                                     r0*((t-2.)*t+1.)*t + \
  11.                                     r1*(t-1.)*t*t )
  12.  
  13. #define dhermite(p0,p1,r0,r1,t) (    p0*(t-1.)*6.*t + \
  14.                                     p1*(-t+1.)*6.*t + \
  15.                                     r0*((3.*t-4.)*t + 1.) + \
  16.                                     r1*(3.*t-2.)*t )
  17. #else
  18. #define hermite(p0,p1,r0,r1,t) (    (p0-p1)*((2.*t-3.)*t*t) + p0 \
  19.                                     r0*((t-2.)*t+1.)*t + \
  20.                                     r1*(t-1.)*t*t )
  21.  
  22. #define dhermite(p0,p1,r0,r1,t) (    (p0-p1)*(t-1.)*6.*t + \
  23.                                     r0*((3.*t-4.)*t + 1.) + \
  24.                                     r1*(3.*t-2.)*t )
  25. #endif
  26.  
  27. #define rand3a(x,y,z)    frand(67*(x)+59*(y)+71*(z))
  28. #define rand3b(x,y,z)    frand(73*(x)+79*(y)+83*(z))
  29. #define rand3c(x,y,z)    frand(89*(x)+97*(y)+101*(z))
  30. #define rand3d(x,y,z)    frand(103*(x)+107*(y)+109*(z))
  31.  
  32. short                    xlim[3][2];                /* integer bounds for point        */
  33. short                    xint[3];                /* fractional part                */
  34. double                    xarg[3];                /* fractional part                */
  35.  
  36. #if 0
  37. double frand(long s)
  38. {    s = s<<13 ^ s;
  39.     return (1.0 - ((s*(s*s*15731+789221)+1376312589)&0x7fffffff)
  40.                     / 1073741824.0 );
  41. }
  42. #else
  43. double frand(long s)
  44. {    return (float)(rand24(s) - 128) * .0078125;    /* / 128 */
  45. }
  46. #endif
  47.  
  48. double calc_hermite(double p0,double p1,double r0,double r1,double t)
  49. {    return (p0*((2.*t-3.)*t*t+1.) +
  50.             p1*(-2.*t+3.)*t*t +
  51.             r0*((t-2.)*t+1.)*t +
  52.             r1*(t-1.)*t*t );
  53. }
  54.  
  55. double calc_dhermite(double p0,double p1,double r0,double r1,double t)
  56. {    return (p0*(t-1.)*6.*t +
  57.             p1*(-t+1.)*6.*t +
  58.             r0*((3.*t-4.)*t + 1.) +
  59.             r1*(3.*t-2.)*t );
  60. }
  61.  
  62. extern char                rtable0[256],
  63.                         rtable1[256],
  64.                         rtable2[256],
  65.                         rtable3[256];
  66.  
  67. long                    xyzrand16t(short,short,short,char *);
  68.  
  69. #define HERMITE_TABSIZ    1024
  70. #define HERMITE_SCALE    0x0400
  71. #define HERMITE_UNITY    0x8000
  72.  
  73. struct HermiteTable {
  74.     short                cp0,cp1,cr0,cr1,        /* hermite coefficients            */
  75.                         dp0,dp1,dr0,dr1;        /* hermite derivative coeff.    */
  76.     float                ep0,ep1,er0,er1,
  77.                         fp0,fp1,fr0,fr1;
  78.     short                t;
  79. };
  80.  
  81. static struct HermiteTable h_table[HERMITE_TABSIZ];
  82. struct HermiteTable        *h_arg[4];
  83.  
  84. init_hermite()
  85. {    float                t;
  86.     short                i;
  87.     struct HermiteTable    *ht;
  88.  
  89.     for (i=0; i<HERMITE_TABSIZ; i++)
  90.     {    ht = &h_table[i];
  91.  
  92.         t = (float)i / (float)HERMITE_TABSIZ;
  93.  
  94.         ht->cp0 = ((2.*t-3.)*t*t+1.) * (float)HERMITE_SCALE;
  95.         ht->cp1 = (-2.*t+3.)*t*t * (float)HERMITE_SCALE;
  96.         ht->cr0 = ((t-2.)*t+1.)*t * (float)HERMITE_SCALE;
  97.         ht->cr1 = (t-1.)*t*t * (float)HERMITE_SCALE;
  98.  
  99.         ht->dp0 = (t-1.)*6.*t * (float)HERMITE_SCALE;
  100.         ht->dp1 = (-t+1.)*6.*t * (float)HERMITE_SCALE;
  101.         ht->dr0 = ((3.*t-4.)*t + 1.) * (float)HERMITE_SCALE;
  102.         ht->dr1 = (3.*t-2.)*t * (float)HERMITE_SCALE;
  103. #if 0
  104.         ht->ep0 = ((2.*t-3.)*t*t+1.);
  105.         ht->ep1 = (-2.*t+3.)*t*t;
  106.         ht->er0 = ((t-2.)*t+1.)*t;
  107.         ht->er1 = (t-1.)*t*t;
  108.  
  109.         ht->fp0 = (t-1.)*6.*t;
  110.         ht->fp1 = (-t+1.)*6.*t;
  111.         ht->fr0 = ((3.*t-4.)*t);
  112.         ht->fr1 = (3.*t-2.)*t;
  113. #endif
  114.         ht->t    = (int)(t * 4096.0);
  115.     }
  116. }
  117.  
  118. long lerp_i(short p1, short p2, short t)
  119. {    return p1 + (((p2 - p1) * t) >> 12);
  120. }
  121.  
  122. #if 0
  123. interpolate(double f[4], short e[4], int i, int n)
  124. {
  125.     if (n == 0)
  126.     {    long            ix = xlim[0][i&1],
  127.                         iy = xlim[1][i>>1&1],
  128.                         iz = xlim[2][i>>2];
  129.  
  130. #if 0
  131.         f[0] = (float)(xyzrand16t(ix,iy,iz,rtable0) - 128);
  132.         f[1] = (float)(xyzrand16t(ix,iy,iz,rtable1) - 128);
  133.         f[2] = (float)(xyzrand16t(ix,iy,iz,rtable2) - 128);
  134.         f[3] = (float)(xyzrand16t(ix,iy,iz,rtable3) - 128);
  135. #endif
  136.  
  137.         e[0] = ((xyzrand16t(ix,iy,iz,rtable0) - 128) << 5) + (1<<4);
  138.         e[1] = ((xyzrand16t(ix,iy,iz,rtable1) - 128) << 5) + (1<<4);
  139.         e[2] = ((xyzrand16t(ix,iy,iz,rtable2) - 128) << 5) + (1<<4);
  140.         e[3] = ((xyzrand16t(ix,iy,iz,rtable3) - 128) << 5) + (1<<4);
  141.  
  142. #if 0
  143.         f[0] = rand3a(ix,iy,iz);
  144.         f[1] = rand3b(ix,iy,iz);
  145.         f[2] = rand3c(ix,iy,iz);
  146.         f[3] = rand3d(ix,iy,iz);
  147. #endif
  148.         return;
  149.     }
  150.     else
  151.     {    double            f0[4],                    /* result for first half        */
  152.                         f1[4];                    /* result for second half        */
  153.         double            t,
  154.                         nt;
  155.         short            e0[4],
  156.                         e1[4];
  157.         struct HermiteTable    *ht;
  158.  
  159.         n--;
  160.         interpolate(f0,e0,i,n);
  161.         interpolate(f1,e1,i|1<<n,n);
  162.  
  163.         ht = h_arg[n];
  164.  
  165.         e[0] = lerp_i( e0[0], e1[0], ht->t );
  166.         e[1] = lerp_i( e0[1], e1[1], ht->t );
  167.         e[2] = lerp_i( e0[2], e1[2], ht->t );
  168.  
  169.         e[3] = (e0[3] * ht->cp0 +
  170.                 e1[3] * ht->cp1 +
  171.                 e0[n] * ht->cr0 +
  172.                 e1[n] * ht->cr1) >> 10;
  173.  
  174. /*        e[3] = (int)(f[3] * 32.0); */
  175.  
  176. #if 0
  177.         t = xarg[n];
  178.         nt = 1.0 - t;
  179.  
  180.         f[0] = nt*f0[0] + t*f1[0];
  181.         f[1] = nt*f0[1] + t*f1[1];
  182.         f[2] = nt*f0[2] + t*f1[2];
  183.  
  184. /*        f[3] = calc_hermite(f0[3], f1[3], f0[n], f1[n], t); */
  185. /*        f[n] = calc_dhermite(f0[3], f1[3], f0[n], f1[n], t); */
  186. #endif
  187.     }
  188. }
  189. #else
  190. interpolate(short e[4], int i, int n)
  191. {
  192.     if (n == 0)
  193.     {    long            ix = xlim[0][i&1],
  194.                         iy = xlim[1][i>>1&1],
  195.                         iz = xlim[2][i>>2];
  196.  
  197.         e[0] = ((xyzrand16t(ix,iy,iz,rtable0) - 128) << 5) + (1<<4);
  198.         e[1] = ((xyzrand16t(ix,iy,iz,rtable1) - 128) << 5) + (1<<4);
  199.         e[2] = ((xyzrand16t(ix,iy,iz,rtable2) - 128) << 5) + (1<<4);
  200.         e[3] = ((xyzrand16t(ix,iy,iz,rtable3) - 128) << 5) + (1<<4);
  201.         return;
  202.     }
  203.     else
  204.     {    short            e0[4],
  205.                         e1[4];
  206.         struct HermiteTable    *ht;
  207.  
  208.         n--;
  209.         interpolate(e0,i,n);
  210.         interpolate(e1,i|1<<n,n);
  211.  
  212.         ht = h_arg[n];
  213.  
  214.         e[0] = lerp_i( e0[0], e1[0], ht->t );
  215.         e[1] = lerp_i( e0[1], e1[1], ht->t );
  216.         e[2] = lerp_i( e0[2], e1[2], ht->t );
  217.  
  218.         e[3] = (e0[3] * ht->cp0 +
  219.                 e1[3] * ht->cp1 +
  220.                 e0[n] * ht->cr0 +
  221.                 e1[n] * ht->cr1) >> 10;
  222.     }
  223. }
  224. #endif
  225.  
  226. double *noise3(double x[3])
  227. {    static double        f[4];
  228.     short                e[4];
  229.     long                l;
  230.  
  231.     l = x[0] * (float)HERMITE_TABSIZ;
  232.     xlim[0][0] = l >> 10; xlim[0][1] = xlim[0][0] + 1;
  233.     h_arg[0] = &h_table[l & (HERMITE_TABSIZ - 1)];
  234.  
  235.     l = x[1] * (float)HERMITE_TABSIZ;
  236.     xlim[1][0] = l >> 10; xlim[1][1] = xlim[1][0] + 1;
  237.     h_arg[1] = &h_table[l & (HERMITE_TABSIZ - 1)];
  238.  
  239.     l = x[2] * (float)HERMITE_TABSIZ;
  240.     xlim[2][0] = l >> 10; xlim[2][1] = xlim[2][0] + 1;
  241.     h_arg[2] = &h_table[l & (HERMITE_TABSIZ - 1)];
  242.  
  243. #if 0
  244.     xlim[0][0] = floor(x[0]); xlim[0][1] = xlim[0][0] + 1;
  245.     xlim[1][0] = floor(x[1]); xlim[1][1] = xlim[1][0] + 1;
  246.     xlim[2][0] = floor(x[2]); xlim[2][1] = xlim[2][0] + 1;
  247.     xarg[0] = x[0] - xlim[0][0];
  248.     xarg[1] = x[1] - xlim[1][0];
  249.     xarg[2] = x[2] - xlim[2][0];
  250.     h_arg[0] = &h_table[ (int)(xarg[0] * HERMITE_TABSIZ) ];
  251.     h_arg[1] = &h_table[ (int)(xarg[1] * HERMITE_TABSIZ) ];
  252.     h_arg[2] = &h_table[ (int)(xarg[2] * HERMITE_TABSIZ) ];
  253. #endif
  254.     interpolate(e,0,3);
  255. /*    f[0] /= 128.0;
  256.     f[1] /= 128.0;
  257.     f[2] /= 128.0;
  258.     f[3] /= 128.0;
  259. */
  260.     f[3] = (float)e[3] / (32.0 * 128.0);
  261.     return (f);
  262. }
  263.  
  264. double *noise2(double x[2])
  265. {    static double        f[4];
  266.     short                e[4];
  267.  
  268.     xlim[0][0] = floor(x[0]); xlim[0][1] = xlim[0][0] + 1;
  269.     xlim[1][0] = floor(x[1]); xlim[1][1] = xlim[1][0] + 1;
  270. /*    xlim[2][0] = floor(x[2]); xlim[2][1] = xlim[2][0] + 1; */
  271.     xarg[0] = x[0] - xlim[0][0];
  272.     xarg[1] = x[1] - xlim[1][0];
  273. /*    xarg[2] = x[2] - xlim[2][0]; */
  274.     h_arg[0] = &h_table[ (int)(xarg[0] * HERMITE_TABSIZ) ];
  275.     h_arg[1] = &h_table[ (int)(xarg[1] * HERMITE_TABSIZ) ];
  276.     h_arg[2] = &h_table[ (int)(xarg[2] * HERMITE_TABSIZ) ];
  277.     interpolate(e,0,2);
  278.     f[0] /= 128.0;
  279.     f[1] /= 128.0;
  280.     f[2] /= 128.0;
  281.     f[3] /= 128.0;
  282.     return (f);
  283. }
  284.