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.6 by greg, Sat Feb 22 02:07:29 2003 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 < /* ====================================================================
12 < * The Radiance Software License, Version 1.0
13 < *
14 < * Copyright (c) 1990 - 2002 The Regents of the University of California,
15 < * through Lawrence Berkeley National Laboratory.   All rights reserved.
16 < *
17 < * Redistribution and use in source and binary forms, with or without
18 < * modification, are permitted provided that the following conditions
19 < * are met:
20 < *
21 < * 1. Redistributions of source code must retain the above copyright
22 < *         notice, this list of conditions and the following disclaimer.
23 < *
24 < * 2. Redistributions in binary form must reproduce the above copyright
25 < *       notice, this list of conditions and the following disclaimer in
26 < *       the documentation and/or other materials provided with the
27 < *       distribution.
28 < *
29 < * 3. The end-user documentation included with the redistribution,
30 < *           if any, must include the following acknowledgment:
31 < *             "This product includes Radiance software
32 < *                 (http://radsite.lbl.gov/)
33 < *                 developed by the Lawrence Berkeley National Laboratory
34 < *               (http://www.lbl.gov/)."
35 < *       Alternately, this acknowledgment may appear in the software itself,
36 < *       if and wherever such third-party acknowledgments normally appear.
37 < *
38 < * 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
39 < *       and "The Regents of the University of California" must
40 < *       not be used to endorse or promote products derived from this
41 < *       software without prior written permission. For written
42 < *       permission, please contact [email protected].
43 < *
44 < * 5. Products derived from this software may not be called "Radiance",
45 < *       nor may "Radiance" appear in their name, without prior written
46 < *       permission of Lawrence Berkeley National Laboratory.
47 < *
48 < * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
49 < * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
50 < * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
51 < * DISCLAIMED.   IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR
52 < * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
53 < * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
54 < * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
55 < * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
56 < * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
57 < * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
58 < * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
59 < * SUCH DAMAGE.
60 < * ====================================================================
61 < *
62 < * This software consists of voluntary contributions made by many
63 < * individuals on behalf of Lawrence Berkeley National Laboratory.   For more
64 < * information on Lawrence Berkeley National Laboratory, please see
65 < * <http://www.lbl.gov/>.
66 < */
12 > #include "copyright.h"
13  
14 < #include  <math.h>
14 > #include  "ray.h"
15 > #include  "func.h"
16  
70 – #define  A              0
71 – #define  B              1
72 – #define  C              2
73 – #define  D              3
74 –
75 – #define  rand3a(x,y,z)  frand(67*(x)+59*(y)+71*(z))
76 – #define  rand3b(x,y,z)  frand(73*(x)+79*(y)+83*(z))
77 – #define  rand3c(x,y,z)  frand(89*(x)+97*(y)+101*(z))
78 – #define  rand3d(x,y,z)  frand(103*(x)+107*(y)+109*(z))
79 –
80 – #define  hpoly1(t)      ((2.0*t-3.0)*t*t+1.0)
81 – #define  hpoly2(t)      (-2.0*t+3.0)*t*t
82 – #define  hpoly3(t)      ((t-2.0)*t+1.0)*t
83 – #define  hpoly4(t)      (t-1.0)*t*t
84 –
85 – #define  hermite(p0,p1,r0,r1,t)  (      p0*hpoly1(t) + \
86 –                                        p1*hpoly2(t) + \
87 –                                        r0*hpoly3(t) + \
88 –                                        r1*hpoly4(t) )
89 –
17   static char  noise_name[4][8] = {"noise3x", "noise3y", "noise3z", "noise3"};
18   static char  fnoise_name[] = "fnoise3";
92 – static char  hermite_name[] = "hermite";
19  
20 < double  *noise3(), fnoise3(), argument(), frand();
95 < static  interpolate();
20 > #define  EPSILON        .0005           /* error allowed in fractal */
21  
97 – static long  xlim[3][2];
98 – static double  xarg[3];
99 –
100 – #define  EPSILON        .001            /* error allowed in fractal */
101 –
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 118 | 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;
132 <        
133 <        t = argument(5);
134 <        return( hermite(argument(1), argument(2),
135 <                        argument(3), argument(4), t) );
136 < }
60 >        int  i;
61  
138 –
139 – setnoisefuncs()                 /* add noise functions to library */
140 – {
141 –        register int  i;
142 –
143 –        funset(hermite_name, 5, ':', l_hermite);
62          funset(fnoise_name, 3, ':', l_noise3);
63          i = 4;
64          while (i--)
# Line 148 | 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   {
155 –        static double  x[3] = {-100000.0, -100000.0, -100000.0};
156 –        static double  f[4];
157 –
158 –        if (x[0]==xnew[0] && x[1]==xnew[1] && x[2]==xnew[2])
159 –                return(f);
160 –        x[0] = xnew[0]; x[1] = xnew[1]; x[2] = xnew[2];
161 –        xlim[0][0] = floor(x[0]); xlim[0][1] = xlim[0][0] + 1;
162 –        xlim[1][0] = floor(x[1]); xlim[1][1] = xlim[1][0] + 1;
163 –        xlim[2][0] = floor(x[2]); xlim[2][1] = xlim[2][0] + 1;
164 –        xarg[0] = x[0] - xlim[0][0];
165 –        xarg[1] = x[1] - xlim[1][0];
166 –        xarg[2] = x[2] - xlim[2][0];
167 –        interpolate(f, 0, 3);
168 –        return(f);
169 – }
170 –
171 –
172 – static
173 – interpolate(f, i, n)
174 – double  f[4];
175 – register int  i, n;
176 – {
177 –        double  f0[4], f1[4], hp1, hp2;
178 –
179 –        if (n == 0) {
180 –                f[A] = rand3a(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
181 –                f[B] = rand3b(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
182 –                f[C] = rand3c(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
183 –                f[D] = rand3d(xlim[0][i&1],xlim[1][i>>1&1],xlim[2][i>>2]);
184 –        } else {
185 –                n--;
186 –                interpolate(f0, i, n);
187 –                interpolate(f1, i | 1<<n, n);
188 –                hp1 = hpoly1(xarg[n]); hp2 = hpoly2(xarg[n]);
189 –                f[A] = f0[A]*hp1 + f1[A]*hp2;
190 –                f[B] = f0[B]*hp1 + f1[B]*hp2;
191 –                f[C] = f0[C]*hp1 + f1[C]*hp2;
192 –                f[D] = f0[D]*hp1 + f1[D]*hp2 +
193 –                                f0[n]*hpoly3(xarg[n]) + f1[n]*hpoly4(xarg[n]);
194 –        }
195 – }
196 –
197 –
198 – double
199 – frand(s)                        /* get random number from seed */
200 – register long  s;
201 – {
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 258 | 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 272 | 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 && 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