Article 32788 of comp.graphics: Newsgroups: comp.graphics Path: kpc!news.kpc.com!decwrl!elroy.jpl.nasa.gov!usc!howland.reston.ans.net!zaphod.mps.ohio-state.edu!cs.utexas.edu!uunet!ftpbox!mothost!lmpsbbs!setzer From: setzer@ssd.comm.mot.com (Thomas Setzer) Subject: Summary on noise(Re: Noise) Organization: Motorola Date: Tue, 23 Mar 1993 17:35:43 GMT Message-ID: <1993Mar23.173543.12932@lmpsbbs.comm.mot.com> Sender: news@lmpsbbs.comm.mot.com (Net News) Nntp-Posting-Host: 145.1.56.126 Lines: 249 Due to several requests, I am posting a summary of the initial responses I got about the noise function. Other references on noise: Ken Perlin, An Image Synthesizer, Siggraph '85 p 287-296: Perlin, K. "Hypertexture," Computer Graphics 19(3), July 1989, pp. 253-262. In the same issue is another paper of interest: Lewis, J.P., "Algorithms for Solid Noise Synthesis," Computer Graphics 19(3), July 1989, pp. 263-270. Thanks to Larry Gritz and others for those(others may have given those references as well, but I saved Larry's copy) And now for the good stuff(read code), thanks to Bernie Kirbys mail. --------------------------------------------------------------- EXCERPTED FROM SIGGRAPH 92, COURSE 23 PROCEDURAL MODELING Ken Perlin New York University 3.6 TURBULENCE AND NOISE 3.6.1 The turbulence function The turbulence function, which you use to make marble, clouds, explosions, etc., is just a simple fractal gen- erating loop built on top of the noise function. It is not a real turbulence model at all. The key trick is the use of the fabs() function, which makes the function have gradient discontinuity "fault lines" at all scales. This fools the eye into thinking it is seeing the results of turbulent flow. The turbulence() function gives the best results when used as a phase shift, as in the familiar marble trick: sin(point + turbulence(point) * point.x); Note the second argument below, lofreq, which sets the lowest desired frequency component of the turbulence. The third, hifreq, argument is used by the function to ensure that the turbulence effect reaches down to the single pixel level, but no further. I usually set this argument equal to the image resolution. float turbulence(point, lofreq, hifreq) float point[3], freq, resolution; { float noise3(), freq, t, p[3]; p[0] = point[0] + 123.456; p[1] = point[1]; p[2] = point[2]; t = 0; for (freq = lofreq ; freq < hifreq ; freq *= 2.) { t += fabs(noise3(p)) / freq; p[0] *= 2.; p[1] *= 2.; p[2] *= 2.; } return t - 0.3; /* readjust to make mean value = 0.0 */ } 3.6.2 The noise function noise3 is a rough approximation to "pink" (band-limited) noise, implemented by a pseudorandom tricubic spline. Given a vector in 3-space, it returns a value between -1.0 and 1.0. There are two principal tricks to make it run fast: - Precompute an array of pseudo-random unit length gra- dients g[n]. - Precompute a permutation array p[] of the first n integers. Given the above two arrays, any integer lattice point (i,j,k) can be quickly mapped to a pseudorandom gradient vector by: g[ (p[ (p[i] + j) % n ] + k) % n] By extending the g[] and p[] arrays, so that g[n+i]=g[i] and p[n+i]=p[i], the above lookup can be replaced by the (somewhat faster): g[ p[ p[i] + j ] + k ] Now for any point in 3-space, we just have to do the following two steps: (1) Get the gradient for each of its surrounding 8 integer lattice points as above. (2) Do a tricubic hermite spline interpolation, giving each lattice point the value 0.0. The second step above is just an evaluation of the her- mite derivative basis function 3 * t * t - 2 * t * t * t in each by a dot product of the gradient at the lattice. Here is my implementation in C of the noise function. Feel free to use it, as long as you reference where you got it. :-) /* noise function over R3 - implemented by a pseudorandom tricubic spline */ #include #include #define DOT(a,b) (a[0] * b[0] + a[1] * b[1] + a[2] * b[2]) #define B 256 static p[B + B + 2]; static float g[B + B + 2][3]; static start = 1; #define setup(i,b0,b1,r0,r1) \ t = vec[i] + 10000.; \ b0 = ((int)t) & (B-1); \ b1 = (b0+1) & (B-1); \ r0 = t - (int)t; \ r1 = r0 - 1.; float noise3(vec) float vec[3]; { int bx0, bx1, by0, by1, bz0, bz1, b00, b10, b01, b11; float rx0, rx1, ry0, ry1, rz0, rz1, *q, sy, sz, a, b, c, d, t, u, v; register i, j; if (start) { start = 0; init(); } setup(0, bx0,bx1, rx0,rx1); setup(1, by0,by1, ry0,ry1); setup(2, bz0,bz1, rz0,rz1); i = p[ bx0 ]; j = p[ bx1 ]; b00 = p[ i + by0 ]; b10 = p[ j + by0 ]; b01 = p[ i + by1 ]; b11 = p[ j + by1 ]; #define at(rx,ry,rz) ( rx * q[0] + ry * q[1] + rz * q[2] ) #define surve(t) ( t * t * (3. - 2. * t) ) #define lerp(t, a, b) ( a + t * (b - a) ) sx = surve(rx0); sy = surve(ry0); sz = surve(rz0); q = g[ b00 + bz0 ] ; u = at(rx0,ry0,rz0); q = g[ b10 + bz0 ] ; v = at(rx1,ry0,rz0); a = lerp(sx, u, v); q = g[ b01 + bz0 ] ; u = at(rx0,ry1,rz0); q = g[ b11 + bz0 ] ; v = at(rx1,ry1,rz0); b = lerp(sx, u, v); c = lerp(sy, a, b); /* interpolate in y at lo x */ q = g[ b00 + bz1 ] ; u = at(rx0,ry0,rz1); q = g[ b10 + bz1 ] ; v = at(rx1,ry0,rz1); a = lerp(sx, u, v); q = g[ b01 + bz1 ] ; u = at(rx0,ry1,rz1); q = g[ b11 + bz1 ] ; v = at(rx1,ry1,rz1); b = lerp(sx, u, v); d = lerp(sy, a, b); /* interpolate in y at hi x */ return 1.5 * lerp(sz, c, d); /* interpolate in z */ } static init() { long random(); int i, j, k; float v[3], s; /* Create an array of random gradient vectors uniformly on the unit sphere */ srandom(1); for (i = 0 ; i < B ; i++) { do { /* Choose uniformly in a cube */ for (j=0 ; j<3 ; j++) v[j] = (float)((random() % (B + B)) - B) / B; s = DOT(v,v); } while (s > 1.0); /* If not in sphere try again */ s = sqrt(s); for (j = 0 ; j < 3 ; j++) /* Else normalize */ g[i][j] = v[j] / s; } /* Create a pseudorandom permutation of [1..B] */ for (i = 0 ; i < B ; i++) p[i] = i; for (i = B ; i > 0 ; i -= 2) { k = p[i]; p[i] = p[j = random() % B]; p[j] = k; } /* Extend g and p arrays to allow for faster indexing */ for (i = 0 ; i < B + 2 ; i++) { p[B + i] = p[i]; for (j = 0 ; j < 3 ; j++) g[B + i][j] = g[i][j]; } }