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.12 by greg, Tue Feb 22 16:45:12 2011 UTC vs.
Revision 2.13 by greg, Tue Oct 8 18:59:44 2013 UTC

# Line 4 | Line 4 | static const char      RCSid[] = "$Id$";
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)
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 "copyright.h"
# Line 13 | Line 14 | static const char      RCSid[] = "$Id$";
14   #include  "ray.h"
15   #include  "func.h"
16  
16 – #define  A              0
17 – #define  B              1
18 – #define  C              2
19 – #define  D              3
20 –
21 – #define  rand3a(x,y,z)  frand(67*(x)+59*(y)+71*(z))
22 – #define  rand3b(x,y,z)  frand(73*(x)+79*(y)+83*(z))
23 – #define  rand3c(x,y,z)  frand(89*(x)+97*(y)+101*(z))
24 – #define  rand3d(x,y,z)  frand(103*(x)+107*(y)+109*(z))
25 –
26 – #define  hpoly1(t)      ((2.0*t-3.0)*t*t+1.0)
27 – #define  hpoly2(t)      (-2.0*t+3.0)*t*t
28 – #define  hpoly3(t)      ((t-2.0)*t+1.0)*t
29 – #define  hpoly4(t)      (t-1.0)*t*t
30 –
31 – #define  hermite(p0,p1,r0,r1,t)  (      p0*hpoly1(t) + \
32 –                                        p1*hpoly2(t) + \
33 –                                        r0*hpoly3(t) + \
34 –                                        r1*hpoly4(t) )
35 –
17   static char  noise_name[4][8] = {"noise3x", "noise3y", "noise3z", "noise3"};
18   static char  fnoise_name[] = "fnoise3";
38 – static char  hermite_name[] = "hermite";
19  
20 < static long  xlim[3][2];
41 < static double  xarg[3];
20 > #define  EPSILON        .0005           /* error allowed in fractal */
21  
43 – #define  EPSILON        .001            /* error allowed in fractal */
44 –
22   #define  frand3(x,y,z)  frand(17*(x)+23*(y)+29*(z))
23  
24 < static double l_noise3(char  *nam);
25 < static double l_hermite(char *nm);
26 < static double * noise3(double  xnew[3]);
27 < static void interpolate(double  f[4], int  i, int  n);
28 < static double frand(long  s);
29 < static double fnoise3(double  p[3]);
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(                       /* compute a noise function */
34 <        register char  *nam
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 69 | Line 46 | l_noise3(                      /* compute a noise function */
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);
# Line 77 | Line 54 | l_noise3(                      /* compute a noise function */
54   }
55  
56  
57 < static double
81 < l_hermite(char *nm)             /* library call for hermite interpolation */
82 < {
83 <        double  t;
84 <        
85 <        t = argument(5);
86 <        return( hermite(argument(1), argument(2),
87 <                        argument(3), argument(4), t) );
88 < }
89 <
90 <
91 < extern void
57 > void
58   setnoisefuncs(void)                     /* add noise functions to library */
59   {
60 <        register int  i;
60 >        int  i;
61  
96 –        funset(hermite_name, 5, ':', l_hermite);
62          funset(fnoise_name, 3, ':', l_noise3);
63          i = 4;
64          while (i--)
# Line 101 | Line 66 | setnoisefuncs(void)                    /* add noise functions to librar
66   }
67  
68  
104 – static double *
105 – noise3(                 /* compute the noise function */
106 –        register double  xnew[3]
107 – )
108 – {
109 –        static double  x[3] = {-100000.0, -100000.0, -100000.0};
110 –        static double  f[4];
111 –
112 –        if (x[0]==xnew[0] && x[1]==xnew[1] && x[2]==xnew[2])
113 –                return(f);
114 –        x[0] = xnew[0]; x[1] = xnew[1]; x[2] = xnew[2];
115 –        xlim[0][0] = floor(x[0]); xlim[0][1] = xlim[0][0] + 1;
116 –        xlim[1][0] = floor(x[1]); xlim[1][1] = xlim[1][0] + 1;
117 –        xlim[2][0] = floor(x[2]); xlim[2][1] = xlim[2][0] + 1;
118 –        xarg[0] = x[0] - xlim[0][0];
119 –        xarg[1] = x[1] - xlim[1][0];
120 –        xarg[2] = x[2] - xlim[2][0];
121 –        interpolate(f, 0, 3);
122 –        return(f);
123 – }
124 –
125 –
126 – static void
127 – interpolate(
128 –        double  f[4],
129 –        register int  i,
130 –        register int  n
131 – )
132 – {
133 –        double  f0[4], f1[4], hp1, hp2;
134 –
135 –        if (n == 0) {
136 –                f[A] = rand3a(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
137 –                f[B] = rand3b(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
138 –                f[C] = rand3c(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
139 –                f[D] = rand3d(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
140 –        } else {
141 –                n--;
142 –                interpolate(f0, i, n);
143 –                interpolate(f1, i | 1<<n, n);
144 –                hp1 = hpoly1(xarg[n]); hp2 = hpoly2(xarg[n]);
145 –                f[A] = f0[A]*hp1 + f1[A]*hp2;
146 –                f[B] = f0[B]*hp1 + f1[B]*hp2;
147 –                f[C] = f0[C]*hp1 + f1[C]*hp2;
148 –                f[D] = f0[D]*hp1 + f1[D]*hp2 +
149 –                                f0[n]*hpoly3(xarg[n]) + f1[n]*hpoly4(xarg[n]);
150 –        }
151 – }
152 –
153 –
69   static double
70   frand(                  /* get random number from seed */
71 <        register long  s
71 >        long  s
72   )
73   {
74          s = s<<13 ^ s;
# Line 163 | Line 78 | frand(                 /* get random number from seed */
78  
79   static double
80   fnoise3(                        /* compute fractal noise function */
81 <        double  p[3]
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 239 | Line 154 | fnoise3(                       /* compute fractal noise function */
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 && x[0]==xnew[0] && (x[1]==xnew[1]) & (x[2]==xnew[2])) {
170 +                if (!(gotV>>i & 1)) {
171 +                        f[i] = noise3partial(f[3], x, i);
172 +                        gotV |= 1<<i;
173 +                }
174 +                return(f[i]);
175 +        }
176 +        gotV = 0x8;
177 +        return(f[3] = perlin_noise(x[0]=xnew[0], x[1]=xnew[1], x[2]=xnew[2]));
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