initial commit
[noise] / simplex.c
1 #include <math.h>
2
3 static double 
4 dot3d(int g[], double x, double y, double z)
5 {
6         return g[0] * x + g[1] * y + g[2] * z;
7 }
8
9 static double 
10 dot2d(int g[], double x, double y)
11 {
12         return g[0] * x + g[1] * y;
13 }
14
15 static double 
16 dot4d(int g[], double x, double y, double z, double w)
17 {
18         return g[0] * x + g[1] * y + g[2] * z + g[3] * w;
19 }
20
21 static double 
22 fade(double t)
23 {
24         return t * t * t * (t * (t * 6 - 15) + 10);
25 }
26
27 static double 
28 mix(double a, double b, double t)
29 {
30         return (1 - t) * a + t * b;
31 }
32
33 static int      grad3[][3] = {
34         {1, 1, 0}, {-1, 1, 0}, {1, -1, 0}, {-1, -1, 0},
35         {1, 0, 1}, {-1, 0, 1}, {1, 0, -1}, {-1, 0, -1},
36         {0, 1, 1}, {0, -1, 1}, {0, 1, -1}, {0, -1, -1}
37 };
38
39 static int      grad4[][4] = {{0, 1, 1, 1}, {0, 1, 1, -1}, {0, 1, -1, 1}, {0, 1, -1, -1},
40 {0, -1, 1, 1}, {0, -1, 1, -1}, {0, -1, -1, 1}, {0, -1, -1, -1},
41 {1, 0, 1, 1}, {1, 0, 1, -1}, {1, 0, -1, 1}, {1, 0, -1, -1},
42 {-1, 0, 1, 1}, {-1, 0, 1, -1}, {-1, 0, -1, 1}, {-1, 0, -1, -1},
43 {1, 1, 0, 1}, {1, 1, 0, -1}, {1, -1, 0, 1}, {1, -1, 0, -1},
44 {-1, 1, 0, 1}, {-1, 1, 0, -1}, {-1, -1, 0, 1}, {-1, -1, 0, -1},
45 {1, 1, 1, 0}, {1, 1, -1, 0}, {1, -1, 1, 0}, {1, -1, -1, 0},
46 {-1, 1, 1, 0}, {-1, 1, -1, 0}, {-1, -1, 1, 0}, {-1, -1, -1, 0}
47 };
48
49 /* Classic Perlin noise in 3D, for comparison */
50
51 static int      perm[] = {
52         151, 160, 137, 91, 90, 15,
53         131, 13, 201, 95, 96, 53, 194, 233, 7, 225, 140, 36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23,
54         190, 6, 148, 247, 120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 11, 32, 57, 177, 33,
55         88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 68, 175, 74, 165, 71, 134, 139, 48, 27, 166,
56         77, 146, 158, 231, 83, 111, 229, 122, 60, 211, 133, 230, 220, 105, 92, 41, 55, 46, 245, 40, 244,
57         102, 143, 54, 65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 187, 208, 89, 18, 169, 200, 196,
58         135, 130, 116, 188, 159, 86, 164, 100, 109, 198, 173, 186, 3, 64, 52, 217, 226, 250, 124, 123,
59         5, 202, 38, 147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58, 17, 182, 189, 28, 42,
60         223, 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101, 155, 167, 43, 172, 9,
61         129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, 112, 104, 218, 246, 97, 228,
62         251, 34, 242, 193, 238, 210, 144, 12, 191, 179, 162, 241, 81, 51, 145, 235, 249, 14, 239, 107,
63         49, 192, 214, 31, 181, 199, 106, 157, 184, 84, 204, 176, 115, 121, 50, 45, 127, 4, 150, 254,
64         138, 236, 205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156, 180,
65
66         /*
67          * To remove the need for index wrapping, double the permutation
68          * table length
69          */
70
71         151, 160, 137, 91, 90, 15,
72         131, 13, 201, 95, 96, 53, 194, 233, 7, 225, 140, 36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23,
73         190, 6, 148, 247, 120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 11, 32, 57, 177, 33,
74         88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 68, 175, 74, 165, 71, 134, 139, 48, 27, 166,
75         77, 146, 158, 231, 83, 111, 229, 122, 60, 211, 133, 230, 220, 105, 92, 41, 55, 46, 245, 40, 244,
76         102, 143, 54, 65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 187, 208, 89, 18, 169, 200, 196,
77         135, 130, 116, 188, 159, 86, 164, 100, 109, 198, 173, 186, 3, 64, 52, 217, 226, 250, 124, 123,
78         5, 202, 38, 147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58, 17, 182, 189, 28, 42,
79         223, 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101, 155, 167, 43, 172, 9,
80         129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, 112, 104, 218, 246, 97, 228,
81         251, 34, 242, 193, 238, 210, 144, 12, 191, 179, 162, 241, 81, 51, 145, 235, 249, 14, 239, 107,
82         49, 192, 214, 31, 181, 199, 106, 157, 184, 84, 204, 176, 115, 121, 50, 45, 127, 4, 150, 254,
83         138, 236, 205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156, 180
84 };
85
86 /* Classic Perlin noise, 3D version */
87 static double 
88 noise3dp(double x, double y, double z)
89 {
90         /* Find unit grid cell containing point */
91         int             X = floor(x);
92         int             Y = floor(y);
93         int             Z = floor(z);
94
95         /* Get relative xyz coordinates of point within that cell */
96         x = x - X;
97         y = y - Y;
98         z = z - Z;
99
100         /*
101          * Wrap the integer cells at 255 (smaller integer period can be
102          * introduced here)
103          */
104         X = X & 255;
105         Y = Y & 255;
106         Z = Z & 255;
107
108         /* Calculate a set of eight hashed gradient indices */
109         int             gi000 = perm[X + perm[Y + perm[Z]]] % 12;
110         int             gi001 = perm[X + perm[Y + perm[Z + 1]]] % 12;
111         int             gi010 = perm[X + perm[Y + 1 + perm[Z]]] % 12;
112         int             gi011 = perm[X + perm[Y + 1 + perm[Z + 1]]] % 12;
113         int             gi100 = perm[X + 1 + perm[Y + perm[Z]]] % 12;
114         int             gi101 = perm[X + 1 + perm[Y + perm[Z + 1]]] % 12;
115         int             gi110 = perm[X + 1 + perm[Y + 1 + perm[Z]]] % 12;
116         int             gi111 = perm[X + 1 + perm[Y + 1 + perm[Z + 1]]] % 12;
117
118         /*
119          * The gradients of each corner are now:
120          * g000 = grad3[gi000];
121          * g001 = grad3[gi001];
122          * g010 = grad3[gi010];
123          * g011 = grad3[gi011];
124          * g100 = grad3[gi100];
125          * g101 = grad3[gi101];
126          * g110 = grad3[gi110];
127          * g111 = grad3[gi111];
128          */
129
130         /* Calculate noise contributions from each of the eight corners */
131         double          n000 = dot3d(grad3[gi000], x, y, z);
132         double          n100 = dot3d(grad3[gi100], x - 1, y, z);
133         double          n010 = dot3d(grad3[gi010], x, y - 1, z);
134         double          n110 = dot3d(grad3[gi110], x - 1, y - 1, z);
135         double          n001 = dot3d(grad3[gi001], x, y, z - 1);
136         double          n101 = dot3d(grad3[gi101], x - 1, y, z - 1);
137         double          n011 = dot3d(grad3[gi011], x, y - 1, z - 1);
138         double          n111 = dot3d(grad3[gi111], x - 1, y - 1, z - 1);
139
140         /* Compute the fade curve value for each of x, y, z */
141         double          u = fade(x);
142         double          v = fade(y);
143         double          w = fade(z);
144
145         /* Interpolate along x the contributions from each of the corners */
146         double          nx00 = mix(n000, n100, u);
147         double          nx01 = mix(n001, n101, u);
148         double          nx10 = mix(n010, n110, u);
149         double          nx11 = mix(n011, n111, u);
150
151         /* Interpolate the four results along y */
152         double          nxy0 = mix(nx00, nx10, v);
153         double          nxy1 = mix(nx01, nx11, v);
154
155         /* Interpolate the two last results along z */
156         double          nxyz = mix(nxy0, nxy1, w);
157
158         return nxyz;
159 }
160
161 /*
162  * A lookup table to traverse the simplex around a given point in 4D.
163  * Details can be found where this table is used, in the 4D noise method.
164  */
165 int             simplex[][4] = {
166         {0, 1, 2, 3}, {0, 1, 3, 2}, {0, 0, 0, 0}, {0, 2, 3, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 2, 3, 0},
167         {0, 2, 1, 3}, {0, 0, 0, 0}, {0, 3, 1, 2}, {0, 3, 2, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 3, 2, 0},
168         {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
169         {1, 2, 0, 3}, {0, 0, 0, 0}, {1, 3, 0, 2}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {2, 3, 0, 1}, {2, 3, 1, 0},
170         {1, 0, 2, 3}, {1, 0, 3, 2}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {2, 0, 3, 1}, {0, 0, 0, 0}, {2, 1, 3, 0},
171         {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
172         {2, 0, 1, 3}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {3, 0, 1, 2}, {3, 0, 2, 1}, {0, 0, 0, 0}, {3, 1, 2, 0},
173 {2, 1, 0, 3}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {3, 1, 0, 2}, {0, 0, 0, 0}, {3, 2, 0, 1}, {3, 2, 1, 0}};
174
175 double grad1(int hash) {
176         double g;
177         int h = hash & 15;
178         g = 1.0 + (h & 7);   /* Gradient value is one of 1.0, 2.0, ..., 8.0 */
179         if (h&8) g = -g;   /* Make half of the gradients negative */
180         return g;
181 }
182
183 double noise1d(double x, double *dx)
184 {
185         int i0 = floor(x);
186         int i1 = i0 + 1;
187         double x0 = x - i0;
188         double x1 = x0 - 1.0;
189
190         double gx0, gx1;
191         double n0, n1;
192         double t20, t40, t21, t41;
193
194         double x20 = x0*x0;
195         double t0 = 1.0 - x20;
196
197         /*  if(t0 < 0.0f) t0 = 0.0f; /* Never happens for 1D: x0<=1 always */
198         t20 = t0 * t0;
199         t40 = t20 * t20;
200         gx0 = grad1(perm[i0 & 0xff]);
201         n0 = t40 * gx0 * x0;
202
203         double x21 = x1*x1;
204         double t1 = 1.0 - x21;
205         /*  if(t1 < 0.0f) t1 = 0.0f; /* Never happens for 1D: |x1|<=1 always */
206         t21 = t1 * t1;
207         t41 = t21 * t21;
208         gx1 = grad1(perm[i1 & 0xff]);
209         n1 = t41 * gx1 * x1;
210
211         /*
212          * Compute derivative according to:
213          * *dnoise_dx = -8.0f * t20 * t0 * x0 * (gx0 * x0) + t40 * gx0;
214          * *dnoise_dx += -8.0f * t21 * t1 * x1 * (gx1 * x1) + t41 * gx1;
215          */
216
217         if (dx) {
218                 *dx = t20 * t0 * gx0 * x20;
219                 *dx += t21 * t1 * gx1 * x21;
220                 *dx *= -8.0;
221                 *dx += t40 * gx0 + t41 * gx1;
222                 *dx *= 0.25; /* Scale derivative to match the noise scaling */
223         }
224
225         /*
226          * The maximum value of this noise is 8*(3/4)^4 = 2.53125
227          * A factor of 0.395 would scale to fit exactly within [-1,1], but
228          * to better match classic Perlin noise, we scale it down some more.
229          * scale = 1.0 / 2.53125
230          */
231         return (n0 + n1)/2.53125;
232 }
233
234 /* 2D simplex noise */
235 static double 
236 noise2d(double xin, double yin)
237 {
238         double          n0, n1, n2;     /* Noise contributions from the three
239                                          * corners */
240
241         /* Skew the input space to determine which simplex cell we're in */
242         double          F2 = 0.5 * (sqrt(3.0) - 1.0);
243         double          s = (xin + yin) * F2;   /* Hairy factor for 2D */
244
245         int             i = floor(xin + s);
246         int             j = floor(yin + s);
247
248         double          G2 = (3.0 - sqrt(3.0)) / 6.0;
249         double          t = (i + j) * G2;
250         double          X0 = i - t;     /* Unskew the cell origin back to
251                                          * (x,y) space */
252         double          Y0 = j - t;
253         double          x0 = xin - X0;  /* The x,y distances from the cell
254                                          * origin */
255         double          y0 = yin - Y0;
256
257         /* For the 2D case, the simplex shape is an equilateral triangle. */
258
259         /* Determine which simplex we are in. */
260         int             i1, j1; /* Offsets for second (middle) corner of
261                                  * simplex in (i,j) coords */
262         if (x0 > y0) {
263                 /* lower triangle, XY order: (0,0)->(1,0)->(1,1) */
264                 i1 = 1;
265                 j1 = 0;
266         } else {
267                 /* upper triangle, YX order: (0,0)->(0,1)->(1,1) */
268                 i1 = 0;
269                 j1 = 1;
270         }
271         /*
272          * A step of (1,0) in (i,j) means a step of (1-c,-c) in (x,y), and
273          * a step of (0,1) in (i,j) means a step of (-c,1-c) in (x,y), where
274          * c = (3-sqrt(3))/6
275          */
276         double          x1 = x0 - i1 + G2;      /* Offsets for middle corner
277                                                  * in (x,y) unskewed coords */
278         double          y1 = y0 - j1 + G2;
279         double          x2 = x0 - 1.0 + 2.0 * G2;       /* Offsets for last
280                                                          * corner in (x,y)
281                                                          * unskewed coords */
282         double          y2 = y0 - 1.0 + 2.0 * G2;
283
284         /* Work out the hashed gradient indices of the three simplex corners */
285         int             ii = i & 255;
286         int             jj = j & 255;
287         int             gi0 = perm[ii + perm[jj]] % 12;
288         int             gi1 = perm[ii + i1 + perm[jj + j1]] % 12;
289         int             gi2 = perm[ii + 1 + perm[jj + 1]] % 12;
290
291         /* Calculate the contribution from the three corners */
292         double          t0 = 0.5 - x0 * x0 - y0 * y0;
293         if (t0 < 0) {
294                 n0 = 0.0;
295         } else {
296                 t0 *= t0;
297                 n0 = t0 * t0 * dot2d(grad3[gi0], x0, y0);       /* (x,y) of grad3 used
298                                                                  * for 2D gradient */
299         }
300
301         double          t1 = 0.5 - x1 * x1 - y1 * y1;
302         if (t1 < 0) {
303                 n1 = 0.0;
304         } else {
305                 t1 *= t1;
306                 n1 = t1 * t1 * dot2d(grad3[gi1], x1, y1);
307         }
308
309         double          t2 = 0.5 - x2 * x2 - y2 * y2;
310         if (t2 < 0) {
311                 n2 = 0.0;
312         } else {
313                 t2 *= t2;
314                 n2 = t2 * t2 * dot2d(grad3[gi2], x2, y2);
315         }
316
317         /*
318          * Add contributions from each corner to get the final noise value.
319          * The result is scaled to return values in the interval [-1,1].
320          */
321         return 70.0 * (n0 + n1 + n2);
322 }
323
324 /* 3D simplex noise */
325 static double 
326 noise3d(double xin, double yin, double zin)
327 {
328         double          n0, n1, n2, n3; /* Noise contributions from the four
329                                          * corners */
330         /* Skew the input space to determine which simplex cell we're in */
331         double          F3 = 1.0 / 3.0;
332         double          s = (xin + yin + zin) * F3;     /* Very nice and simple
333                                                          * skew factor for 3D */
334         int             i = floor(xin + s);
335         int             j = floor(yin + s);
336         int             k = floor(zin + s);
337         double          G3 = 1.0 / 6.0; /* Very nice and simple unskew
338                                          * factor, too */
339         double          t = (i + j + k) * G3;
340         double          X0 = i - t;     /* Unskew the cell origin back to
341                                          * (x,y,z) space */
342         double          Y0 = j - t;
343         double          Z0 = k - t;
344         double          x0 = xin - X0;  /* The x,y,z distances from the cell
345                                          * origin */
346         double          y0 = yin - Y0;
347         double          z0 = zin - Z0;
348
349         /*
350          * For the 3D case, the simplex shape is a slightly irregular
351          * tetrahedron.
352          */
353
354         /* Determine which simplex we are in. */
355
356         int             i1, j1, k1;     /* Offsets for second corner of
357                                          * simplex in (i,j,k) coords */
358
359         int             i2, j2, k2;     /* Offsets for third corner of
360                                          * simplex in (i,j,k) coords */
361         if (x0 >= y0) {
362                 if (y0 >= z0) {
363                         i1 = 1;
364                         j1 = 0;
365                         k1 = 0;
366                         i2 = 1;
367                         j2 = 1;
368                         k2 = 0;
369                 }
370                  /* X Y Z order */ 
371                 else if (x0 >= z0) {
372                         i1 = 1;
373                         j1 = 0;
374                         k1 = 0;
375                         i2 = 1;
376                         j2 = 0;
377                         k2 = 1;
378                 }
379                  /* X Z Y order */ 
380                 else {
381                         i1 = 0;
382                         j1 = 0;
383                         k1 = 1;
384                         i2 = 1;
385                         j2 = 0;
386                         k2 = 1;
387                 }               /* Z X Y order */
388         } else {                /* x0<y0 */
389                 if (y0 < z0) {
390                         i1 = 0;
391                         j1 = 0;
392                         k1 = 1;
393                         i2 = 0;
394                         j2 = 1;
395                         k2 = 1;
396                 }
397                  /* Z Y X order */ 
398                 else if (x0 < z0) {
399                         i1 = 0;
400                         j1 = 1;
401                         k1 = 0;
402                         i2 = 0;
403                         j2 = 1;
404                         k2 = 1;
405                 }
406                  /* Y Z X order */ 
407                 else {
408                         i1 = 0;
409                         j1 = 1;
410                         k1 = 0;
411                         i2 = 1;
412                         j2 = 1;
413                         k2 = 0;
414                 }               /* Y X Z order */
415         }
416
417         /*
418          * A step of (1,0,0) in (i,j,k) means a step of (1-c,-c,-c) in
419          * (x,y,z),
420          */
421
422         /*
423          * a step of (0,1,0) in (i,j,k) means a step of (-c,1-c,-c) in
424          * (x,y,z), and
425          */
426
427         /*
428          * a step of (0,0,1) in (i,j,k) means a step of (-c,-c,1-c) in
429          * (x,y,z), where
430          */
431
432         /* c = 1/6. */
433         double          x1 = x0 - i1 + G3;      /* Offsets for second corner
434                                                  * in (x,y,z) coords */
435         double          y1 = y0 - j1 + G3;
436         double          z1 = z0 - k1 + G3;
437         double          x2 = x0 - i2 + 2.0 * G3;        /* Offsets for third
438                                                          * corner in (x,y,z)
439                                                          * coords */
440         double          y2 = y0 - j2 + 2.0 * G3;
441         double          z2 = z0 - k2 + 2.0 * G3;
442         double          x3 = x0 - 1.0 + 3.0 * G3;       /* Offsets for last
443                                                          * corner in (x,y,z)
444                                                          * coords */
445         double          y3 = y0 - 1.0 + 3.0 * G3;
446         double          z3 = z0 - 1.0 + 3.0 * G3;
447
448         /* Work out the hashed gradient indices of the four simplex corners */
449         int             ii = i & 255;
450         int             jj = j & 255;
451         int             kk = k & 255;
452         int             gi0 = perm[ii + perm[jj + perm[kk]]] % 12;
453         int             gi1 = perm[ii + i1 + perm[jj + j1 + perm[kk + k1]]] % 12;
454         int             gi2 = perm[ii + i2 + perm[jj + j2 + perm[kk + k2]]] % 12;
455         int             gi3 = perm[ii + 1 + perm[jj + 1 + perm[kk + 1]]] % 12;
456
457         /* Calculate the contribution from the four corners */
458         double          t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0;
459         if (t0 < 0)
460                 n0 = 0.0;
461         else {
462                 t0 *= t0;
463                 n0 = t0 * t0 * dot3d(grad3[gi0], x0, y0, z0);
464         }
465
466         double          t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1;
467         if (t1 < 0) {
468                 n1 = 0.0;
469         }
470         else {
471                 t1 *= t1;
472                 n1 = t1 * t1 * dot3d(grad3[gi1], x1, y1, z1);
473         }
474
475         double          t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2;
476         if (t2 < 0) {
477                 n2 = 0.0;
478         }
479         else {
480                 t2 *= t2;
481                 n2 = t2 * t2 * dot3d(grad3[gi2], x2, y2, z2);
482         }
483
484         double          t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3;
485         if (t3 < 0) {
486                 n3 = 0.0;
487         }
488         else {
489                 t3 *= t3;
490                 n3 = t3 * t3 * dot3d(grad3[gi3], x3, y3, z3);
491         }
492         /* Add contributions from each corner to get the final noise value. */
493         /* The result is scaled to stay just inside [-1,1] */
494         return 32.0 * (n0 + n1 + n2 + n3);
495 }
496
497 /* 4D simplex noise */
498 double 
499 noise4d(double x, double y, double z, double w)
500 {
501
502         /* The skewing and unskewing factors are hairy again for the 4D case */
503         double          F4 = (sqrt(5.0) - 1.0) / 4.0;
504         double          G4 = (5.0 - sqrt(5.0)) / 20.0;
505         double          n0, n1, n2, n3, n4;     /* Noise contributions from
506                                                  * the five corners */
507         /*
508          * Skew the (x,y,z,w) space to determine which cell of 24 simplices
509          * we're in
510          */
511         double          s = (x + y + z + w) * F4;       /* Factor for 4D skewing */
512         int             i = floor(x + s);
513         int             j = floor(y + s);
514         int             k = floor(z + s);
515         int             l = floor(w + s);
516         double          t = (i + j + k + l) * G4;       /* Factor for 4D
517                                                          * unskewing */
518         double          X0 = i - t;     /* Unskew the cell origin back to
519                                          * (x,y,z,w) space */
520         double          Y0 = j - t;
521         double          Z0 = k - t;
522         double          W0 = l - t;
523         double          x0 = x - X0;    /* The x,y,z,w distances from the
524                                          * cell origin */
525         double          y0 = y - Y0;
526         double          z0 = z - Z0;
527         double          w0 = w - W0;
528         /*
529          * For the 4D case, the simplex is a 4D shape I won't even try to
530          * describe.
531          */
532         /*
533          * To find out which of the 24 possible simplices we're in, we need
534          * to
535          */
536         /* determine the magnitude ordering of x0, y0, z0 and w0. */
537         /*
538          * The method below is a good way of finding the ordering of x,y,z,w
539          * and
540          */
541         /* then find the correct traversal order for the simplex we’re in. */
542         /*
543          * First, six pair-wise comparisons are performed between each
544          * possible pair
545          */
546         /*
547          * of the four coordinates, and the results are used to add up binary
548          * bits
549          */
550         /* for an integer index. */
551         int             c1 = (x0 > y0) ? 32 : 0;
552         int             c2 = (x0 > z0) ? 16 : 0;
553         int             c3 = (y0 > z0) ? 8 : 0;
554         int             c4 = (x0 > w0) ? 4 : 0;
555         int             c5 = (y0 > w0) ? 2 : 0;
556         int             c6 = (z0 > w0) ? 1 : 0;
557         int             c = c1 + c2 + c3 + c4 + c5 + c6;
558         int             i1, j1, k1, l1; /* The integer offsets for the second
559                                          * simplex corner */
560         int             i2, j2, k2, l2; /* The integer offsets for the third
561                                          * simplex corner */
562         int             i3, j3, k3, l3; /* The integer offsets for the fourth
563                                          * simplex corner */
564         /*
565          * simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some
566          * order.
567          */
568         /*
569          * Many values of c will never occur, since e.g. x>y>z>w makes x<z,
570          * y<w and x<w
571          */
572         /*
573          * impossible. Only the 24 indices which have non-zero entries make
574          * any sense.
575          */
576         /*
577          * We use a thresholding to set the coordinates in turn from the
578          * largest magnitude.
579          */
580         /*
581          * The number 3 in the "simplex" array is at the position of the
582          * largest coordinate.
583          */
584         i1 = simplex[c][0] >= 3 ? 1 : 0;
585         j1 = simplex[c][1] >= 3 ? 1 : 0;
586         k1 = simplex[c][2] >= 3 ? 1 : 0;
587         l1 = simplex[c][3] >= 3 ? 1 : 0;
588         /*
589          * The number 2 in the "simplex" array is at the second largest
590          * coordinate.
591          */
592         i2 = simplex[c][0] >= 2 ? 1 : 0;
593         j2 = simplex[c][1] >= 2 ? 1 : 0;
594         k2 = simplex[c][2] >= 2 ? 1 : 0;
595         l2 = simplex[c][3] >= 2 ? 1 : 0;
596         /*
597          * The number 1 in the "simplex" array is at the second smallest
598          * coordinate.
599          */
600         i3 = simplex[c][0] >= 1 ? 1 : 0;
601         j3 = simplex[c][1] >= 1 ? 1 : 0;
602         k3 = simplex[c][2] >= 1 ? 1 : 0;
603         l3 = simplex[c][3] >= 1 ? 1 : 0;
604         /*
605          * The fifth corner has all coordinate offsets = 1, so no need to
606          * look that up.
607          */
608         double          x1 = x0 - i1 + G4;      /* Offsets for second corner
609                                                  * in (x,y,z,w) coords */
610         double          y1 = y0 - j1 + G4;
611         double          z1 = z0 - k1 + G4;
612         double          w1 = w0 - l1 + G4;
613         double          x2 = x0 - i2 + 2.0 * G4;        /* Offsets for third
614                                                          * corner in (x,y,z,w)
615                                                          * coords */
616         double          y2 = y0 - j2 + 2.0 * G4;
617         double          z2 = z0 - k2 + 2.0 * G4;
618         double          w2 = w0 - l2 + 2.0 * G4;
619         double          x3 = x0 - i3 + 3.0 * G4;        /* Offsets for fourth
620                                                          * corner in (x,y,z,w)
621                                                          * coords */
622         double          y3 = y0 - j3 + 3.0 * G4;
623         double          z3 = z0 - k3 + 3.0 * G4;
624         double          w3 = w0 - l3 + 3.0 * G4;
625         double          x4 = x0 - 1.0 + 4.0 * G4;       /* Offsets for last
626                                                          * corner in (x,y,z,w)
627                                                          * coords */
628         double          y4 = y0 - 1.0 + 4.0 * G4;
629         double          z4 = z0 - 1.0 + 4.0 * G4;
630         double          w4 = w0 - 1.0 + 4.0 * G4;
631         /* Work out the hashed gradient indices of the five simplex corners */
632         int             ii = i & 255;
633         int             jj = j & 255;
634         int             kk = k & 255;
635         int             ll = l & 255;
636         int             gi0 = perm[ii + perm[jj + perm[kk + perm[ll]]]] % 32;
637         int             gi1 = perm[ii + i1 + perm[jj + j1 + perm[kk + k1 + perm[ll + l1]]]] % 32;
638         int             gi2 = perm[ii + i2 + perm[jj + j2 + perm[kk + k2 + perm[ll + l2]]]] % 32;
639         int             gi3 = perm[ii + i3 + perm[jj + j3 + perm[kk + k3 + perm[ll + l3]]]] % 32;
640         int             gi4 = perm[ii + 1 + perm[jj + 1 + perm[kk + 1 + perm[ll + 1]]]] % 32;
641         /* Calculate the contribution from the five corners */
642         double          t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0 - w0 * w0;
643         if (t0 < 0) {
644                 n0 = 0.0;
645         }
646         else {
647                 t0 *= t0;
648                 n0 = t0 * t0 * dot4d(grad4[gi0], x0, y0, z0, w0);
649         }
650         double          t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1 - w1 * w1;
651         if (t1 < 0) {
652                 n1 = 0.0;
653         }
654         else {
655                 t1 *= t1;
656                 n1 = t1 * t1 * dot4d(grad4[gi1], x1, y1, z1, w1);
657         }
658         double          t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2 - w2 * w2;
659         if (t2 < 0) {
660                 n2 = 0.0;
661         }
662         else {
663                 t2 *= t2;
664                 n2 = t2 * t2 * dot4d(grad4[gi2], x2, y2, z2, w2);
665         }
666         double          t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3 - w3 * w3;
667         if (t3 < 0) {
668                 n3 = 0.0;
669         }
670         else {
671                 t3 *= t3;
672                 n3 = t3 * t3 * dot4d(grad4[gi3], x3, y3, z3, w3);
673         }
674         double          t4 = 0.6 - x4 * x4 - y4 * y4 - z4 * z4 - w4 * w4;
675         if (t4 < 0) {
676                 n4 = 0.0;
677         }
678         else {
679                 t4 *= t4;
680                 n4 = t4 * t4 * dot4d(grad4[gi4], x4, y4, z4, w4);
681         }
682
683         /* Sum up and scale the result to cover the range [-1,1] */
684         return 27.0 * (n0 + n1 + n2 + n3 + n4);
685 }
686
687 double noise_simplex(int dim, double scale, double *pt) {
688         switch (dim) {
689                 case 1:
690                         return noise1d(pt[0]*scale,0);
691                         break;
692                 case 2:
693                         return noise2d(pt[0]*scale,pt[1]*scale);
694                         break;
695                 case 3:
696                         return noise3d(pt[0]*scale,pt[1]*scale,pt[2]*scale);
697                         break;
698                 case 4:
699                         return noise4d(pt[0]*scale,pt[1]*scale,pt[2]*scale,pt[3]*scale);
700                         break;
701                 default:
702                         break;
703         }
704         return 0.0;
705 }
706
707 double noise_fractal(int dim, double *pt, int octaves, double lacunarity, double persistence) {
708         double noise,scale,lac;
709         noise = 0.0;
710         scale = 1.0;
711         lac = 1.0;
712
713         for (scale=1.0; octaves--;scale*=persistence) {
714                 noise += noise_simplex(dim, lac, pt) * scale;
715                 scale *= persistence;
716                 lac *= lacunarity;
717         }
718         return noise;
719 }
720
721 #ifdef MAIN
722
723 #include <stdio.h>
724 #include <stdlib.h>
725 #include <float.h>
726
727 int 
728 main(int argc, char *argv[])
729 {
730         double  r;
731         double  base;
732         double  s[2];
733
734         int x,y;
735
736         int octaves;
737         int pos = 0, neg = 0;
738         double max, min;
739 #ifdef INFINITY
740         max = -INFINITY;
741         min = +INFINITY;
742 #else
743         max = DBL_MIN;
744         min = DBL_MAX;
745 #endif
746
747         octaves = argc > 1 ? atoi(argv[1]) : 1;
748         base = 40.0;
749
750         for (x=0;x<256;x++) {
751                 s[0] = x/base;
752                 for (y=0;y<256;y++) {
753                         s[1] = y/base;
754                         r = noise_fractal(2, s, octaves, 2.0, 0.7);
755                         r /= 1.5;
756
757                         /* grain */
758                         r = noise_simplex(2,0.20,s)*10;
759                         r = (r - (int)r);
760                         /* fine */
761                         r += noise_simplex(2, 75.0, s) * 0.20;
762 #if 1
763                         /* streak */
764                         r += noise2d(s[0]*0.45, s[1]*200.0) * 0.3;
765 #endif
766
767                         r /= 2.5;
768                         printf("%f\n", r);
769                         if (r > 1.0) {
770                                 pos++;
771                         } else if (r < -1.0) {
772                                 neg++;
773                         }
774                         if (r > max) max = r;
775                         if (r < min) min = r;
776
777                 }
778         }
779         fprintf(stderr, "max value: %f\n", max);
780         fprintf(stderr, "min value: %f\n", min);
781         if (pos) fprintf(stderr, "%d values > +1.0\n", pos);
782         if (neg) fprintf(stderr, "%d values < -1.0\n", neg);
783         
784         return 0;
785
786         switch (--argc) {
787                 case 1:
788                 r = noise1d(atof(argv[1]),0);
789                 break;
790
791         case 2:
792                 r = noise2d(atof(argv[1]), atof(argv[2]));
793                 break;
794         case 3:
795                 r = noise3d(atof(argv[1]), atof(argv[2]), atof(argv[3]));
796                 break;
797         case 4:
798                 r = noise4d(atof(argv[1]), atof(argv[2]), atof(argv[3]), atof(argv[4]));
799                 break;
800         default:
801                 fprintf(stderr, "Wrong number of arguments: %d, must be 2, 3, or 4\n", argc);
802                 exit(EXIT_FAILURE);
803                 break;
804         };
805
806         printf("%f\n", r);
807
808         return 0;
809 }
810
811 #endif