001package squidpony.squidmath; 002 003/** 004 * This is Ken Perlin's third revision of his noise function. It is sometimes 005 * referred to as "Simplex Noise". Results are bound by (-1, 1) inclusive. 006 * 007 * 008 * It is significantly faster than his earlier versions. This particular version 009 * was originally from Stefan Gustavson. This is much preferred to the earlier 010 * versions of Perlin Noise due to the reasons noted in the articles: 011 * <ul> 012 * <li>http://www.itn.liu.se/~stegu/simplexnoise/simplexnoise.pdf</li> 013 * <li>http://webstaff.itn.liu.se/~stegu/TNM022-2005/perlinnoiselinks/ch02.pdf</li> 014 * </ul> 015 * 016 */ 017public class PerlinNoise { 018 019 private static final int grad3[][] = {{1, 1, 0}, {-1, 1, 0}, {1, -1, 0}, 020 {-1, -1, 0}, {1, 0, 1}, {-1, 0, 1}, 021 {1, 0, -1}, {-1, 0, -1}, {0, 1, 1}, 022 {0, -1, 1}, {0, 1, -1}, {0, -1, -1}}; 023 private static final int grad4[][] = {{0, 1, 1, 1}, {0, 1, 1, -1}, 024 {0, 1, -1, 1}, {0, 1, -1, -1}, 025 {0, -1, 1, 1}, {0, -1, 1, -1}, 026 {0, -1, -1, 1}, {0, -1, -1, -1}, 027 {1, 0, 1, 1}, {1, 0, 1, -1}, 028 {1, 0, -1, 1}, {1, 0, -1, -1}, 029 {-1, 0, 1, 1}, {-1, 0, 1, -1}, 030 {-1, 0, -1, 1}, {-1, 0, -1, -1}, 031 {1, 1, 0, 1}, {1, 1, 0, -1}, 032 {1, -1, 0, 1}, {1, -1, 0, -1}, 033 {-1, 1, 0, 1}, {-1, 1, 0, -1}, 034 {-1, -1, 0, 1}, {-1, -1, 0, -1}, 035 {1, 1, 1, 0}, {1, 1, -1, 0}, 036 {1, -1, 1, 0}, {1, -1, -1, 0}, 037 {-1, 1, 1, 0}, {-1, 1, -1, 0}, 038 {-1, -1, 1, 0}, {-1, -1, -1, 0}}; 039 private static final int p[] = {151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96, 040 53, 194, 233, 7, 225, 140, 36, 103, 30, 69, 142, 041 8, 99, 37, 240, 21, 10, 23, 190, 6, 148, 247, 042 120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 043 117, 35, 11, 32, 57, 177, 33, 88, 237, 149, 56, 044 87, 174, 20, 125, 136, 171, 168, 68, 175, 74, 045 165, 71, 134, 139, 48, 27, 166, 77, 146, 158, 046 231, 83, 111, 229, 122, 60, 211, 133, 230, 220, 047 105, 92, 41, 55, 46, 245, 40, 244, 102, 143, 54, 048 65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 049 187, 208, 89, 18, 169, 200, 196, 135, 130, 116, 050 188, 159, 86, 164, 100, 109, 198, 173, 186, 3, 051 64, 52, 217, 226, 250, 124, 123, 5, 202, 38, 147, 052 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, 053 47, 16, 58, 17, 182, 189, 28, 42, 223, 183, 170, 054 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 055 153, 101, 155, 167, 43, 172, 9, 129, 22, 39, 253, 056 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, 057 112, 104, 218, 246, 97, 228, 251, 34, 242, 193, 058 238, 210, 144, 12, 191, 179, 162, 241, 81, 51, 059 145, 235, 249, 14, 239, 107, 49, 192, 214, 31, 060 181, 199, 106, 157, 184, 84, 204, 176, 115, 121, 061 50, 45, 127, 4, 150, 254, 138, 236, 205, 93, 222, 062 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 063 215, 61, 156, 180}; 064 private static final double F2 = 0.5 * (Math.sqrt(3.0) - 1.0); 065 private static final double G2 = (3.0 - Math.sqrt(3.0)) / 6.0; 066 private static final double F3 = 1.0 / 3.0; 067 private static final double G3 = 1.0 / 6.0; 068 private static final double F4 = (Math.sqrt(5.0) - 1.0) / 4.0; 069 private static final double G4 = (5.0 - Math.sqrt(5.0)) / 20.0; 070 // To remove the need for index wrapping, double the permutation table 071 // length 072 private static final int perm[] = new int[512]; 073 074 static { 075 for (int i = 0; i < 512; i++) { 076 perm[i] = p[i & 255]; 077 } 078 } 079 private PerlinNoise() 080 { 081 082 } 083 // A lookup table to traverse the simplex around a given point in 4D. 084 // Details can be found where this table is used, in the 4D noise method. 085 private static final int simplex[][] 086 = {{0, 1, 2, 3}, {0, 1, 3, 2}, {0, 0, 0, 0}, {0, 2, 3, 1}, 087 {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 2, 3, 0}, 088 {0, 2, 1, 3}, {0, 0, 0, 0}, {0, 3, 1, 2}, {0, 3, 2, 1}, 089 {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 3, 2, 0}, 090 {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, 091 {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, 092 {1, 2, 0, 3}, {0, 0, 0, 0}, {1, 3, 0, 2}, {0, 0, 0, 0}, 093 {0, 0, 0, 0}, {0, 0, 0, 0}, {2, 3, 0, 1}, {2, 3, 1, 0}, 094 {1, 0, 2, 3}, {1, 0, 3, 2}, {0, 0, 0, 0}, {0, 0, 0, 0}, 095 {0, 0, 0, 0}, {2, 0, 3, 1}, {0, 0, 0, 0}, {2, 1, 3, 0}, 096 {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, 097 {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, 098 {2, 0, 1, 3}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, 099 {3, 0, 1, 2}, {3, 0, 2, 1}, {0, 0, 0, 0}, {3, 1, 2, 0}, 100 {2, 1, 0, 3}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, 101 {3, 1, 0, 2}, {0, 0, 0, 0}, {3, 2, 0, 1}, {3, 2, 1, 0}}; 102 103 private static double dot(int g[], double x, double y) { 104 return g[0] * x + g[1] * y; 105 } 106 107 private static double dot(int g[], double x, double y, double z) { 108 return g[0] * x + g[1] * y + g[2] * z; 109 } 110 111 private static double dot(int g[], double x, double y, double z, double w) { 112 return g[0] * x + g[1] * y + g[2] * z + g[3] * w; 113 } 114 115 /** 116 * 2D simplex noise. 117 * 118 * @param xin x input 119 * @param yin y input 120 * @return noise from -1.0 to 1.0, inclusive 121 */ 122 public static double noise(double xin, double yin) { 123 double noise0, noise1, noise2; // from the three corners 124 // Skew the input space to determine which simplex cell we're in 125 double skew = (xin + yin) * F2; // Hairy factor for 2D 126 int i = (int) Math.floor(xin + skew); 127 int j = (int) Math.floor(yin + skew); 128 double t = (i + j) * G2; 129 double X0 = i - t; // Unskew the cell origin back to (x,y) space 130 double Y0 = j - t; 131 double x0 = xin - X0; // The x,y distances from the cell origin 132 double y0 = yin - Y0; 133 // For the 2D case, the simplex shape is an equilateral triangle. 134 // Determine which simplex we are in. 135 int i1, j1; // Offsets for second (middle) corner of simplex in (i,j) 136 // coords 137 if (x0 > y0) { 138 i1 = 1; 139 j1 = 0; 140 } // lower triangle, XY order: (0,0)->(1,0)->(1,1) 141 else { 142 i1 = 0; 143 j1 = 1; 144 } // upper triangle, YX order: (0,0)->(0,1)->(1,1) 145 // A step of (1,0) in (i,j) means a step of (1-c,-c) in (x,y), and 146 // a step of (0,1) in (i,j) means a step of (-c,1-c) in (x,y), 147 // where 148 // c = (3-sqrt(3))/6 149 double x1 = x0 - i1 + G2; // Offsets for middle corner in (x,y) 150 // unskewed coords 151 double y1 = y0 - j1 + G2; 152 double x2 = x0 - 1.0 + 2.0 * G2; // Offsets for last corner in (x,y) 153 // unskewed coords 154 double y2 = y0 - 1.0 + 2.0 * G2; 155 // Work out the hashed gradient indices of the three simplex corners 156 int ii = i & 255; 157 int jj = j & 255; 158 int gi0 = perm[ii + perm[jj]] % 12; 159 int gi1 = perm[ii + i1 + perm[jj + j1]] % 12; 160 int gi2 = perm[ii + 1 + perm[jj + 1]] % 12; 161 // Calculate the contribution from the three corners 162 double t0 = 0.5 - x0 * x0 - y0 * y0; 163 if (t0 < 0) { 164 noise0 = 0.0; 165 } else { 166 t0 *= t0; 167 noise0 = t0 * t0 * dot(grad3[gi0], x0, y0); // (x,y) of grad3 used 168 // for 2D gradient 169 } 170 double t1 = 0.5 - x1 * x1 - y1 * y1; 171 if (t1 < 0) { 172 noise1 = 0.0; 173 } else { 174 t1 *= t1; 175 noise1 = t1 * t1 * dot(grad3[gi1], x1, y1); 176 } 177 double t2 = 0.5 - x2 * x2 - y2 * y2; 178 if (t2 < 0) { 179 noise2 = 0.0; 180 } else { 181 t2 *= t2; 182 noise2 = t2 * t2 * dot(grad3[gi2], x2, y2); 183 } 184 // Add contributions from each corner to get the final noise value. 185 // The result is scaled to return values in the interval [-1,1]. 186 return 70.0 * (noise0 + noise1 + noise2); 187 } 188 189 /** 190 * 3D simplex noise. 191 * 192 * @param xin X input 193 * @param yin Y input 194 * @param zin Z input 195 * @return noise from -1.0 to 1.0, inclusive 196 */ 197 public static double noise(double xin, double yin, double zin) { 198 double n0, n1, n2, n3; // Noise contributions from the four corners 199 // Skew the input space to determine which simplex cell we're in 200 double s = (xin + yin + zin) * F3; // Very nice and simple skew 201 // factor for 3D 202 int i = (int) Math.floor(xin + s); 203 int j = (int) Math.floor(yin + s); 204 int k = (int) Math.floor(zin + s); 205 double t = (i + j + k) * G3; 206 double X0 = i - t; // Unskew the cell origin back to (x,y,z) space 207 double Y0 = j - t; 208 double Z0 = k - t; 209 double x0 = xin - X0; // The x,y,z distances from the cell origin 210 double y0 = yin - Y0; 211 double z0 = zin - Z0; 212 // For the 3D case, the simplex shape is a slightly irregular 213 // tetrahedron. 214 // Determine which simplex we are in. 215 int i1, j1, k1; // Offsets for second corner of simplex in (i,j,k) 216 // coords 217 int i2, j2, k2; // Offsets for third corner of simplex in (i,j,k) 218 // coords 219 if (x0 >= y0) { 220 if (y0 >= z0) { 221 i1 = 1; 222 j1 = 0; 223 k1 = 0; 224 i2 = 1; 225 j2 = 1; 226 k2 = 0; 227 } // X Y Z order 228 else if (x0 >= z0) { 229 i1 = 1; 230 j1 = 0; 231 k1 = 0; 232 i2 = 1; 233 j2 = 0; 234 k2 = 1; 235 } // X Z Y order 236 else { 237 i1 = 0; 238 j1 = 0; 239 k1 = 1; 240 i2 = 1; 241 j2 = 0; 242 k2 = 1; 243 } // Z X Y order 244 } else { // x0<y0 245 if (y0 < z0) { 246 i1 = 0; 247 j1 = 0; 248 k1 = 1; 249 i2 = 0; 250 j2 = 1; 251 k2 = 1; 252 } // Z Y X order 253 else if (x0 < z0) { 254 i1 = 0; 255 j1 = 1; 256 k1 = 0; 257 i2 = 0; 258 j2 = 1; 259 k2 = 1; 260 } // Y Z X order 261 else { 262 i1 = 0; 263 j1 = 1; 264 k1 = 0; 265 i2 = 1; 266 j2 = 1; 267 k2 = 0; 268 } // Y X Z order 269 } 270 // A step of (1,0,0) in (i,j,k) means a step of (1-c,-c,-c) in 271 // (x,y,z), 272 // a step of (0,1,0) in (i,j,k) means a step of (-c,1-c,-c) in 273 // (x,y,z), and 274 // a step of (0,0,1) in (i,j,k) means a step of (-c,-c,1-c) in 275 // (x,y,z), where 276 // c = 1/6. 277 double x1 = x0 - i1 + G3; // Offsets for second corner in (x,y,z) 278 // coords 279 double y1 = y0 - j1 + G3; 280 double z1 = z0 - k1 + G3; 281 double x2 = x0 - i2 + 2.0 * G3; // Offsets for third corner in 282 // (x,y,z) coords 283 double y2 = y0 - j2 + 2.0 * G3; 284 double z2 = z0 - k2 + 2.0 * G3; 285 double x3 = x0 - 1.0 + 3.0 * G3; // Offsets for last corner in 286 // (x,y,z) coords 287 double y3 = y0 - 1.0 + 3.0 * G3; 288 double z3 = z0 - 1.0 + 3.0 * G3; 289 // Work out the hashed gradient indices of the four simplex corners 290 int ii = i & 255; 291 int jj = j & 255; 292 int kk = k & 255; 293 int gi0 = perm[ii + perm[jj + perm[kk]]] % 12; 294 int gi1 = perm[ii + i1 + perm[jj + j1 + perm[kk + k1]]] % 12; 295 int gi2 = perm[ii + i2 + perm[jj + j2 + perm[kk + k2]]] % 12; 296 int gi3 = perm[ii + 1 + perm[jj + 1 + perm[kk + 1]]] % 12; 297 // Calculate the contribution from the four corners 298 double t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0; 299 if (t0 < 0) { 300 n0 = 0.0; 301 } else { 302 t0 *= t0; 303 n0 = t0 * t0 * dot(grad3[gi0], x0, y0, z0); 304 } 305 double t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1; 306 if (t1 < 0) { 307 n1 = 0.0; 308 } else { 309 t1 *= t1; 310 n1 = t1 * t1 * dot(grad3[gi1], x1, y1, z1); 311 } 312 double t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2; 313 if (t2 < 0) { 314 n2 = 0.0; 315 } else { 316 t2 *= t2; 317 n2 = t2 * t2 * dot(grad3[gi2], x2, y2, z2); 318 } 319 double t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3; 320 if (t3 < 0) { 321 n3 = 0.0; 322 } else { 323 t3 *= t3; 324 n3 = t3 * t3 * dot(grad3[gi3], x3, y3, z3); 325 } 326 // Add contributions from each corner to get the final noise value. 327 // The result is scaled to stay just inside [-1,1] 328 return 32.0 * (n0 + n1 + n2 + n3); 329 } 330 331 /** 332 * 4D simplex noise. 333 * 334 * @param x X position 335 * @param y Y position 336 * @param z Z position 337 * @param w Fourth-dimensional position. It is I, the Fourth-Dimensional Ziltoid the Omniscient! 338 * @return noise from -1.0 to 1.0, inclusive 339 */ 340 public static double noise(double x, double y, double z, double w) { 341 // The skewing and unskewing factors are hairy again for the 4D case 342 double n0, n1, n2, n3, n4; // Noise contributions from the five 343 // corners 344 // Skew the (x,y,z,w) space to determine which cell of 24 simplices 345 // we're in 346 double s = (x + y + z + w) * F4; // Factor for 4D skewing 347 int i = (int) Math.floor(x + s); 348 int j = (int) Math.floor(y + s); 349 int k = (int) Math.floor(z + s); 350 int l = (int) Math.floor(w + s); 351 double t = (i + j + k + l) * G4; // Factor for 4D unskewing 352 double X0 = i - t; // Unskew the cell origin back to (x,y,z,w) space 353 double Y0 = j - t; 354 double Z0 = k - t; 355 double W0 = l - t; 356 double x0 = x - X0; // The x,y,z,w distances from the cell origin 357 double y0 = y - Y0; 358 double z0 = z - Z0; 359 double w0 = w - W0; 360 // For the 4D case, the simplex is a 4D shape I won't even try to 361 // describe. 362 // To find out which of the 24 possible simplices we're in, we need 363 // to 364 // determine the magnitude ordering of x0, y0, z0 and w0. 365 // The method below is a good way of finding the ordering of x,y,z,w 366 // and 367 // then find the correct traversal order for the simplex we’re in. 368 // First, six pair-wise comparisons are performed between each 369 // possible pair 370 // of the four coordinates, and the results are used to add up binary 371 // bits 372 // for an integer index. 373 int c1 = (x0 > y0) ? 32 : 0; 374 int c2 = (x0 > z0) ? 16 : 0; 375 int c3 = (y0 > z0) ? 8 : 0; 376 int c4 = (x0 > w0) ? 4 : 0; 377 int c5 = (y0 > w0) ? 2 : 0; 378 int c6 = (z0 > w0) ? 1 : 0; 379 int c = c1 + c2 + c3 + c4 + c5 + c6; 380 int i1, j1, k1, l1; // The integer offsets for the second simplex 381 // corner 382 int i2, j2, k2, l2; // The integer offsets for the third simplex 383 // corner 384 int i3, j3, k3, l3; // The integer offsets for the fourth simplex 385 // corner 386 // simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some 387 // order. 388 // Many values of c will never occur, since e.g. x>y>z>w makes x<z, 389 // y<w and x<w 390 // impossible. Only the 24 indices which have non-zero entries make 391 // any sense. 392 // We use a thresholding to set the coordinates in turn from the 393 // largest magnitude. 394 // The number 3 in the "simplex" array is at the position of the 395 // largest coordinate. 396 i1 = simplex[c][0] >= 3 ? 1 : 0; 397 j1 = simplex[c][1] >= 3 ? 1 : 0; 398 k1 = simplex[c][2] >= 3 ? 1 : 0; 399 l1 = simplex[c][3] >= 3 ? 1 : 0; 400 // The number 2 in the "simplex" array is at the second largest 401 // coordinate. 402 i2 = simplex[c][0] >= 2 ? 1 : 0; 403 j2 = simplex[c][1] >= 2 ? 1 : 0; 404 k2 = simplex[c][2] >= 2 ? 1 : 0; 405 l2 = simplex[c][3] >= 2 ? 1 : 0; 406 // The number 1 in the "simplex" array is at the second smallest 407 // coordinate. 408 i3 = simplex[c][0] >= 1 ? 1 : 0; 409 j3 = simplex[c][1] >= 1 ? 1 : 0; 410 k3 = simplex[c][2] >= 1 ? 1 : 0; 411 l3 = simplex[c][3] >= 1 ? 1 : 0; 412 // The fifth corner has all coordinate offsets = 1, so no need to 413 // look that up. 414 double x1 = x0 - i1 + G4; // Offsets for second corner in (x,y,z,w) 415 // coords 416 double y1 = y0 - j1 + G4; 417 double z1 = z0 - k1 + G4; 418 double w1 = w0 - l1 + G4; 419 double x2 = x0 - i2 + 2.0 * G4; // Offsets for third corner in 420 // (x,y,z,w) coords 421 double y2 = y0 - j2 + 2.0 * G4; 422 double z2 = z0 - k2 + 2.0 * G4; 423 double w2 = w0 - l2 + 2.0 * G4; 424 double x3 = x0 - i3 + 3.0 * G4; // Offsets for fourth corner in 425 // (x,y,z,w) coords 426 double y3 = y0 - j3 + 3.0 * G4; 427 double z3 = z0 - k3 + 3.0 * G4; 428 double w3 = w0 - l3 + 3.0 * G4; 429 double x4 = x0 - 1.0 + 4.0 * G4; // Offsets for last corner in 430 // (x,y,z,w) coords 431 double y4 = y0 - 1.0 + 4.0 * G4; 432 double z4 = z0 - 1.0 + 4.0 * G4; 433 double w4 = w0 - 1.0 + 4.0 * G4; 434 // Work out the hashed gradient indices of the five simplex corners 435 int ii = i & 255; 436 int jj = j & 255; 437 int kk = k & 255; 438 int ll = l & 255; 439 int gi0 = perm[ii + perm[jj + perm[kk + perm[ll]]]] % 32; 440 int gi1 441 = perm[ii + i1 + perm[jj + j1 + perm[kk + k1 + perm[ll + l1]]]] % 32; 442 int gi2 443 = perm[ii + i2 + perm[jj + j2 + perm[kk + k2 + perm[ll + l2]]]] % 32; 444 int gi3 445 = perm[ii + i3 + perm[jj + j3 + perm[kk + k3 + perm[ll + l3]]]] % 32; 446 int gi4 447 = perm[ii + 1 + perm[jj + 1 + perm[kk + 1 + perm[ll + 1]]]] % 32; 448 // Calculate the contribution from the five corners 449 double t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0 - w0 * w0; 450 if (t0 < 0) { 451 n0 = 0.0; 452 } else { 453 t0 *= t0; 454 n0 = t0 * t0 * dot(grad4[gi0], x0, y0, z0, w0); 455 } 456 double t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1 - w1 * w1; 457 if (t1 < 0) { 458 n1 = 0.0; 459 } else { 460 t1 *= t1; 461 n1 = t1 * t1 * dot(grad4[gi1], x1, y1, z1, w1); 462 } 463 double t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2 - w2 * w2; 464 if (t2 < 0) { 465 n2 = 0.0; 466 } else { 467 t2 *= t2; 468 n2 = t2 * t2 * dot(grad4[gi2], x2, y2, z2, w2); 469 } 470 double t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3 - w3 * w3; 471 if (t3 < 0) { 472 n3 = 0.0; 473 } else { 474 t3 *= t3; 475 n3 = t3 * t3 * dot(grad4[gi3], x3, y3, z3, w3); 476 } 477 double t4 = 0.6 - x4 * x4 - y4 * y4 - z4 * z4 - w4 * w4; 478 if (t4 < 0) { 479 n4 = 0.0; 480 } else { 481 t4 *= t4; 482 n4 = t4 * t4 * dot(grad4[gi4], x4, y4, z4, w4); 483 } 484 // Sum up and scale the result to cover the range [-1,1] 485 return 27.0 * (n0 + n1 + n2 + n3 + n4); 486 } 487}