Home | History | Annotate | Download | only in program
      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 (at) 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 (at) 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 }
    639