src/mesa/program/prog_noise.c

   1 /*
   2  * Mesa 3-D graphics library
   3  * Version:  6.5
   4  *
   5  * Copyright (C) 2006  Brian Paul   All Rights Reserved.
   6  *
   7  * Permission is hereby granted, free of charge, to any person obtaining a
   8  * copy of this software and associated documentation files (the "Software"),
   9  * to deal in the Software without restriction, including without limitation
  10  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
  11  * and/or sell copies of the Software, and to permit persons to whom the
  12  * Software is furnished to do so, subject to the following conditions:
  13  *
  14  * The above copyright notice and this permission notice shall be included
  15  * in all copies or substantial portions of the Software.
  16  *
  17  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  18  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  19  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
  20  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
  21  * OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
  22  * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  23  * OTHER DEALINGS IN THE SOFTWARE.
  24  */
  25
  26 /*
  27  * SimplexNoise1234
  28  * Copyright (c) 2003-2005, Stefan Gustavson
  29  *
  30  * Contact: stegu@itn.liu.se
  31  */
  32
  33 /**
  34  * \file
  35  * \brief C implementation of Perlin Simplex Noise over 1, 2, 3 and 4 dims.
  36  * \author Stefan Gustavson (stegu@itn.liu.se)
  37  *
  38  *
  39  * This implementation is "Simplex Noise" as presented by
  40  * Ken Perlin at a relatively obscure and not often cited course
  41  * session "Real-Time Shading" at Siggraph 2001 (before real
  42  * time shading actually took on), under the title "hardware noise".
  43  * The 3D function is numerically equivalent to his Java reference
  44  * code available in the PDF course notes, although I re-implemented
  45  * it from scratch to get more readable code. The 1D, 2D and 4D cases
  46  * were implemented from scratch by me from Ken Perlin's text.
  47  *
  48  * This file has no dependencies on any other file, not even its own
  49  * header file. The header file is made for use by external code only.
  50  */
  51
  52
  53 #include "main/imports.h"
  54 #include "prog_noise.h"
  55
  56 #define FASTFLOOR(x) ( ((x)>0) ? ((int)x) : (((int)x)-1) )
  57
  58 /*
  59  * ---------------------------------------------------------------------
  60  * Static data
  61  */
  62
  63 /**
  64  * Permutation table. This is just a random jumble of all numbers 0-255,
  65  * repeated twice to avoid wrapping the index at 255 for each lookup.
  66  * This needs to be exactly the same for all instances on all platforms,
  67  * so it's easiest to just keep it as static explicit data.
  68  * This also removes the need for any initialisation of this class.
  69  *
  70  * Note that making this an int[] instead of a char[] might make the
  71  * code run faster on platforms with a high penalty for unaligned single
  72  * byte addressing. Intel x86 is generally single-byte-friendly, but
  73  * some other CPUs are faster with 4-aligned reads.
  74  * However, a char[] is smaller, which avoids cache trashing, and that
  75  * is probably the most important aspect on most architectures.
  76  * This array is accessed a *lot* by the noise functions.
  77  * A vector-valued noise over 3D accesses it 96 times, and a
  78  * float-valued 4D noise 64 times. We want this to fit in the cache!
  79  */
  80 static const unsigned char perm[512] = { 151, 160, 137, 91, 90, 15,
  81    131, 13, 201, 95, 96, 53, 194, 233, 7, 225, 140, 36, 103, 30, 69, 142, 8,
  82       99, 37, 240, 21, 10, 23,
  83    190, 6, 148, 247, 120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35,
  84       11, 32, 57, 177, 33,
  85    88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 68, 175, 74, 165, 71,
  86       134, 139, 48, 27, 166,
  87    77, 146, 158, 231, 83, 111, 229, 122, 60, 211, 133, 230, 220, 105, 92, 41,
  88       55, 46, 245, 40, 244,
  89    102, 143, 54, 65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 187, 208, 89,
  90       18, 169, 200, 196,
  91    135, 130, 116, 188, 159, 86, 164, 100, 109, 198, 173, 186, 3, 64, 52, 217,
  92       226, 250, 124, 123,
  93    5, 202, 38, 147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58,
  94       17, 182, 189, 28, 42,
  95    223, 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101, 155,
  96       167, 43, 172, 9,
  97    129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, 112, 104,
  98       218, 246, 97, 228,
  99    251, 34, 242, 193, 238, 210, 144, 12, 191, 179, 162, 241, 81, 51, 145, 235,
 100       249, 14, 239, 107,
 101    49, 192, 214, 31, 181, 199, 106, 157, 184, 84, 204, 176, 115, 121, 50, 45,
 102       127, 4, 150, 254,
 103    138, 236, 205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66,
 104       215, 61, 156, 180,
 105    151, 160, 137, 91, 90, 15,
 106    131, 13, 201, 95, 96, 53, 194, 233, 7, 225, 140, 36, 103, 30, 69, 142, 8,
 107       99, 37, 240, 21, 10, 23,
 108    190, 6, 148, 247, 120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35,
 109       11, 32, 57, 177, 33,
 110    88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 68, 175, 74, 165, 71,
 111       134, 139, 48, 27, 166,
 112    77, 146, 158, 231, 83, 111, 229, 122, 60, 211, 133, 230, 220, 105, 92, 41,
 113       55, 46, 245, 40, 244,
 114    102, 143, 54, 65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 187, 208, 89,
 115       18, 169, 200, 196,
 116    135, 130, 116, 188, 159, 86, 164, 100, 109, 198, 173, 186, 3, 64, 52, 217,
 117       226, 250, 124, 123,
 118    5, 202, 38, 147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58,
 119       17, 182, 189, 28, 42,
 120    223, 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101, 155,
 121       167, 43, 172, 9,
 122    129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, 112, 104,
 123       218, 246, 97, 228,
 124    251, 34, 242, 193, 238, 210, 144, 12, 191, 179, 162, 241, 81, 51, 145, 235,
 125       249, 14, 239, 107,
 126    49, 192, 214, 31, 181, 199, 106, 157, 184, 84, 204, 176, 115, 121, 50, 45,
 127       127, 4, 150, 254,
 128    138, 236, 205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66,
 129       215, 61, 156, 180
 130 };
 131
 132 /*
 133  * ---------------------------------------------------------------------
 134  */
 135
 136 /*
 137  * Helper functions to compute gradients-dot-residualvectors (1D to 4D)
 138  * Note that these generate gradients of more than unit length. To make
 139  * a close match with the value range of classic Perlin noise, the final
 140  * noise values need to be rescaled to fit nicely within [-1,1].
 141  * (The simplex noise functions as such also have different scaling.)
 142  * Note also that these noise functions are the most practical and useful
 143  * signed version of Perlin noise. To return values according to the
 144  * RenderMan specification from the SL noise() and pnoise() functions,
 145  * the noise values need to be scaled and offset to [0,1], like this:
 146  * float SLnoise = (SimplexNoise1234::noise(x,y,z) + 1.0) * 0.5;
 147  */
 148
 149 static float
 150 grad1(int hash, float x)
 151 {
 152    int h = hash & 15;
 153    float grad = 1.0f + (h & 7); /* Gradient value 1.0, 2.0, ..., 8.0 */
 154    if (h & 8)
 155       grad = -grad;             /* Set a random sign for the gradient */
 156    return (grad * x);           /* Multiply the gradient with the distance */
 157 }
 158
 159 static float
 160 grad2(int hash, float x, float y)
 161 {
 162    int h = hash & 7;            /* Convert low 3 bits of hash code */
 163    float u = h < 4 ? x : y;     /* into 8 simple gradient directions, */
 164    float v = h < 4 ? y : x;     /* and compute the dot product with (x,y). */
 165    return ((h & 1) ? -u : u) + ((h & 2) ? -2.0f * v : 2.0f * v);
 166 }
 167
 168 static float
 169 grad3(int hash, float x, float y, float z)
 170 {
 171    int h = hash & 15;           /* Convert low 4 bits of hash code into 12 simple */
 172    float u = h < 8 ? x : y;     /* gradient directions, and compute dot product. */
 173    float v = h < 4 ? y : h == 12 || h == 14 ? x : z;    /* Fix repeats at h = 12 to 15 */
 174    return ((h & 1) ? -u : u) + ((h & 2) ? -v : v);
 175 }
 176
 177 static float
 178 grad4(int hash, float x, float y, float z, float t)
 179 {
 180    int h = hash & 31;           /* Convert low 5 bits of hash code into 32 simple */
 181    float u = h < 24 ? x : y;    /* gradient directions, and compute dot product. */
 182    float v = h < 16 ? y : z;
 183    float w = h < 8 ? z : t;
 184    return ((h & 1) ? -u : u) + ((h & 2) ? -v : v) + ((h & 4) ? -w : w);
 185 }
 186
 187 /**
 188  * A lookup table to traverse the simplex around a given point in 4D.
 189  * Details can be found where this table is used, in the 4D noise method.
 190  * TODO: This should not be required, backport it from Bill's GLSL code!
 191  */
 192 static unsigned char simplex[64][4] = {
 193    {0, 1, 2, 3}, {0, 1, 3, 2}, {0, 0, 0, 0}, {0, 2, 3, 1},
 194    {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 2, 3, 0},
 195    {0, 2, 1, 3}, {0, 0, 0, 0}, {0, 3, 1, 2}, {0, 3, 2, 1},
 196    {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 3, 2, 0},
 197    {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
 198    {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
 199    {1, 2, 0, 3}, {0, 0, 0, 0}, {1, 3, 0, 2}, {0, 0, 0, 0},
 200    {0, 0, 0, 0}, {0, 0, 0, 0}, {2, 3, 0, 1}, {2, 3, 1, 0},
 201    {1, 0, 2, 3}, {1, 0, 3, 2}, {0, 0, 0, 0}, {0, 0, 0, 0},
 202    {0, 0, 0, 0}, {2, 0, 3, 1}, {0, 0, 0, 0}, {2, 1, 3, 0},
 203    {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
 204    {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
 205    {2, 0, 1, 3}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
 206    {3, 0, 1, 2}, {3, 0, 2, 1}, {0, 0, 0, 0}, {3, 1, 2, 0},
 207    {2, 1, 0, 3}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
 208    {3, 1, 0, 2}, {0, 0, 0, 0}, {3, 2, 0, 1}, {3, 2, 1, 0}
 209 };
 210
 211
 212 /** 1D simplex noise */
 213 GLfloat
 214 _mesa_noise1(GLfloat x)
 215 {
 216    int i0 = FASTFLOOR(x);
 217    int i1 = i0 + 1;
 218    float x0 = x - i0;
 219    float x1 = x0 - 1.0f;
 220    float t1 = 1.0f - x1 * x1;
 221    float n0, n1;
 222
 223    float t0 = 1.0f - x0 * x0;
 224 /*  if(t0 < 0.0f) t0 = 0.0f; // this never happens for the 1D case */
 225    t0 *= t0;
 226    n0 = t0 * t0 * grad1(perm[i0 & 0xff], x0);
 227
 228 /*  if(t1 < 0.0f) t1 = 0.0f; // this never happens for the 1D case */
 229    t1 *= t1;
 230    n1 = t1 * t1 * grad1(perm[i1 & 0xff], x1);
 231    /* The maximum value of this noise is 8*(3/4)^4 = 2.53125 */
 232    /* A factor of 0.395 would scale to fit exactly within [-1,1], but */
 233    /* we want to match PRMan's 1D noise, so we scale it down some more. */
 234    return 0.25f * (n0 + n1);
 235 }
 236
 237
 238 /** 2D simplex noise */
 239 GLfloat
 240 _mesa_noise2(GLfloat x, GLfloat y)
 241 {
 242 #define F2 0.366025403f         /* F2 = 0.5*(sqrt(3.0)-1.0) */
 243 #define G2 0.211324865f         /* G2 = (3.0-Math.sqrt(3.0))/6.0 */
 244
 245    float n0, n1, n2;            /* Noise contributions from the three corners */
 246
 247    /* Skew the input space to determine which simplex cell we're in */
 248    float s = (x + y) * F2;      /* Hairy factor for 2D */
 249    float xs = x + s;
 250    float ys = y + s;
 251    int i = FASTFLOOR(xs);
 252    int j = FASTFLOOR(ys);
 253
 254    float t = (float) (i + j) * G2;
 255    float X0 = i - t;            /* Unskew the cell origin back to (x,y) space */
 256    float Y0 = j - t;
 257    float x0 = x - X0;           /* The x,y distances from the cell origin */
 258    float y0 = y - Y0;
 259
 260    float x1, y1, x2, y2;
 261    int ii, jj;
 262    float t0, t1, t2;
 263
 264    /* For the 2D case, the simplex shape is an equilateral triangle. */
 265    /* Determine which simplex we are in. */
 266    int i1, j1;                  /* Offsets for second (middle) corner of simplex in (i,j) coords */
 267    if (x0 > y0) {
 268       i1 = 1;
 269       j1 = 0;
 270    }                            /* lower triangle, XY order: (0,0)->(1,0)->(1,1) */
 271    else {
 272       i1 = 0;
 273       j1 = 1;
 274    }                            /* upper triangle, YX order: (0,0)->(0,1)->(1,1) */
 275
 276    /* A step of (1,0) in (i,j) means a step of (1-c,-c) in (x,y), and */
 277    /* a step of (0,1) in (i,j) means a step of (-c,1-c) in (x,y), where */
 278    /* c = (3-sqrt(3))/6 */
 279
 280    x1 = x0 - i1 + G2;           /* Offsets for middle corner in (x,y) unskewed coords */
 281    y1 = y0 - j1 + G2;
 282    x2 = x0 - 1.0f + 2.0f * G2;  /* Offsets for last corner in (x,y) unskewed coords */
 283    y2 = y0 - 1.0f + 2.0f * G2;
 284
 285    /* Wrap the integer indices at 256, to avoid indexing perm[] out of bounds */
 286    ii = i % 256;
 287    jj = j % 256;
 288
 289    /* Calculate the contribution from the three corners */
 290    t0 = 0.5f - x0 * x0 - y0 * y0;
 291    if (t0 < 0.0f)
 292       n0 = 0.0f;
 293    else {
 294       t0 *= t0;
 295       n0 = t0 * t0 * grad2(perm[ii + perm[jj]], x0, y0);
 296    }
 297
 298    t1 = 0.5f - x1 * x1 - y1 * y1;
 299    if (t1 < 0.0f)
 300       n1 = 0.0f;
 301    else {
 302       t1 *= t1;
 303       n1 = t1 * t1 * grad2(perm[ii + i1 + perm[jj + j1]], x1, y1);
 304    }
 305
 306    t2 = 0.5f - x2 * x2 - y2 * y2;
 307    if (t2 < 0.0f)
 308       n2 = 0.0f;
 309    else {
 310       t2 *= t2;
 311       n2 = t2 * t2 * grad2(perm[ii + 1 + perm[jj + 1]], x2, y2);
 312    }
 313
 314    /* Add contributions from each corner to get the final noise value. */
 315    /* The result is scaled to return values in the interval [-1,1]. */
 316    return 40.0f * (n0 + n1 + n2);       /* TODO: The scale factor is preliminary! */
 317 }
 318
 319
 320 /** 3D simplex noise */
 321 GLfloat
 322 _mesa_noise3(GLfloat x, GLfloat y, GLfloat z)
 323 {
 324 /* Simple skewing factors for the 3D case */
 325 #define F3 0.333333333f
 326 #define G3 0.166666667f
 327
 328    float n0, n1, n2, n3;        /* Noise contributions from the four corners */
 329
 330    /* Skew the input space to determine which simplex cell we're in */
 331    float s = (x + y + z) * F3;  /* Very nice and simple skew factor for 3D */
 332    float xs = x + s;
 333    float ys = y + s;
 334    float zs = z + s;
 335    int i = FASTFLOOR(xs);
 336    int j = FASTFLOOR(ys);
 337    int k = FASTFLOOR(zs);
 338
 339    float t = (float) (i + j + k) * G3;
 340    float X0 = i - t;            /* Unskew the cell origin back to (x,y,z) space */
 341    float Y0 = j - t;
 342    float Z0 = k - t;
 343    float x0 = x - X0;           /* The x,y,z distances from the cell origin */
 344    float y0 = y - Y0;
 345    float z0 = z - Z0;
 346
 347    float x1, y1, z1, x2, y2, z2, x3, y3, z3;
 348    int ii, jj, kk;
 349    float t0, t1, t2, t3;
 350
 351    /* For the 3D case, the simplex shape is a slightly irregular tetrahedron. */
 352    /* Determine which simplex we are in. */
 353    int i1, j1, k1;              /* Offsets for second corner of simplex in (i,j,k) coords */
 354    int i2, j2, k2;              /* Offsets for third corner of simplex in (i,j,k) coords */
 355
 356 /* This code would benefit from a backport from the GLSL version! */
 357    if (x0 >= y0) {
 358       if (y0 >= z0) {
 359          i1 = 1;
 360          j1 = 0;
 361          k1 = 0;
 362          i2 = 1;
 363          j2 = 1;
 364          k2 = 0;
 365       }                         /* X Y Z order */
 366       else if (x0 >= z0) {
 367          i1 = 1;
 368          j1 = 0;
 369          k1 = 0;
 370          i2 = 1;
 371          j2 = 0;
 372          k2 = 1;
 373       }                         /* X Z Y order */
 374       else {
 375          i1 = 0;
 376          j1 = 0;
 377          k1 = 1;
 378          i2 = 1;
 379          j2 = 0;
 380          k2 = 1;
 381       }                         /* Z X Y order */
 382    }
 383    else {                       /* x0<y0 */
 384       if (y0 < z0) {
 385          i1 = 0;
 386          j1 = 0;
 387          k1 = 1;
 388          i2 = 0;
 389          j2 = 1;
 390          k2 = 1;
 391       }                         /* Z Y X order */
 392       else if (x0 < z0) {
 393          i1 = 0;
 394          j1 = 1;
 395          k1 = 0;
 396          i2 = 0;
 397          j2 = 1;
 398          k2 = 1;
 399       }                         /* Y Z X order */
 400       else {
 401          i1 = 0;
 402          j1 = 1;
 403          k1 = 0;
 404          i2 = 1;
 405          j2 = 1;
 406          k2 = 0;
 407       }                         /* Y X Z order */
 408    }
 409
 410    /* A step of (1,0,0) in (i,j,k) means a step of (1-c,-c,-c) in
 411     * (x,y,z), a step of (0,1,0) in (i,j,k) means a step of
 412     * (-c,1-c,-c) in (x,y,z), and a step of (0,0,1) in (i,j,k) means a
 413     * step of (-c,-c,1-c) in (x,y,z), where c = 1/6.
 414     */
 415
 416    x1 = x0 - i1 + G3;         /* Offsets for second corner in (x,y,z) coords */
 417    y1 = y0 - j1 + G3;
 418    z1 = z0 - k1 + G3;
 419    x2 = x0 - i2 + 2.0f * G3;  /* Offsets for third corner in (x,y,z) coords */
 420    y2 = y0 - j2 + 2.0f * G3;
 421    z2 = z0 - k2 + 2.0f * G3;
 422    x3 = x0 - 1.0f + 3.0f * G3;/* Offsets for last corner in (x,y,z) coords */
 423    y3 = y0 - 1.0f + 3.0f * G3;
 424    z3 = z0 - 1.0f + 3.0f * G3;
 425
 426    /* Wrap the integer indices at 256 to avoid indexing perm[] out of bounds */
 427    ii = i % 256;
 428    jj = j % 256;
 429    kk = k % 256;
 430
 431    /* Calculate the contribution from the four corners */
 432    t0 = 0.6f - x0 * x0 - y0 * y0 - z0 * z0;
 433    if (t0 < 0.0f)
 434       n0 = 0.0f;
 435    else {
 436       t0 *= t0;
 437       n0 = t0 * t0 * grad3(perm[ii + perm[jj + perm[kk]]], x0, y0, z0);
 438    }
 439
 440    t1 = 0.6f - x1 * x1 - y1 * y1 - z1 * z1;
 441    if (t1 < 0.0f)
 442       n1 = 0.0f;
 443    else {
 444       t1 *= t1;
 445       n1 =
 446          t1 * t1 * grad3(perm[ii + i1 + perm[jj + j1 + perm[kk + k1]]], x1,
 447                          y1, z1);
 448    }
 449
 450    t2 = 0.6f - x2 * x2 - y2 * y2 - z2 * z2;
 451    if (t2 < 0.0f)
 452       n2 = 0.0f;
 453    else {
 454       t2 *= t2;
 455       n2 =
 456          t2 * t2 * grad3(perm[ii + i2 + perm[jj + j2 + perm[kk + k2]]], x2,
 457                          y2, z2);
 458    }
 459
 460    t3 = 0.6f - x3 * x3 - y3 * y3 - z3 * z3;
 461    if (t3 < 0.0f)
 462       n3 = 0.0f;
 463    else {
 464       t3 *= t3;
 465       n3 =
 466          t3 * t3 * grad3(perm[ii + 1 + perm[jj + 1 + perm[kk + 1]]], x3, y3,
 467                          z3);
 468    }
 469
 470    /* Add contributions from each corner to get the final noise value.
 471     * The result is scaled to stay just inside [-1,1]
 472     */
 473    return 32.0f * (n0 + n1 + n2 + n3);  /* TODO: The scale factor is preliminary! */
 474 }
 475
 476
 477 /** 4D simplex noise */
 478 GLfloat
 479 _mesa_noise4(GLfloat x, GLfloat y, GLfloat z, GLfloat w)
 480 {
 481    /* The skewing and unskewing factors are hairy again for the 4D case */
 482 #define F4 0.309016994f         /* F4 = (Math.sqrt(5.0)-1.0)/4.0 */
 483 #define G4 0.138196601f         /* G4 = (5.0-Math.sqrt(5.0))/20.0 */
 484
 485    float n0, n1, n2, n3, n4;    /* Noise contributions from the five corners */
 486
 487    /* Skew the (x,y,z,w) space to determine which cell of 24 simplices we're in */
 488    float s = (x + y + z + w) * F4;      /* Factor for 4D skewing */
 489    float xs = x + s;
 490    float ys = y + s;
 491    float zs = z + s;
 492    float ws = w + s;
 493    int i = FASTFLOOR(xs);
 494    int j = FASTFLOOR(ys);
 495    int k = FASTFLOOR(zs);
 496    int l = FASTFLOOR(ws);
 497
 498    float t = (i + j + k + l) * G4;      /* Factor for 4D unskewing */
 499    float X0 = i - t;            /* Unskew the cell origin back to (x,y,z,w) space */
 500    float Y0 = j - t;
 501    float Z0 = k - t;
 502    float W0 = l - t;
 503
 504    float x0 = x - X0;           /* The x,y,z,w distances from the cell origin */
 505    float y0 = y - Y0;
 506    float z0 = z - Z0;
 507    float w0 = w - W0;
 508
 509    /* For the 4D case, the simplex is a 4D shape I won't even try to describe.
 510     * To find out which of the 24 possible simplices we're in, we need to
 511     * determine the magnitude ordering of x0, y0, z0 and w0.
 512     * The method below is a good way of finding the ordering of x,y,z,w and
 513     * then find the correct traversal order for the simplex we're in.
 514     * First, six pair-wise comparisons are performed between each possible pair
 515     * of the four coordinates, and the results are used to add up binary bits
 516     * for an integer index.
 517     */
 518    int c1 = (x0 > y0) ? 32 : 0;
 519    int c2 = (x0 > z0) ? 16 : 0;
 520    int c3 = (y0 > z0) ? 8 : 0;
 521    int c4 = (x0 > w0) ? 4 : 0;
 522    int c5 = (y0 > w0) ? 2 : 0;
 523    int c6 = (z0 > w0) ? 1 : 0;
 524    int c = c1 + c2 + c3 + c4 + c5 + c6;
 525
 526    int i1, j1, k1, l1;  /* The integer offsets for the second simplex corner */
 527    int i2, j2, k2, l2;  /* The integer offsets for the third simplex corner */
 528    int i3, j3, k3, l3;  /* The integer offsets for the fourth simplex corner */
 529
 530    float x1, y1, z1, w1, x2, y2, z2, w2, x3, y3, z3, w3, x4, y4, z4, w4;
 531    int ii, jj, kk, ll;
 532    float t0, t1, t2, t3, t4;
 533
 534    /*
 535     * simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some
 536     * order.  Many values of c will never occur, since e.g. x>y>z>w
 537     * makes x<z, y<w and x<w impossible. Only the 24 indices which
 538     * have non-zero entries make any sense.  We use a thresholding to
 539     * set the coordinates in turn from the largest magnitude.  The
 540     * number 3 in the "simplex" array is at the position of the
 541     * largest coordinate.
 542     */
 543    i1 = simplex[c][0] >= 3 ? 1 : 0;
 544    j1 = simplex[c][1] >= 3 ? 1 : 0;
 545    k1 = simplex[c][2] >= 3 ? 1 : 0;
 546    l1 = simplex[c][3] >= 3 ? 1 : 0;
 547    /* The number 2 in the "simplex" array is at the second largest coordinate. */
 548    i2 = simplex[c][0] >= 2 ? 1 : 0;
 549    j2 = simplex[c][1] >= 2 ? 1 : 0;
 550    k2 = simplex[c][2] >= 2 ? 1 : 0;
 551    l2 = simplex[c][3] >= 2 ? 1 : 0;
 552    /* The number 1 in the "simplex" array is at the second smallest coordinate. */
 553    i3 = simplex[c][0] >= 1 ? 1 : 0;
 554    j3 = simplex[c][1] >= 1 ? 1 : 0;
 555    k3 = simplex[c][2] >= 1 ? 1 : 0;
 556    l3 = simplex[c][3] >= 1 ? 1 : 0;
 557    /* The fifth corner has all coordinate offsets = 1, so no need to look that up. */
 558
 559    x1 = x0 - i1 + G4;           /* Offsets for second corner in (x,y,z,w) coords */
 560    y1 = y0 - j1 + G4;
 561    z1 = z0 - k1 + G4;
 562    w1 = w0 - l1 + G4;
 563    x2 = x0 - i2 + 2.0f * G4;    /* Offsets for third corner in (x,y,z,w) coords */
 564    y2 = y0 - j2 + 2.0f * G4;
 565    z2 = z0 - k2 + 2.0f * G4;
 566    w2 = w0 - l2 + 2.0f * G4;
 567    x3 = x0 - i3 + 3.0f * G4;    /* Offsets for fourth corner in (x,y,z,w) coords */
 568    y3 = y0 - j3 + 3.0f * G4;
 569    z3 = z0 - k3 + 3.0f * G4;
 570    w3 = w0 - l3 + 3.0f * G4;
 571    x4 = x0 - 1.0f + 4.0f * G4;  /* Offsets for last corner in (x,y,z,w) coords */
 572    y4 = y0 - 1.0f + 4.0f * G4;
 573    z4 = z0 - 1.0f + 4.0f * G4;
 574    w4 = w0 - 1.0f + 4.0f * G4;
 575
 576    /* Wrap the integer indices at 256, to avoid indexing perm[] out of bounds */
 577    ii = i % 256;
 578    jj = j % 256;
 579    kk = k % 256;
 580    ll = l % 256;
 581
 582    /* Calculate the contribution from the five corners */
 583    t0 = 0.6f - x0 * x0 - y0 * y0 - z0 * z0 - w0 * w0;
 584    if (t0 < 0.0f)
 585       n0 = 0.0f;
 586    else {
 587       t0 *= t0;
 588       n0 =
 589          t0 * t0 * grad4(perm[ii + perm[jj + perm[kk + perm[ll]]]], x0, y0,
 590                          z0, w0);
 591    }
 592
 593    t1 = 0.6f - x1 * x1 - y1 * y1 - z1 * z1 - w1 * w1;
 594    if (t1 < 0.0f)
 595       n1 = 0.0f;
 596    else {
 597       t1 *= t1;
 598       n1 =
 599          t1 * t1 *
 600          grad4(perm[ii + i1 + perm[jj + j1 + perm[kk + k1 + perm[ll + l1]]]],
 601                x1, y1, z1, w1);
 602    }
 603
 604    t2 = 0.6f - x2 * x2 - y2 * y2 - z2 * z2 - w2 * w2;
 605    if (t2 < 0.0f)
 606       n2 = 0.0f;
 607    else {
 608       t2 *= t2;
 609       n2 =
 610          t2 * t2 *
 611          grad4(perm[ii + i2 + perm[jj + j2 + perm[kk + k2 + perm[ll + l2]]]],
 612                x2, y2, z2, w2);
 613    }
 614
 615    t3 = 0.6f - x3 * x3 - y3 * y3 - z3 * z3 - w3 * w3;
 616    if (t3 < 0.0f)
 617       n3 = 0.0f;
 618    else {
 619       t3 *= t3;
 620       n3 =
 621          t3 * t3 *
 622          grad4(perm[ii + i3 + perm[jj + j3 + perm[kk + k3 + perm[ll + l3]]]],
 623                x3, y3, z3, w3);
 624    }
 625
 626    t4 = 0.6f - x4 * x4 - y4 * y4 - z4 * z4 - w4 * w4;
 627    if (t4 < 0.0f)
 628       n4 = 0.0f;
 629    else {
 630       t4 *= t4;
 631       n4 =
 632          t4 * t4 *
 633          grad4(perm[ii + 1 + perm[jj + 1 + perm[kk + 1 + perm[ll + 1]]]], x4,
 634                y4, z4, w4);
 635    }
 636
 637    /* Sum up and scale the result to cover the range [-1,1] */
 638    return 27.0f * (n0 + n1 + n2 + n3 + n4);     /* TODO: The scale factor is preliminary! */
 639 }