1 // noisegen.cpp 2 // 3 // Copyright (C) 2003, 2004 Jason Bevins 4 // 5 // This library is free software; you can redistribute it and/or modify it 6 // under the terms of the GNU Lesser General Public License as published by 7 // the Free Software Foundation; either version 2.1 of the License, or (at 8 // your option) any later version. 9 // 10 // This library is distributed in the hope that it will be useful, but WITHOUT 11 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 12 // FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 13 // License (COPYING.txt) for more details. 14 // 15 // You should have received a copy of the GNU Lesser General Public License 16 // along with this library; if not, write to the Free Software Foundation, 17 // Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 18 // 19 // The developer's email is jlbezigvins@gmzigail.com (for great email, take 20 // off every 'zig'.) 21 // 22 module noise.noisegen; 23 24 import noise.interp; 25 import noise.vectortable; 26 import noise.basictypes; 27 import std.math; 28 29 // Specifies the version of the coherent-noise functions to use. 30 // - Set to 2 to use the current version. 31 // - Set to 1 to use the flawed version from the original version of libnoise. 32 // If your application requires coherent-noise values that were generated by 33 // an earlier version of libnoise, change this constant to the appropriate 34 // value and recompile libnoise. 35 const int NOISE_VERSION = 2; 36 37 // These constants control certain parameters that all coherent-noise 38 // functions require. 39 static if (NOISE_VERSION == 1) { 40 // Constants used by the original version of libnoise. 41 // Because X_NOISE_GEN is not relatively prime to the other values, and 42 // Z_NOISE_GEN is close to 256 (the number of random gradient vectors), 43 // patterns show up in high-frequency coherent noise. 44 const int X_NOISE_GEN = 1; 45 const int Y_NOISE_GEN = 31337; 46 const int Z_NOISE_GEN = 263; 47 const int SEED_NOISE_GEN = 1013; 48 const int SHIFT_NOISE_GEN = 13; 49 } else { 50 // Constants used by the current version of libnoise. 51 const int X_NOISE_GEN = 1619; 52 const int Y_NOISE_GEN = 31337; 53 const int Z_NOISE_GEN = 6971; 54 const int SEED_NOISE_GEN = 1013; 55 const int SHIFT_NOISE_GEN = 8; 56 } 57 58 /// Enumerates the noise quality. 59 enum NoiseQuality 60 { 61 62 /// Generates coherent noise quickly. When a coherent-noise function with 63 /// this quality setting is used to generate a bump-map image, there are 64 /// noticeable "creasing" artifacts in the resulting image. This is 65 /// because the derivative of that function is discontinuous at integer 66 /// boundaries. 67 QUALITY_FAST = 0, 68 69 /// Generates standard-quality coherent noise. When a coherent-noise 70 /// function with this quality setting is used to generate a bump-map 71 /// image, there are some minor "creasing" artifacts in the resulting 72 /// image. This is because the second derivative of that function is 73 /// discontinuous at integer boundaries. 74 QUALITY_STD = 1, 75 76 /// Generates the best-quality coherent noise. When a coherent-noise 77 /// function with this quality setting is used to generate a bump-map 78 /// image, there are no "creasing" artifacts in the resulting image. This 79 /// is because the first and second derivatives of that function are 80 /// continuous at integer boundaries. 81 QUALITY_BEST = 2 82 83 }; 84 85 86 /// Generates a gradient-coherent-noise value from the coordinates of a 87 /// three-dimensional input value. 88 /// 89 /// @param x The @a x coordinate of the input value. 90 /// @param y The @a y coordinate of the input value. 91 /// @param z The @a z coordinate of the input value. 92 /// @param seed The random number seed. 93 /// @param noiseQuality The quality of the coherent-noise. 94 /// 95 /// @returns The generated gradient-coherent-noise value. 96 /// 97 /// The return value ranges from -1.0 to +1.0. 98 /// 99 /// For an explanation of the difference between <i>gradient</i> noise and 100 /// <i>value</i> noise, see the comments for the GradientNoise3D() function. 101 double GradientCoherentNoise3D (double x, double y, double z, int seed = 0, 102 NoiseQuality noiseQuality = NoiseQuality.QUALITY_STD) 103 { 104 // Create a unit-length cube aligned along an integer boundary. This cube 105 // surrounds the input point. 106 int x0 = (x > 0.0? cast(int)x: cast(int)x - 1); 107 int x1 = x0 + 1; 108 int y0 = (y > 0.0? cast(int)y: cast(int)y - 1); 109 int y1 = y0 + 1; 110 int z0 = (z > 0.0? cast(int)z: cast(int)z - 1); 111 int z1 = z0 + 1; 112 113 // Map the difference between the coordinates of the input value and the 114 // coordinates of the cube's outer-lower-left vertex onto an S-curve. 115 double xs = 0, ys = 0, zs = 0; 116 final switch (noiseQuality) { 117 case NoiseQuality.QUALITY_FAST: 118 xs = (x - cast(double)x0); 119 ys = (y - cast(double)y0); 120 zs = (z - cast(double)z0); 121 break; 122 case NoiseQuality.QUALITY_STD: 123 xs = SCurve3 (x - cast(double)x0); 124 ys = SCurve3 (y - cast(double)y0); 125 zs = SCurve3 (z - cast(double)z0); 126 break; 127 case NoiseQuality.QUALITY_BEST: 128 xs = SCurve5 (x - cast(double)x0); 129 ys = SCurve5 (y - cast(double)y0); 130 zs = SCurve5 (z - cast(double)z0); 131 break; 132 } 133 134 // Now calculate the noise values at each vertex of the cube. To generate 135 // the coherent-noise value at the input point, interpolate these eight 136 // noise values using the S-curve value as the interpolant (trilinear 137 // interpolation.) 138 double n0, n1, ix0, ix1, iy0, iy1; 139 n0 = GradientNoise3D (x, y, z, x0, y0, z0, seed); 140 n1 = GradientNoise3D (x, y, z, x1, y0, z0, seed); 141 ix0 = LinearInterp (n0, n1, xs); 142 n0 = GradientNoise3D (x, y, z, x0, y1, z0, seed); 143 n1 = GradientNoise3D (x, y, z, x1, y1, z0, seed); 144 ix1 = LinearInterp (n0, n1, xs); 145 iy0 = LinearInterp (ix0, ix1, ys); 146 n0 = GradientNoise3D (x, y, z, x0, y0, z1, seed); 147 n1 = GradientNoise3D (x, y, z, x1, y0, z1, seed); 148 ix0 = LinearInterp (n0, n1, xs); 149 n0 = GradientNoise3D (x, y, z, x0, y1, z1, seed); 150 n1 = GradientNoise3D (x, y, z, x1, y1, z1, seed); 151 ix1 = LinearInterp (n0, n1, xs); 152 iy1 = LinearInterp (ix0, ix1, ys); 153 154 return LinearInterp (iy0, iy1, zs); 155 } 156 157 /// Generates a gradient-noise value from the coordinates of a 158 /// three-dimensional input value and the integer coordinates of a 159 /// nearby three-dimensional value. 160 /// 161 /// @param fx The floating-point @a x coordinate of the input value. 162 /// @param fy The floating-point @a y coordinate of the input value. 163 /// @param fz The floating-point @a z coordinate of the input value. 164 /// @param ix The integer @a x coordinate of a nearby value. 165 /// @param iy The integer @a y coordinate of a nearby value. 166 /// @param iz The integer @a z coordinate of a nearby value. 167 /// @param seed The random number seed. 168 /// 169 /// @returns The generated gradient-noise value. 170 /// 171 /// @pre The difference between @a fx and @a ix must be less than or equal 172 /// to one. 173 /// 174 /// @pre The difference between @a fy and @a iy must be less than or equal 175 /// to one. 176 /// 177 /// @pre The difference between @a fz and @a iz must be less than or equal 178 /// to one. 179 /// 180 /// A <i>gradient</i>-noise function generates better-quality noise than a 181 /// <i>value</i>-noise function. Most noise modules use gradient noise for 182 /// this reason, although it takes much longer to calculate. 183 /// 184 /// The return value ranges from -1.0 to +1.0. 185 /// 186 /// This function generates a gradient-noise value by performing the 187 /// following steps: 188 /// - It first calculates a random normalized vector based on the 189 /// nearby integer value passed to this function. 190 /// - It then calculates a new value by adding this vector to the 191 /// nearby integer value passed to this function. 192 /// - It then calculates the dot product of the above-generated value 193 /// and the floating-point input value passed to this function. 194 /// 195 /// A noise function differs from a random-number generator because it 196 /// always returns the same output value if the same input value is passed 197 /// to it. 198 double GradientNoise3D (double fx, double fy, double fz, int ix, int iy, 199 int iz, int seed = 0) 200 { 201 // Randomly generate a gradient vector given the integer coordinates of the 202 // input value. This implementation generates a random number and uses it 203 // as an index into a normalized-vector lookup table. 204 int vectorIndex = ( 205 X_NOISE_GEN * ix 206 + Y_NOISE_GEN * iy 207 + Z_NOISE_GEN * iz 208 + SEED_NOISE_GEN * seed) 209 & 0xffffffff; 210 vectorIndex ^= (vectorIndex >> SHIFT_NOISE_GEN); 211 vectorIndex &= 0xff; 212 213 double xvGradient = g_randomVectors[(vectorIndex << 2) ]; 214 double yvGradient = g_randomVectors[(vectorIndex << 2) + 1]; 215 double zvGradient = g_randomVectors[(vectorIndex << 2) + 2]; 216 217 // Set up us another vector equal to the distance between the two vectors 218 // passed to this function. 219 double xvPoint = (fx - cast(double)ix); 220 double yvPoint = (fy - cast(double)iy); 221 double zvPoint = (fz - cast(double)iz); 222 223 // Now compute the dot product of the gradient vector with the distance 224 // vector. The resulting value is gradient noise. Apply a scaling value 225 // so that this noise value ranges from -1.0 to 1.0. 226 return ((xvGradient * xvPoint) 227 + (yvGradient * yvPoint) 228 + (zvGradient * zvPoint)) * 2.12; 229 } 230 231 /// Generates an integer-noise value from the coordinates of a 232 /// three-dimensional input value. 233 /// 234 /// @param x The integer @a x coordinate of the input value. 235 /// @param y The integer @a y coordinate of the input value. 236 /// @param z The integer @a z coordinate of the input value. 237 /// @param seed A random number seed. 238 /// 239 /// @returns The generated integer-noise value. 240 /// 241 /// The return value ranges from 0 to 2147483647. 242 /// 243 /// A noise function differs from a random-number generator because it 244 /// always returns the same output value if the same input value is passed 245 /// to it. 246 int IntValueNoise3D (int x, int y, int z, int seed = 0) 247 { 248 // All constants are primes and must remain prime in order for this noise 249 // function to work correctly. 250 int n = ( 251 X_NOISE_GEN * x 252 + Y_NOISE_GEN * y 253 + Z_NOISE_GEN * z 254 + SEED_NOISE_GEN * seed) 255 & 0x7fffffff; 256 n = (n >> 13) ^ n; 257 return (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff; 258 } 259 260 /// Modifies a floating-point value so that it can be stored in a 261 /// noise::int32 variable. 262 /// 263 /// @param n A floating-point number. 264 /// 265 /// @returns The modified floating-point number. 266 /// 267 /// This function does not modify @a n. 268 /// 269 /// In libnoise, the noise-generating algorithms are all integer-based; 270 /// they use variables of type noise::int32. Before calling a noise 271 /// function, pass the @a x, @a y, and @a z coordinates to this function to 272 /// ensure that these coordinates can be cast to a noise::int32 value. 273 /// 274 /// Although you could do a straight cast from double to noise::int32, the 275 /// resulting value may differ between platforms. By using this function, 276 /// you ensure that the resulting value is identical between platforms. 277 double MakeInt32Range (double n) 278 { 279 if (n >= 1073741824.0) { 280 return (2.0 * fmod (n, 1073741824.0)) - 1073741824.0; 281 } else if (n <= -1073741824.0) { 282 return (2.0 * fmod (n, 1073741824.0)) + 1073741824.0; 283 } else { 284 return n; 285 } 286 } 287 288 /// Generates a value-coherent-noise value from the coordinates of a 289 /// three-dimensional input value. 290 /// 291 /// @param x The @a x coordinate of the input value. 292 /// @param y The @a y coordinate of the input value. 293 /// @param z The @a z coordinate of the input value. 294 /// @param seed The random number seed. 295 /// @param noiseQuality The quality of the coherent-noise. 296 /// 297 /// @returns The generated value-coherent-noise value. 298 /// 299 /// The return value ranges from -1.0 to +1.0. 300 /// 301 /// For an explanation of the difference between <i>gradient</i> noise and 302 /// <i>value</i> noise, see the comments for the GradientNoise3D() function. 303 double ValueCoherentNoise3D (double x, double y, double z, int seed = 0, 304 NoiseQuality noiseQuality = NoiseQuality.QUALITY_STD) 305 { 306 // Create a unit-length cube aligned along an integer boundary. This cube 307 // surrounds the input point. 308 int x0 = (x > 0.0? cast(int)x: cast(int)x - 1); 309 int x1 = x0 + 1; 310 int y0 = (y > 0.0? cast(int)y: cast(int)y - 1); 311 int y1 = y0 + 1; 312 int z0 = (z > 0.0? cast(int)z: cast(int)z - 1); 313 int z1 = z0 + 1; 314 315 // Map the difference between the coordinates of the input value and the 316 // coordinates of the cube's outer-lower-left vertex onto an S-curve. 317 double xs = 0, ys = 0, zs = 0; 318 final switch (noiseQuality) { 319 case NoiseQuality.QUALITY_FAST: 320 xs = (x - cast(double)x0); 321 ys = (y - cast(double)y0); 322 zs = (z - cast(double)z0); 323 break; 324 case NoiseQuality.QUALITY_STD: 325 xs = SCurve3 (x - cast(double)x0); 326 ys = SCurve3 (y - cast(double)y0); 327 zs = SCurve3 (z - cast(double)z0); 328 break; 329 case NoiseQuality.QUALITY_BEST: 330 xs = SCurve5 (x - cast(double)x0); 331 ys = SCurve5 (y - cast(double)y0); 332 zs = SCurve5 (z - cast(double)z0); 333 break; 334 } 335 336 // Now calculate the noise values at each vertex of the cube. To generate 337 // the coherent-noise value at the input point, interpolate these eight 338 // noise values using the S-curve value as the interpolant (trilinear 339 // interpolation.) 340 double n0, n1, ix0, ix1, iy0, iy1; 341 n0 = ValueNoise3D (x0, y0, z0, seed); 342 n1 = ValueNoise3D (x1, y0, z0, seed); 343 ix0 = LinearInterp (n0, n1, xs); 344 n0 = ValueNoise3D (x0, y1, z0, seed); 345 n1 = ValueNoise3D (x1, y1, z0, seed); 346 ix1 = LinearInterp (n0, n1, xs); 347 iy0 = LinearInterp (ix0, ix1, ys); 348 n0 = ValueNoise3D (x0, y0, z1, seed); 349 n1 = ValueNoise3D (x1, y0, z1, seed); 350 ix0 = LinearInterp (n0, n1, xs); 351 n0 = ValueNoise3D (x0, y1, z1, seed); 352 n1 = ValueNoise3D (x1, y1, z1, seed); 353 ix1 = LinearInterp (n0, n1, xs); 354 iy1 = LinearInterp (ix0, ix1, ys); 355 return LinearInterp (iy0, iy1, zs); 356 } 357 358 /// Generates a value-noise value from the coordinates of a 359 /// three-dimensional input value. 360 /// 361 /// @param x The @a x coordinate of the input value. 362 /// @param y The @a y coordinate of the input value. 363 /// @param z The @a z coordinate of the input value. 364 /// @param seed A random number seed. 365 /// 366 /// @returns The generated value-noise value. 367 /// 368 /// The return value ranges from -1.0 to +1.0. 369 /// 370 /// A noise function differs from a random-number generator because it 371 /// always returns the same output value if the same input value is passed 372 /// to it. 373 double ValueNoise3D (int x, int y, int z, int seed = 0) 374 { 375 return 1.0 - (cast(double)IntValueNoise3D (x, y, z, seed) / 1073741824.0); 376 } 377 378 379 // 2D simplex noise 380 double SimplexNoise(double xin, double yin) { //727.20, 131.32 381 double n0, n1, n2; // Noise contributions from the three corners 382 // Skew the input space to determine which simplex cell we're in 383 double s = (xin+yin)*F2; // Hairy factor for 2D //s = 314.196 384 int i = fastfloor(xin+s); //1041 385 int j = fastfloor(yin+s); //445 386 double t = (i+j)*G2; //314.028 387 double X0 = i-t; // Unskew the cell origin back to (x,y) space //726.972 388 double Y0 = j-t; //130.972 389 double x0 = xin-X0; // The x,y distances from the cell origin //0.228 390 double y0 = yin-Y0; // 0.348 391 // For the 2D case, the simplex shape is an equilateral triangle. 392 // Determine which simplex we are in. 393 int i1, j1; // Offsets for second (middle) corner of simplex in (i,j) coords 394 if(x0>y0) {i1=1; j1=0;} // lower triangle, XY order: (0,0)->(1,0)->(1,1) 395 else {i1=0; j1=1;} // upper triangle, YX order: (0,0)->(0,1)->(1,1) 396 // A step of (1,0) in (i,j) means a step of (1-c,-c) in (x,y), and 397 // a step of (0,1) in (i,j) means a step of (-c,1-c) in (x,y), where 398 // c = (3-sqrt(3))/6 399 double x1 = x0 - i1 + G2; // Offsets for middle corner in (x,y) unskewed coords 400 double y1 = y0 - j1 + G2; 401 double x2 = x0 - 1.0 + 2.0 * G2; // Offsets for last corner in (x,y) unskewed coords 402 double y2 = y0 - 1.0 + 2.0 * G2; 403 // Work out the hashed gradient indices of the three simplex corners 404 int ii = i & 255; 405 int jj = j & 255; 406 int gi0 = permMod12[ii+perm[jj]]; 407 int gi1 = permMod12[ii+i1+perm[jj+j1]]; 408 int gi2 = permMod12[ii+1+perm[jj+1]]; 409 // Calculate the contribution from the three corners 410 double t0 = 0.5 - (x0*x0) - (y0*y0); 411 if(t0<0) n0 = 0.0; 412 else { 413 t0 *= t0; 414 n0 = t0 * t0 * dot(grad3[gi0], x0, y0); // (x,y) of grad3 used for 2D gradient 415 } 416 double t1 = 0.5 - (x1*x1) - (y1*y1); 417 if(t1<0) n1 = 0.0; 418 else { 419 t1 *= t1; 420 n1 = t1 * t1 * dot(grad3[gi1], x1, y1); 421 } 422 double t2 = 0.5 - (x2*x2) - (y2*y2); 423 if(t2<0) n2 = 0.0; 424 else { 425 t2 *= t2; 426 n2 = t2 * t2 * dot(grad3[gi2], x2, y2); 427 } 428 // Add contributions from each corner to get the final noise value. 429 // The result is scaled to return values in the interval [-1,1]. 430 return 70.1043345 * (n0 + n1 + n2); 431 } 432 433 // 3D simplex noise 434 double SimplexNoise(double xin, double yin, double zin) { 435 double n0, n1, n2, n3; // Noise contributions from the four corners 436 // Skew the input space to determine which simplex cell we're in 437 double s = (xin+yin+zin)*F3; // Very nice and simple skew factor for 3D 438 int i = fastfloor(xin+s); 439 int j = fastfloor(yin+s); 440 int k = fastfloor(zin+s); 441 double t = (i+j+k)*G3; 442 double X0 = i-t; // Unskew the cell origin back to (x,y,z) space 443 double Y0 = j-t; 444 double Z0 = k-t; 445 double x0 = xin-X0; // The x,y,z distances from the cell origin 446 double y0 = yin-Y0; 447 double z0 = zin-Z0; 448 // For the 3D case, the simplex shape is a slightly irregular tetrahedron. 449 // Determine which simplex we are in. 450 int i1, j1, k1; // Offsets for second corner of simplex in (i,j,k) coords 451 int i2, j2, k2; // Offsets for third corner of simplex in (i,j,k) coords 452 if(x0>=y0) { 453 if(y0>=z0) { i1=1; j1=0; k1=0; i2=1; j2=1; k2=0; } // X Y Z order 454 else if(x0>=z0) { i1=1; j1=0; k1=0; i2=1; j2=0; k2=1; } // X Z Y order 455 else { i1=0; j1=0; k1=1; i2=1; j2=0; k2=1; } // Z X Y order 456 } 457 else { // x0<y0 458 if(y0<z0) { i1=0; j1=0; k1=1; i2=0; j2=1; k2=1; } // Z Y X order 459 else if(x0<z0) { i1=0; j1=1; k1=0; i2=0; j2=1; k2=1; } // Y Z X order 460 else { i1=0; j1=1; k1=0; i2=1; j2=1; k2=0; } // Y X Z order 461 } 462 // A step of (1,0,0) in (i,j,k) means a step of (1-c,-c,-c) in (x,y,z), 463 // a step of (0,1,0) in (i,j,k) means a step of (-c,1-c,-c) in (x,y,z), and 464 // a step of (0,0,1) in (i,j,k) means a step of (-c,-c,1-c) in (x,y,z), where 465 // c = 1/6. 466 double x1 = x0 - i1 + G3; // Offsets for second corner in (x,y,z) coords 467 double y1 = y0 - j1 + G3; 468 double z1 = z0 - k1 + G3; 469 double x2 = x0 - i2 + 2.0*G3; // Offsets for third corner in (x,y,z) coords 470 double y2 = y0 - j2 + 2.0*G3; 471 double z2 = z0 - k2 + 2.0*G3; 472 double x3 = x0 - 1.0 + 3.0*G3; // Offsets for last corner in (x,y,z) coords 473 double y3 = y0 - 1.0 + 3.0*G3; 474 double z3 = z0 - 1.0 + 3.0*G3; 475 // Work out the hashed gradient indices of the four simplex corners 476 int ii = i & 255; 477 int jj = j & 255; 478 int kk = k & 255; 479 int gi0 = permMod12[ii+perm[jj+perm[kk]]]; 480 int gi1 = permMod12[ii+i1+perm[jj+j1+perm[kk+k1]]]; 481 int gi2 = permMod12[ii+i2+perm[jj+j2+perm[kk+k2]]]; 482 int gi3 = permMod12[ii+1+perm[jj+1+perm[kk+1]]]; 483 // Calculate the contribution from the four corners, the 0.5 is changed from most 484 // (all?) implementations i've seen (it was 0.6), as per recommendation in this pdf 485 // by Stefan Gustavsen, the writer of the implementation this is based off: 486 // http://staffwww.itn.liu.se/~stegu/simplexnoise/simplexnoise.pdf 487 double t0 = 0.5 - x0*x0 - y0*y0 - z0*z0; 488 if(t0<0) n0 = 0.0; 489 else { 490 t0 *= t0; 491 n0 = t0 * t0 * dot(grad3[gi0], x0, y0, z0); 492 } 493 double t1 = 0.5 - x1*x1 - y1*y1 - z1*z1; 494 if(t1<0) n1 = 0.0; 495 else { 496 t1 *= t1; 497 n1 = t1 * t1 * dot(grad3[gi1], x1, y1, z1); 498 } 499 double t2 = 0.5 - x2*x2 - y2*y2 - z2*z2; 500 if(t2<0) n2 = 0.0; 501 else { 502 t2 *= t2; 503 n2 = t2 * t2 * dot(grad3[gi2], x2, y2, z2); 504 } 505 double t3 = 0.5 - x3*x3 - y3*y3 - z3*z3; 506 if(t3<0) n3 = 0.0; 507 else { 508 t3 *= t3; 509 n3 = t3 * t3 * dot(grad3[gi3], x3, y3, z3); 510 } 511 // Add contributions from each corner to get the final noise value. 512 // The result is scaled to stay just inside [-1,1] 513 return 76.8*(n0 + n1 + n2 + n3); 514 } 515 516 // 4D simplex noise, better simplex rank ordering method 2012-03-09 517 double SimplexNoise(double x, double y, double z, double w) { 518 519 double n0, n1, n2, n3, n4; // Noise contributions from the five corners 520 // Skew the (x,y,z,w) space to determine which cell of 24 simplices we're in 521 double s = (x + y + z + w) * F4; // Factor for 4D skewing 522 int i = fastfloor(x + s); 523 int j = fastfloor(y + s); 524 int k = fastfloor(z + s); 525 int l = fastfloor(w + s); 526 double t = (i + j + k + l) * G4; // Factor for 4D unskewing 527 double X0 = i - t; // Unskew the cell origin back to (x,y,z,w) space 528 double Y0 = j - t; 529 double Z0 = k - t; 530 double W0 = l - t; 531 double x0 = x - X0; // The x,y,z,w distances from the cell origin 532 double y0 = y - Y0; 533 double z0 = z - Z0; 534 double w0 = w - W0; 535 // For the 4D case, the simplex is a 4D shape I won't even try to describe. 536 // To find out which of the 24 possible simplices we're in, we need to 537 // determine the magnitude ordering of x0, y0, z0 and w0. 538 // Six pair-wise comparisons are performed between each possible pair 539 // of the four coordinates, and the results are used to rank the numbers. 540 int rankx = 0; 541 int ranky = 0; 542 int rankz = 0; 543 int rankw = 0; 544 if(x0 > y0) rankx++; else ranky++; 545 if(x0 > z0) rankx++; else rankz++; 546 if(x0 > w0) rankx++; else rankw++; 547 if(y0 > z0) ranky++; else rankz++; 548 if(y0 > w0) ranky++; else rankw++; 549 if(z0 > w0) rankz++; else rankw++; 550 int i1, j1, k1, l1; // The integer offsets for the second simplex corner 551 int i2, j2, k2, l2; // The integer offsets for the third simplex corner 552 int i3, j3, k3, l3; // The integer offsets for the fourth simplex corner 553 // simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some order. 554 // Many values of c will never occur, since e.g. x>y>z>w makes x<z, y<w and x<w 555 // impossible. Only the 24 indices which have non-zero entries make any sense. 556 // We use a thresholding to set the coordinates in turn from the largest magnitude. 557 // Rank 3 denotes the largest coordinate. 558 i1 = rankx >= 3 ? 1 : 0; 559 j1 = ranky >= 3 ? 1 : 0; 560 k1 = rankz >= 3 ? 1 : 0; 561 l1 = rankw >= 3 ? 1 : 0; 562 // Rank 2 denotes the second largest coordinate. 563 i2 = rankx >= 2 ? 1 : 0; 564 j2 = ranky >= 2 ? 1 : 0; 565 k2 = rankz >= 2 ? 1 : 0; 566 l2 = rankw >= 2 ? 1 : 0; 567 // Rank 1 denotes the second smallest coordinate. 568 i3 = rankx >= 1 ? 1 : 0; 569 j3 = ranky >= 1 ? 1 : 0; 570 k3 = rankz >= 1 ? 1 : 0; 571 l3 = rankw >= 1 ? 1 : 0; 572 // The fifth corner has all coordinate offsets = 1, so no need to compute that. 573 double x1 = x0 - i1 + G4; // Offsets for second corner in (x,y,z,w) coords 574 double y1 = y0 - j1 + G4; 575 double z1 = z0 - k1 + G4; 576 double w1 = w0 - l1 + G4; 577 double x2 = x0 - i2 + 2.0*G4; // Offsets for third corner in (x,y,z,w) coords 578 double y2 = y0 - j2 + 2.0*G4; 579 double z2 = z0 - k2 + 2.0*G4; 580 double w2 = w0 - l2 + 2.0*G4; 581 double x3 = x0 - i3 + 3.0*G4; // Offsets for fourth corner in (x,y,z,w) coords 582 double y3 = y0 - j3 + 3.0*G4; 583 double z3 = z0 - k3 + 3.0*G4; 584 double w3 = w0 - l3 + 3.0*G4; 585 double x4 = x0 - 1.0 + 4.0*G4; // Offsets for last corner in (x,y,z,w) coords 586 double y4 = y0 - 1.0 + 4.0*G4; 587 double z4 = z0 - 1.0 + 4.0*G4; 588 double w4 = w0 - 1.0 + 4.0*G4; 589 // Work out the hashed gradient indices of the five simplex corners 590 int ii = i & 255; 591 int jj = j & 255; 592 int kk = k & 255; 593 int ll = l & 255; 594 int gi0 = permMod32[ii+perm[jj+perm[kk+perm[ll]]]]; 595 int gi1 = permMod32[ii+i1+perm[jj+j1+perm[kk+k1+perm[ll+l1]]]]; 596 int gi2 = permMod32[ii+i2+perm[jj+j2+perm[kk+k2+perm[ll+l2]]]]; 597 int gi3 = permMod32[ii+i3+perm[jj+j3+perm[kk+k3+perm[ll+l3]]]]; 598 int gi4 = permMod32[ii+1+perm[jj+1+perm[kk+1+perm[ll+1]]]]; 599 // Calculate the contribution from the five corners 600 double t0 = 0.5 - x0*x0 - y0*y0 - z0*z0 - w0*w0; 601 if(t0<0) n0 = 0.0; 602 else { 603 t0 *= t0; 604 n0 = t0 * t0 * dot(grad4[gi0], x0, y0, z0, w0); 605 } 606 double t1 = 0.5 - x1*x1 - y1*y1 - z1*z1 - w1*w1; 607 if(t1<0) n1 = 0.0; 608 else { 609 t1 *= t1; 610 n1 = t1 * t1 * dot(grad4[gi1], x1, y1, z1, w1); 611 } 612 double t2 = 0.5 - x2*x2 - y2*y2 - z2*z2 - w2*w2; 613 if(t2<0) n2 = 0.0; 614 else { 615 t2 *= t2; 616 n2 = t2 * t2 * dot(grad4[gi2], x2, y2, z2, w2); 617 } 618 double t3 = 0.5 - x3*x3 - y3*y3 - z3*z3 - w3*w3; 619 if(t3<0) n3 = 0.0; 620 else { 621 t3 *= t3; 622 n3 = t3 * t3 * dot(grad4[gi3], x3, y3, z3, w3); 623 } 624 double t4 = 0.5 - x4*x4 - y4*y4 - z4*z4 - w4*w4; 625 if(t4<0) n4 = 0.0; 626 else { 627 t4 *= t4; 628 n4 = t4 * t4 * dot(grad4[gi4], x4, y4, z4, w4); 629 } 630 // Sum up and scale the result to cover the range [-1,1] 631 return 62.7*(n0 + n1 + n2 + n3 + n4); 632 } 633 634 double dot(Grad g, double x, double y) { 635 return g.x*x + g.y*y; 636 } 637 638 double dot(Grad g, double x, double y, double z) { 639 return g.x*x + g.y*y + g.z*z; 640 } 641 642 double dot(Grad g, double x, double y, double z, double w) { 643 return g.x*x + g.y*y + g.z*z + g.w*w; 644 } 645 646 int fastfloor(double x) { 647 int xi = cast(int)x; 648 return x<xi ? xi-1 : xi; 649 } 650 651 short[512] modulo_lookup(immutable int modulo) { 652 short[512] temp; 653 for(int i = 0; i < perm.length; i++) { 654 temp[i] = perm[i] % modulo; 655 } 656 return temp; 657 } 658 659 //lookup table for selecting gradients, doubled to avoid having to deal 660 //with index wrapping. 661 immutable short[512] perm = [151,160,137,91,90,15, 662 131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23, 663 190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33, 664 88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166, 665 77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244, 666 102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196, 667 135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123, 668 5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42, 669 223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9, 670 129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228, 671 251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107, 672 49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254, 673 138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180, 674 151,160,137,91,90,15, 675 131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23, 676 190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33, 677 88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166, 678 77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244, 679 102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196, 680 135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123, 681 5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42, 682 223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9, 683 129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228, 684 251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107, 685 49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254, 686 138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180]; 687 688 immutable short[] permMod12 = modulo_lookup(12); 689 immutable short[] permMod32 = modulo_lookup(32); 690 691 //gradient tables 692 immutable Grad[12] grad3 = [ Grad(1,1,0),Grad(-1,1,0),Grad(1,-1,0),Grad(-1,-1,0), 693 Grad(1,0,1),Grad(-1,0,1),Grad(1,0,-1),Grad(-1,0,-1), 694 Grad(0,1,1),Grad(0,-1,1),Grad(0,1,-1),Grad(0,-1,-1)]; 695 696 immutable Grad[32] grad4 = [ Grad(0,1,1,1),Grad(0,1,1,-1),Grad(0,1,-1,1),Grad(0,1,-1,-1), 697 Grad(0,-1,1,1),Grad(0,-1,1,-1),Grad(0,-1,-1,1),Grad(0,-1,-1,-1), 698 Grad(1,0,1,1),Grad(1,0,1,-1),Grad(1,0,-1,1),Grad(1,0,-1,-1), 699 Grad(-1,0,1,1),Grad(-1,0,1,-1),Grad(-1,0,-1,1),Grad(-1,0,-1,-1), 700 Grad(1,1,0,1),Grad(1,1,0,-1),Grad(1,-1,0,1),Grad(1,-1,0,-1), 701 Grad(-1,1,0,1),Grad(-1,1,0,-1),Grad(-1,-1,0,1),Grad(-1,-1,0,-1), 702 Grad(1,1,1,0),Grad(1,1,-1,0),Grad(1,-1,1,0),Grad(1,-1,-1,0), 703 Grad(-1,1,1,0),Grad(-1,1,-1,0),Grad(-1,-1,1,0),Grad(-1,-1,-1,0)]; 704 705 // Skewing and unskewing factors for 2, 3, and 4 dimensions 706 const double F2 = 0.5*(sqrt(cast(real)3.0)-1.0); 707 const double G2 = (3.0-sqrt(cast(real)3.0))/6.0; 708 const double F3 = 1.0/3.0; 709 const double G3 = 1.0/6.0; 710 const double F4 = (sqrt(cast(real)5.0)-1.0)/4.0; 711 const double G4 = (5.0-sqrt(cast(real)5.0))/20.0; 712 713 private { 714 struct Grad { 715 int x = 0; 716 int y = 0; 717 int z = 0; 718 int w = 0; 719 this(int x, int y, int z) { 720 this.x = x; 721 this.y = y; 722 this.z = z; 723 } 724 725 this(int x, int y, int z, int w) { 726 this.x = x; 727 this.y = y; 728 this.z = z; 729 this.w = w; 730 } 731 } 732 }