pd.if.org Git - noise/blob - 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