ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/noise3.c
(Generate patch)

Comparing ray/src/rt/noise3.c (file contents):
Revision 2.4 by greg, Mon Aug 9 10:03:45 1993 UTC vs.
Revision 2.14 by greg, Wed Nov 21 20:23:15 2018 UTC

# Line 1 | Line 1
1 – /* Copyright (c) 1988 Regents of the University of California */
2 –
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ LBL";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6 –
4   /*
5   *  noise3.c - noise functions for random textures.
6   *
7 < *     Credit for the smooth algorithm goes to Ken Perlin.
8 < *     (ref. SIGGRAPH Vol 19, No 3, pp 287-96)
9 < *
13 < *     4/15/86
14 < *     5/19/88  Added fractal noise function
7 > *     Credit for the smooth algorithm goes to Ken Perlin, and code
8 > *      translation/implementation to Rahul Narain.
9 > *      (ref. Improving Noise, Computer Graphics; Vol. 35 No. 3., 2002)
10   */
11  
12 < #include  <math.h>
12 > #include "copyright.h"
13  
14 < #define  A              0
15 < #define  B              1
21 < #define  C              2
22 < #define  D              3
14 > #include  "ray.h"
15 > #include  "func.h"
16  
17 < #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"};
17 > static char  noise_name[4][8] = {"noise3x", "noise3y", "noise3z", "noise3"};
18   static char  fnoise_name[] = "fnoise3";
41 – static char  hermite_name[] = "hermite";
19  
20 < double  *noise3(), fnoise3(), argument(), frand();
44 < static  interpolate();
20 > #define  EPSILON        .0005           /* error allowed in fractal */
21  
46 – static long  xlim[3][2];
47 – static double  xarg[3];
48 –
49 – #define  EPSILON        .001            /* error allowed in fractal */
50 –
22   #define  frand3(x,y,z)  frand(17*(x)+23*(y)+29*(z))
23  
24 + static double l_noise3(char *nam);
25 + static double noise3(double xnew[3], int i);
26 + static double noise3partial(double f3, double x[3], int i);
27 + static double perlin_noise (double x, double y, double z);
28 + static double frand(long s);
29 + static double fnoise3(double x[3]);
30  
31 +
32   static double
33 < l_noise3(nam)                   /* compute a noise function */
34 < register char  *nam;
33 > l_noise3(                       /* compute a noise function */
34 >        char  *nam
35 > )
36   {
37 <        register int  i;
37 >        int  i;
38          double  x[3];
39                                          /* get point */
40          x[0] = argument(1);
# Line 67 | Line 46 | register char  *nam;
46          i = 4;
47          while (i--)
48                  if (nam == noise_name[i])
49 <                        return(noise3(x)[i]);
49 >                        return(noise3(x,i));
50          eputs(nam);
51          eputs(": called l_noise3!\n");
52          quit(1);
53 +        return 1; /* pro forma return */
54   }
55  
56  
57 < double
58 < l_hermite()                     /* library call for hermite interpolation */
57 > void
58 > setnoisefuncs(void)                     /* add noise functions to library */
59   {
60 <        double  t;
81 <        
82 <        t = argument(5);
83 <        return( hermite(argument(1), argument(2),
84 <                        argument(3), argument(4), t) );
85 < }
60 >        int  i;
61  
87 –
88 – setnoisefuncs()                 /* add noise functions to library */
89 – {
90 –        register int  i;
91 –
92 –        funset(hermite_name, 5, ':', l_hermite);
62          funset(fnoise_name, 3, ':', l_noise3);
63          i = 4;
64          while (i--)
# Line 97 | Line 66 | setnoisefuncs()                        /* add noise functions to library */
66   }
67  
68  
69 < double *
70 < noise3(xnew)                    /* compute the noise function */
71 < register double  xnew[3];
69 > static double
70 > frand(                  /* get random number from seed */
71 >        long  s
72 > )
73   {
104 –        static double  x[3] = {-100000.0, -100000.0, -100000.0};
105 –        static double  f[4];
106 –
107 –        if (x[0]==xnew[0] && x[1]==xnew[1] && x[2]==xnew[2])
108 –                return(f);
109 –        x[0] = xnew[0]; x[1] = xnew[1]; x[2] = xnew[2];
110 –        xlim[0][0] = floor(x[0]); xlim[0][1] = xlim[0][0] + 1;
111 –        xlim[1][0] = floor(x[1]); xlim[1][1] = xlim[1][0] + 1;
112 –        xlim[2][0] = floor(x[2]); xlim[2][1] = xlim[2][0] + 1;
113 –        xarg[0] = x[0] - xlim[0][0];
114 –        xarg[1] = x[1] - xlim[1][0];
115 –        xarg[2] = x[2] - xlim[2][0];
116 –        interpolate(f, 0, 3);
117 –        return(f);
118 – }
119 –
120 –
121 – static
122 – interpolate(f, i, n)
123 – double  f[4];
124 – register int  i, n;
125 – {
126 –        double  f0[4], f1[4], hp1, hp2;
127 –
128 –        if (n == 0) {
129 –                f[A] = rand3a(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
130 –                f[B] = rand3b(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
131 –                f[C] = rand3c(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
132 –                f[D] = rand3d(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
133 –        } else {
134 –                n--;
135 –                interpolate(f0, i, n);
136 –                interpolate(f1, i | 1<<n, n);
137 –                hp1 = hpoly1(xarg[n]); hp2 = hpoly2(xarg[n]);
138 –                f[A] = f0[A]*hp1 + f1[A]*hp2;
139 –                f[B] = f0[B]*hp1 + f1[B]*hp2;
140 –                f[C] = f0[C]*hp1 + f1[C]*hp2;
141 –                f[D] = f0[D]*hp1 + f1[D]*hp2 +
142 –                                f0[n]*hpoly3(xarg[n]) + f1[n]*hpoly4(xarg[n]);
143 –        }
144 – }
145 –
146 –
147 – double
148 – frand(s)                        /* get random number from seed */
149 – register long  s;
150 – {
74          s = s<<13 ^ s;
75          return(1.0-((s*(s*s*15731+789221)+1376312589)&0x7fffffff)/1073741824.0);
76   }
77  
78  
79 < double
80 < fnoise3(p)                      /* compute fractal noise function */
81 < double  p[3];
79 > static double
80 > fnoise3(                        /* compute fractal noise function */
81 >        double  x[3]
82 > )
83   {
84          long  t[3], v[3], beg[3];
85          double  fval[8], fc;
86          int  branch;
87 <        register long  s;
88 <        register int  i, j;
87 >        long  s;
88 >        int  i, j;
89                                                  /* get starting cube */
90          s = (long)(1.0/EPSILON);
91          for (i = 0; i < 3; i++) {
92 <                t[i] = s*p[i];
93 <                beg[i] = s*floor(p[i]);
92 >                t[i] = s*x[i];
93 >                beg[i] = s*floor(x[i]);
94          }
95          for (j = 0; j < 8; j++) {
96                  for (i = 0; i < 3; i++) {
# Line 207 | Line 131 | double  p[3];
131                          v[i] = beg[i] + s;
132                  }
133                  for (i = 0; i < 3; i++) {       /* do edges */
134 <                        j = (i+1)%3;
134 >                        if ((j = i+1) == 3) j = 0;
135                          if (branch & 1<<j)
136                                  v[j] += s;
137                          else
138                                  v[j] -= s;
139 <                        j = (i+2)%3;
139 >                        if (++j == 3) j = 0;
140                          if (branch & 1<<j)
141                                  v[j] += s;
142                          else
# Line 221 | Line 145 | double  p[3];
145                          fc += fval[branch | 1<<i];
146                          fc = 0.5*fc + s*EPSILON*frand3(v[0],v[1],v[2]);
147                          fval[branch^1<<i] = fc;
148 <                        j = (i+1)%3;
148 >                        if ((j = i+1) == 3) j = 0;
149                          v[j] = beg[j] + s;
150 <                        j = (i+2)%3;
150 >                        if (++j == 3) j = 0;
151                          v[j] = beg[j] + s;
152                  }
153                  for (i = 0; i < 3; i++)         /* new cube */
154                          if (branch & 1<<i)
155                                  beg[i] += s;
156          }
157 + }
158 +
159 +
160 + static double
161 + noise3(                 /* compute the revised Perlin noise function */
162 +        double  xnew[3], int i
163 + )
164 + {
165 +        static int     gotV;
166 +        static double  x[3];
167 +        static double  f[4];
168 +
169 +        if (!gotV || xnew[0] != x[0] || (xnew[1] != x[1]) | (xnew[2] != x[2])) {
170 +                f[3] = perlin_noise(x[0]=xnew[0], x[1]=xnew[1], x[2]=xnew[2]);
171 +                gotV = 0x8;
172 +        }
173 +        if (!(gotV>>i & 1)) {
174 +                f[i] = noise3partial(f[3], x, i);
175 +                gotV |= 1<<i;
176 +        }
177 +        return(f[i]);
178 + }
179 +
180 + static double
181 + noise3partial(          /* compute partial derivative for ith coordinate */
182 +        double f3, double x[3], int i
183 + )
184 + {
185 +        double  fc;
186 +
187 +        switch (i) {
188 +        case 0:
189 +                fc = perlin_noise(x[0]-EPSILON, x[1], x[2]);
190 +                break;
191 +        case 1:
192 +                fc = perlin_noise(x[0], x[1]-EPSILON, x[2]);
193 +                break;
194 +        case 2:
195 +                fc = perlin_noise(x[0], x[1], x[2]-EPSILON);
196 +                break;
197 +        default:
198 +                return(.0);
199 +        }
200 +        return((f3 - fc)/EPSILON);
201 + }
202 +
203 + #define fade(t) ((t)*(t)*(t)*((t)*((t)*6. - 15.) + 10.))
204 +
205 + static double lerp(double t, double a, double b) {return a + t*(b - a);}
206 +
207 + static double
208 + grad(int hash, double x, double y, double z)
209 + {
210 +    int h = hash & 15;                      // CONVERT LO 4 BITS OF HASH CODE
211 +    double u = h<8 ? x : y,                 // INTO 12 GRADIENT DIRECTIONS.
212 +                 v = h<4 ? y : h==12|h==14 ? x : z;
213 +    return (!(h&1) ? u : -u) + (!(h&2) ? v : -v);
214 + }
215 +
216 + static const int permutation[256] = {151,160,137,91,90,15,
217 +    131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,
218 +    190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,
219 +    88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,
220 +    77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,
221 +    102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,
222 +    135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,
223 +    5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,
224 +    223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,
225 +    129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,
226 +    251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,
227 +    49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
228 +    138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180
229 + };
230 +
231 + #define p(i)    permutation[(i)&0xff]
232 +
233 + static double
234 + perlin_noise(double x, double y, double z)
235 + {
236 +    int         X, Y, Z;
237 +    double      u, v, w;
238 +    int         A, AA, AB, B, BA, BB;
239 +
240 +    X = (int)x-(x<0);                    // FIND UNIT CUBE THAT
241 +    Y = (int)y-(y<0);                    // CONTAINS POINT.
242 +    Z = (int)z-(z<0);
243 +    x -= (double)X;                          // FIND RELATIVE X,Y,Z
244 +    y -= (double)Y;                          // OF POINT IN CUBE.
245 +    z -= (double)Z;
246 +    X &= 0xff; Y &= 0xff; Z &= 0xff;
247 +
248 +    u = fade(x);                                  // COMPUTE FADE CURVES
249 +    v = fade(y);                                  // FOR EACH OF X,Y,Z.
250 +    w = fade(z);
251 +
252 +    A = p(X  )+Y; AA = p(A)+Z; AB = p(A+1)+Z;          // HASH COORDINATES OF
253 +    B = p(X+1)+Y; BA = p(B)+Z; BB = p(B+1)+Z;          // THE 8 CUBE CORNERS,
254 +
255 +    return lerp(w, lerp(v, lerp(u, grad(p(AA  ), x  , y  , z  ),  // AND ADD
256 +                                   grad(p(BA  ), x-1, y  , z  )), // BLENDED
257 +                           lerp(u, grad(p(AB  ), x  , y-1, z  ),  // RESULTS
258 +                                   grad(p(BB  ), x-1, y-1, z  ))),// FROM  8
259 +                   lerp(v, lerp(u, grad(p(AA+1), x  , y  , z-1),  // CORNERS
260 +                                   grad(p(BA+1), x-1, y  , z-1)), // OF CUBE
261 +                           lerp(u, grad(p(AB+1), x  , y-1, z-1),
262 +                                   grad(p(BB+1), x-1, y-1, z-1))));
263   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines