1 /* 2 * This is a C++ port of version Stefan Gustavson's public domain 3 * implementation of simplex noise (Version 2012-03-09), which can be 4 * found at <http://webstaff.itn.liu.se/~stegu/simplexnoise/>. 5 * 6 * (Simplex Noise is a new (2001) algorithm created by Ken Perlin to 7 * replace his classic "Perlin" noise algorithm.) 8 * 9 * It was ported by Brendan Hickey (brendan@bhickey.net) and released on 10 * 2012-09-16. 11 * 12 * It is made available under the Creative Commons CC0 license. 13 * 14 * A speed-improved simplex noise algorithm for 2D, 3D and 4D in C++. 15 * 16 * Based on example code by Stefan Gustavson (stegu@itn.liu.se). 17 * Optimisations by Peter Eastman (peastman@drizzle.stanford.edu). 18 * Better rank ordering method by Stefan Gustavson in 2012. 19 * This could be speeded up even further, but it's useful as it is. 20 * 21 * Clumsily ported to some horrendous C/C++ mix by 22 * Brendan Hickey (brendan@bhickey.net) 23 * 24 * Version 2012-09-16 25 * 26 * This code was placed in the public domain by its original author, 27 * Stefan Gustavson. You may use it as you see fit, but 28 * attribution is appreciated. 29 * 30 */ 31 32 #include "AppHdr.h" 33 34 #include "perlin.h" 35 36 #include <cmath> 37 #include <stdint.h> 38 39 namespace perlin 40 { 41 // Inner class to speed up gradient computations 42 // ([in Java,] array access is a lot slower than member access) 43 class Grad 44 { 45 public: 46 const double x, y, z, w; Grad(double _x,double _y,double _z)47 Grad(double _x, double _y, double _z) : x(_x), y(_y), z(_z), w(0) {} Grad(double _x,double _y,double _z,double _w)48 Grad(double _x, double _y, double _z, double _w) : x(_x), y(_y), z(_z), w(_w) {} 49 }; 50 51 static const Grad grad3[] = 52 { 53 Grad(1,1,0), Grad(-1,1,0), Grad(1,-1,0), Grad(-1,-1,0), 54 Grad(1,0,1), Grad(-1,0,1), Grad(1,0,-1), Grad(-1,0,-1), 55 Grad(0,1,1), Grad(0,-1,1), Grad(0,1,-1), Grad(0,-1,-1) 56 }; 57 58 static const Grad grad4[] = 59 { 60 Grad(0,1,1,1), Grad(0,1,1,-1), Grad(0,1,-1,1), Grad(0,1,-1,-1), 61 Grad(0,-1,1,1), Grad(0,-1,1,-1), Grad(0,-1,-1,1), Grad(0,-1,-1,-1), 62 Grad(1,0,1,1), Grad(1,0,1,-1), Grad(1,0,-1,1), Grad(1,0,-1,-1), 63 Grad(-1,0,1,1), Grad(-1,0,1,-1), Grad(-1,0,-1,1), Grad(-1,0,-1,-1), 64 Grad(1,1,0,1), Grad(1,1,0,-1), Grad(1,-1,0,1), Grad(1,-1,0,-1), 65 Grad(-1,1,0,1), Grad(-1,1,0,-1), Grad(-1,-1,0,1), Grad(-1,-1,0,-1), 66 Grad(1,1,1,0), Grad(1,1,-1,0), Grad(1,-1,1,0), Grad(1,-1,-1,0), 67 Grad(-1,1,1,0), Grad(-1,1,-1,0), Grad(-1,-1,1,0), Grad(-1,-1,-1,0) 68 }; 69 70 static const uint8_t perm[] = {151,160,137,91,90,15, 71 131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23, 72 190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33, 73 88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166, 74 77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244, 75 102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196, 76 135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123, 77 5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42, 78 223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9, 79 129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228, 80 251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107, 81 49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254, 82 138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180, 83 // wrap 84 151,160,137,91,90,15, 85 131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23, 86 190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33, 87 88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166, 88 77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244, 89 102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196, 90 135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123, 91 5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42, 92 223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9, 93 129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228, 94 251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107, 95 49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254, 96 138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180}; 97 permMod12(const uint32_t x)98 static IMMUTABLE uint8_t permMod12(const uint32_t x) 99 { 100 return perm[x] % 12; 101 } 102 103 // Skewing and unskewing factors for 2, 3, and 4 dimensions 104 static const double F2 = 0.5 * (sqrt(3.0) - 1.0); 105 static const double G2 = (3.0 - sqrt(3.0)) / 6.0; 106 static const double F3 = 1.0 / 3.0; 107 static const double G3 = 1.0 / 6.0; 108 static const double F4 = (sqrt(5.0) - 1.0) / 4.0; 109 static const double G4 = (5.0 - sqrt(5.0)) / 20.0; 110 111 // Use uint64_t so that noise() can work sensibly for 112 // coordinates from the full range of uint32_t; otherwise scaling, 113 // signedness, and skew will give us considerably less than that. fastfloor(const double x)114 static uint64_t fastfloor(const double x) 115 { 116 uint64_t xi = (uint64_t) x; 117 return x < xi ? xi-1 : xi; 118 } 119 dot(Grad g,double x,double y)120 static double dot(Grad g, double x, double y) 121 { 122 return g.x*x + g.y*y; 123 } dot(Grad g,double x,double y,double z)124 static double dot(Grad g, double x, double y, double z) 125 { 126 return g.x*x + g.y*y + g.z*z; 127 } dot(Grad g,double x,double y,double z,double w)128 static double dot(Grad g, double x, double y, double z, double w) 129 { 130 return g.x*x + g.y*y + g.z*z + g.w*w; 131 } 132 133 134 // 2D simplex noise noise(double xin,double yin)135 double noise(double xin, double yin) 136 { 137 double n0, n1, n2; // Noise contributions from the three corners 138 // Skew the input space to determine which simplex cell we're in 139 double s = (xin+yin)*F2; // Hairy factor for 2D 140 uint64_t i = fastfloor(xin+s); 141 uint64_t j = fastfloor(yin+s); 142 double t = (i+j)*G2; 143 double X0 = i-t; // Unskew the cell origin back to (x,y) space 144 double Y0 = j-t; 145 double x0 = xin-X0; // The x,y distances from the cell origin 146 double y0 = yin-Y0; 147 // For the 2D case, the simplex shape is an equilateral triangle. 148 // Determine which simplex we are in. 149 int i1, j1; // Offsets for second (middle) corner of simplex in (i,j) coords 150 if (x0 > y0) 151 i1=1, j1=0; // lower triangle, XY order: (0,0)->(1,0)->(1,1) 152 else 153 i1=0, j1=1; // upper triangle, YX order: (0,0)->(0,1)->(1,1) 154 // A step of (1,0) in (i,j) means a step of (1-c,-c) in (x,y), and 155 // a step of (0,1) in (i,j) means a step of (-c,1-c) in (x,y), where 156 // c = (3-sqrt(3))/6 157 double x1 = x0 - i1 + G2; // Offsets for middle corner in (x,y) unskewed coords 158 double y1 = y0 - j1 + G2; 159 double x2 = x0 - 1.0 + 2.0 * G2; // Offsets for last corner in (x,y) unskewed coords 160 double y2 = y0 - 1.0 + 2.0 * G2; 161 // Work out the hashed gradient indices of the three simplex corners 162 int ii = i & 255; 163 int jj = j & 255; 164 int gi0 = permMod12(ii+perm[jj]); 165 int gi1 = permMod12(ii+i1+perm[jj+j1]); 166 int gi2 = permMod12(ii+1+perm[jj+1]); 167 // Calculate the contribution from the three corners 168 double t0 = 0.5 - x0*x0-y0*y0; 169 if (t0 < 0) 170 n0 = 0.0; 171 else 172 { 173 t0 *= t0; 174 n0 = t0 * t0 * dot(grad3[gi0], x0, y0); // (x,y) of grad3 used for 2D gradient 175 } 176 double t1 = 0.5 - x1*x1-y1*y1; 177 if (t1 < 0) 178 n1 = 0.0; 179 else 180 { 181 t1 *= t1; 182 n1 = t1 * t1 * dot(grad3[gi1], x1, y1); 183 } 184 double t2 = 0.5 - x2*x2-y2*y2; 185 if (t2 < 0) 186 n2 = 0.0; 187 else 188 { 189 t2 *= t2; 190 n2 = t2 * t2 * dot(grad3[gi2], x2, y2); 191 } 192 // Add contributions from each corner to get the final noise value. 193 // The result is scaled to return values in the interval [-1,1]. 194 return 70.0 * (n0 + n1 + n2); 195 } 196 197 // 3D simplex noise noise(double xin,double yin,double zin)198 double noise(double xin, double yin, double zin) 199 { 200 double n0, n1, n2, n3; // Noise contributions from the four corners 201 // Skew the input space to determine which simplex cell we're in 202 double s = (xin+yin+zin)*F3; // Very nice and simple skew factor for 3D 203 uint64_t i = fastfloor(xin+s); 204 uint64_t j = fastfloor(yin+s); 205 uint64_t k = fastfloor(zin+s); 206 double t = (i+j+k)*G3; 207 double X0 = i-t; // Unskew the cell origin back to (x,y,z) space 208 double Y0 = j-t; 209 double Z0 = k-t; 210 double x0 = xin-X0; // The x,y,z distances from the cell origin 211 double y0 = yin-Y0; 212 double z0 = zin-Z0; 213 // For the 3D case, the simplex shape is a slightly irregular tetrahedron. 214 // Determine which simplex we are in. 215 int i1, j1, k1; // Offsets for second corner of simplex in (i,j,k) coords 216 int i2, j2, k2; // Offsets for third corner of simplex in (i,j,k) coords 217 if (x0 >= y0) 218 { 219 if (y0 >= z0) 220 i1=1, j1=0, k1=0, i2=1, j2=1, k2=0; // X Y Z order 221 else if (x0 >= z0) 222 i1=1, j1=0, k1=0, i2=1, j2=0, k2=1; // X Z Y order 223 else 224 i1=0, j1=0, k1=1, i2=1, j2=0, k2=1; // Z X Y order 225 } 226 else 227 { // x0 < y0 228 if (y0 < z0) 229 i1=0, j1=0, k1=1, i2=0, j2=1, k2=1; // Z Y X order 230 else if (x0 < z0) 231 i1=0, j1=1, k1=0, i2=0, j2=1, k2=1; // Y Z X order 232 else 233 i1=0, j1=1, k1=0, i2=1, j2=1, k2=0; // Y X Z order 234 } 235 // A step of (1,0,0) in (i,j,k) means a step of (1-c,-c,-c) in (x,y,z), 236 // a step of (0,1,0) in (i,j,k) means a step of (-c,1-c,-c) in (x,y,z), and 237 // a step of (0,0,1) in (i,j,k) means a step of (-c,-c,1-c) in (x,y,z), where 238 // c = 1/6. 239 double x1 = x0 - i1 + G3; // Offsets for second corner in (x,y,z) coords 240 double y1 = y0 - j1 + G3; 241 double z1 = z0 - k1 + G3; 242 double x2 = x0 - i2 + 2.0*G3; // Offsets for third corner in (x,y,z) coords 243 double y2 = y0 - j2 + 2.0*G3; 244 double z2 = z0 - k2 + 2.0*G3; 245 double x3 = x0 - 1.0 + 3.0*G3; // Offsets for last corner in (x,y,z) coords 246 double y3 = y0 - 1.0 + 3.0*G3; 247 double z3 = z0 - 1.0 + 3.0*G3; 248 // Work out the hashed gradient indices of the four simplex corners 249 int ii = i & 255; 250 int jj = j & 255; 251 int kk = k & 255; 252 int gi0 = permMod12(ii+perm[jj+perm[kk]]); 253 int gi1 = permMod12(ii+i1+perm[jj+j1+perm[kk+k1]]); 254 int gi2 = permMod12(ii+i2+perm[jj+j2+perm[kk+k2]]); 255 int gi3 = permMod12(ii+1+perm[jj+1+perm[kk+1]]); 256 // Calculate the contribution from the four corners 257 double t0 = 0.6 - x0*x0 - y0*y0 - z0*z0; 258 if (t0 < 0) 259 n0 = 0.0; 260 else 261 { 262 t0 *= t0; 263 n0 = t0 * t0 * dot(grad3[gi0], x0, y0, z0); 264 } 265 double t1 = 0.6 - x1*x1 - y1*y1 - z1*z1; 266 if (t1 < 0) 267 n1 = 0.0; 268 else 269 { 270 t1 *= t1; 271 n1 = t1 * t1 * dot(grad3[gi1], x1, y1, z1); 272 } 273 double t2 = 0.6 - x2*x2 - y2*y2 - z2*z2; 274 if (t2 < 0) 275 n2 = 0.0; 276 else 277 { 278 t2 *= t2; 279 n2 = t2 * t2 * dot(grad3[gi2], x2, y2, z2); 280 } 281 double t3 = 0.6 - x3*x3 - y3*y3 - z3*z3; 282 if (t3<0) 283 n3 = 0.0; 284 else 285 { 286 t3 *= t3; 287 n3 = t3 * t3 * dot(grad3[gi3], x3, y3, z3); 288 } 289 // Add contributions from each corner to get the final noise value. 290 // The result is scaled to stay just inside [-1,1] 291 return 32.0 * (n0 + n1 + n2 + n3); 292 } 293 294 295 // 4D simplex noise, better simplex rank ordering method 2012-03-09 noise(double x,double y,double z,double w)296 double noise(double x, double y, double z, double w) 297 { 298 double n0, n1, n2, n3, n4; // Noise contributions from the five corners 299 // Skew the (x,y,z,w) space to determine which cell of 24 simplices we're in 300 double s = (x + y + z + w) * F4; // Factor for 4D skewing 301 uint64_t i = fastfloor(x + s); 302 uint64_t j = fastfloor(y + s); 303 uint64_t k = fastfloor(z + s); 304 uint64_t l = fastfloor(w + s); 305 double t = (i + j + k + l) * G4; // Factor for 4D unskewing 306 double X0 = i - t; // Unskew the cell origin back to (x,y,z,w) space 307 double Y0 = j - t; 308 double Z0 = k - t; 309 double W0 = l - t; 310 double x0 = x - X0; // The x,y,z,w distances from the cell origin 311 double y0 = y - Y0; 312 double z0 = z - Z0; 313 double w0 = w - W0; 314 // For the 4D case, the simplex is a 4D shape I won't even try to describe. 315 // To find out which of the 24 possible simplices we're in, we need to 316 // determine the magnitude ordering of x0, y0, z0 and w0. 317 // Six pair-wise comparisons are performed between each possible pair 318 // of the four coordinates, and the results are used to rank the numbers. 319 int rankx = 0; 320 int ranky = 0; 321 int rankz = 0; 322 int rankw = 0; 323 ++(x0 > y0 ? rankx : ranky); 324 ++(x0 > z0 ? rankx : rankz); 325 ++(x0 > w0 ? rankx : rankw); 326 ++(y0 > z0 ? ranky : rankz); 327 ++(y0 > w0 ? ranky : rankw); 328 ++(z0 > w0 ? rankz : rankw); 329 int i1, j1, k1, l1; // The integer offsets for the second simplex corner 330 int i2, j2, k2, l2; // The integer offsets for the third simplex corner 331 int i3, j3, k3, l3; // The integer offsets for the fourth simplex corner 332 // simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some order. 333 // Many values of c will never occur, since e.g. x>y>z>w makes x<z, y<w and x<w 334 // impossible. Only the 24 indices which have non-zero entries make any sense. 335 // We use a thresholding to set the coordinates in turn from the largest magnitude. 336 // Rank 3 denotes the largest coordinate. 337 i1 = rankx >= 3 ? 1 : 0; 338 j1 = ranky >= 3 ? 1 : 0; 339 k1 = rankz >= 3 ? 1 : 0; 340 l1 = rankw >= 3 ? 1 : 0; 341 // Rank 2 denotes the second largest coordinate. 342 i2 = rankx >= 2 ? 1 : 0; 343 j2 = ranky >= 2 ? 1 : 0; 344 k2 = rankz >= 2 ? 1 : 0; 345 l2 = rankw >= 2 ? 1 : 0; 346 // Rank 1 denotes the second smallest coordinate. 347 i3 = rankx >= 1 ? 1 : 0; 348 j3 = ranky >= 1 ? 1 : 0; 349 k3 = rankz >= 1 ? 1 : 0; 350 l3 = rankw >= 1 ? 1 : 0; 351 // The fifth corner has all coordinate offsets = 1, so no need to compute that. 352 double x1 = x0 - i1 + G4; // Offsets for second corner in (x,y,z,w) coords 353 double y1 = y0 - j1 + G4; 354 double z1 = z0 - k1 + G4; 355 double w1 = w0 - l1 + G4; 356 double x2 = x0 - i2 + 2.0*G4; // Offsets for third corner in (x,y,z,w) coords 357 double y2 = y0 - j2 + 2.0*G4; 358 double z2 = z0 - k2 + 2.0*G4; 359 double w2 = w0 - l2 + 2.0*G4; 360 double x3 = x0 - i3 + 3.0*G4; // Offsets for fourth corner in (x,y,z,w) coords 361 double y3 = y0 - j3 + 3.0*G4; 362 double z3 = z0 - k3 + 3.0*G4; 363 double w3 = w0 - l3 + 3.0*G4; 364 double x4 = x0 - 1.0 + 4.0*G4; // Offsets for last corner in (x,y,z,w) coords 365 double y4 = y0 - 1.0 + 4.0*G4; 366 double z4 = z0 - 1.0 + 4.0*G4; 367 double w4 = w0 - 1.0 + 4.0*G4; 368 // Work out the hashed gradient indices of the five simplex corners 369 int ii = i & 255; 370 int jj = j & 255; 371 int kk = k & 255; 372 int ll = l & 255; 373 int gi0 = perm[ii+perm[jj+perm[kk+perm[ll]]]] % 32; 374 int gi1 = perm[ii+i1+perm[jj+j1+perm[kk+k1+perm[ll+l1]]]] % 32; 375 int gi2 = perm[ii+i2+perm[jj+j2+perm[kk+k2+perm[ll+l2]]]] % 32; 376 int gi3 = perm[ii+i3+perm[jj+j3+perm[kk+k3+perm[ll+l3]]]] % 32; 377 int gi4 = perm[ii+1+perm[jj+1+perm[kk+1+perm[ll+1]]]] % 32; 378 // Calculate the contribution from the five corners 379 double t0 = 0.6 - x0*x0 - y0*y0 - z0*z0 - w0*w0; 380 if (t0 < 0) 381 n0 = 0.0; 382 else 383 { 384 t0 *= t0; 385 n0 = t0 * t0 * dot(grad4[gi0], x0, y0, z0, w0); 386 } 387 double t1 = 0.6 - x1*x1 - y1*y1 - z1*z1 - w1*w1; 388 if (t1 < 0) 389 n1 = 0.0; 390 else 391 { 392 t1 *= t1; 393 n1 = t1 * t1 * dot(grad4[gi1], x1, y1, z1, w1); 394 } 395 double t2 = 0.6 - x2*x2 - y2*y2 - z2*z2 - w2*w2; 396 if (t2 < 0) 397 n2 = 0.0; 398 else 399 { 400 t2 *= t2; 401 n2 = t2 * t2 * dot(grad4[gi2], x2, y2, z2, w2); 402 } 403 double t3 = 0.6 - x3*x3 - y3*y3 - z3*z3 - w3*w3; 404 if (t3 < 0) 405 n3 = 0.0; 406 else 407 { 408 t3 *= t3; 409 n3 = t3 * t3 * dot(grad4[gi3], x3, y3, z3, w3); 410 } 411 double t4 = 0.6 - x4*x4 - y4*y4 - z4*z4 - w4*w4; 412 if (t4 < 0) 413 n4 = 0.0; 414 else 415 { 416 t4 *= t4; 417 n4 = t4 * t4 * dot(grad4[gi4], x4, y4, z4, w4); 418 } 419 // Sum up and scale the result to cover the range [-1,1] 420 return 27.0 * (n0 + n1 + n2 + n3 + n4); 421 } 422 423 // This is *not* in Stefan Gustavson's Java original 424 // FIXME: what does it do? fBM(double x,double y,double z,uint32_t octaves)425 double fBM(double x, double y, double z, uint32_t octaves) 426 { 427 if (octaves < 1) 428 return 0.0; 429 if (octaves == 1) 430 return noise(x, y, z); 431 432 uint32_t divisor = 1; 433 double norm = 0.0; 434 double value = 0; 435 double xi = x; 436 double yi = y; 437 double zi = z; 438 for (uint32_t octave = 0; octave < octaves; ++octave) 439 { 440 value += noise(xi / divisor, yi / divisor, zi / divisor) / divisor; 441 norm += 1 / divisor; 442 divisor *= 2; 443 double xt = yi * sin(1.41421356) + cos(1.41421356); 444 yi = yi * cos(1.41421356) + sin(1.41421356); 445 xi = xt; 446 zi += 1.7; 447 } 448 return value / norm; 449 } 450 } 451