Back to index

radiance  4R0+20100331
noise3.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: noise3.c,v 2.10 2006/03/02 17:16:56 greg Exp $";
00003 #endif
00004 /*
00005  *  noise3.c - noise functions for random textures.
00006  *
00007  *     Credit for the smooth algorithm goes to Ken Perlin.
00008  *     (ref. SIGGRAPH Vol 19, No 3, pp 287-96)
00009  */
00010 
00011 #include "copyright.h"
00012 
00013 #include  <math.h>
00014 
00015 #include  "calcomp.h"
00016 #include  "func.h"
00017 
00018 #define  A           0
00019 #define  B           1
00020 #define  C           2
00021 #define  D           3
00022 
00023 #define  rand3a(x,y,z)      frand(67*(x)+59*(y)+71*(z))
00024 #define  rand3b(x,y,z)      frand(73*(x)+79*(y)+83*(z))
00025 #define  rand3c(x,y,z)      frand(89*(x)+97*(y)+101*(z))
00026 #define  rand3d(x,y,z)      frand(103*(x)+107*(y)+109*(z))
00027 
00028 #define  hpoly1(t)   ((2.0*t-3.0)*t*t+1.0)
00029 #define  hpoly2(t)   (-2.0*t+3.0)*t*t
00030 #define  hpoly3(t)   ((t-2.0)*t+1.0)*t
00031 #define  hpoly4(t)   (t-1.0)*t*t
00032 
00033 #define  hermite(p0,p1,r0,r1,t)  ( p0*hpoly1(t) + \
00034                                    p1*hpoly2(t) + \
00035                                    r0*hpoly3(t) + \
00036                                    r1*hpoly4(t) )
00037 
00038 static char  noise_name[4][8] = {"noise3x", "noise3y", "noise3z", "noise3"};
00039 static char  fnoise_name[] = "fnoise3";
00040 static char  hermite_name[] = "hermite";
00041 
00042 //double  *noise3(), fnoise3(), frand();
00043 //static  interpolate();
00044 
00045 static long  xlim[3][2];
00046 static double  xarg[3];
00047 
00048 #define  EPSILON     .001          /* error allowed in fractal */
00049 
00050 #define  frand3(x,y,z)      frand(17*(x)+23*(y)+29*(z))
00051 
00052 static double l_noise3(char  *nam);
00053 static double l_hermite(char *nm);
00054 static double * noise3(double  xnew[3]);
00055 static void interpolate(double  f[4], int  i, int  n);
00056 static double frand(long  s);
00057 static double fnoise3(double  p[3]);
00058 
00059 
00060 static double
00061 l_noise3(                   /* compute a noise function */
00062        register char  *nam
00063 )
00064 {
00065        register int  i;
00066        double  x[3];
00067                                    /* get point */
00068        x[0] = argument(1);
00069        x[1] = argument(2);
00070        x[2] = argument(3);
00071                                    /* make appropriate call */
00072        if (nam == fnoise_name)
00073               return(fnoise3(x));
00074        i = 4;
00075        while (i--)
00076               if (nam == noise_name[i])
00077                      return(noise3(x)[i]);
00078        eputs(nam);
00079        eputs(": called l_noise3!\n");
00080        quit(1);
00081        return 1; /* pro forma return */
00082 }
00083 
00084 
00085 static double
00086 l_hermite(char *nm)         /* library call for hermite interpolation */
00087 {
00088        double  t;
00089        
00090        t = argument(5);
00091        return( hermite(argument(1), argument(2),
00092                      argument(3), argument(4), t) );
00093 }
00094 
00095 
00096 extern void
00097 setnoisefuncs(void)                /* add noise functions to library */
00098 {
00099        register int  i;
00100 
00101        funset(hermite_name, 5, ':', l_hermite);
00102        funset(fnoise_name, 3, ':', l_noise3);
00103        i = 4;
00104        while (i--)
00105               funset(noise_name[i], 3, ':', l_noise3);
00106 }
00107 
00108 
00109 static double *
00110 noise3(                     /* compute the noise function */
00111        register double  xnew[3]
00112 )
00113 {
00114        static double  x[3] = {-100000.0, -100000.0, -100000.0};
00115        static double  f[4];
00116 
00117        if (x[0]==xnew[0] && x[1]==xnew[1] && x[2]==xnew[2])
00118               return(f);
00119        x[0] = xnew[0]; x[1] = xnew[1]; x[2] = xnew[2];
00120        xlim[0][0] = floor(x[0]); xlim[0][1] = xlim[0][0] + 1;
00121        xlim[1][0] = floor(x[1]); xlim[1][1] = xlim[1][0] + 1;
00122        xlim[2][0] = floor(x[2]); xlim[2][1] = xlim[2][0] + 1;
00123        xarg[0] = x[0] - xlim[0][0];
00124        xarg[1] = x[1] - xlim[1][0];
00125        xarg[2] = x[2] - xlim[2][0];
00126        interpolate(f, 0, 3);
00127        return(f);
00128 }
00129 
00130 
00131 static void
00132 interpolate(
00133        double  f[4],
00134        register int  i,
00135        register int  n
00136 )
00137 {
00138        double  f0[4], f1[4], hp1, hp2;
00139 
00140        if (n == 0) {
00141               f[A] = rand3a(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
00142               f[B] = rand3b(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
00143               f[C] = rand3c(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
00144               f[D] = rand3d(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
00145        } else {
00146               n--;
00147               interpolate(f0, i, n);
00148               interpolate(f1, i | 1<<n, n);
00149               hp1 = hpoly1(xarg[n]); hp2 = hpoly2(xarg[n]);
00150               f[A] = f0[A]*hp1 + f1[A]*hp2;
00151               f[B] = f0[B]*hp1 + f1[B]*hp2;
00152               f[C] = f0[C]*hp1 + f1[C]*hp2;
00153               f[D] = f0[D]*hp1 + f1[D]*hp2 +
00154                             f0[n]*hpoly3(xarg[n]) + f1[n]*hpoly4(xarg[n]);
00155        }
00156 }
00157 
00158 
00159 static double
00160 frand(               /* get random number from seed */
00161        register long  s
00162 )
00163 {
00164        s = s<<13 ^ s;
00165        return(1.0-((s*(s*s*15731+789221)+1376312589)&0x7fffffff)/1073741824.0);
00166 }
00167 
00168 
00169 static double
00170 fnoise3(                    /* compute fractal noise function */
00171        double  p[3]
00172 )
00173 {
00174        long  t[3], v[3], beg[3];
00175        double  fval[8], fc;
00176        int  branch;
00177        register long  s;
00178        register int  i, j;
00179                                           /* get starting cube */
00180        s = (long)(1.0/EPSILON);
00181        for (i = 0; i < 3; i++) {
00182               t[i] = s*p[i];
00183               beg[i] = s*floor(p[i]);
00184        }
00185        for (j = 0; j < 8; j++) {
00186               for (i = 0; i < 3; i++) {
00187                      v[i] = beg[i];
00188                      if (j & 1<<i)
00189                             v[i] += s;
00190               }
00191               fval[j] = frand3(v[0],v[1],v[2]);
00192        }
00193                                           /* compute fractal */
00194        for ( ; ; ) {
00195               fc = 0.0;
00196               for (j = 0; j < 8; j++)
00197                      fc += fval[j];
00198               fc *= 0.125;
00199               if ((s >>= 1) == 0)
00200                      return(fc);          /* close enough */
00201               branch = 0;
00202               for (i = 0; i < 3; i++) {   /* do center */
00203                      v[i] = beg[i] + s;
00204                      if (t[i] > v[i]) {
00205                             branch |= 1<<i;
00206                      }
00207               }
00208               fc += s*EPSILON*frand3(v[0],v[1],v[2]);
00209               fval[~branch & 7] = fc;
00210               for (i = 0; i < 3; i++) {   /* do faces */
00211                      if (branch & 1<<i)
00212                             v[i] += s;
00213                      else
00214                             v[i] -= s;
00215                      fc = 0.0;
00216                      for (j = 0; j < 8; j++)
00217                             if (~(j^branch) & 1<<i)
00218                                    fc += fval[j];
00219                      fc = 0.25*fc + s*EPSILON*frand3(v[0],v[1],v[2]);
00220                      fval[~(branch^1<<i) & 7] = fc;
00221                      v[i] = beg[i] + s;
00222               }
00223               for (i = 0; i < 3; i++) {   /* do edges */
00224                      if ((j = i+1) == 3) j = 0;
00225                      if (branch & 1<<j)
00226                             v[j] += s;
00227                      else
00228                             v[j] -= s;
00229                      if (++j == 3) j = 0;
00230                      if (branch & 1<<j)
00231                             v[j] += s;
00232                      else
00233                             v[j] -= s;
00234                      fc = fval[branch & ~(1<<i)];
00235                      fc += fval[branch | 1<<i];
00236                      fc = 0.5*fc + s*EPSILON*frand3(v[0],v[1],v[2]);
00237                      fval[branch^1<<i] = fc;
00238                      if ((j = i+1) == 3) j = 0;
00239                      v[j] = beg[j] + s;
00240                      if (++j == 3) j = 0;
00241                      v[j] = beg[j] + s;
00242               }
00243               for (i = 0; i < 3; i++)            /* new cube */
00244                      if (branch & 1<<i)
00245                             beg[i] += s;
00246        }
00247 }