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 }