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