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