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