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