1 module voxelman.math.opensimplex; 2 3 /* 4 * OpenSimplex Noise in Java. 5 * by Kurt Spencer 6 * 7 * v1.1 (October 5, 2014) 8 * - Added 2D and 4D implementations. 9 * - Proper gradient sets for all dimensions, from a 10 * dimensionally-generalizable scheme with an actual 11 * rhyme and reason behind it. 12 * - Removed default permutation array in favor of 13 * default seed. 14 * - Changed seed-based constructor to be independent 15 * of any particular randomization library, so results 16 * will be the same when ported to other languages. 17 */ 18 19 public struct OpenSimplexNoise { 20 21 private enum double STRETCH_CONSTANT_2D = -0.211324865405187; //(1/Math.sqrt(2+1)-1)/2; 22 private enum double SQUISH_CONSTANT_2D = 0.366025403784439; //(Math.sqrt(2+1)-1)/2; 23 private enum double STRETCH_CONSTANT_3D = -1.0 / 6; //(1/Math.sqrt(3+1)-1)/3; 24 private enum double SQUISH_CONSTANT_3D = 1.0 / 3; //(Math.sqrt(3+1)-1)/3; 25 private enum double STRETCH_CONSTANT_4D = -0.138196601125011; //(1/Math.sqrt(4+1)-1)/4; 26 private enum double SQUISH_CONSTANT_4D = 0.309016994374947; //(Math.sqrt(4+1)-1)/4; 27 28 private enum double NORM_CONSTANT_2D = 47; 29 private enum double NORM_CONSTANT_3D = 103; 30 private enum double NORM_CONSTANT_4D = 30; 31 32 private enum long DEFAULT_SEED = 0; 33 34 private short[256] perm; 35 private short[256] permGradIndex3D; 36 37 //public this() { 38 // this(DEFAULT_SEED); 39 //} 40 41 public this(short[] perm) { 42 this.perm = perm; 43 44 for (int i = 0; i < 256; i++) { 45 //Since 3D has 24 gradients, simple bitmask won't work, so precompute modulo array. 46 permGradIndex3D[i] = cast(short)((perm[i] % (gradients3D.length / 3)) * 3); 47 } 48 } 49 50 //Initializes the class using a permutation array generated from a 64-bit seed. 51 //Generates a proper permutation (i.e. doesn't merely perform N successive pair swaps on a base array) 52 //Uses a simple 64-bit LCG. 53 public this(long seed) { 54 short[256] source; 55 for (short i = 0; i < 256; i++) 56 source[i] = i; 57 seed = seed * 6364136223846793005L + 1442695040888963407L; 58 seed = seed * 6364136223846793005L + 1442695040888963407L; 59 seed = seed * 6364136223846793005L + 1442695040888963407L; 60 for (int i = 255; i >= 0; i--) { 61 seed = seed * 6364136223846793005L + 1442695040888963407L; 62 int r = cast(int)((seed + 31) % (i + 1)); 63 if (r < 0) 64 r += (i + 1); 65 perm[i] = source[r]; 66 permGradIndex3D[i] = cast(short)((perm[i] % (gradients3D.length / 3)) * 3); 67 source[r] = source[i]; 68 } 69 } 70 71 //2D OpenSimplex Noise. 72 public double noise(double x, double y) { 73 74 //Place input coordinates onto grid. 75 double stretchOffset = (x + y) * STRETCH_CONSTANT_2D; 76 double xs = x + stretchOffset; 77 double ys = y + stretchOffset; 78 79 //Floor to get grid coordinates of rhombus (stretched square) super-cell origin. 80 int xsb = fastFloor(xs); 81 int ysb = fastFloor(ys); 82 83 //Skew out to get actual coordinates of rhombus origin. We'll need these later. 84 double squishOffset = (xsb + ysb) * SQUISH_CONSTANT_2D; 85 double xb = xsb + squishOffset; 86 double yb = ysb + squishOffset; 87 88 //Compute grid coordinates relative to rhombus origin. 89 double xins = xs - xsb; 90 double yins = ys - ysb; 91 92 //Sum those together to get a value that determines which region we're in. 93 double inSum = xins + yins; 94 95 //Positions relative to origin point. 96 double dx0 = x - xb; 97 double dy0 = y - yb; 98 99 //We'll be defining these inside the next block and using them afterwards. 100 double dx_ext, dy_ext; 101 int xsv_ext, ysv_ext; 102 103 double value = 0; 104 105 //Contribution (1,0) 106 double dx1 = dx0 - 1 - SQUISH_CONSTANT_2D; 107 double dy1 = dy0 - 0 - SQUISH_CONSTANT_2D; 108 double attn1 = 2 - dx1 * dx1 - dy1 * dy1; 109 if (attn1 > 0) { 110 attn1 *= attn1; 111 value += attn1 * attn1 * extrapolate(xsb + 1, ysb + 0, dx1, dy1); 112 } 113 114 //Contribution (0,1) 115 double dx2 = dx0 - 0 - SQUISH_CONSTANT_2D; 116 double dy2 = dy0 - 1 - SQUISH_CONSTANT_2D; 117 double attn2 = 2 - dx2 * dx2 - dy2 * dy2; 118 if (attn2 > 0) { 119 attn2 *= attn2; 120 value += attn2 * attn2 * extrapolate(xsb + 0, ysb + 1, dx2, dy2); 121 } 122 123 if (inSum <= 1) { //We're inside the triangle (2-Simplex) at (0,0) 124 double zins = 1 - inSum; 125 if (zins > xins || zins > yins) { //(0,0) is one of the closest two triangular vertices 126 if (xins > yins) { 127 xsv_ext = xsb + 1; 128 ysv_ext = ysb - 1; 129 dx_ext = dx0 - 1; 130 dy_ext = dy0 + 1; 131 } else { 132 xsv_ext = xsb - 1; 133 ysv_ext = ysb + 1; 134 dx_ext = dx0 + 1; 135 dy_ext = dy0 - 1; 136 } 137 } else { //(1,0) and (0,1) are the closest two vertices. 138 xsv_ext = xsb + 1; 139 ysv_ext = ysb + 1; 140 dx_ext = dx0 - 1 - 2 * SQUISH_CONSTANT_2D; 141 dy_ext = dy0 - 1 - 2 * SQUISH_CONSTANT_2D; 142 } 143 } else { //We're inside the triangle (2-Simplex) at (1,1) 144 double zins = 2 - inSum; 145 if (zins < xins || zins < yins) { //(0,0) is one of the closest two triangular vertices 146 if (xins > yins) { 147 xsv_ext = xsb + 2; 148 ysv_ext = ysb + 0; 149 dx_ext = dx0 - 2 - 2 * SQUISH_CONSTANT_2D; 150 dy_ext = dy0 + 0 - 2 * SQUISH_CONSTANT_2D; 151 } else { 152 xsv_ext = xsb + 0; 153 ysv_ext = ysb + 2; 154 dx_ext = dx0 + 0 - 2 * SQUISH_CONSTANT_2D; 155 dy_ext = dy0 - 2 - 2 * SQUISH_CONSTANT_2D; 156 } 157 } else { //(1,0) and (0,1) are the closest two vertices. 158 dx_ext = dx0; 159 dy_ext = dy0; 160 xsv_ext = xsb; 161 ysv_ext = ysb; 162 } 163 xsb += 1; 164 ysb += 1; 165 dx0 = dx0 - 1 - 2 * SQUISH_CONSTANT_2D; 166 dy0 = dy0 - 1 - 2 * SQUISH_CONSTANT_2D; 167 } 168 169 //Contribution (0,0) or (1,1) 170 double attn0 = 2 - dx0 * dx0 - dy0 * dy0; 171 if (attn0 > 0) { 172 attn0 *= attn0; 173 value += attn0 * attn0 * extrapolate(xsb, ysb, dx0, dy0); 174 } 175 176 //Extra Vertex 177 double attn_ext = 2 - dx_ext * dx_ext - dy_ext * dy_ext; 178 if (attn_ext > 0) { 179 attn_ext *= attn_ext; 180 value += attn_ext * attn_ext * extrapolate(xsv_ext, ysv_ext, dx_ext, dy_ext); 181 } 182 183 return value / NORM_CONSTANT_2D; 184 } 185 186 //3D OpenSimplex Noise. 187 public double noise(double x, double y, double z) { 188 189 //Place input coordinates on simplectic honeycomb. 190 double stretchOffset = (x + y + z) * STRETCH_CONSTANT_3D; 191 double xs = x + stretchOffset; 192 double ys = y + stretchOffset; 193 double zs = z + stretchOffset; 194 195 //Floor to get simplectic honeycomb coordinates of rhombohedron (stretched cube) super-cell origin. 196 int xsb = fastFloor(xs); 197 int ysb = fastFloor(ys); 198 int zsb = fastFloor(zs); 199 200 //Skew out to get actual coordinates of rhombohedron origin. We'll need these later. 201 double squishOffset = (xsb + ysb + zsb) * SQUISH_CONSTANT_3D; 202 double xb = xsb + squishOffset; 203 double yb = ysb + squishOffset; 204 double zb = zsb + squishOffset; 205 206 //Compute simplectic honeycomb coordinates relative to rhombohedral origin. 207 double xins = xs - xsb; 208 double yins = ys - ysb; 209 double zins = zs - zsb; 210 211 //Sum those together to get a value that determines which region we're in. 212 double inSum = xins + yins + zins; 213 214 //Positions relative to origin point. 215 double dx0 = x - xb; 216 double dy0 = y - yb; 217 double dz0 = z - zb; 218 219 //We'll be defining these inside the next block and using them afterwards. 220 double dx_ext0, dy_ext0, dz_ext0; 221 double dx_ext1, dy_ext1, dz_ext1; 222 int xsv_ext0, ysv_ext0, zsv_ext0; 223 int xsv_ext1, ysv_ext1, zsv_ext1; 224 225 double value = 0; 226 if (inSum <= 1) { //We're inside the tetrahedron (3-Simplex) at (0,0,0) 227 228 //Determine which two of (0,0,1), (0,1,0), (1,0,0) are closest. 229 byte aPoint = 0x01; 230 double aScore = xins; 231 byte bPoint = 0x02; 232 double bScore = yins; 233 if (aScore >= bScore && zins > bScore) { 234 bScore = zins; 235 bPoint = 0x04; 236 } else if (aScore < bScore && zins > aScore) { 237 aScore = zins; 238 aPoint = 0x04; 239 } 240 241 //Now we determine the two lattice points not part of the tetrahedron that may contribute. 242 //This depends on the closest two tetrahedral vertices, including (0,0,0) 243 double wins = 1 - inSum; 244 if (wins > aScore || wins > bScore) { //(0,0,0) is one of the closest two tetrahedral vertices. 245 byte c = (bScore > aScore ? bPoint : aPoint); //Our other closest vertex is the closest out of a and b. 246 247 if ((c & 0x01) == 0) { 248 xsv_ext0 = xsb - 1; 249 xsv_ext1 = xsb; 250 dx_ext0 = dx0 + 1; 251 dx_ext1 = dx0; 252 } else { 253 xsv_ext0 = xsv_ext1 = xsb + 1; 254 dx_ext0 = dx_ext1 = dx0 - 1; 255 } 256 257 if ((c & 0x02) == 0) { 258 ysv_ext0 = ysv_ext1 = ysb; 259 dy_ext0 = dy_ext1 = dy0; 260 if ((c & 0x01) == 0) { 261 ysv_ext1 -= 1; 262 dy_ext1 += 1; 263 } else { 264 ysv_ext0 -= 1; 265 dy_ext0 += 1; 266 } 267 } else { 268 ysv_ext0 = ysv_ext1 = ysb + 1; 269 dy_ext0 = dy_ext1 = dy0 - 1; 270 } 271 272 if ((c & 0x04) == 0) { 273 zsv_ext0 = zsb; 274 zsv_ext1 = zsb - 1; 275 dz_ext0 = dz0; 276 dz_ext1 = dz0 + 1; 277 } else { 278 zsv_ext0 = zsv_ext1 = zsb + 1; 279 dz_ext0 = dz_ext1 = dz0 - 1; 280 } 281 } else { //(0,0,0) is not one of the closest two tetrahedral vertices. 282 byte c = cast(byte)(aPoint | bPoint); //Our two extra vertices are determined by the closest two. 283 284 if ((c & 0x01) == 0) { 285 xsv_ext0 = xsb; 286 xsv_ext1 = xsb - 1; 287 dx_ext0 = dx0 - 2 * SQUISH_CONSTANT_3D; 288 dx_ext1 = dx0 + 1 - SQUISH_CONSTANT_3D; 289 } else { 290 xsv_ext0 = xsv_ext1 = xsb + 1; 291 dx_ext0 = dx0 - 1 - 2 * SQUISH_CONSTANT_3D; 292 dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_3D; 293 } 294 295 if ((c & 0x02) == 0) { 296 ysv_ext0 = ysb; 297 ysv_ext1 = ysb - 1; 298 dy_ext0 = dy0 - 2 * SQUISH_CONSTANT_3D; 299 dy_ext1 = dy0 + 1 - SQUISH_CONSTANT_3D; 300 } else { 301 ysv_ext0 = ysv_ext1 = ysb + 1; 302 dy_ext0 = dy0 - 1 - 2 * SQUISH_CONSTANT_3D; 303 dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_3D; 304 } 305 306 if ((c & 0x04) == 0) { 307 zsv_ext0 = zsb; 308 zsv_ext1 = zsb - 1; 309 dz_ext0 = dz0 - 2 * SQUISH_CONSTANT_3D; 310 dz_ext1 = dz0 + 1 - SQUISH_CONSTANT_3D; 311 } else { 312 zsv_ext0 = zsv_ext1 = zsb + 1; 313 dz_ext0 = dz0 - 1 - 2 * SQUISH_CONSTANT_3D; 314 dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_3D; 315 } 316 } 317 318 //Contribution (0,0,0) 319 double attn0 = 2 - dx0 * dx0 - dy0 * dy0 - dz0 * dz0; 320 if (attn0 > 0) { 321 attn0 *= attn0; 322 value += attn0 * attn0 * extrapolate(xsb + 0, ysb + 0, zsb + 0, dx0, dy0, dz0); 323 } 324 325 //Contribution (1,0,0) 326 double dx1 = dx0 - 1 - SQUISH_CONSTANT_3D; 327 double dy1 = dy0 - 0 - SQUISH_CONSTANT_3D; 328 double dz1 = dz0 - 0 - SQUISH_CONSTANT_3D; 329 double attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1; 330 if (attn1 > 0) { 331 attn1 *= attn1; 332 value += attn1 * attn1 * extrapolate(xsb + 1, ysb + 0, zsb + 0, dx1, dy1, dz1); 333 } 334 335 //Contribution (0,1,0) 336 double dx2 = dx0 - 0 - SQUISH_CONSTANT_3D; 337 double dy2 = dy0 - 1 - SQUISH_CONSTANT_3D; 338 double dz2 = dz1; 339 double attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2; 340 if (attn2 > 0) { 341 attn2 *= attn2; 342 value += attn2 * attn2 * extrapolate(xsb + 0, ysb + 1, zsb + 0, dx2, dy2, dz2); 343 } 344 345 //Contribution (0,0,1) 346 double dx3 = dx2; 347 double dy3 = dy1; 348 double dz3 = dz0 - 1 - SQUISH_CONSTANT_3D; 349 double attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3; 350 if (attn3 > 0) { 351 attn3 *= attn3; 352 value += attn3 * attn3 * extrapolate(xsb + 0, ysb + 0, zsb + 1, dx3, dy3, dz3); 353 } 354 } else if (inSum >= 2) { //We're inside the tetrahedron (3-Simplex) at (1,1,1) 355 356 //Determine which two tetrahedral vertices are the closest, out of (1,1,0), (1,0,1), (0,1,1) but not (1,1,1). 357 byte aPoint = 0x06; 358 double aScore = xins; 359 byte bPoint = 0x05; 360 double bScore = yins; 361 if (aScore <= bScore && zins < bScore) { 362 bScore = zins; 363 bPoint = 0x03; 364 } else if (aScore > bScore && zins < aScore) { 365 aScore = zins; 366 aPoint = 0x03; 367 } 368 369 //Now we determine the two lattice points not part of the tetrahedron that may contribute. 370 //This depends on the closest two tetrahedral vertices, including (1,1,1) 371 double wins = 3 - inSum; 372 if (wins < aScore || wins < bScore) { //(1,1,1) is one of the closest two tetrahedral vertices. 373 byte c = (bScore < aScore ? bPoint : aPoint); //Our other closest vertex is the closest out of a and b. 374 375 if ((c & 0x01) != 0) { 376 xsv_ext0 = xsb + 2; 377 xsv_ext1 = xsb + 1; 378 dx_ext0 = dx0 - 2 - 3 * SQUISH_CONSTANT_3D; 379 dx_ext1 = dx0 - 1 - 3 * SQUISH_CONSTANT_3D; 380 } else { 381 xsv_ext0 = xsv_ext1 = xsb; 382 dx_ext0 = dx_ext1 = dx0 - 3 * SQUISH_CONSTANT_3D; 383 } 384 385 if ((c & 0x02) != 0) { 386 ysv_ext0 = ysv_ext1 = ysb + 1; 387 dy_ext0 = dy_ext1 = dy0 - 1 - 3 * SQUISH_CONSTANT_3D; 388 if ((c & 0x01) != 0) { 389 ysv_ext1 += 1; 390 dy_ext1 -= 1; 391 } else { 392 ysv_ext0 += 1; 393 dy_ext0 -= 1; 394 } 395 } else { 396 ysv_ext0 = ysv_ext1 = ysb; 397 dy_ext0 = dy_ext1 = dy0 - 3 * SQUISH_CONSTANT_3D; 398 } 399 400 if ((c & 0x04) != 0) { 401 zsv_ext0 = zsb + 1; 402 zsv_ext1 = zsb + 2; 403 dz_ext0 = dz0 - 1 - 3 * SQUISH_CONSTANT_3D; 404 dz_ext1 = dz0 - 2 - 3 * SQUISH_CONSTANT_3D; 405 } else { 406 zsv_ext0 = zsv_ext1 = zsb; 407 dz_ext0 = dz_ext1 = dz0 - 3 * SQUISH_CONSTANT_3D; 408 } 409 } else { //(1,1,1) is not one of the closest two tetrahedral vertices. 410 byte c = cast(byte)(aPoint & bPoint); //Our two extra vertices are determined by the closest two. 411 412 if ((c & 0x01) != 0) { 413 xsv_ext0 = xsb + 1; 414 xsv_ext1 = xsb + 2; 415 dx_ext0 = dx0 - 1 - SQUISH_CONSTANT_3D; 416 dx_ext1 = dx0 - 2 - 2 * SQUISH_CONSTANT_3D; 417 } else { 418 xsv_ext0 = xsv_ext1 = xsb; 419 dx_ext0 = dx0 - SQUISH_CONSTANT_3D; 420 dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D; 421 } 422 423 if ((c & 0x02) != 0) { 424 ysv_ext0 = ysb + 1; 425 ysv_ext1 = ysb + 2; 426 dy_ext0 = dy0 - 1 - SQUISH_CONSTANT_3D; 427 dy_ext1 = dy0 - 2 - 2 * SQUISH_CONSTANT_3D; 428 } else { 429 ysv_ext0 = ysv_ext1 = ysb; 430 dy_ext0 = dy0 - SQUISH_CONSTANT_3D; 431 dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D; 432 } 433 434 if ((c & 0x04) != 0) { 435 zsv_ext0 = zsb + 1; 436 zsv_ext1 = zsb + 2; 437 dz_ext0 = dz0 - 1 - SQUISH_CONSTANT_3D; 438 dz_ext1 = dz0 - 2 - 2 * SQUISH_CONSTANT_3D; 439 } else { 440 zsv_ext0 = zsv_ext1 = zsb; 441 dz_ext0 = dz0 - SQUISH_CONSTANT_3D; 442 dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D; 443 } 444 } 445 446 //Contribution (1,1,0) 447 double dx3 = dx0 - 1 - 2 * SQUISH_CONSTANT_3D; 448 double dy3 = dy0 - 1 - 2 * SQUISH_CONSTANT_3D; 449 double dz3 = dz0 - 0 - 2 * SQUISH_CONSTANT_3D; 450 double attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3; 451 if (attn3 > 0) { 452 attn3 *= attn3; 453 value += attn3 * attn3 * extrapolate(xsb + 1, ysb + 1, zsb + 0, dx3, dy3, dz3); 454 } 455 456 //Contribution (1,0,1) 457 double dx2 = dx3; 458 double dy2 = dy0 - 0 - 2 * SQUISH_CONSTANT_3D; 459 double dz2 = dz0 - 1 - 2 * SQUISH_CONSTANT_3D; 460 double attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2; 461 if (attn2 > 0) { 462 attn2 *= attn2; 463 value += attn2 * attn2 * extrapolate(xsb + 1, ysb + 0, zsb + 1, dx2, dy2, dz2); 464 } 465 466 //Contribution (0,1,1) 467 double dx1 = dx0 - 0 - 2 * SQUISH_CONSTANT_3D; 468 double dy1 = dy3; 469 double dz1 = dz2; 470 double attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1; 471 if (attn1 > 0) { 472 attn1 *= attn1; 473 value += attn1 * attn1 * extrapolate(xsb + 0, ysb + 1, zsb + 1, dx1, dy1, dz1); 474 } 475 476 //Contribution (1,1,1) 477 dx0 = dx0 - 1 - 3 * SQUISH_CONSTANT_3D; 478 dy0 = dy0 - 1 - 3 * SQUISH_CONSTANT_3D; 479 dz0 = dz0 - 1 - 3 * SQUISH_CONSTANT_3D; 480 double attn0 = 2 - dx0 * dx0 - dy0 * dy0 - dz0 * dz0; 481 if (attn0 > 0) { 482 attn0 *= attn0; 483 value += attn0 * attn0 * extrapolate(xsb + 1, ysb + 1, zsb + 1, dx0, dy0, dz0); 484 } 485 } else { //We're inside the octahedron (Rectified 3-Simplex) in between. 486 double aScore; 487 byte aPoint; 488 bool aIsFurtherSide; 489 double bScore; 490 byte bPoint; 491 bool bIsFurtherSide; 492 493 //Decide between point (0,0,1) and (1,1,0) as closest 494 double p1 = xins + yins; 495 if (p1 > 1) { 496 aScore = p1 - 1; 497 aPoint = 0x03; 498 aIsFurtherSide = true; 499 } else { 500 aScore = 1 - p1; 501 aPoint = 0x04; 502 aIsFurtherSide = false; 503 } 504 505 //Decide between point (0,1,0) and (1,0,1) as closest 506 double p2 = xins + zins; 507 if (p2 > 1) { 508 bScore = p2 - 1; 509 bPoint = 0x05; 510 bIsFurtherSide = true; 511 } else { 512 bScore = 1 - p2; 513 bPoint = 0x02; 514 bIsFurtherSide = false; 515 } 516 517 //The closest out of the two (1,0,0) and (0,1,1) will replace the furthest out of the two decided above, if closer. 518 double p3 = yins + zins; 519 if (p3 > 1) { 520 double score = p3 - 1; 521 if (aScore <= bScore && aScore < score) { 522 aScore = score; 523 aPoint = 0x06; 524 aIsFurtherSide = true; 525 } else if (aScore > bScore && bScore < score) { 526 bScore = score; 527 bPoint = 0x06; 528 bIsFurtherSide = true; 529 } 530 } else { 531 double score = 1 - p3; 532 if (aScore <= bScore && aScore < score) { 533 aScore = score; 534 aPoint = 0x01; 535 aIsFurtherSide = false; 536 } else if (aScore > bScore && bScore < score) { 537 bScore = score; 538 bPoint = 0x01; 539 bIsFurtherSide = false; 540 } 541 } 542 543 //Where each of the two closest points are determines how the extra two vertices are calculated. 544 if (aIsFurtherSide == bIsFurtherSide) { 545 if (aIsFurtherSide) { //Both closest points on (1,1,1) side 546 547 //One of the two extra points is (1,1,1) 548 dx_ext0 = dx0 - 1 - 3 * SQUISH_CONSTANT_3D; 549 dy_ext0 = dy0 - 1 - 3 * SQUISH_CONSTANT_3D; 550 dz_ext0 = dz0 - 1 - 3 * SQUISH_CONSTANT_3D; 551 xsv_ext0 = xsb + 1; 552 ysv_ext0 = ysb + 1; 553 zsv_ext0 = zsb + 1; 554 555 //Other extra point is based on the shared axis. 556 byte c = cast(byte)(aPoint & bPoint); 557 if ((c & 0x01) != 0) { 558 dx_ext1 = dx0 - 2 - 2 * SQUISH_CONSTANT_3D; 559 dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D; 560 dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D; 561 xsv_ext1 = xsb + 2; 562 ysv_ext1 = ysb; 563 zsv_ext1 = zsb; 564 } else if ((c & 0x02) != 0) { 565 dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D; 566 dy_ext1 = dy0 - 2 - 2 * SQUISH_CONSTANT_3D; 567 dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D; 568 xsv_ext1 = xsb; 569 ysv_ext1 = ysb + 2; 570 zsv_ext1 = zsb; 571 } else { 572 dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D; 573 dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D; 574 dz_ext1 = dz0 - 2 - 2 * SQUISH_CONSTANT_3D; 575 xsv_ext1 = xsb; 576 ysv_ext1 = ysb; 577 zsv_ext1 = zsb + 2; 578 } 579 } else {//Both closest points on (0,0,0) side 580 581 //One of the two extra points is (0,0,0) 582 dx_ext0 = dx0; 583 dy_ext0 = dy0; 584 dz_ext0 = dz0; 585 xsv_ext0 = xsb; 586 ysv_ext0 = ysb; 587 zsv_ext0 = zsb; 588 589 //Other extra point is based on the omitted axis. 590 byte c = cast(byte)(aPoint | bPoint); 591 if ((c & 0x01) == 0) { 592 dx_ext1 = dx0 + 1 - SQUISH_CONSTANT_3D; 593 dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_3D; 594 dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_3D; 595 xsv_ext1 = xsb - 1; 596 ysv_ext1 = ysb + 1; 597 zsv_ext1 = zsb + 1; 598 } else if ((c & 0x02) == 0) { 599 dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_3D; 600 dy_ext1 = dy0 + 1 - SQUISH_CONSTANT_3D; 601 dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_3D; 602 xsv_ext1 = xsb + 1; 603 ysv_ext1 = ysb - 1; 604 zsv_ext1 = zsb + 1; 605 } else { 606 dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_3D; 607 dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_3D; 608 dz_ext1 = dz0 + 1 - SQUISH_CONSTANT_3D; 609 xsv_ext1 = xsb + 1; 610 ysv_ext1 = ysb + 1; 611 zsv_ext1 = zsb - 1; 612 } 613 } 614 } else { //One point on (0,0,0) side, one point on (1,1,1) side 615 byte c1, c2; 616 if (aIsFurtherSide) { 617 c1 = aPoint; 618 c2 = bPoint; 619 } else { 620 c1 = bPoint; 621 c2 = aPoint; 622 } 623 624 //One contribution is a permutation of (1,1,-1) 625 if ((c1 & 0x01) == 0) { 626 dx_ext0 = dx0 + 1 - SQUISH_CONSTANT_3D; 627 dy_ext0 = dy0 - 1 - SQUISH_CONSTANT_3D; 628 dz_ext0 = dz0 - 1 - SQUISH_CONSTANT_3D; 629 xsv_ext0 = xsb - 1; 630 ysv_ext0 = ysb + 1; 631 zsv_ext0 = zsb + 1; 632 } else if ((c1 & 0x02) == 0) { 633 dx_ext0 = dx0 - 1 - SQUISH_CONSTANT_3D; 634 dy_ext0 = dy0 + 1 - SQUISH_CONSTANT_3D; 635 dz_ext0 = dz0 - 1 - SQUISH_CONSTANT_3D; 636 xsv_ext0 = xsb + 1; 637 ysv_ext0 = ysb - 1; 638 zsv_ext0 = zsb + 1; 639 } else { 640 dx_ext0 = dx0 - 1 - SQUISH_CONSTANT_3D; 641 dy_ext0 = dy0 - 1 - SQUISH_CONSTANT_3D; 642 dz_ext0 = dz0 + 1 - SQUISH_CONSTANT_3D; 643 xsv_ext0 = xsb + 1; 644 ysv_ext0 = ysb + 1; 645 zsv_ext0 = zsb - 1; 646 } 647 648 //One contribution is a permutation of (0,0,2) 649 dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D; 650 dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D; 651 dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D; 652 xsv_ext1 = xsb; 653 ysv_ext1 = ysb; 654 zsv_ext1 = zsb; 655 if ((c2 & 0x01) != 0) { 656 dx_ext1 -= 2; 657 xsv_ext1 += 2; 658 } else if ((c2 & 0x02) != 0) { 659 dy_ext1 -= 2; 660 ysv_ext1 += 2; 661 } else { 662 dz_ext1 -= 2; 663 zsv_ext1 += 2; 664 } 665 } 666 667 //Contribution (1,0,0) 668 double dx1 = dx0 - 1 - SQUISH_CONSTANT_3D; 669 double dy1 = dy0 - 0 - SQUISH_CONSTANT_3D; 670 double dz1 = dz0 - 0 - SQUISH_CONSTANT_3D; 671 double attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1; 672 if (attn1 > 0) { 673 attn1 *= attn1; 674 value += attn1 * attn1 * extrapolate(xsb + 1, ysb + 0, zsb + 0, dx1, dy1, dz1); 675 } 676 677 //Contribution (0,1,0) 678 double dx2 = dx0 - 0 - SQUISH_CONSTANT_3D; 679 double dy2 = dy0 - 1 - SQUISH_CONSTANT_3D; 680 double dz2 = dz1; 681 double attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2; 682 if (attn2 > 0) { 683 attn2 *= attn2; 684 value += attn2 * attn2 * extrapolate(xsb + 0, ysb + 1, zsb + 0, dx2, dy2, dz2); 685 } 686 687 //Contribution (0,0,1) 688 double dx3 = dx2; 689 double dy3 = dy1; 690 double dz3 = dz0 - 1 - SQUISH_CONSTANT_3D; 691 double attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3; 692 if (attn3 > 0) { 693 attn3 *= attn3; 694 value += attn3 * attn3 * extrapolate(xsb + 0, ysb + 0, zsb + 1, dx3, dy3, dz3); 695 } 696 697 //Contribution (1,1,0) 698 double dx4 = dx0 - 1 - 2 * SQUISH_CONSTANT_3D; 699 double dy4 = dy0 - 1 - 2 * SQUISH_CONSTANT_3D; 700 double dz4 = dz0 - 0 - 2 * SQUISH_CONSTANT_3D; 701 double attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4; 702 if (attn4 > 0) { 703 attn4 *= attn4; 704 value += attn4 * attn4 * extrapolate(xsb + 1, ysb + 1, zsb + 0, dx4, dy4, dz4); 705 } 706 707 //Contribution (1,0,1) 708 double dx5 = dx4; 709 double dy5 = dy0 - 0 - 2 * SQUISH_CONSTANT_3D; 710 double dz5 = dz0 - 1 - 2 * SQUISH_CONSTANT_3D; 711 double attn5 = 2 - dx5 * dx5 - dy5 * dy5 - dz5 * dz5; 712 if (attn5 > 0) { 713 attn5 *= attn5; 714 value += attn5 * attn5 * extrapolate(xsb + 1, ysb + 0, zsb + 1, dx5, dy5, dz5); 715 } 716 717 //Contribution (0,1,1) 718 double dx6 = dx0 - 0 - 2 * SQUISH_CONSTANT_3D; 719 double dy6 = dy4; 720 double dz6 = dz5; 721 double attn6 = 2 - dx6 * dx6 - dy6 * dy6 - dz6 * dz6; 722 if (attn6 > 0) { 723 attn6 *= attn6; 724 value += attn6 * attn6 * extrapolate(xsb + 0, ysb + 1, zsb + 1, dx6, dy6, dz6); 725 } 726 } 727 728 //First extra vertex 729 double attn_ext0 = 2 - dx_ext0 * dx_ext0 - dy_ext0 * dy_ext0 - dz_ext0 * dz_ext0; 730 if (attn_ext0 > 0) 731 { 732 attn_ext0 *= attn_ext0; 733 value += attn_ext0 * attn_ext0 * extrapolate(xsv_ext0, ysv_ext0, zsv_ext0, dx_ext0, dy_ext0, dz_ext0); 734 } 735 736 //Second extra vertex 737 double attn_ext1 = 2 - dx_ext1 * dx_ext1 - dy_ext1 * dy_ext1 - dz_ext1 * dz_ext1; 738 if (attn_ext1 > 0) 739 { 740 attn_ext1 *= attn_ext1; 741 value += attn_ext1 * attn_ext1 * extrapolate(xsv_ext1, ysv_ext1, zsv_ext1, dx_ext1, dy_ext1, dz_ext1); 742 } 743 744 return value / NORM_CONSTANT_3D; 745 } 746 747 //4D OpenSimplex Noise. 748 public double noise(double x, double y, double z, double w) { 749 750 //Place input coordinates on simplectic honeycomb. 751 double stretchOffset = (x + y + z + w) * STRETCH_CONSTANT_4D; 752 double xs = x + stretchOffset; 753 double ys = y + stretchOffset; 754 double zs = z + stretchOffset; 755 double ws = w + stretchOffset; 756 757 //Floor to get simplectic honeycomb coordinates of rhombo-hypercube super-cell origin. 758 int xsb = fastFloor(xs); 759 int ysb = fastFloor(ys); 760 int zsb = fastFloor(zs); 761 int wsb = fastFloor(ws); 762 763 //Skew out to get actual coordinates of stretched rhombo-hypercube origin. We'll need these later. 764 double squishOffset = (xsb + ysb + zsb + wsb) * SQUISH_CONSTANT_4D; 765 double xb = xsb + squishOffset; 766 double yb = ysb + squishOffset; 767 double zb = zsb + squishOffset; 768 double wb = wsb + squishOffset; 769 770 //Compute simplectic honeycomb coordinates relative to rhombo-hypercube origin. 771 double xins = xs - xsb; 772 double yins = ys - ysb; 773 double zins = zs - zsb; 774 double wins = ws - wsb; 775 776 //Sum those together to get a value that determines which region we're in. 777 double inSum = xins + yins + zins + wins; 778 779 //Positions relative to origin point. 780 double dx0 = x - xb; 781 double dy0 = y - yb; 782 double dz0 = z - zb; 783 double dw0 = w - wb; 784 785 //We'll be defining these inside the next block and using them afterwards. 786 double dx_ext0, dy_ext0, dz_ext0, dw_ext0; 787 double dx_ext1, dy_ext1, dz_ext1, dw_ext1; 788 double dx_ext2, dy_ext2, dz_ext2, dw_ext2; 789 int xsv_ext0, ysv_ext0, zsv_ext0, wsv_ext0; 790 int xsv_ext1, ysv_ext1, zsv_ext1, wsv_ext1; 791 int xsv_ext2, ysv_ext2, zsv_ext2, wsv_ext2; 792 793 double value = 0; 794 if (inSum <= 1) { //We're inside the pentachoron (4-Simplex) at (0,0,0,0) 795 796 //Determine which two of (0,0,0,1), (0,0,1,0), (0,1,0,0), (1,0,0,0) are closest. 797 byte aPoint = 0x01; 798 double aScore = xins; 799 byte bPoint = 0x02; 800 double bScore = yins; 801 if (aScore >= bScore && zins > bScore) { 802 bScore = zins; 803 bPoint = 0x04; 804 } else if (aScore < bScore && zins > aScore) { 805 aScore = zins; 806 aPoint = 0x04; 807 } 808 if (aScore >= bScore && wins > bScore) { 809 bScore = wins; 810 bPoint = 0x08; 811 } else if (aScore < bScore && wins > aScore) { 812 aScore = wins; 813 aPoint = 0x08; 814 } 815 816 //Now we determine the three lattice points not part of the pentachoron that may contribute. 817 //This depends on the closest two pentachoron vertices, including (0,0,0,0) 818 double uins = 1 - inSum; 819 if (uins > aScore || uins > bScore) { //(0,0,0,0) is one of the closest two pentachoron vertices. 820 byte c = (bScore > aScore ? bPoint : aPoint); //Our other closest vertex is the closest out of a and b. 821 if ((c & 0x01) == 0) { 822 xsv_ext0 = xsb - 1; 823 xsv_ext1 = xsv_ext2 = xsb; 824 dx_ext0 = dx0 + 1; 825 dx_ext1 = dx_ext2 = dx0; 826 } else { 827 xsv_ext0 = xsv_ext1 = xsv_ext2 = xsb + 1; 828 dx_ext0 = dx_ext1 = dx_ext2 = dx0 - 1; 829 } 830 831 if ((c & 0x02) == 0) { 832 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb; 833 dy_ext0 = dy_ext1 = dy_ext2 = dy0; 834 if ((c & 0x01) == 0x01) { 835 ysv_ext0 -= 1; 836 dy_ext0 += 1; 837 } else { 838 ysv_ext1 -= 1; 839 dy_ext1 += 1; 840 } 841 } else { 842 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb + 1; 843 dy_ext0 = dy_ext1 = dy_ext2 = dy0 - 1; 844 } 845 846 if ((c & 0x04) == 0) { 847 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb; 848 dz_ext0 = dz_ext1 = dz_ext2 = dz0; 849 if ((c & 0x03) != 0) { 850 if ((c & 0x03) == 0x03) { 851 zsv_ext0 -= 1; 852 dz_ext0 += 1; 853 } else { 854 zsv_ext1 -= 1; 855 dz_ext1 += 1; 856 } 857 } else { 858 zsv_ext2 -= 1; 859 dz_ext2 += 1; 860 } 861 } else { 862 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb + 1; 863 dz_ext0 = dz_ext1 = dz_ext2 = dz0 - 1; 864 } 865 866 if ((c & 0x08) == 0) { 867 wsv_ext0 = wsv_ext1 = wsb; 868 wsv_ext2 = wsb - 1; 869 dw_ext0 = dw_ext1 = dw0; 870 dw_ext2 = dw0 + 1; 871 } else { 872 wsv_ext0 = wsv_ext1 = wsv_ext2 = wsb + 1; 873 dw_ext0 = dw_ext1 = dw_ext2 = dw0 - 1; 874 } 875 } else { //(0,0,0,0) is not one of the closest two pentachoron vertices. 876 byte c = cast(byte)(aPoint | bPoint); //Our three extra vertices are determined by the closest two. 877 878 if ((c & 0x01) == 0) { 879 xsv_ext0 = xsv_ext2 = xsb; 880 xsv_ext1 = xsb - 1; 881 dx_ext0 = dx0 - 2 * SQUISH_CONSTANT_4D; 882 dx_ext1 = dx0 + 1 - SQUISH_CONSTANT_4D; 883 dx_ext2 = dx0 - SQUISH_CONSTANT_4D; 884 } else { 885 xsv_ext0 = xsv_ext1 = xsv_ext2 = xsb + 1; 886 dx_ext0 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 887 dx_ext1 = dx_ext2 = dx0 - 1 - SQUISH_CONSTANT_4D; 888 } 889 890 if ((c & 0x02) == 0) { 891 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb; 892 dy_ext0 = dy0 - 2 * SQUISH_CONSTANT_4D; 893 dy_ext1 = dy_ext2 = dy0 - SQUISH_CONSTANT_4D; 894 if ((c & 0x01) == 0x01) { 895 ysv_ext1 -= 1; 896 dy_ext1 += 1; 897 } else { 898 ysv_ext2 -= 1; 899 dy_ext2 += 1; 900 } 901 } else { 902 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb + 1; 903 dy_ext0 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 904 dy_ext1 = dy_ext2 = dy0 - 1 - SQUISH_CONSTANT_4D; 905 } 906 907 if ((c & 0x04) == 0) { 908 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb; 909 dz_ext0 = dz0 - 2 * SQUISH_CONSTANT_4D; 910 dz_ext1 = dz_ext2 = dz0 - SQUISH_CONSTANT_4D; 911 if ((c & 0x03) == 0x03) { 912 zsv_ext1 -= 1; 913 dz_ext1 += 1; 914 } else { 915 zsv_ext2 -= 1; 916 dz_ext2 += 1; 917 } 918 } else { 919 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb + 1; 920 dz_ext0 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 921 dz_ext1 = dz_ext2 = dz0 - 1 - SQUISH_CONSTANT_4D; 922 } 923 924 if ((c & 0x08) == 0) { 925 wsv_ext0 = wsv_ext1 = wsb; 926 wsv_ext2 = wsb - 1; 927 dw_ext0 = dw0 - 2 * SQUISH_CONSTANT_4D; 928 dw_ext1 = dw0 - SQUISH_CONSTANT_4D; 929 dw_ext2 = dw0 + 1 - SQUISH_CONSTANT_4D; 930 } else { 931 wsv_ext0 = wsv_ext1 = wsv_ext2 = wsb + 1; 932 dw_ext0 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 933 dw_ext1 = dw_ext2 = dw0 - 1 - SQUISH_CONSTANT_4D; 934 } 935 } 936 937 //Contribution (0,0,0,0) 938 double attn0 = 2 - dx0 * dx0 - dy0 * dy0 - dz0 * dz0 - dw0 * dw0; 939 if (attn0 > 0) { 940 attn0 *= attn0; 941 value += attn0 * attn0 * extrapolate(xsb + 0, ysb + 0, zsb + 0, wsb + 0, dx0, dy0, dz0, dw0); 942 } 943 944 //Contribution (1,0,0,0) 945 double dx1 = dx0 - 1 - SQUISH_CONSTANT_4D; 946 double dy1 = dy0 - 0 - SQUISH_CONSTANT_4D; 947 double dz1 = dz0 - 0 - SQUISH_CONSTANT_4D; 948 double dw1 = dw0 - 0 - SQUISH_CONSTANT_4D; 949 double attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1 - dw1 * dw1; 950 if (attn1 > 0) { 951 attn1 *= attn1; 952 value += attn1 * attn1 * extrapolate(xsb + 1, ysb + 0, zsb + 0, wsb + 0, dx1, dy1, dz1, dw1); 953 } 954 955 //Contribution (0,1,0,0) 956 double dx2 = dx0 - 0 - SQUISH_CONSTANT_4D; 957 double dy2 = dy0 - 1 - SQUISH_CONSTANT_4D; 958 double dz2 = dz1; 959 double dw2 = dw1; 960 double attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2 - dw2 * dw2; 961 if (attn2 > 0) { 962 attn2 *= attn2; 963 value += attn2 * attn2 * extrapolate(xsb + 0, ysb + 1, zsb + 0, wsb + 0, dx2, dy2, dz2, dw2); 964 } 965 966 //Contribution (0,0,1,0) 967 double dx3 = dx2; 968 double dy3 = dy1; 969 double dz3 = dz0 - 1 - SQUISH_CONSTANT_4D; 970 double dw3 = dw1; 971 double attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3 - dw3 * dw3; 972 if (attn3 > 0) { 973 attn3 *= attn3; 974 value += attn3 * attn3 * extrapolate(xsb + 0, ysb + 0, zsb + 1, wsb + 0, dx3, dy3, dz3, dw3); 975 } 976 977 //Contribution (0,0,0,1) 978 double dx4 = dx2; 979 double dy4 = dy1; 980 double dz4 = dz1; 981 double dw4 = dw0 - 1 - SQUISH_CONSTANT_4D; 982 double attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4 - dw4 * dw4; 983 if (attn4 > 0) { 984 attn4 *= attn4; 985 value += attn4 * attn4 * extrapolate(xsb + 0, ysb + 0, zsb + 0, wsb + 1, dx4, dy4, dz4, dw4); 986 } 987 } else if (inSum >= 3) { //We're inside the pentachoron (4-Simplex) at (1,1,1,1) 988 //Determine which two of (1,1,1,0), (1,1,0,1), (1,0,1,1), (0,1,1,1) are closest. 989 byte aPoint = 0x0E; 990 double aScore = xins; 991 byte bPoint = 0x0D; 992 double bScore = yins; 993 if (aScore <= bScore && zins < bScore) { 994 bScore = zins; 995 bPoint = 0x0B; 996 } else if (aScore > bScore && zins < aScore) { 997 aScore = zins; 998 aPoint = 0x0B; 999 } 1000 if (aScore <= bScore && wins < bScore) { 1001 bScore = wins; 1002 bPoint = 0x07; 1003 } else if (aScore > bScore && wins < aScore) { 1004 aScore = wins; 1005 aPoint = 0x07; 1006 } 1007 1008 //Now we determine the three lattice points not part of the pentachoron that may contribute. 1009 //This depends on the closest two pentachoron vertices, including (0,0,0,0) 1010 double uins = 4 - inSum; 1011 if (uins < aScore || uins < bScore) { //(1,1,1,1) is one of the closest two pentachoron vertices. 1012 byte c = (bScore < aScore ? bPoint : aPoint); //Our other closest vertex is the closest out of a and b. 1013 1014 if ((c & 0x01) != 0) { 1015 xsv_ext0 = xsb + 2; 1016 xsv_ext1 = xsv_ext2 = xsb + 1; 1017 dx_ext0 = dx0 - 2 - 4 * SQUISH_CONSTANT_4D; 1018 dx_ext1 = dx_ext2 = dx0 - 1 - 4 * SQUISH_CONSTANT_4D; 1019 } else { 1020 xsv_ext0 = xsv_ext1 = xsv_ext2 = xsb; 1021 dx_ext0 = dx_ext1 = dx_ext2 = dx0 - 4 * SQUISH_CONSTANT_4D; 1022 } 1023 1024 if ((c & 0x02) != 0) { 1025 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb + 1; 1026 dy_ext0 = dy_ext1 = dy_ext2 = dy0 - 1 - 4 * SQUISH_CONSTANT_4D; 1027 if ((c & 0x01) != 0) { 1028 ysv_ext1 += 1; 1029 dy_ext1 -= 1; 1030 } else { 1031 ysv_ext0 += 1; 1032 dy_ext0 -= 1; 1033 } 1034 } else { 1035 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb; 1036 dy_ext0 = dy_ext1 = dy_ext2 = dy0 - 4 * SQUISH_CONSTANT_4D; 1037 } 1038 1039 if ((c & 0x04) != 0) { 1040 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb + 1; 1041 dz_ext0 = dz_ext1 = dz_ext2 = dz0 - 1 - 4 * SQUISH_CONSTANT_4D; 1042 if ((c & 0x03) != 0x03) { 1043 if ((c & 0x03) == 0) { 1044 zsv_ext0 += 1; 1045 dz_ext0 -= 1; 1046 } else { 1047 zsv_ext1 += 1; 1048 dz_ext1 -= 1; 1049 } 1050 } else { 1051 zsv_ext2 += 1; 1052 dz_ext2 -= 1; 1053 } 1054 } else { 1055 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb; 1056 dz_ext0 = dz_ext1 = dz_ext2 = dz0 - 4 * SQUISH_CONSTANT_4D; 1057 } 1058 1059 if ((c & 0x08) != 0) { 1060 wsv_ext0 = wsv_ext1 = wsb + 1; 1061 wsv_ext2 = wsb + 2; 1062 dw_ext0 = dw_ext1 = dw0 - 1 - 4 * SQUISH_CONSTANT_4D; 1063 dw_ext2 = dw0 - 2 - 4 * SQUISH_CONSTANT_4D; 1064 } else { 1065 wsv_ext0 = wsv_ext1 = wsv_ext2 = wsb; 1066 dw_ext0 = dw_ext1 = dw_ext2 = dw0 - 4 * SQUISH_CONSTANT_4D; 1067 } 1068 } else { //(1,1,1,1) is not one of the closest two pentachoron vertices. 1069 byte c = cast(byte)(aPoint & bPoint); //Our three extra vertices are determined by the closest two. 1070 1071 if ((c & 0x01) != 0) { 1072 xsv_ext0 = xsv_ext2 = xsb + 1; 1073 xsv_ext1 = xsb + 2; 1074 dx_ext0 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1075 dx_ext1 = dx0 - 2 - 3 * SQUISH_CONSTANT_4D; 1076 dx_ext2 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; 1077 } else { 1078 xsv_ext0 = xsv_ext1 = xsv_ext2 = xsb; 1079 dx_ext0 = dx0 - 2 * SQUISH_CONSTANT_4D; 1080 dx_ext1 = dx_ext2 = dx0 - 3 * SQUISH_CONSTANT_4D; 1081 } 1082 1083 if ((c & 0x02) != 0) { 1084 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb + 1; 1085 dy_ext0 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1086 dy_ext1 = dy_ext2 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; 1087 if ((c & 0x01) != 0) { 1088 ysv_ext2 += 1; 1089 dy_ext2 -= 1; 1090 } else { 1091 ysv_ext1 += 1; 1092 dy_ext1 -= 1; 1093 } 1094 } else { 1095 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb; 1096 dy_ext0 = dy0 - 2 * SQUISH_CONSTANT_4D; 1097 dy_ext1 = dy_ext2 = dy0 - 3 * SQUISH_CONSTANT_4D; 1098 } 1099 1100 if ((c & 0x04) != 0) { 1101 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb + 1; 1102 dz_ext0 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1103 dz_ext1 = dz_ext2 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; 1104 if ((c & 0x03) != 0) { 1105 zsv_ext2 += 1; 1106 dz_ext2 -= 1; 1107 } else { 1108 zsv_ext1 += 1; 1109 dz_ext1 -= 1; 1110 } 1111 } else { 1112 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb; 1113 dz_ext0 = dz0 - 2 * SQUISH_CONSTANT_4D; 1114 dz_ext1 = dz_ext2 = dz0 - 3 * SQUISH_CONSTANT_4D; 1115 } 1116 1117 if ((c & 0x08) != 0) { 1118 wsv_ext0 = wsv_ext1 = wsb + 1; 1119 wsv_ext2 = wsb + 2; 1120 dw_ext0 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1121 dw_ext1 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; 1122 dw_ext2 = dw0 - 2 - 3 * SQUISH_CONSTANT_4D; 1123 } else { 1124 wsv_ext0 = wsv_ext1 = wsv_ext2 = wsb; 1125 dw_ext0 = dw0 - 2 * SQUISH_CONSTANT_4D; 1126 dw_ext1 = dw_ext2 = dw0 - 3 * SQUISH_CONSTANT_4D; 1127 } 1128 } 1129 1130 //Contribution (1,1,1,0) 1131 double dx4 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; 1132 double dy4 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; 1133 double dz4 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; 1134 double dw4 = dw0 - 3 * SQUISH_CONSTANT_4D; 1135 double attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4 - dw4 * dw4; 1136 if (attn4 > 0) { 1137 attn4 *= attn4; 1138 value += attn4 * attn4 * extrapolate(xsb + 1, ysb + 1, zsb + 1, wsb + 0, dx4, dy4, dz4, dw4); 1139 } 1140 1141 //Contribution (1,1,0,1) 1142 double dx3 = dx4; 1143 double dy3 = dy4; 1144 double dz3 = dz0 - 3 * SQUISH_CONSTANT_4D; 1145 double dw3 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; 1146 double attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3 - dw3 * dw3; 1147 if (attn3 > 0) { 1148 attn3 *= attn3; 1149 value += attn3 * attn3 * extrapolate(xsb + 1, ysb + 1, zsb + 0, wsb + 1, dx3, dy3, dz3, dw3); 1150 } 1151 1152 //Contribution (1,0,1,1) 1153 double dx2 = dx4; 1154 double dy2 = dy0 - 3 * SQUISH_CONSTANT_4D; 1155 double dz2 = dz4; 1156 double dw2 = dw3; 1157 double attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2 - dw2 * dw2; 1158 if (attn2 > 0) { 1159 attn2 *= attn2; 1160 value += attn2 * attn2 * extrapolate(xsb + 1, ysb + 0, zsb + 1, wsb + 1, dx2, dy2, dz2, dw2); 1161 } 1162 1163 //Contribution (0,1,1,1) 1164 double dx1 = dx0 - 3 * SQUISH_CONSTANT_4D; 1165 double dz1 = dz4; 1166 double dy1 = dy4; 1167 double dw1 = dw3; 1168 double attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1 - dw1 * dw1; 1169 if (attn1 > 0) { 1170 attn1 *= attn1; 1171 value += attn1 * attn1 * extrapolate(xsb + 0, ysb + 1, zsb + 1, wsb + 1, dx1, dy1, dz1, dw1); 1172 } 1173 1174 //Contribution (1,1,1,1) 1175 dx0 = dx0 - 1 - 4 * SQUISH_CONSTANT_4D; 1176 dy0 = dy0 - 1 - 4 * SQUISH_CONSTANT_4D; 1177 dz0 = dz0 - 1 - 4 * SQUISH_CONSTANT_4D; 1178 dw0 = dw0 - 1 - 4 * SQUISH_CONSTANT_4D; 1179 double attn0 = 2 - dx0 * dx0 - dy0 * dy0 - dz0 * dz0 - dw0 * dw0; 1180 if (attn0 > 0) { 1181 attn0 *= attn0; 1182 value += attn0 * attn0 * extrapolate(xsb + 1, ysb + 1, zsb + 1, wsb + 1, dx0, dy0, dz0, dw0); 1183 } 1184 } else if (inSum <= 2) { //We're inside the first dispentachoron (Rectified 4-Simplex) 1185 double aScore; 1186 byte aPoint; 1187 bool aIsBiggerSide = true; 1188 double bScore; 1189 byte bPoint; 1190 bool bIsBiggerSide = true; 1191 1192 //Decide between (1,1,0,0) and (0,0,1,1) 1193 if (xins + yins > zins + wins) { 1194 aScore = xins + yins; 1195 aPoint = 0x03; 1196 } else { 1197 aScore = zins + wins; 1198 aPoint = 0x0C; 1199 } 1200 1201 //Decide between (1,0,1,0) and (0,1,0,1) 1202 if (xins + zins > yins + wins) { 1203 bScore = xins + zins; 1204 bPoint = 0x05; 1205 } else { 1206 bScore = yins + wins; 1207 bPoint = 0x0A; 1208 } 1209 1210 //Closer between (1,0,0,1) and (0,1,1,0) will replace the further of a and b, if closer. 1211 if (xins + wins > yins + zins) { 1212 double score = xins + wins; 1213 if (aScore >= bScore && score > bScore) { 1214 bScore = score; 1215 bPoint = 0x09; 1216 } else if (aScore < bScore && score > aScore) { 1217 aScore = score; 1218 aPoint = 0x09; 1219 } 1220 } else { 1221 double score = yins + zins; 1222 if (aScore >= bScore && score > bScore) { 1223 bScore = score; 1224 bPoint = 0x06; 1225 } else if (aScore < bScore && score > aScore) { 1226 aScore = score; 1227 aPoint = 0x06; 1228 } 1229 } 1230 1231 //Decide if (1,0,0,0) is closer. 1232 double p1 = 2 - inSum + xins; 1233 if (aScore >= bScore && p1 > bScore) { 1234 bScore = p1; 1235 bPoint = 0x01; 1236 bIsBiggerSide = false; 1237 } else if (aScore < bScore && p1 > aScore) { 1238 aScore = p1; 1239 aPoint = 0x01; 1240 aIsBiggerSide = false; 1241 } 1242 1243 //Decide if (0,1,0,0) is closer. 1244 double p2 = 2 - inSum + yins; 1245 if (aScore >= bScore && p2 > bScore) { 1246 bScore = p2; 1247 bPoint = 0x02; 1248 bIsBiggerSide = false; 1249 } else if (aScore < bScore && p2 > aScore) { 1250 aScore = p2; 1251 aPoint = 0x02; 1252 aIsBiggerSide = false; 1253 } 1254 1255 //Decide if (0,0,1,0) is closer. 1256 double p3 = 2 - inSum + zins; 1257 if (aScore >= bScore && p3 > bScore) { 1258 bScore = p3; 1259 bPoint = 0x04; 1260 bIsBiggerSide = false; 1261 } else if (aScore < bScore && p3 > aScore) { 1262 aScore = p3; 1263 aPoint = 0x04; 1264 aIsBiggerSide = false; 1265 } 1266 1267 //Decide if (0,0,0,1) is closer. 1268 double p4 = 2 - inSum + wins; 1269 if (aScore >= bScore && p4 > bScore) { 1270 bScore = p4; 1271 bPoint = 0x08; 1272 bIsBiggerSide = false; 1273 } else if (aScore < bScore && p4 > aScore) { 1274 aScore = p4; 1275 aPoint = 0x08; 1276 aIsBiggerSide = false; 1277 } 1278 1279 //Where each of the two closest points are determines how the extra three vertices are calculated. 1280 if (aIsBiggerSide == bIsBiggerSide) { 1281 if (aIsBiggerSide) { //Both closest points on the bigger side 1282 byte c1 = cast(byte)(aPoint | bPoint); 1283 byte c2 = cast(byte)(aPoint & bPoint); 1284 if ((c1 & 0x01) == 0) { 1285 xsv_ext0 = xsb; 1286 xsv_ext1 = xsb - 1; 1287 dx_ext0 = dx0 - 3 * SQUISH_CONSTANT_4D; 1288 dx_ext1 = dx0 + 1 - 2 * SQUISH_CONSTANT_4D; 1289 } else { 1290 xsv_ext0 = xsv_ext1 = xsb + 1; 1291 dx_ext0 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; 1292 dx_ext1 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1293 } 1294 1295 if ((c1 & 0x02) == 0) { 1296 ysv_ext0 = ysb; 1297 ysv_ext1 = ysb - 1; 1298 dy_ext0 = dy0 - 3 * SQUISH_CONSTANT_4D; 1299 dy_ext1 = dy0 + 1 - 2 * SQUISH_CONSTANT_4D; 1300 } else { 1301 ysv_ext0 = ysv_ext1 = ysb + 1; 1302 dy_ext0 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; 1303 dy_ext1 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1304 } 1305 1306 if ((c1 & 0x04) == 0) { 1307 zsv_ext0 = zsb; 1308 zsv_ext1 = zsb - 1; 1309 dz_ext0 = dz0 - 3 * SQUISH_CONSTANT_4D; 1310 dz_ext1 = dz0 + 1 - 2 * SQUISH_CONSTANT_4D; 1311 } else { 1312 zsv_ext0 = zsv_ext1 = zsb + 1; 1313 dz_ext0 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; 1314 dz_ext1 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1315 } 1316 1317 if ((c1 & 0x08) == 0) { 1318 wsv_ext0 = wsb; 1319 wsv_ext1 = wsb - 1; 1320 dw_ext0 = dw0 - 3 * SQUISH_CONSTANT_4D; 1321 dw_ext1 = dw0 + 1 - 2 * SQUISH_CONSTANT_4D; 1322 } else { 1323 wsv_ext0 = wsv_ext1 = wsb + 1; 1324 dw_ext0 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; 1325 dw_ext1 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1326 } 1327 1328 //One combination is a permutation of (0,0,0,2) based on c2 1329 xsv_ext2 = xsb; 1330 ysv_ext2 = ysb; 1331 zsv_ext2 = zsb; 1332 wsv_ext2 = wsb; 1333 dx_ext2 = dx0 - 2 * SQUISH_CONSTANT_4D; 1334 dy_ext2 = dy0 - 2 * SQUISH_CONSTANT_4D; 1335 dz_ext2 = dz0 - 2 * SQUISH_CONSTANT_4D; 1336 dw_ext2 = dw0 - 2 * SQUISH_CONSTANT_4D; 1337 if ((c2 & 0x01) != 0) { 1338 xsv_ext2 += 2; 1339 dx_ext2 -= 2; 1340 } else if ((c2 & 0x02) != 0) { 1341 ysv_ext2 += 2; 1342 dy_ext2 -= 2; 1343 } else if ((c2 & 0x04) != 0) { 1344 zsv_ext2 += 2; 1345 dz_ext2 -= 2; 1346 } else { 1347 wsv_ext2 += 2; 1348 dw_ext2 -= 2; 1349 } 1350 1351 } else { //Both closest points on the smaller side 1352 //One of the two extra points is (0,0,0,0) 1353 xsv_ext2 = xsb; 1354 ysv_ext2 = ysb; 1355 zsv_ext2 = zsb; 1356 wsv_ext2 = wsb; 1357 dx_ext2 = dx0; 1358 dy_ext2 = dy0; 1359 dz_ext2 = dz0; 1360 dw_ext2 = dw0; 1361 1362 //Other two points are based on the omitted axes. 1363 byte c = cast(byte)(aPoint | bPoint); 1364 1365 if ((c & 0x01) == 0) { 1366 xsv_ext0 = xsb - 1; 1367 xsv_ext1 = xsb; 1368 dx_ext0 = dx0 + 1 - SQUISH_CONSTANT_4D; 1369 dx_ext1 = dx0 - SQUISH_CONSTANT_4D; 1370 } else { 1371 xsv_ext0 = xsv_ext1 = xsb + 1; 1372 dx_ext0 = dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_4D; 1373 } 1374 1375 if ((c & 0x02) == 0) { 1376 ysv_ext0 = ysv_ext1 = ysb; 1377 dy_ext0 = dy_ext1 = dy0 - SQUISH_CONSTANT_4D; 1378 if ((c & 0x01) == 0x01) 1379 { 1380 ysv_ext0 -= 1; 1381 dy_ext0 += 1; 1382 } else { 1383 ysv_ext1 -= 1; 1384 dy_ext1 += 1; 1385 } 1386 } else { 1387 ysv_ext0 = ysv_ext1 = ysb + 1; 1388 dy_ext0 = dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_4D; 1389 } 1390 1391 if ((c & 0x04) == 0) { 1392 zsv_ext0 = zsv_ext1 = zsb; 1393 dz_ext0 = dz_ext1 = dz0 - SQUISH_CONSTANT_4D; 1394 if ((c & 0x03) == 0x03) 1395 { 1396 zsv_ext0 -= 1; 1397 dz_ext0 += 1; 1398 } else { 1399 zsv_ext1 -= 1; 1400 dz_ext1 += 1; 1401 } 1402 } else { 1403 zsv_ext0 = zsv_ext1 = zsb + 1; 1404 dz_ext0 = dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_4D; 1405 } 1406 1407 if ((c & 0x08) == 0) 1408 { 1409 wsv_ext0 = wsb; 1410 wsv_ext1 = wsb - 1; 1411 dw_ext0 = dw0 - SQUISH_CONSTANT_4D; 1412 dw_ext1 = dw0 + 1 - SQUISH_CONSTANT_4D; 1413 } else { 1414 wsv_ext0 = wsv_ext1 = wsb + 1; 1415 dw_ext0 = dw_ext1 = dw0 - 1 - SQUISH_CONSTANT_4D; 1416 } 1417 1418 } 1419 } else { //One point on each "side" 1420 byte c1, c2; 1421 if (aIsBiggerSide) { 1422 c1 = aPoint; 1423 c2 = bPoint; 1424 } else { 1425 c1 = bPoint; 1426 c2 = aPoint; 1427 } 1428 1429 //Two contributions are the bigger-sided point with each 0 replaced with -1. 1430 if ((c1 & 0x01) == 0) { 1431 xsv_ext0 = xsb - 1; 1432 xsv_ext1 = xsb; 1433 dx_ext0 = dx0 + 1 - SQUISH_CONSTANT_4D; 1434 dx_ext1 = dx0 - SQUISH_CONSTANT_4D; 1435 } else { 1436 xsv_ext0 = xsv_ext1 = xsb + 1; 1437 dx_ext0 = dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_4D; 1438 } 1439 1440 if ((c1 & 0x02) == 0) { 1441 ysv_ext0 = ysv_ext1 = ysb; 1442 dy_ext0 = dy_ext1 = dy0 - SQUISH_CONSTANT_4D; 1443 if ((c1 & 0x01) == 0x01) { 1444 ysv_ext0 -= 1; 1445 dy_ext0 += 1; 1446 } else { 1447 ysv_ext1 -= 1; 1448 dy_ext1 += 1; 1449 } 1450 } else { 1451 ysv_ext0 = ysv_ext1 = ysb + 1; 1452 dy_ext0 = dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_4D; 1453 } 1454 1455 if ((c1 & 0x04) == 0) { 1456 zsv_ext0 = zsv_ext1 = zsb; 1457 dz_ext0 = dz_ext1 = dz0 - SQUISH_CONSTANT_4D; 1458 if ((c1 & 0x03) == 0x03) { 1459 zsv_ext0 -= 1; 1460 dz_ext0 += 1; 1461 } else { 1462 zsv_ext1 -= 1; 1463 dz_ext1 += 1; 1464 } 1465 } else { 1466 zsv_ext0 = zsv_ext1 = zsb + 1; 1467 dz_ext0 = dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_4D; 1468 } 1469 1470 if ((c1 & 0x08) == 0) { 1471 wsv_ext0 = wsb; 1472 wsv_ext1 = wsb - 1; 1473 dw_ext0 = dw0 - SQUISH_CONSTANT_4D; 1474 dw_ext1 = dw0 + 1 - SQUISH_CONSTANT_4D; 1475 } else { 1476 wsv_ext0 = wsv_ext1 = wsb + 1; 1477 dw_ext0 = dw_ext1 = dw0 - 1 - SQUISH_CONSTANT_4D; 1478 } 1479 1480 //One contribution is a permutation of (0,0,0,2) based on the smaller-sided point 1481 xsv_ext2 = xsb; 1482 ysv_ext2 = ysb; 1483 zsv_ext2 = zsb; 1484 wsv_ext2 = wsb; 1485 dx_ext2 = dx0 - 2 * SQUISH_CONSTANT_4D; 1486 dy_ext2 = dy0 - 2 * SQUISH_CONSTANT_4D; 1487 dz_ext2 = dz0 - 2 * SQUISH_CONSTANT_4D; 1488 dw_ext2 = dw0 - 2 * SQUISH_CONSTANT_4D; 1489 if ((c2 & 0x01) != 0) { 1490 xsv_ext2 += 2; 1491 dx_ext2 -= 2; 1492 } else if ((c2 & 0x02) != 0) { 1493 ysv_ext2 += 2; 1494 dy_ext2 -= 2; 1495 } else if ((c2 & 0x04) != 0) { 1496 zsv_ext2 += 2; 1497 dz_ext2 -= 2; 1498 } else { 1499 wsv_ext2 += 2; 1500 dw_ext2 -= 2; 1501 } 1502 } 1503 1504 //Contribution (1,0,0,0) 1505 double dx1 = dx0 - 1 - SQUISH_CONSTANT_4D; 1506 double dy1 = dy0 - 0 - SQUISH_CONSTANT_4D; 1507 double dz1 = dz0 - 0 - SQUISH_CONSTANT_4D; 1508 double dw1 = dw0 - 0 - SQUISH_CONSTANT_4D; 1509 double attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1 - dw1 * dw1; 1510 if (attn1 > 0) { 1511 attn1 *= attn1; 1512 value += attn1 * attn1 * extrapolate(xsb + 1, ysb + 0, zsb + 0, wsb + 0, dx1, dy1, dz1, dw1); 1513 } 1514 1515 //Contribution (0,1,0,0) 1516 double dx2 = dx0 - 0 - SQUISH_CONSTANT_4D; 1517 double dy2 = dy0 - 1 - SQUISH_CONSTANT_4D; 1518 double dz2 = dz1; 1519 double dw2 = dw1; 1520 double attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2 - dw2 * dw2; 1521 if (attn2 > 0) { 1522 attn2 *= attn2; 1523 value += attn2 * attn2 * extrapolate(xsb + 0, ysb + 1, zsb + 0, wsb + 0, dx2, dy2, dz2, dw2); 1524 } 1525 1526 //Contribution (0,0,1,0) 1527 double dx3 = dx2; 1528 double dy3 = dy1; 1529 double dz3 = dz0 - 1 - SQUISH_CONSTANT_4D; 1530 double dw3 = dw1; 1531 double attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3 - dw3 * dw3; 1532 if (attn3 > 0) { 1533 attn3 *= attn3; 1534 value += attn3 * attn3 * extrapolate(xsb + 0, ysb + 0, zsb + 1, wsb + 0, dx3, dy3, dz3, dw3); 1535 } 1536 1537 //Contribution (0,0,0,1) 1538 double dx4 = dx2; 1539 double dy4 = dy1; 1540 double dz4 = dz1; 1541 double dw4 = dw0 - 1 - SQUISH_CONSTANT_4D; 1542 double attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4 - dw4 * dw4; 1543 if (attn4 > 0) { 1544 attn4 *= attn4; 1545 value += attn4 * attn4 * extrapolate(xsb + 0, ysb + 0, zsb + 0, wsb + 1, dx4, dy4, dz4, dw4); 1546 } 1547 1548 //Contribution (1,1,0,0) 1549 double dx5 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1550 double dy5 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1551 double dz5 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; 1552 double dw5 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; 1553 double attn5 = 2 - dx5 * dx5 - dy5 * dy5 - dz5 * dz5 - dw5 * dw5; 1554 if (attn5 > 0) { 1555 attn5 *= attn5; 1556 value += attn5 * attn5 * extrapolate(xsb + 1, ysb + 1, zsb + 0, wsb + 0, dx5, dy5, dz5, dw5); 1557 } 1558 1559 //Contribution (1,0,1,0) 1560 double dx6 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1561 double dy6 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; 1562 double dz6 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1563 double dw6 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; 1564 double attn6 = 2 - dx6 * dx6 - dy6 * dy6 - dz6 * dz6 - dw6 * dw6; 1565 if (attn6 > 0) { 1566 attn6 *= attn6; 1567 value += attn6 * attn6 * extrapolate(xsb + 1, ysb + 0, zsb + 1, wsb + 0, dx6, dy6, dz6, dw6); 1568 } 1569 1570 //Contribution (1,0,0,1) 1571 double dx7 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1572 double dy7 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; 1573 double dz7 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; 1574 double dw7 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1575 double attn7 = 2 - dx7 * dx7 - dy7 * dy7 - dz7 * dz7 - dw7 * dw7; 1576 if (attn7 > 0) { 1577 attn7 *= attn7; 1578 value += attn7 * attn7 * extrapolate(xsb + 1, ysb + 0, zsb + 0, wsb + 1, dx7, dy7, dz7, dw7); 1579 } 1580 1581 //Contribution (0,1,1,0) 1582 double dx8 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; 1583 double dy8 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1584 double dz8 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1585 double dw8 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; 1586 double attn8 = 2 - dx8 * dx8 - dy8 * dy8 - dz8 * dz8 - dw8 * dw8; 1587 if (attn8 > 0) { 1588 attn8 *= attn8; 1589 value += attn8 * attn8 * extrapolate(xsb + 0, ysb + 1, zsb + 1, wsb + 0, dx8, dy8, dz8, dw8); 1590 } 1591 1592 //Contribution (0,1,0,1) 1593 double dx9 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; 1594 double dy9 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1595 double dz9 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; 1596 double dw9 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1597 double attn9 = 2 - dx9 * dx9 - dy9 * dy9 - dz9 * dz9 - dw9 * dw9; 1598 if (attn9 > 0) { 1599 attn9 *= attn9; 1600 value += attn9 * attn9 * extrapolate(xsb + 0, ysb + 1, zsb + 0, wsb + 1, dx9, dy9, dz9, dw9); 1601 } 1602 1603 //Contribution (0,0,1,1) 1604 double dx10 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; 1605 double dy10 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; 1606 double dz10 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1607 double dw10 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1608 double attn10 = 2 - dx10 * dx10 - dy10 * dy10 - dz10 * dz10 - dw10 * dw10; 1609 if (attn10 > 0) { 1610 attn10 *= attn10; 1611 value += attn10 * attn10 * extrapolate(xsb + 0, ysb + 0, zsb + 1, wsb + 1, dx10, dy10, dz10, dw10); 1612 } 1613 } else { //We're inside the second dispentachoron (Rectified 4-Simplex) 1614 double aScore; 1615 byte aPoint; 1616 bool aIsBiggerSide = true; 1617 double bScore; 1618 byte bPoint; 1619 bool bIsBiggerSide = true; 1620 1621 //Decide between (0,0,1,1) and (1,1,0,0) 1622 if (xins + yins < zins + wins) { 1623 aScore = xins + yins; 1624 aPoint = 0x0C; 1625 } else { 1626 aScore = zins + wins; 1627 aPoint = 0x03; 1628 } 1629 1630 //Decide between (0,1,0,1) and (1,0,1,0) 1631 if (xins + zins < yins + wins) { 1632 bScore = xins + zins; 1633 bPoint = 0x0A; 1634 } else { 1635 bScore = yins + wins; 1636 bPoint = 0x05; 1637 } 1638 1639 //Closer between (0,1,1,0) and (1,0,0,1) will replace the further of a and b, if closer. 1640 if (xins + wins < yins + zins) { 1641 double score = xins + wins; 1642 if (aScore <= bScore && score < bScore) { 1643 bScore = score; 1644 bPoint = 0x06; 1645 } else if (aScore > bScore && score < aScore) { 1646 aScore = score; 1647 aPoint = 0x06; 1648 } 1649 } else { 1650 double score = yins + zins; 1651 if (aScore <= bScore && score < bScore) { 1652 bScore = score; 1653 bPoint = 0x09; 1654 } else if (aScore > bScore && score < aScore) { 1655 aScore = score; 1656 aPoint = 0x09; 1657 } 1658 } 1659 1660 //Decide if (0,1,1,1) is closer. 1661 double p1 = 3 - inSum + xins; 1662 if (aScore <= bScore && p1 < bScore) { 1663 bScore = p1; 1664 bPoint = 0x0E; 1665 bIsBiggerSide = false; 1666 } else if (aScore > bScore && p1 < aScore) { 1667 aScore = p1; 1668 aPoint = 0x0E; 1669 aIsBiggerSide = false; 1670 } 1671 1672 //Decide if (1,0,1,1) is closer. 1673 double p2 = 3 - inSum + yins; 1674 if (aScore <= bScore && p2 < bScore) { 1675 bScore = p2; 1676 bPoint = 0x0D; 1677 bIsBiggerSide = false; 1678 } else if (aScore > bScore && p2 < aScore) { 1679 aScore = p2; 1680 aPoint = 0x0D; 1681 aIsBiggerSide = false; 1682 } 1683 1684 //Decide if (1,1,0,1) is closer. 1685 double p3 = 3 - inSum + zins; 1686 if (aScore <= bScore && p3 < bScore) { 1687 bScore = p3; 1688 bPoint = 0x0B; 1689 bIsBiggerSide = false; 1690 } else if (aScore > bScore && p3 < aScore) { 1691 aScore = p3; 1692 aPoint = 0x0B; 1693 aIsBiggerSide = false; 1694 } 1695 1696 //Decide if (1,1,1,0) is closer. 1697 double p4 = 3 - inSum + wins; 1698 if (aScore <= bScore && p4 < bScore) { 1699 bScore = p4; 1700 bPoint = 0x07; 1701 bIsBiggerSide = false; 1702 } else if (aScore > bScore && p4 < aScore) { 1703 aScore = p4; 1704 aPoint = 0x07; 1705 aIsBiggerSide = false; 1706 } 1707 1708 //Where each of the two closest points are determines how the extra three vertices are calculated. 1709 if (aIsBiggerSide == bIsBiggerSide) { 1710 if (aIsBiggerSide) { //Both closest points on the bigger side 1711 byte c1 = cast(byte)(aPoint & bPoint); 1712 byte c2 = cast(byte)(aPoint | bPoint); 1713 1714 //Two contributions are permutations of (0,0,0,1) and (0,0,0,2) based on c1 1715 xsv_ext0 = xsv_ext1 = xsb; 1716 ysv_ext0 = ysv_ext1 = ysb; 1717 zsv_ext0 = zsv_ext1 = zsb; 1718 wsv_ext0 = wsv_ext1 = wsb; 1719 dx_ext0 = dx0 - SQUISH_CONSTANT_4D; 1720 dy_ext0 = dy0 - SQUISH_CONSTANT_4D; 1721 dz_ext0 = dz0 - SQUISH_CONSTANT_4D; 1722 dw_ext0 = dw0 - SQUISH_CONSTANT_4D; 1723 dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_4D; 1724 dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_4D; 1725 dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_4D; 1726 dw_ext1 = dw0 - 2 * SQUISH_CONSTANT_4D; 1727 if ((c1 & 0x01) != 0) { 1728 xsv_ext0 += 1; 1729 dx_ext0 -= 1; 1730 xsv_ext1 += 2; 1731 dx_ext1 -= 2; 1732 } else if ((c1 & 0x02) != 0) { 1733 ysv_ext0 += 1; 1734 dy_ext0 -= 1; 1735 ysv_ext1 += 2; 1736 dy_ext1 -= 2; 1737 } else if ((c1 & 0x04) != 0) { 1738 zsv_ext0 += 1; 1739 dz_ext0 -= 1; 1740 zsv_ext1 += 2; 1741 dz_ext1 -= 2; 1742 } else { 1743 wsv_ext0 += 1; 1744 dw_ext0 -= 1; 1745 wsv_ext1 += 2; 1746 dw_ext1 -= 2; 1747 } 1748 1749 //One contribution is a permutation of (1,1,1,-1) based on c2 1750 xsv_ext2 = xsb + 1; 1751 ysv_ext2 = ysb + 1; 1752 zsv_ext2 = zsb + 1; 1753 wsv_ext2 = wsb + 1; 1754 dx_ext2 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1755 dy_ext2 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1756 dz_ext2 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1757 dw_ext2 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1758 if ((c2 & 0x01) == 0) { 1759 xsv_ext2 -= 2; 1760 dx_ext2 += 2; 1761 } else if ((c2 & 0x02) == 0) { 1762 ysv_ext2 -= 2; 1763 dy_ext2 += 2; 1764 } else if ((c2 & 0x04) == 0) { 1765 zsv_ext2 -= 2; 1766 dz_ext2 += 2; 1767 } else { 1768 wsv_ext2 -= 2; 1769 dw_ext2 += 2; 1770 } 1771 } else { //Both closest points on the smaller side 1772 //One of the two extra points is (1,1,1,1) 1773 xsv_ext2 = xsb + 1; 1774 ysv_ext2 = ysb + 1; 1775 zsv_ext2 = zsb + 1; 1776 wsv_ext2 = wsb + 1; 1777 dx_ext2 = dx0 - 1 - 4 * SQUISH_CONSTANT_4D; 1778 dy_ext2 = dy0 - 1 - 4 * SQUISH_CONSTANT_4D; 1779 dz_ext2 = dz0 - 1 - 4 * SQUISH_CONSTANT_4D; 1780 dw_ext2 = dw0 - 1 - 4 * SQUISH_CONSTANT_4D; 1781 1782 //Other two points are based on the shared axes. 1783 byte c = cast(byte)(aPoint & bPoint); 1784 1785 if ((c & 0x01) != 0) { 1786 xsv_ext0 = xsb + 2; 1787 xsv_ext1 = xsb + 1; 1788 dx_ext0 = dx0 - 2 - 3 * SQUISH_CONSTANT_4D; 1789 dx_ext1 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; 1790 } else { 1791 xsv_ext0 = xsv_ext1 = xsb; 1792 dx_ext0 = dx_ext1 = dx0 - 3 * SQUISH_CONSTANT_4D; 1793 } 1794 1795 if ((c & 0x02) != 0) { 1796 ysv_ext0 = ysv_ext1 = ysb + 1; 1797 dy_ext0 = dy_ext1 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; 1798 if ((c & 0x01) == 0) 1799 { 1800 ysv_ext0 += 1; 1801 dy_ext0 -= 1; 1802 } else { 1803 ysv_ext1 += 1; 1804 dy_ext1 -= 1; 1805 } 1806 } else { 1807 ysv_ext0 = ysv_ext1 = ysb; 1808 dy_ext0 = dy_ext1 = dy0 - 3 * SQUISH_CONSTANT_4D; 1809 } 1810 1811 if ((c & 0x04) != 0) { 1812 zsv_ext0 = zsv_ext1 = zsb + 1; 1813 dz_ext0 = dz_ext1 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; 1814 if ((c & 0x03) == 0) 1815 { 1816 zsv_ext0 += 1; 1817 dz_ext0 -= 1; 1818 } else { 1819 zsv_ext1 += 1; 1820 dz_ext1 -= 1; 1821 } 1822 } else { 1823 zsv_ext0 = zsv_ext1 = zsb; 1824 dz_ext0 = dz_ext1 = dz0 - 3 * SQUISH_CONSTANT_4D; 1825 } 1826 1827 if ((c & 0x08) != 0) 1828 { 1829 wsv_ext0 = wsb + 1; 1830 wsv_ext1 = wsb + 2; 1831 dw_ext0 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; 1832 dw_ext1 = dw0 - 2 - 3 * SQUISH_CONSTANT_4D; 1833 } else { 1834 wsv_ext0 = wsv_ext1 = wsb; 1835 dw_ext0 = dw_ext1 = dw0 - 3 * SQUISH_CONSTANT_4D; 1836 } 1837 } 1838 } else { //One point on each "side" 1839 byte c1, c2; 1840 if (aIsBiggerSide) { 1841 c1 = aPoint; 1842 c2 = bPoint; 1843 } else { 1844 c1 = bPoint; 1845 c2 = aPoint; 1846 } 1847 1848 //Two contributions are the bigger-sided point with each 1 replaced with 2. 1849 if ((c1 & 0x01) != 0) { 1850 xsv_ext0 = xsb + 2; 1851 xsv_ext1 = xsb + 1; 1852 dx_ext0 = dx0 - 2 - 3 * SQUISH_CONSTANT_4D; 1853 dx_ext1 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; 1854 } else { 1855 xsv_ext0 = xsv_ext1 = xsb; 1856 dx_ext0 = dx_ext1 = dx0 - 3 * SQUISH_CONSTANT_4D; 1857 } 1858 1859 if ((c1 & 0x02) != 0) { 1860 ysv_ext0 = ysv_ext1 = ysb + 1; 1861 dy_ext0 = dy_ext1 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; 1862 if ((c1 & 0x01) == 0) { 1863 ysv_ext0 += 1; 1864 dy_ext0 -= 1; 1865 } else { 1866 ysv_ext1 += 1; 1867 dy_ext1 -= 1; 1868 } 1869 } else { 1870 ysv_ext0 = ysv_ext1 = ysb; 1871 dy_ext0 = dy_ext1 = dy0 - 3 * SQUISH_CONSTANT_4D; 1872 } 1873 1874 if ((c1 & 0x04) != 0) { 1875 zsv_ext0 = zsv_ext1 = zsb + 1; 1876 dz_ext0 = dz_ext1 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; 1877 if ((c1 & 0x03) == 0) { 1878 zsv_ext0 += 1; 1879 dz_ext0 -= 1; 1880 } else { 1881 zsv_ext1 += 1; 1882 dz_ext1 -= 1; 1883 } 1884 } else { 1885 zsv_ext0 = zsv_ext1 = zsb; 1886 dz_ext0 = dz_ext1 = dz0 - 3 * SQUISH_CONSTANT_4D; 1887 } 1888 1889 if ((c1 & 0x08) != 0) { 1890 wsv_ext0 = wsb + 1; 1891 wsv_ext1 = wsb + 2; 1892 dw_ext0 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; 1893 dw_ext1 = dw0 - 2 - 3 * SQUISH_CONSTANT_4D; 1894 } else { 1895 wsv_ext0 = wsv_ext1 = wsb; 1896 dw_ext0 = dw_ext1 = dw0 - 3 * SQUISH_CONSTANT_4D; 1897 } 1898 1899 //One contribution is a permutation of (1,1,1,-1) based on the smaller-sided point 1900 xsv_ext2 = xsb + 1; 1901 ysv_ext2 = ysb + 1; 1902 zsv_ext2 = zsb + 1; 1903 wsv_ext2 = wsb + 1; 1904 dx_ext2 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1905 dy_ext2 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1906 dz_ext2 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1907 dw_ext2 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1908 if ((c2 & 0x01) == 0) { 1909 xsv_ext2 -= 2; 1910 dx_ext2 += 2; 1911 } else if ((c2 & 0x02) == 0) { 1912 ysv_ext2 -= 2; 1913 dy_ext2 += 2; 1914 } else if ((c2 & 0x04) == 0) { 1915 zsv_ext2 -= 2; 1916 dz_ext2 += 2; 1917 } else { 1918 wsv_ext2 -= 2; 1919 dw_ext2 += 2; 1920 } 1921 } 1922 1923 //Contribution (1,1,1,0) 1924 double dx4 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; 1925 double dy4 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; 1926 double dz4 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; 1927 double dw4 = dw0 - 3 * SQUISH_CONSTANT_4D; 1928 double attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4 - dw4 * dw4; 1929 if (attn4 > 0) { 1930 attn4 *= attn4; 1931 value += attn4 * attn4 * extrapolate(xsb + 1, ysb + 1, zsb + 1, wsb + 0, dx4, dy4, dz4, dw4); 1932 } 1933 1934 //Contribution (1,1,0,1) 1935 double dx3 = dx4; 1936 double dy3 = dy4; 1937 double dz3 = dz0 - 3 * SQUISH_CONSTANT_4D; 1938 double dw3 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; 1939 double attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3 - dw3 * dw3; 1940 if (attn3 > 0) { 1941 attn3 *= attn3; 1942 value += attn3 * attn3 * extrapolate(xsb + 1, ysb + 1, zsb + 0, wsb + 1, dx3, dy3, dz3, dw3); 1943 } 1944 1945 //Contribution (1,0,1,1) 1946 double dx2 = dx4; 1947 double dy2 = dy0 - 3 * SQUISH_CONSTANT_4D; 1948 double dz2 = dz4; 1949 double dw2 = dw3; 1950 double attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2 - dw2 * dw2; 1951 if (attn2 > 0) { 1952 attn2 *= attn2; 1953 value += attn2 * attn2 * extrapolate(xsb + 1, ysb + 0, zsb + 1, wsb + 1, dx2, dy2, dz2, dw2); 1954 } 1955 1956 //Contribution (0,1,1,1) 1957 double dx1 = dx0 - 3 * SQUISH_CONSTANT_4D; 1958 double dz1 = dz4; 1959 double dy1 = dy4; 1960 double dw1 = dw3; 1961 double attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1 - dw1 * dw1; 1962 if (attn1 > 0) { 1963 attn1 *= attn1; 1964 value += attn1 * attn1 * extrapolate(xsb + 0, ysb + 1, zsb + 1, wsb + 1, dx1, dy1, dz1, dw1); 1965 } 1966 1967 //Contribution (1,1,0,0) 1968 double dx5 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1969 double dy5 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1970 double dz5 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; 1971 double dw5 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; 1972 double attn5 = 2 - dx5 * dx5 - dy5 * dy5 - dz5 * dz5 - dw5 * dw5; 1973 if (attn5 > 0) { 1974 attn5 *= attn5; 1975 value += attn5 * attn5 * extrapolate(xsb + 1, ysb + 1, zsb + 0, wsb + 0, dx5, dy5, dz5, dw5); 1976 } 1977 1978 //Contribution (1,0,1,0) 1979 double dx6 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1980 double dy6 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; 1981 double dz6 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1982 double dw6 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; 1983 double attn6 = 2 - dx6 * dx6 - dy6 * dy6 - dz6 * dz6 - dw6 * dw6; 1984 if (attn6 > 0) { 1985 attn6 *= attn6; 1986 value += attn6 * attn6 * extrapolate(xsb + 1, ysb + 0, zsb + 1, wsb + 0, dx6, dy6, dz6, dw6); 1987 } 1988 1989 //Contribution (1,0,0,1) 1990 double dx7 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1991 double dy7 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; 1992 double dz7 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; 1993 double dw7 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1994 double attn7 = 2 - dx7 * dx7 - dy7 * dy7 - dz7 * dz7 - dw7 * dw7; 1995 if (attn7 > 0) { 1996 attn7 *= attn7; 1997 value += attn7 * attn7 * extrapolate(xsb + 1, ysb + 0, zsb + 0, wsb + 1, dx7, dy7, dz7, dw7); 1998 } 1999 2000 //Contribution (0,1,1,0) 2001 double dx8 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; 2002 double dy8 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 2003 double dz8 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 2004 double dw8 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; 2005 double attn8 = 2 - dx8 * dx8 - dy8 * dy8 - dz8 * dz8 - dw8 * dw8; 2006 if (attn8 > 0) { 2007 attn8 *= attn8; 2008 value += attn8 * attn8 * extrapolate(xsb + 0, ysb + 1, zsb + 1, wsb + 0, dx8, dy8, dz8, dw8); 2009 } 2010 2011 //Contribution (0,1,0,1) 2012 double dx9 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; 2013 double dy9 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 2014 double dz9 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; 2015 double dw9 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 2016 double attn9 = 2 - dx9 * dx9 - dy9 * dy9 - dz9 * dz9 - dw9 * dw9; 2017 if (attn9 > 0) { 2018 attn9 *= attn9; 2019 value += attn9 * attn9 * extrapolate(xsb + 0, ysb + 1, zsb + 0, wsb + 1, dx9, dy9, dz9, dw9); 2020 } 2021 2022 //Contribution (0,0,1,1) 2023 double dx10 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; 2024 double dy10 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; 2025 double dz10 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 2026 double dw10 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 2027 double attn10 = 2 - dx10 * dx10 - dy10 * dy10 - dz10 * dz10 - dw10 * dw10; 2028 if (attn10 > 0) { 2029 attn10 *= attn10; 2030 value += attn10 * attn10 * extrapolate(xsb + 0, ysb + 0, zsb + 1, wsb + 1, dx10, dy10, dz10, dw10); 2031 } 2032 } 2033 2034 //First extra vertex 2035 double attn_ext0 = 2 - dx_ext0 * dx_ext0 - dy_ext0 * dy_ext0 - dz_ext0 * dz_ext0 - dw_ext0 * dw_ext0; 2036 if (attn_ext0 > 0) 2037 { 2038 attn_ext0 *= attn_ext0; 2039 value += attn_ext0 * attn_ext0 * extrapolate(xsv_ext0, ysv_ext0, zsv_ext0, wsv_ext0, dx_ext0, dy_ext0, dz_ext0, dw_ext0); 2040 } 2041 2042 //Second extra vertex 2043 double attn_ext1 = 2 - dx_ext1 * dx_ext1 - dy_ext1 * dy_ext1 - dz_ext1 * dz_ext1 - dw_ext1 * dw_ext1; 2044 if (attn_ext1 > 0) 2045 { 2046 attn_ext1 *= attn_ext1; 2047 value += attn_ext1 * attn_ext1 * extrapolate(xsv_ext1, ysv_ext1, zsv_ext1, wsv_ext1, dx_ext1, dy_ext1, dz_ext1, dw_ext1); 2048 } 2049 2050 //Third extra vertex 2051 double attn_ext2 = 2 - dx_ext2 * dx_ext2 - dy_ext2 * dy_ext2 - dz_ext2 * dz_ext2 - dw_ext2 * dw_ext2; 2052 if (attn_ext2 > 0) 2053 { 2054 attn_ext2 *= attn_ext2; 2055 value += attn_ext2 * attn_ext2 * extrapolate(xsv_ext2, ysv_ext2, zsv_ext2, wsv_ext2, dx_ext2, dy_ext2, dz_ext2, dw_ext2); 2056 } 2057 2058 return value / NORM_CONSTANT_4D; 2059 } 2060 2061 private double extrapolate(int xsb, int ysb, double dx, double dy) 2062 { 2063 int index = perm[(perm[xsb & 0xFF] + ysb) & 0xFF] & 0x0E; 2064 return gradients2D[index] * dx 2065 + gradients2D[index + 1] * dy; 2066 } 2067 2068 private double extrapolate(int xsb, int ysb, int zsb, double dx, double dy, double dz) 2069 { 2070 int index = permGradIndex3D[(perm[(perm[xsb & 0xFF] + ysb) & 0xFF] + zsb) & 0xFF]; 2071 return gradients3D[index] * dx 2072 + gradients3D[index + 1] * dy 2073 + gradients3D[index + 2] * dz; 2074 } 2075 2076 private double extrapolate(int xsb, int ysb, int zsb, int wsb, double dx, double dy, double dz, double dw) 2077 { 2078 int index = perm[(perm[(perm[(perm[xsb & 0xFF] + ysb) & 0xFF] + zsb) & 0xFF] + wsb) & 0xFF] & 0xFC; 2079 return gradients4D[index] * dx 2080 + gradients4D[index + 1] * dy 2081 + gradients4D[index + 2] * dz 2082 + gradients4D[index + 3] * dw; 2083 } 2084 2085 private static int fastFloor(double x) { 2086 int xi = cast(int)x; 2087 return x < xi ? xi - 1 : xi; 2088 } 2089 2090 //Gradients for 2D. They approximate the directions to the 2091 //vertices of an octagon from the center. 2092 private static byte[16] gradients2D = [ 2093 5, 2, 2, 5, 2094 -5, 2, -2, 5, 2095 5, -2, 2, -5, 2096 -5, -2, -2, -5, 2097 ]; 2098 2099 //Gradients for 3D. They approximate the directions to the 2100 //vertices of a rhombicuboctahedron from the center, skewed so 2101 //that the triangular and square facets can be inscribed inside 2102 //circles of the same radius. 2103 private static byte[72] gradients3D = [ 2104 -11, 4, 4, -4, 11, 4, -4, 4, 11, 2105 11, 4, 4, 4, 11, 4, 4, 4, 11, 2106 -11, -4, 4, -4, -11, 4, -4, -4, 11, 2107 11, -4, 4, 4, -11, 4, 4, -4, 11, 2108 -11, 4, -4, -4, 11, -4, -4, 4, -11, 2109 11, 4, -4, 4, 11, -4, 4, 4, -11, 2110 -11, -4, -4, -4, -11, -4, -4, -4, -11, 2111 11, -4, -4, 4, -11, -4, 4, -4, -11, 2112 ]; 2113 2114 //Gradients for 4D. They approximate the directions to the 2115 //vertices of a disprismatotesseractihexadecachoron from the center, 2116 //skewed so that the tetrahedral and cubic facets can be inscribed inside 2117 //spheres of the same radius. 2118 private static byte[256] gradients4D = [ 2119 3, 1, 1, 1, 1, 3, 1, 1, 1, 1, 3, 1, 1, 1, 1, 3, 2120 -3, 1, 1, 1, -1, 3, 1, 1, -1, 1, 3, 1, -1, 1, 1, 3, 2121 3, -1, 1, 1, 1, -3, 1, 1, 1, -1, 3, 1, 1, -1, 1, 3, 2122 -3, -1, 1, 1, -1, -3, 1, 1, -1, -1, 3, 1, -1, -1, 1, 3, 2123 3, 1, -1, 1, 1, 3, -1, 1, 1, 1, -3, 1, 1, 1, -1, 3, 2124 -3, 1, -1, 1, -1, 3, -1, 1, -1, 1, -3, 1, -1, 1, -1, 3, 2125 3, -1, -1, 1, 1, -3, -1, 1, 1, -1, -3, 1, 1, -1, -1, 3, 2126 -3, -1, -1, 1, -1, -3, -1, 1, -1, -1, -3, 1, -1, -1, -1, 3, 2127 3, 1, 1, -1, 1, 3, 1, -1, 1, 1, 3, -1, 1, 1, 1, -3, 2128 -3, 1, 1, -1, -1, 3, 1, -1, -1, 1, 3, -1, -1, 1, 1, -3, 2129 3, -1, 1, -1, 1, -3, 1, -1, 1, -1, 3, -1, 1, -1, 1, -3, 2130 -3, -1, 1, -1, -1, -3, 1, -1, -1, -1, 3, -1, -1, -1, 1, -3, 2131 3, 1, -1, -1, 1, 3, -1, -1, 1, 1, -3, -1, 1, 1, -1, -3, 2132 -3, 1, -1, -1, -1, 3, -1, -1, -1, 1, -3, -1, -1, 1, -1, -3, 2133 3, -1, -1, -1, 1, -3, -1, -1, 1, -1, -3, -1, 1, -1, -1, -3, 2134 -3, -1, -1, -1, -1, -3, -1, -1, -1, -1, -3, -1, -1, -1, -1, -3, 2135 ]; 2136 }