1 /* 2 * Mesa 3-D graphics library 3 * 4 * Copyright (C) 2006 Brian Paul All Rights Reserved. 5 * 6 * Permission is hereby granted, free of charge, to any person obtaining a 7 * copy of this software and associated documentation files (the "Software"), 8 * to deal in the Software without restriction, including without limitation 9 * the rights to use, copy, modify, merge, publish, distribute, sublicense, 10 * and/or sell copies of the Software, and to permit persons to whom the 11 * Software is furnished to do so, subject to the following conditions: 12 * 13 * The above copyright notice and this permission notice shall be included 14 * in all copies or substantial portions of the Software. 15 * 16 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 17 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 18 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 19 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR 20 * OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 21 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 22 * 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 static const 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 const 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 unsigned 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 unsigned 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 & 0xff; 286 jj = j & 0xff; 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 unsigned 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 unsigned int i1, j1, k1; /* Offsets for second corner of simplex in (i,j,k) coords */ 353 unsigned 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 & 0xff; 427 jj = j & 0xff; 428 kk = k & 0xff; 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 unsigned int i1, j1, k1, l1; /* The integer offsets for the second simplex corner */ 526 unsigned int i2, j2, k2, l2; /* The integer offsets for the third simplex corner */ 527 unsigned 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 unsigned 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 & 0xff; 577 jj = j & 0xff; 578 kk = k & 0xff; 579 ll = l & 0xff; 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
