home *** CD-ROM | disk | FTP | other *** search
/ Encyclopedia of Graphics File Formats Companion / GFF_CD.ISO / formats / ttddd / code / noise.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-06-20  |  15.6 KB  |  438 lines
  1. /*---------------------------------------------------------------\
  2. |  Name:          NOISE.C                                        |
  3. |  Author:        Charles Congdon                                |
  4. |  Created:       09/23/92                                       |
  5. |  Last Updated:  11/04/92                                       |
  6. |                                                                |
  7. |  Purpose:                                                      |
  8. |       Implements 3D Perlin noise functions, support for the    |
  9. |  Worley 3D noise interface, etc.  The noise function itself is |
  10. |  from Graphics Gems II, page 396 - the rest guesswork and      |
  11. |  intentional changes of my own.                                |
  12. |  Greg Ward is the author of the Graphics Gems II code          |
  13. |  reproduced here.                                              |
  14. \---------------------------------------------------------------*/
  15.  
  16. #define INITNOISE 1
  17.  
  18. #include "t3dlib.h"
  19. #include "noise.h"
  20.  
  21. /*--------------------------------------\
  22. |  3D Noise function supporting cast    |
  23. \--------------------------------------*/
  24. #define HERMITE(p0,p1,r0,r1,t,t5,var32,tt)  ((p0)*(1.0 - (var32)) +   \
  25.                                              (p1)*(var32) +           \
  26.                           (r0)*((t5)-(tt)+(t)) +   \
  27.                          (r1)*(t5) )
  28.  
  29. /*-----------------------------------------------------------------\
  30. |  The following assumes a 32-bit two's compliment architecture    |
  31. \-----------------------------------------------------------------*/
  32. #define FRANDM(s,s1) (s) =  ((s)<<13) ^ (s);                                      \
  33.                      (s) =  ((s)*((s)*(s)*15731+789221) + 1376312589)&0x7fffffff; \
  34.                      (s1) = ((double)1.0 - (double)(s) / (double)1073741824.0)
  35.  
  36. double frand (register SLONG s)
  37. {
  38.    double s1;
  39.    
  40.    s  = (s<<13)^s;
  41.    s  = (s*(s*s*15731 + 789221) + 1376312589) & 0x7fffffff;
  42.    s1 = (double)1.0 - (double)((ULNG)s) / (double)1073741824.0;
  43.    return (s1);
  44. }
  45.  
  46. /*-----------------------------------------------------------\
  47. |  Noise function interpolator - optimized as best as I can  |
  48. |  without completely re-writing (tempting).                 |
  49. \-----------------------------------------------------------*/
  50. void ninterp(double f[4], short int i, register short int n)
  51. {
  52.    register SLONG      x, y, z;
  53.    register SLONG      s1;
  54.    register double     t;    /* Stack those floating point registers ! */
  55.    register double     tt;
  56.    register double     t5;
  57.    register double     var32;
  58.    register double     ttt;
  59.    register double     f0[4], f1[4];
  60.           
  61.    if (n == 0)
  62.    {
  63.       /*----------------------------------------\
  64.       |  At zero, just return lattice value.    |
  65.       \----------------------------------------*/
  66.       x = xlim[0][i&1];
  67.       y = xlim[1][(i>>1)&1];
  68.       z = xlim[2][i>>2];
  69.  
  70.       s1 = (SLONG) (((SLONG)67)*x + ((SLONG)59)*y + ((SLONG)71)*z);
  71.       FRANDM(s1, t);
  72.       
  73.       s1 = (SLONG) (((SLONG)73)*x + ((SLONG)79)*y + ((SLONG)83)*z);
  74.       FRANDM(s1, tt);
  75.       
  76.       s1 = (SLONG) (((SLONG)89)*x + ((SLONG)97)*y + ((SLONG)101)*z);
  77.       FRANDM(s1, t5);
  78.       
  79.       s1 = (SLONG) (((SLONG)103)*x + ((SLONG)107)*y + ((SLONG)109)*z);
  80.       FRANDM(s1, var32);
  81.  
  82.       f[0] = t;
  83.       f[1] = tt;
  84.       f[2] = t5;
  85.       f[3] = var32;
  86.       return;
  87.    }
  88.    
  89.    n--;                        /* Decrease order */
  90.    ninterp(f0, i, n);                    /* Compute first half */
  91.    ninterp(f1, (short int) (i | (1<<n)), n);   /* Compute second half */
  92.  
  93. #ifdef CONSUME_MY_CPU
  94.    /*----------------------------------------------------------------\
  95.    |  Use Hermite interpolation for tangent vector direction and     |
  96.    |  its magnitude.                                                 |
  97.    \----------------------------------------------------------------*/
  98.    t  = xarg[n];
  99.    tt    = t * t;
  100.    ttt   = tt * t;
  101.    t5    = ttt - tt;
  102.    var32 = tt + tt + tt - ttt - ttt;
  103.    
  104.    f[0] = HERMITE(f0[0], f1[0], f0[n], f1[n], t, t5, var32, tt);
  105.    f[1] = HERMITE(f0[1], f1[1], f0[n], f1[n], t, t5, var32, tt);
  106.    f[2] = HERMITE(f0[2], f1[2], f0[n], f1[n], t, t5, var32, tt);
  107.    f[3] = HERMITE(f0[3], f1[3], f0[n], f1[n], t, t5, var32, tt);
  108.       
  109. #else  /* CONSUME_MY_CPU */
  110.    /*----------------------------------------------------\
  111.    |  Use linear interpolation to find tangent vector    |
  112.    \----------------------------------------------------*/
  113.    t  = xarg[n];
  114.    tt = (double)1.0 - t;
  115.    f[0] = tt * f0[0] + t * f1[0];
  116.    f[1] = tt * f0[1] + t * f1[1];
  117.    f[2] = tt * f0[2] + t * f1[2];
  118.    
  119.    tt    = t * t;
  120.    ttt   = tt * t;
  121.    t5    = ttt - tt;
  122.    var32 = tt + tt + tt - ttt - ttt;
  123.    
  124.    /*--------------------------------------------------------\
  125.    |  And Hermite interpolation for the vector magnitude.    |
  126.    \--------------------------------------------------------*/
  127.    
  128.    f[3] = HERMITE(f0[3], f1[3], f0[n], f1[n], t, t5, var32, tt);
  129. #endif  /* CONSUME_MY_CPU */
  130.  
  131. }
  132.  
  133. /*--------------------------------------------------------------------\
  134. |  Here is a recursive implementation of the Perlin Noise Function    |
  135. |  as found in Graphics Gems II, page 396.  It varies randomly        |
  136. |  between 1 and -1 with an autocorollation distance of about         |
  137. |  1.  It also appears to have no pronounced harmonics.               |
  138. \--------------------------------------------------------------------*/
  139. double Noise3D(vector point, vector noisev)
  140. {
  141.    static double   f[4];
  142.    register double floored;
  143.    register double pnt;
  144.  
  145.  
  146.    /*------------------------\
  147.    |  Define the unit cube.  |
  148.    \------------------------*/
  149.    pnt = point[0];  /* Stuff the floating point registers */
  150.    floored = floor(pnt);
  151. #ifdef FLOOR_ERROR
  152.    if (pnt < (double) 0.0)
  153.       floored = floored - (double)1.0;
  154. #endif
  155.    xlim[0][0] = (SLONG)floored;
  156.        xlim[0][1] = xlim[0][0] + (SLONG)1;
  157.        xarg[0] = pnt - floored;
  158.               
  159.    pnt = point[1];
  160.    floored = floor(pnt);
  161. #ifdef FLOOR_ERROR
  162.    if (pnt < (double) 0.0)
  163.       floored = floored - (double)1.0;
  164. #endif  
  165.    xlim[1][0] = (SLONG)floored;
  166.        xlim[1][1] = xlim[1][0] + (SLONG)1;
  167.        xarg[1] = pnt - floored;
  168.  
  169.    pnt = point[2];
  170.    floored = floor(pnt);
  171. #ifdef FLOOR_ERROR
  172.    if (pnt < (double) 0.0)
  173.       floored = floored - (double)1.0;
  174. #endif
  175.    xlim[2][0] = (SLONG)floored;
  176.        xlim[2][1] = xlim[2][0] + (SLONG)1;
  177.        xarg[2] = pnt - floored;
  178.  
  179.    ninterp(f, 0, 3);   /* Now interpolate the between the random values
  180.                           at the unit cube's lattice points */
  181.    
  182.    noisev[0] = f[0];   /* Return the noise vector */
  183.    noisev[1] = f[1];
  184.    noisev[2] = f[2];
  185.    
  186.    return(f[3]);       /* And the noise value */
  187. }
  188.  
  189. /*------------------------------------------------------------\
  190. |  My partial implementation of the Worley Fractal Noise      |
  191. |  function.  Steve uses float[6] for input and output        |
  192. |  vectors, for reasons I cannot fathom.                      |
  193. |  I return the vector magnitude as well, where his seems     |
  194. |  not to.  No big deal.  If it serves the purpose, who needs |
  195. |  anything more...                                           |
  196. \------------------------------------------------------------*/
  197. double wfractalval3(vector in, double Initial_scale, double NumScales,
  198.                   double Scale_Ratio, double Amplitude_Ratio, 
  199.            double Time_Ratio, double Time, vector out)
  200. {     /* wfractval3 */
  201.    register double       sratio1;
  202.    register double       scaleaccum;
  203.    register double       ampaccum;
  204.    register double       ampl;
  205.    register vector       naccum;
  206.    register vector       inval;
  207.    register vector       outval;
  208.    register short int    i;
  209.    register double       valaccum;
  210.       
  211.             double       randseed;
  212.         unsigned long int    seedint;
  213.         short int    seed[3];
  214.    static   short int    *pseed;
  215.             short int    axisrot;
  216.         short int    iterations;
  217.         double       dummynoise;
  218.             double       noisecent[IIMAX_SCALES][3];
  219.    static   double       noisecent0[IIMAX_SCALES][3];
  220.    static   double       axisgoal[IIMAX_SCALES][3];
  221.         
  222.            
  223.    /*---------------------------------\
  224.    |  Do a little pre-calculation.    |
  225.    \---------------------------------*/
  226.    ampl = Amplitude_Ratio;
  227.    sratio1 = 1.0 / Scale_Ratio;
  228.    
  229.    scaleaccum = 1.0 / Initial_scale;  /* Fit initial scale into unit cube */
  230.    ampaccum = 1.0;
  231.    
  232.    /*-----------------------------------\
  233.    |  Get info from NumScales value.    |
  234.    \-----------------------------------*/
  235.    if (NumScales < 0.0)
  236.    {
  237.       axisrot = 1;
  238.       NumScales = fabs(NumScales);
  239.    }
  240.    else
  241.    {
  242.       axisrot = 0;
  243.    }
  244.    
  245.    /*--------------------------------------------------\
  246.    |  Get the random number generator's seed and the   |
  247.    |  number of iterations.                            |
  248.    \--------------------------------------------------*/
  249.    
  250.    randseed = floor(NumScales);
  251.    iterations = (short int) randseed;
  252.  
  253.    if (iterations > IIMAX_SCALES)              /* Enforce scale limits */
  254.       iterations = (short int) IIMAX_SCALES;
  255.    if (iterations < (short int)0)
  256.       iterations = (short int)0;
  257.    
  258.    if (randseed != NumScales)
  259.    {
  260.       randseed = NumScales - randseed;
  261.    }
  262.    else
  263.    {
  264.       randseed = 0.0;
  265.    }
  266.    
  267.    
  268.    /*-----------------------------------------------------------\
  269.    |  Initialize the random number generator and noise values.  |
  270.    \-----------------------------------------------------------*/
  271.    if (Noiseinit == 0)
  272.    {  /* Noiseinit */
  273.       /*--------------------------------------------------------\
  274.       |  If you do not have rand48-type functions (which SAS    |
  275.       |  claim to be from the AT&T Unix System V run-time       |
  276.       |  library), modify the below to play with frand above.   |
  277.       \--------------------------------------------------------*/
  278.       seedint = (unsigned long int) (2032353798 * randseed);
  279.       seed[0] = seedint >> 16;
  280.       seed[1] = seedint - (seed[0] << 16);
  281.       seed[2] = 0x330E;
  282.       pseed = seed48(seed);
  283.    
  284.       /*---------------------------------------\
  285.       |  Randomly displace the axes centers -  |
  286.       |  sorta - starting position of centers  |
  287.       |  is not as important as where they     |
  288.       |  move in subsequent frames.            |
  289.       \---------------------------------------*/
  290.       noisecent0[0][0] = erand48(pseed) * Initial_scale;
  291.       noisecent0[0][1] = erand48(pseed) * Initial_scale;
  292.       noisecent0[0][2] = erand48(pseed) * Initial_scale;
  293.       for (i = 1; i < iterations; i++)
  294.       {
  295.          noisecent0[i][0] = noisecent0[0][0];
  296.          noisecent0[i][1] = noisecent0[0][1];
  297.          noisecent0[i][2] = noisecent0[0][2];
  298.       }
  299.    
  300.       /*-------------------------------------------------------\
  301.       |  Now figure out how the axis for each noise "center"   |
  302.       |  is going to move.                                     |
  303.       \-------------------------------------------------------*/
  304.       if (Time > 0.0)
  305.       {   /* move noise vector */
  306.          /*-----------------------------------------------\
  307.          |  Create a random vector of the proper scale    |
  308.          \-----------------------------------------------*/
  309.          axisgoal[0][0] = erand48(pseed);
  310.          axisgoal[0][1] = erand48(pseed);
  311.          axisgoal[0][2] = erand48(pseed);
  312.          normalze(Initial_scale, axisgoal[0], axisgoal[0]);
  313.          if (axisrot == 1)
  314.          {   /* Move each axis individually */
  315.         dummynoise = Initial_scale * Scale_Ratio;
  316.         for (i = 1; i < iterations; i++)
  317.         {
  318.            axisgoal[i][0] = erand48(pseed);
  319.            axisgoal[i][1] = erand48(pseed);
  320.            axisgoal[i][2] = erand48(pseed);
  321.            normalze(dummynoise, axisgoal[i], axisgoal[i]);
  322.            dummynoise = dummynoise * Scale_Ratio;
  323.         }
  324.        
  325.          }   /* Move each axis individually */
  326.          else
  327.          {   /* All axes move together */
  328.         dummynoise = Initial_scale * Scale_Ratio;
  329.         for (i = 1; i < iterations; i++)
  330.         {
  331.            axisgoal[i][0] = axisgoal[0][0];
  332.            axisgoal[i][1] = axisgoal[0][1];
  333.            axisgoal[i][2] = axisgoal[0][2];
  334.            normalze(dummynoise, axisgoal[i], axisgoal[i]);
  335.            dummynoise = dummynoise * Scale_Ratio;
  336.         }
  337.          }   /* All axis move together */
  338.  
  339.       }   /* move noise vector */
  340.       Noiseinit = 1;
  341.    }  /* Noiseinit */   
  342.  
  343.    /*-----------------------------------------------\
  344.    |  Position the noise centers as appropriate.    |
  345.    \-----------------------------------------------*/
  346.    if (Time > 0.0)
  347.    {   /* Move to proper frame */
  348.       /*-----------------------------------\
  349.       |  Now move us to the proper frame.  |
  350.       |  NOTE that you have to pass FRAME  |
  351.       |  in from some external source.     |
  352.       \-----------------------------------*/
  353.       dummynoise = 1.0;
  354.       for (i = 0; i < iterations; i++)
  355.       {   /* Move to frame */
  356.          noisecent[i][0] = noisecent0[i][0] + 
  357.                        (Frame/Time) * axisgoal[i][0] * dummynoise;
  358.          noisecent[i][1] = noisecent0[i][1] + 
  359.                            (Frame/Time) * axisgoal[i][1] * dummynoise;
  360.          noisecent[i][2] = noisecent0[i][2] + 
  361.                        (Frame/Time) * axisgoal[i][2] * dummynoise;
  362.          dummynoise = dummynoise * Time_Ratio;
  363.       }   /* Move to frame */
  364.    }   /* Move to proper frame */
  365.    else
  366.    {   /* Frame 0 if Time == 0.0 */
  367.       for (i = 0; i < iterations; i++)
  368.       {
  369.          noisecent[i][0] = noisecent0[i][0];
  370.          noisecent[i][1] = noisecent0[i][1];
  371.          noisecent[i][2] = noisecent0[i][2];
  372.       }
  373.    }   /* Frame 0 if Time == 0.0 */
  374.    
  375.    /*---------------------------------\
  376.    |  OK, sum up the noise values.    |
  377.    \---------------------------------*/
  378.    valaccum = 0.0;
  379.    naccum[0] = 0.0;
  380.    naccum[1] = 0.0;
  381.    naccum[2] = 0.0;
  382.    for (i = 0; i < iterations; i++)
  383.    {
  384.       inval[0] = scaleaccum * in[0] + noisecent[i][0]; 
  385.       inval[1] = scaleaccum * in[1] + noisecent[i][1]; 
  386.       inval[2] = scaleaccum * in[2] + noisecent[i][2]; 
  387.       dummynoise = Noise3D(inval, outval);
  388.       valaccum  = valaccum  + ampaccum * dummynoise;
  389.       naccum[0] = naccum[0] + ampaccum * outval[0];
  390.       naccum[1] = naccum[1] + ampaccum * outval[1];
  391.       naccum[2] = naccum[2] + ampaccum * outval[2];
  392.       scaleaccum = scaleaccum * sratio1;  /* Scale scaling factor */
  393.       ampaccum = ampaccum * ampl;         /* Amplitude scaling factor */
  394.    }
  395.    
  396.    out[0] = naccum[0];
  397.    out[1] = naccum[1];
  398.    out[2] = naccum[2];
  399.    return(valaccum);
  400. }     /* wfractval3 */
  401.  
  402. /*-----------------------------------------------\
  403. |  Vector normalization to the desired scale.    |
  404. \-----------------------------------------------*/
  405. void normalze(double scale, vector in, vector out)
  406. {   /* normalze */
  407.     double n;
  408.     
  409.     n = sqrt(in[0]*in[0] + in[1]*in[1] + in[2]*in[2]);
  410.     if (n != 0.0)
  411.     {
  412.        n = scale / n;
  413.        out[0] = n * in[0];
  414.        out[1] = n * in[1];
  415.        out[2] = n * in[2];
  416.     }
  417.     else
  418.     {
  419.        out[0] = in[0];
  420.        out[1] = in[1];
  421.        out[2] = in[2];
  422.     }
  423. }   /* normalze */
  424.  
  425. /*--------------------------------------------------\
  426. |  A shorthand way of using the imitation Worley    |
  427. |  noise function.                                  |
  428. \--------------------------------------------------*/
  429.  
  430. void NiceN3D(vector point, vector noisev, double scale, double NumScales)
  431. {
  432.    wfractalval3(point, scale, NumScales , ISCALE_RATIO, IAMP_RATIO, 
  433.                                           ITIME_RATIO, ITIME, noisev);
  434. }
  435.  
  436.  
  437.  
  438.