1 module gfm.math.simplexnoise; 2 3 import std.math, 4 std.random; 5 6 /// Simplex noise in 2D, 3D and 4D. 7 /// Translated from "Simplex noise demystified", Stefan Gustavson 8 class SimplexNoise(UniformRNG) if (isUniformRNG!UniformRNG) 9 { 10 public 11 { 12 /// Create from a RNG. 13 this(ref UniformRNG rng) 14 { 15 for (int i = 0; i < 256; ++i) 16 { 17 perm_table[i] = i; 18 } 19 20 for (int i = 0; i < 256; ++i) 21 { 22 int which = uniform(0, 256, rng); 23 int temp = perm_table[i]; 24 perm_table[i] = perm_table[which]; 25 perm_table[which] = temp; 26 } 27 } 28 29 /// Returns: 2D simplex noise. 30 double noise(double xin, double yin) 31 { 32 double n0, n1, n2; 33 // Noise contributions from the three corners 34 // Skew the input space to determine which simplex cell we're in 35 immutable double F2 = 0.5*(sqrt(3.0)-1.0); 36 double s = (xin+yin)*F2; 37 // Hairy factor for 2D 38 int i = fastfloor(xin+s); 39 int j = fastfloor(yin+s); 40 immutable double G2 = (3.0-sqrt(3.0))/6.0; 41 double t = (i+j)*G2; 42 double X0 = i-t; 43 // Unskew the cell origin back to (x,y) space 44 double Y0 = j-t; 45 double x0 = xin-X0; 46 // The x,y distances from the cell origin 47 double y0 = yin-Y0; 48 // For the 2D case, the simplex shape is an equilateral triangle. 49 // Determine which simplex we are in. 50 int i1, j1; 51 // Offsets for second (middle) corner of simplex in (i,j) coords 52 if(x0>y0) {i1=1; j1=0;} 53 // lower triangle, XY order: (0,0)->(1,0)->(1,1) 54 else {i1=0; j1=1;} 55 // upper triangle, YX order: (0,0)->(0,1)->(1,1) 56 // A step of (1,0) in (i,j) means a step of (1-c,-c) in (x,y), and 57 // a step of (0,1) in (i,j) means a step of (-c,1-c) in (x,y), where 58 // c = (3-sqrt(3))/6 59 double x1 = x0 - i1 + G2; 60 // Offsets for middle corner in (x,y) unskewed coords 61 double y1 = y0 - j1 + G2; 62 double x2 = x0 - 1.0 + 2.0 * G2; 63 // Offsets for last corner in (x,y) unskewed coords 64 double y2 = y0 - 1.0 + 2.0 * G2; 65 // Work out the hashed gradient indices of the three simplex corners 66 int ii = i & 255; 67 int jj = j & 255; 68 int gi0 = perm(ii+perm(jj)) % 12; 69 int gi1 = perm(ii+i1+perm(jj+j1)) % 12; 70 int gi2 = perm(ii+1+perm(jj+1)) % 12; 71 // Calculate the contribution from the three corners 72 double t0 = 0.5 - x0*x0-y0*y0; 73 if(t0<0) n0 = 0.0; 74 else { 75 t0 *= t0; 76 n0 = t0 * t0 * dot(grad3[gi0], x0, y0); 77 // (x,y) of grad3 used for 2D gradient 78 } 79 double t1 = 0.5 - x1*x1-y1*y1; 80 if(t1<0) n1 = 0.0; 81 else { 82 t1 *= t1; 83 n1 = t1 * t1 * dot(grad3[gi1], x1, y1); 84 } 85 double t2 = 0.5 - x2*x2-y2*y2; 86 if(t2<0) n2 = 0.0; 87 else { 88 t2 *= t2; 89 n2 = t2 * t2 * dot(grad3[gi2], x2, y2); 90 } 91 // Add contributions from each corner to get the final noise value. 92 // The result is scaled to return values in the interval [-1,1]. 93 return 70.0 * (n0 + n1 + n2); 94 } 95 96 97 /// Returns: 3D simplex noise. 98 double noise(double xin, double yin, double zin) 99 { 100 double n0, n1, n2, n3; 101 // Noise contributions from the four corners 102 // Skew the input space to determine which simplex cell we're in 103 immutable double F3 = 1.0/3.0; 104 double s = (xin+yin+zin)*F3; 105 // Very nice and simple skew factor for 3D 106 int i = fastfloor(xin+s); 107 int j = fastfloor(yin+s); 108 int k = fastfloor(zin+s); 109 immutable double G3 = 1.0/6.0; 110 // Very nice and simple unskew factor, too 111 double t = (i+j+k)*G3; 112 double X0 = i-t; 113 // Unskew the cell origin back to (x,y,z) space 114 double Y0 = j-t; 115 double Z0 = k-t; 116 double x0 = xin-X0; 117 // The x,y,z distances from the cell origin 118 double y0 = yin-Y0; 119 double z0 = zin-Z0; 120 // For the 3D case, the simplex shape is a slightly irregular tetrahedron. 121 // Determine which simplex we are in. 122 int i1, j1, k1; 123 // Offsets for second corner of simplex in (i,j,k) coords 124 int i2, j2, k2; 125 // Offsets for third corner of simplex in (i,j,k) coords 126 if(x0>=y0) { 127 if(y0>=z0) 128 { i1=1; j1=0; k1=0; i2=1; j2=1; k2=0; } 129 // X Y Z order 130 else if(x0>=z0) { i1=1; j1=0; k1=0; i2=1; j2=0; k2=1; } 131 // X Z Y order 132 else { i1=0; j1=0; k1=1; i2=1; j2=0; k2=1; } 133 // Z X Y order 134 } 135 else { 136 // x0<y0 137 if(y0<z0) { i1=0; j1=0; k1=1; i2=0; j2=1; k2=1; } 138 // Z Y X order 139 else if(x0<z0) { i1=0; j1=1; k1=0; i2=0; j2=1; k2=1; } 140 // Y Z X order 141 else { i1=0; j1=1; k1=0; i2=1; j2=1; k2=0; } 142 // Y X Z order 143 } 144 // A step of (1,0,0) in (i,j,k) means a step of (1-c,-c,-c) in (x,y,z), 145 // a step of (0,1,0) in (i,j,k) means a step of (-c,1-c,-c) in (x,y,z), and 146 // a step of (0,0,1) in (i,j,k) means a step of (-c,-c,1-c) in (x,y,z), where 147 // c = 1/6. 148 double x1 = x0 - i1 + G3; 149 // Offsets for second corner in (x,y,z) coords 150 double y1 = y0 - j1 + G3; 151 double z1 = z0 - k1 + G3; 152 double x2 = x0 - i2 + 2.0*G3; 153 // Offsets for third corner in (x,y,z) coords 154 double y2 = y0 - j2 + 2.0*G3; 155 double z2 = z0 - k2 + 2.0*G3; 156 double x3 = x0 - 1.0 + 3.0*G3; 157 // Offsets for last corner in (x,y,z) coords 158 double y3 = y0 - 1.0 + 3.0*G3; 159 double z3 = z0 - 1.0 + 3.0*G3; 160 // Work out the hashed gradient indices of the four simplex corners 161 int ii = i & 255; 162 int jj = j & 255; 163 int kk = k & 255; 164 int gi0 = perm(ii+perm(jj+perm(kk))) % 12; 165 int gi1 = perm(ii+i1+perm(jj+j1+perm(kk+k1))) % 12; 166 int gi2 = perm(ii+i2+perm(jj+j2+perm(kk+k2))) % 12; 167 int gi3 = perm(ii+1+perm(jj+1+perm(kk+1))) % 12; 168 // Calculate the contribution from the four corners 169 double t0 = 0.6 - x0*x0 - y0*y0 - z0*z0; 170 if(t0<0) n0 = 0.0; 171 else { 172 t0 *= t0; 173 n0 = t0 * t0 * dot(grad3[gi0], x0, y0, z0); 174 } 175 double t1 = 0.6 - x1*x1 - y1*y1 - z1*z1; 176 if(t1<0) n1 = 0.0; 177 else { 178 t1 *= t1; 179 n1 = t1 * t1 * dot(grad3[gi1], x1, y1, z1); 180 } 181 double t2 = 0.6 - x2*x2 - y2*y2 - z2*z2; 182 if(t2<0) n2 = 0.0; 183 else { 184 t2 *= t2; 185 n2 = t2 * t2 * dot(grad3[gi2], x2, y2, z2); 186 } 187 double t3 = 0.6 - x3*x3 - y3*y3 - z3*z3; 188 if(t3<0) n3 = 0.0; 189 else { 190 t3 *= t3; 191 n3 = t3 * t3 * dot(grad3[gi3], x3, y3, z3); 192 } 193 // Add contributions from each corner to get the final noise value. 194 // The result is scaled to stay just inside [-1,1] 195 return 32.0*(n0 + n1 + n2 + n3); 196 } 197 198 199 /// Returns: 4D simplex noise. 200 double noise(double x, double y, double z, double w) 201 { 202 // The skewing and unskewing factors are hairy again for the 4D case 203 immutable double F4 = (sqrt(5.0)-1.0)/4.0; 204 immutable double G4 = (5.0-sqrt(5.0))/20.0; 205 double n0, n1, n2, n3, n4; 206 // Noise contributions from the five corners 207 // Skew the (x,y,z,w) space to determine which cell of 24 simplices we're in 208 double s = (x + y + z + w) * F4; 209 // Factor for 4D skewing 210 int i = fastfloor(x + s); 211 int j = fastfloor(y + s); 212 int k = fastfloor(z + s); 213 int l = fastfloor(w + s); 214 double t = (i + j + k + l) * G4; 215 // Factor for 4D unskewing 216 double X0 = i - t; 217 // Unskew the cell origin back to (x,y,z,w) space 218 double Y0 = j - t; 219 double Z0 = k - t; 220 double W0 = l - t; 221 double x0 = x - X0; 222 // The x,y,z,w distances from the cell origin 223 double y0 = y - Y0; 224 double z0 = z - Z0; 225 double w0 = w - W0; 226 // For the 4D case, the simplex is a 4D shape I won't even try to describe. 227 // To find out which of the 24 possible simplices we're in, we need to 228 // determine the magnitude ordering of x0, y0, z0 and w0. 229 // The method below is a good way of finding the ordering of x,y,z,w and 230 // then find the correct traversal order for the simplex we’re in. 231 // First, six pair-wise comparisons are performed between each possible pair 232 // of the four coordinates, and the results are used to add up binary bits 233 // for an integer index. 234 int c1 = (x0 > y0) ? 32 : 0; 235 int c2 = (x0 > z0) ? 16 : 0; 236 int c3 = (y0 > z0) ? 8 : 0; 237 int c4 = (x0 > w0) ? 4 : 0; 238 int c5 = (y0 > w0) ? 2 : 0; 239 int c6 = (z0 > w0) ? 1 : 0; 240 int c = c1 + c2 + c3 + c4 + c5 + c6; 241 int i1, j1, k1, l1; 242 // The integer offsets for the second simplex corner 243 int i2, j2, k2, l2; 244 // The integer offsets for the third simplex corner 245 int i3, j3, k3, l3; 246 // The integer offsets for the fourth simplex corner 247 // simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some order. 248 // Many values of c will never occur, since e.g. x>y>z>w makes x<z, y<w and x<w 249 // impossible. Only the 24 indices which have non-zero entries make any sense. 250 // We use a thresholding to set the coordinates in turn from the largest magnitude. 251 // The number 3 in the "simplex" array is at the position of the largest coordinate. 252 i1 = simplex[c][0]>=3 ? 1 : 0; 253 j1 = simplex[c][1]>=3 ? 1 : 0; 254 k1 = simplex[c][2]>=3 ? 1 : 0; 255 l1 = simplex[c][3]>=3 ? 1 : 0; 256 // The number 2 in the "simplex" array is at the second largest coordinate. 257 i2 = simplex[c][0]>=2 ? 1 : 0; 258 j2 = simplex[c][1]>=2 ? 1 : 0; 259 k2 = simplex[c][2]>=2 ? 1 : 0; 260 l2 = simplex[c][3]>=2 ? 1 : 0; 261 // The number 1 in the "simplex" array is at the second smallest coordinate. 262 i3 = simplex[c][0]>=1 ? 1 : 0; 263 j3 = simplex[c][1]>=1 ? 1 : 0; 264 k3 = simplex[c][2]>=1 ? 1 : 0; 265 l3 = simplex[c][3]>=1 ? 1 : 0; 266 // The fifth corner has all coordinate offsets = 1, so no need to look that up. 267 double x1 = x0 - i1 + G4; 268 // Offsets for second corner in (x,y,z,w) coords 269 double y1 = y0 - j1 + G4; 270 double z1 = z0 - k1 + G4; 271 double w1 = w0 - l1 + G4; 272 double x2 = x0 - i2 + 2.0*G4; 273 // Offsets for third corner in (x,y,z,w) coords 274 double y2 = y0 - j2 + 2.0*G4; 275 double z2 = z0 - k2 + 2.0*G4; 276 double w2 = w0 - l2 + 2.0*G4; 277 double x3 = x0 - i3 + 3.0*G4; 278 // Offsets for fourth corner in (x,y,z,w) coords 279 double y3 = y0 - j3 + 3.0*G4; 280 double z3 = z0 - k3 + 3.0*G4; 281 double w3 = w0 - l3 + 3.0*G4; 282 double x4 = x0 - 1.0 + 4.0*G4; 283 // Offsets for last corner in (x,y,z,w) coords 284 double y4 = y0 - 1.0 + 4.0*G4; 285 double z4 = z0 - 1.0 + 4.0*G4; 286 double w4 = w0 - 1.0 + 4.0*G4; 287 // Work out the hashed gradient indices of the five simplex corners 288 int ii = i & 255; 289 int jj = j & 255; 290 int kk = k & 255; 291 int ll = l & 255; 292 int gi0 = perm(ii+perm(jj+perm(kk+perm(ll)))) % 32; 293 int gi1 = perm(ii+i1+perm(jj+j1+perm(kk+k1+perm(ll+l1)))) % 32; 294 int gi2 = perm(ii+i2+perm(jj+j2+perm(kk+k2+perm(ll+l2)))) % 32; 295 int gi3 = perm(ii+i3+perm(jj+j3+perm(kk+k3+perm(ll+l3)))) % 32; 296 int gi4 = perm(ii+1+perm(jj+1+perm(kk+1+perm(ll+1)))) % 32; 297 // Calculate the contribution from the five corners 298 double t0 = 0.6 - x0*x0 - y0*y0 - z0*z0 - w0*w0; 299 if(t0<0) n0 = 0.0; 300 else { 301 t0 *= t0; 302 n0 = t0 * t0 * dot(grad4[gi0], x0, y0, z0, w0); 303 } 304 double t1 = 0.6 - x1*x1 - y1*y1 - z1*z1 - w1*w1; 305 if(t1<0) n1 = 0.0; 306 else { 307 t1 *= t1; 308 n1 = t1 * t1 * dot(grad4[gi1], x1, y1, z1, w1); 309 } 310 double t2 = 0.6 - x2*x2 - y2*y2 - z2*z2 - w2*w2; 311 if(t2<0) n2 = 0.0; 312 else { 313 t2 *= t2; 314 n2 = t2 * t2 * dot(grad4[gi2], x2, y2, z2, w2); 315 } 316 double t3 = 0.6 - x3*x3 - y3*y3 - z3*z3 - w3*w3; 317 if(t3<0) n3 = 0.0; 318 else { 319 t3 *= t3; 320 n3 = t3 * t3 * dot(grad4[gi3], x3, y3, z3, w3); 321 } 322 double t4 = 0.6 - x4*x4 - y4*y4 - z4*z4 - w4*w4; 323 if(t4<0) n4 = 0.0; 324 else { 325 t4 *= t4; 326 n4 = t4 * t4 * dot(grad4[gi4], x4, y4, z4, w4); 327 } 328 // Sum up and scale the result to cover the range [-1,1] 329 return 27.0 * (n0 + n1 + n2 + n3 + n4); 330 } 331 } 332 333 private 334 { 335 int perm(int n) 336 { 337 return perm_table[n & 255]; 338 } 339 340 // permutation table 341 int[256] perm_table; 342 343 static int fastfloor(double x) 344 { 345 int ix = cast(int)x; 346 return x > 0 ? ix : ix - 1; 347 } 348 349 static double dot(immutable int g[3], double x, double y) 350 { 351 return g[0]*x + g[1]*y; 352 } 353 354 static double dot(immutable int g[3], double x, double y, double z) 355 { 356 return g[0]*x + g[1]*y + g[2]*z; 357 } 358 359 static double dot(immutable int g[4], double x, double y, double z, double w) 360 { 361 return g[0]*x + g[1]*y + g[2]*z + g[3]*w; 362 } 363 364 static immutable int[3][12] grad3 = 365 [ 366 [ 1, 1, 0], [-1, 1, 0], [ 1,-1, 0], [-1,-1, 0], 367 [ 1, 0, 1], [-1, 0, 1], [ 1, 0,-1], [-1, 0,-1], 368 [ 0, 1, 1], [ 0,-1, 1], [ 0, 1,-1], [ 0,-1,-1] 369 ]; 370 371 372 static immutable int[4][32] grad4 = 373 [ 374 [ 0, 1, 1, 1], [ 0, 1, 1,-1], [ 0, 1,-1, 1], [ 0, 1,-1,-1], 375 [ 0,-1, 1, 1], [ 0,-1, 1,-1], [ 0,-1,-1, 1], [ 0,-1,-1,-1], 376 [ 1, 0, 1, 1], [ 1, 0, 1,-1], [ 1, 0,-1, 1], [ 1, 0,-1,-1], 377 [-1, 0, 1, 1], [-1, 0, 1,-1], [-1, 0,-1, 1], [-1, 0,-1,-1], 378 [ 1, 1, 0, 1], [ 1, 1, 0,-1], [ 1,-1, 0, 1], [ 1,-1, 0,-1], 379 [-1, 1, 0, 1], [-1, 1, 0,-1], [-1,-1, 0, 1], [-1,-1, 0,-1], 380 [ 1, 1, 1, 0], [ 1, 1,-1, 0], [ 1,-1, 1, 0], [ 1,-1,-1, 0], 381 [-1, 1, 1, 0], [-1, 1,-1, 0], [-1,-1, 1, 0], [-1,-1,-1, 0] 382 ]; 383 384 // A lookup table to traverse the simplex around a given point in 4D. 385 // Details can be found where this table is used, in the 4D noise method. 386 static immutable int[4][64] simplex = 387 [ 388 [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], 389 [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], 390 [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], 391 [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], 392 [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], 393 [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], 394 [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], 395 [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] 396 ]; 397 } 398 } 399 400 401