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