5 /* ---- Ken Perlin's implementation of noise ---- */
9 #define NP 12 /* 2^N */
12 #define s_curve(t) (t * t * (3.0f - 2.0f * t))
14 #define setup(elem, b0, b1, r0, r1) \
23 #define setup_p(elem, b0, b1, r0, r1, p) \
26 b0 = (((int)t) & BM) % p; \
27 b1 = ((b0 + 1) & BM) % p; \
33 static int perm[B + B + 2]; /* permuted index from g_n onto themselves */
34 static float grad3[B + B + 2][3]; /* 3D random gradients */
35 static float grad2[B + B + 2][2]; /* 2D random gradients */
36 static float grad1[B + B + 2]; /* 1D random ... slopes */
37 static int tables_valid;
39 #define init_once() if(!tables_valid) init_noise()
41 #define SQ(x) ((x) * (x))
42 #define LERP(a, b, t) ((a) + ((b) - (a)) * (t))
44 static void init_noise()
49 /* calculate random gradients */
51 perm[i] = i; /* .. and initialize permutation mapping to identity */
53 grad1[i] = (float)((rand() % (B + B)) - B) / B;
55 grad2[i][0] = (float)((rand() % (B + B)) - B) / B;
56 grad2[i][1] = (float)((rand() % (B + B)) - B) / B;
57 len = sqrt(SQ(grad2[i][0]) + SQ(grad2[i][1]));
63 grad3[i][0] = (float)((rand() % (B + B)) - B) / B;
64 grad3[i][1] = (float)((rand() % (B + B)) - B) / B;
65 grad3[i][2] = (float)((rand() % (B + B)) - B) / B;
66 len = sqrt(SQ(grad3[i][0]) + SQ(grad3[i][1]) + SQ(grad3[i][2]));
74 /* permute indices by swapping them randomly */
76 int rand_idx = rand() % B;
79 perm[i] = perm[rand_idx];
83 /* fill up the rest of the arrays by duplicating the existing gradients */
84 /* and permutations */
85 for(i=0; i<B+2; i++) {
86 perm[B + i] = perm[i];
87 grad1[B + i] = grad1[i];
88 grad2[B + i][0] = grad2[i][0];
89 grad2[B + i][1] = grad2[i][1];
90 grad3[B + i][0] = grad3[i][0];
91 grad3[B + i][1] = grad3[i][1];
92 grad3[B + i][2] = grad3[i][2];
102 float rx0, rx1, sx, u, v;
106 setup(x, bx0, bx1, rx0, rx1);
108 u = rx0 * grad1[perm[bx0]];
109 v = rx1 * grad1[perm[bx1]];
110 return LERP(u, v, sx);
113 float noise2(float x, float y)
115 int i, j, b00, b10, b01, b11;
116 int bx0, bx1, by0, by1;
117 float rx0, rx1, ry0, ry1;
118 float sx, sy, u, v, a, b;
122 setup(x, bx0, bx1, rx0, rx1);
123 setup(y, by0, by1, ry0, ry1);
133 /* calculate hermite inteprolating factors */
137 /* interpolate along the left edge */
138 u = grad2[b00][0] * rx0 + grad2[b00][1] * ry0;
139 v = grad2[b10][0] * rx1 + grad2[b10][1] * ry0;
142 /* interpolate along the right edge */
143 u = grad2[b01][0] * rx0 + grad2[b01][1] * ry1;
144 v = grad2[b11][0] * rx1 + grad2[b11][1] * ry1;
147 /* interpolate between them */
148 return LERP(a, b, sy);
151 float noise3(float x, float y, float z)
154 int bx0, bx1, by0, by1, bz0, bz1;
155 int b00, b10, b01, b11;
156 float rx0, rx1, ry0, ry1, rz0, rz1;
158 float u, v, a, b, c, d;
162 setup(x, bx0, bx1, rx0, rx1);
163 setup(y, by0, by1, ry0, ry1);
164 setup(z, bz0, bz1, rz0, rz1);
174 /* calculate hermite interpolating factors */
179 /* interpolate along the top slice of the cell */
180 u = grad3[b00 + bz0][0] * rx0 + grad3[b00 + bz0][1] * ry0 + grad3[b00 + bz0][2] * rz0;
181 v = grad3[b10 + bz0][0] * rx1 + grad3[b10 + bz0][1] * ry0 + grad3[b10 + bz0][2] * rz0;
184 u = grad3[b01 + bz0][0] * rx0 + grad3[b01 + bz0][1] * ry1 + grad3[b01 + bz0][2] * rz0;
185 v = grad3[b11 + bz0][0] * rx1 + grad3[b11 + bz0][1] * ry1 + grad3[b11 + bz0][2] * rz0;
190 /* interpolate along the bottom slice of the cell */
191 u = grad3[b00 + bz0][0] * rx0 + grad3[b00 + bz0][1] * ry0 + grad3[b00 + bz0][2] * rz1;
192 v = grad3[b10 + bz0][0] * rx1 + grad3[b10 + bz0][1] * ry0 + grad3[b10 + bz0][2] * rz1;
195 u = grad3[b01 + bz0][0] * rx0 + grad3[b01 + bz0][1] * ry1 + grad3[b01 + bz0][2] * rz1;
196 v = grad3[b11 + bz0][0] * rx1 + grad3[b11 + bz0][1] * ry1 + grad3[b11 + bz0][2] * rz1;
201 /* interpolate between slices */
202 return LERP(c, d, sz);
206 float pnoise1(float x, int period)
209 float rx0, rx1, sx, u, v;
213 setup_p(x, bx0, bx1, rx0, rx1, period);
215 u = rx0 * grad1[perm[bx0]];
216 v = rx1 * grad1[perm[bx1]];
217 return LERP(u, v, sx);
220 float pnoise2(float x, float y, int per_x, int per_y)
222 int i, j, b00, b10, b01, b11;
223 int bx0, bx1, by0, by1;
224 float rx0, rx1, ry0, ry1;
225 float sx, sy, u, v, a, b;
229 setup_p(x, bx0, bx1, rx0, rx1, per_x);
230 setup_p(y, by0, by1, ry0, ry1, per_y);
240 /* calculate hermite inteprolating factors */
244 /* interpolate along the left edge */
245 u = grad2[b00][0] * rx0 + grad2[b00][1] * ry0;
246 v = grad2[b10][0] * rx1 + grad2[b10][1] * ry0;
249 /* interpolate along the right edge */
250 u = grad2[b01][0] * rx0 + grad2[b01][1] * ry1;
251 v = grad2[b11][0] * rx1 + grad2[b11][1] * ry1;
254 /* interpolate between them */
255 return LERP(a, b, sy);
258 float pnoise3(float x, float y, float z, int per_x, int per_y, int per_z)
261 int bx0, bx1, by0, by1, bz0, bz1;
262 int b00, b10, b01, b11;
263 float rx0, rx1, ry0, ry1, rz0, rz1;
265 float u, v, a, b, c, d;
269 setup_p(x, bx0, bx1, rx0, rx1, per_x);
270 setup_p(y, by0, by1, ry0, ry1, per_y);
271 setup_p(z, bz0, bz1, rz0, rz1, per_z);
281 /* calculate hermite interpolating factors */
286 /* interpolate along the top slice of the cell */
287 u = grad3[b00 + bz0][0] * rx0 + grad3[b00 + bz0][1] * ry0 + grad3[b00 + bz0][2] * rz0;
288 v = grad3[b10 + bz0][0] * rx1 + grad3[b10 + bz0][1] * ry0 + grad3[b10 + bz0][2] * rz0;
291 u = grad3[b01 + bz0][0] * rx0 + grad3[b01 + bz0][1] * ry1 + grad3[b01 + bz0][2] * rz0;
292 v = grad3[b11 + bz0][0] * rx1 + grad3[b11 + bz0][1] * ry1 + grad3[b11 + bz0][2] * rz0;
297 /* interpolate along the bottom slice of the cell */
298 u = grad3[b00 + bz0][0] * rx0 + grad3[b00 + bz0][1] * ry0 + grad3[b00 + bz0][2] * rz1;
299 v = grad3[b10 + bz0][0] * rx1 + grad3[b10 + bz0][1] * ry0 + grad3[b10 + bz0][2] * rz1;
302 u = grad3[b01 + bz0][0] * rx0 + grad3[b01 + bz0][1] * ry1 + grad3[b01 + bz0][2] * rz1;
303 v = grad3[b11 + bz0][0] * rx1 + grad3[b11 + bz0][1] * ry1 + grad3[b11 + bz0][2] * rz1;
308 /* interpolate between slices */
309 return LERP(c, d, sz);
313 float fbm1(float x, int octaves)
316 float res = 0.0f, freq = 1.0f;
317 for(i=0; i<octaves; i++) {
318 res += noise1(x * freq) / freq;
324 float fbm2(float x, float y, int octaves)
327 float res = 0.0f, freq = 1.0f;
328 for(i=0; i<octaves; i++) {
329 res += noise2(x * freq, y * freq) / freq;
335 float fbm3(float x, float y, float z, int octaves)
338 float res = 0.0f, freq = 1.0f;
339 for(i=0; i<octaves; i++) {
340 res += noise3(x * freq, y * freq, z * freq) / freq;
348 float pfbm1(float x, int per, int octaves)
351 float res = 0.0f, freq = 1.0f;
352 for(i=0; i<octaves; i++) {
353 res += pnoise1(x * freq, per) / freq;
360 float pfbm2(float x, float y, int per_x, int per_y, int octaves)
363 float res = 0.0f, freq = 1.0f;
364 for(i=0; i<octaves; i++) {
365 res += pnoise2(x * freq, y * freq, per_x, per_y) / freq;
373 float pfbm3(float x, float y, float z, int per_x, int per_y, int per_z, int octaves)
376 float res = 0.0f, freq = 1.0f;
377 for(int i=0; i<octaves; i++) {
378 res += pnoise3(x * freq, y * freq, z * freq, per_x, per_y, per_z) / freq;
388 float turbulence1(float x, int octaves)
391 float res = 0.0f, freq = 1.0f;
392 for(i=0; i<octaves; i++) {
393 res += fabs(noise1(x * freq) / freq);
399 float turbulence2(float x, float y, int octaves)
402 float res = 0.0f, freq = 1.0f;
403 for(i=0; i<octaves; i++) {
404 res += fabs(noise2(x * freq, y * freq) / freq);
410 float turbulence3(float x, float y, float z, int octaves)
413 float res = 0.0f, freq = 1.0f;
414 for(i=0; i<octaves; i++) {
415 res += fabs(noise3(x * freq, y * freq, z * freq) / freq);
422 float pturbulence1(float x, int per, int octaves)
425 float res = 0.0f, freq = 1.0f;
426 for(i=0; i<octaves; i++) {
427 res += fabs(pnoise1(x * freq, per) / freq);
434 float pturbulence2(float x, float y, int per_x, int per_y, int octaves)
437 float res = 0.0f, freq = 1.0f;
438 for(i=0; i<octaves; i++) {
439 res += fabs(pnoise2(x * freq, y * freq, per_x, per_y) / freq);
447 float pturbulence3(float x, float y, float z, int per_x, int per_y, int per_z, int octaves)
450 float res = 0.0f, freq = 1.0f;
451 for(i=0; i<octaves; i++) {
452 res += fabs(pnoise3(x * freq, y * freq, z * freq, per_x, per_y, per_z) / freq);