1 /** 2 Holds the (mor or less) independet OpenSimplexNoise implementation - original source and implementation: https://gist.github.com/KdotJPG/b1270127455a94ac5d19 and http://uniblock.tumblr.com/ 3 */ 4 module dosimplex.osimplex4d; 5 6 import dosimplex.util; 7 8 /// "Config" enum for the noise generator 9 private enum : double { 10 STRETCH_CONSTANT_4D = -0.138196601125011, ///(1/Math.sqrt(4+1)-1)/4; 11 SQUISH_CONSTANT_4D = 0.309016994374947, ///(Math.sqrt(4+1)-1)/4; 12 NORM_CONSTANT_4D = 30., 13 }; 14 15 @nogc @safe pure double osNoise4D(double x, double y, double z, double w, const ref short[256] perm) { 16 17 //Place input coordinates on simplectic honeycomb. 18 double stretchOffset = (x + y + z + w) * STRETCH_CONSTANT_4D; 19 double xs = x + stretchOffset; 20 double ys = y + stretchOffset; 21 double zs = z + stretchOffset; 22 double ws = w + stretchOffset; 23 24 //Floor to get simplectic honeycomb coordinates of rhombo-hypercube super-cell origin. 25 int xsb = fastFloor(xs); 26 int ysb = fastFloor(ys); 27 int zsb = fastFloor(zs); 28 int wsb = fastFloor(ws); 29 30 //Skew out to get actual coordinates of stretched rhombo-hypercube origin. We'll need these later. 31 double squishOffset = (xsb + ysb + zsb + wsb) * SQUISH_CONSTANT_4D; 32 double xb = xsb + squishOffset; 33 double yb = ysb + squishOffset; 34 double zb = zsb + squishOffset; 35 double wb = wsb + squishOffset; 36 37 //Compute simplectic honeycomb coordinates relative to rhombo-hypercube origin. 38 double xins = xs - xsb; 39 double yins = ys - ysb; 40 double zins = zs - zsb; 41 double wins = ws - wsb; 42 43 //Sum those together to get a value that determines which region we're in. 44 double inSum = xins + yins + zins + wins; 45 46 //Positions relative to origin point. 47 double dx0 = x - xb; 48 double dy0 = y - yb; 49 double dz0 = z - zb; 50 double dw0 = w - wb; 51 52 //We'll be defining these inside the next block and using them afterwards. 53 double dx_ext0, dy_ext0, dz_ext0, dw_ext0; 54 double dx_ext1, dy_ext1, dz_ext1, dw_ext1; 55 double dx_ext2, dy_ext2, dz_ext2, dw_ext2; 56 int xsv_ext0, ysv_ext0, zsv_ext0, wsv_ext0; 57 int xsv_ext1, ysv_ext1, zsv_ext1, wsv_ext1; 58 int xsv_ext2, ysv_ext2, zsv_ext2, wsv_ext2; 59 60 double value = 0; 61 if (inSum <= 1) { //We're inside the pentachoron (4-Simplex) at (0,0,0,0) 62 63 //Determine which two of (0,0,0,1), (0,0,1,0), (0,1,0,0), (1,0,0,0) are closest. 64 byte aPoint = 0x01; 65 double aScore = xins; 66 byte bPoint = 0x02; 67 double bScore = yins; 68 if (aScore >= bScore && zins > bScore) { 69 bScore = zins; 70 bPoint = 0x04; 71 } else if (aScore < bScore && zins > aScore) { 72 aScore = zins; 73 aPoint = 0x04; 74 } 75 if (aScore >= bScore && wins > bScore) { 76 bScore = wins; 77 bPoint = 0x08; 78 } else if (aScore < bScore && wins > aScore) { 79 aScore = wins; 80 aPoint = 0x08; 81 } 82 83 //Now we determine the three lattice points not part of the pentachoron that may contribute. 84 //This depends on the closest two pentachoron vertices, including (0,0,0,0) 85 double uins = 1 - inSum; 86 if (uins > aScore || uins > bScore) { //(0,0,0,0) is one of the closest two pentachoron vertices. 87 byte c = (bScore > aScore ? bPoint : aPoint); //Our other closest vertex is the closest out of a and b. 88 if ((c & 0x01) == 0) { 89 xsv_ext0 = xsb - 1; 90 xsv_ext1 = xsv_ext2 = xsb; 91 dx_ext0 = dx0 + 1; 92 dx_ext1 = dx_ext2 = dx0; 93 } else { 94 xsv_ext0 = xsv_ext1 = xsv_ext2 = xsb + 1; 95 dx_ext0 = dx_ext1 = dx_ext2 = dx0 - 1; 96 } 97 98 if ((c & 0x02) == 0) { 99 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb; 100 dy_ext0 = dy_ext1 = dy_ext2 = dy0; 101 if ((c & 0x01) == 0x01) { 102 ysv_ext0 -= 1; 103 dy_ext0 += 1; 104 } else { 105 ysv_ext1 -= 1; 106 dy_ext1 += 1; 107 } 108 } else { 109 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb + 1; 110 dy_ext0 = dy_ext1 = dy_ext2 = dy0 - 1; 111 } 112 113 if ((c & 0x04) == 0) { 114 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb; 115 dz_ext0 = dz_ext1 = dz_ext2 = dz0; 116 if ((c & 0x03) != 0) { 117 if ((c & 0x03) == 0x03) { 118 zsv_ext0 -= 1; 119 dz_ext0 += 1; 120 } else { 121 zsv_ext1 -= 1; 122 dz_ext1 += 1; 123 } 124 } else { 125 zsv_ext2 -= 1; 126 dz_ext2 += 1; 127 } 128 } else { 129 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb + 1; 130 dz_ext0 = dz_ext1 = dz_ext2 = dz0 - 1; 131 } 132 133 if ((c & 0x08) == 0) { 134 wsv_ext0 = wsv_ext1 = wsb; 135 wsv_ext2 = wsb - 1; 136 dw_ext0 = dw_ext1 = dw0; 137 dw_ext2 = dw0 + 1; 138 } else { 139 wsv_ext0 = wsv_ext1 = wsv_ext2 = wsb + 1; 140 dw_ext0 = dw_ext1 = dw_ext2 = dw0 - 1; 141 } 142 } else { //(0,0,0,0) is not one of the closest two pentachoron vertices. 143 byte c = cast(byte)(aPoint | bPoint); //Our three extra vertices are determined by the closest two. 144 145 if ((c & 0x01) == 0) { 146 xsv_ext0 = xsv_ext2 = xsb; 147 xsv_ext1 = xsb - 1; 148 dx_ext0 = dx0 - 2 * SQUISH_CONSTANT_4D; 149 dx_ext1 = dx0 + 1 - SQUISH_CONSTANT_4D; 150 dx_ext2 = dx0 - SQUISH_CONSTANT_4D; 151 } else { 152 xsv_ext0 = xsv_ext1 = xsv_ext2 = xsb + 1; 153 dx_ext0 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 154 dx_ext1 = dx_ext2 = dx0 - 1 - SQUISH_CONSTANT_4D; 155 } 156 157 if ((c & 0x02) == 0) { 158 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb; 159 dy_ext0 = dy0 - 2 * SQUISH_CONSTANT_4D; 160 dy_ext1 = dy_ext2 = dy0 - SQUISH_CONSTANT_4D; 161 if ((c & 0x01) == 0x01) { 162 ysv_ext1 -= 1; 163 dy_ext1 += 1; 164 } else { 165 ysv_ext2 -= 1; 166 dy_ext2 += 1; 167 } 168 } else { 169 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb + 1; 170 dy_ext0 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 171 dy_ext1 = dy_ext2 = dy0 - 1 - SQUISH_CONSTANT_4D; 172 } 173 174 if ((c & 0x04) == 0) { 175 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb; 176 dz_ext0 = dz0 - 2 * SQUISH_CONSTANT_4D; 177 dz_ext1 = dz_ext2 = dz0 - SQUISH_CONSTANT_4D; 178 if ((c & 0x03) == 0x03) { 179 zsv_ext1 -= 1; 180 dz_ext1 += 1; 181 } else { 182 zsv_ext2 -= 1; 183 dz_ext2 += 1; 184 } 185 } else { 186 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb + 1; 187 dz_ext0 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 188 dz_ext1 = dz_ext2 = dz0 - 1 - SQUISH_CONSTANT_4D; 189 } 190 191 if ((c & 0x08) == 0) { 192 wsv_ext0 = wsv_ext1 = wsb; 193 wsv_ext2 = wsb - 1; 194 dw_ext0 = dw0 - 2 * SQUISH_CONSTANT_4D; 195 dw_ext1 = dw0 - SQUISH_CONSTANT_4D; 196 dw_ext2 = dw0 + 1 - SQUISH_CONSTANT_4D; 197 } else { 198 wsv_ext0 = wsv_ext1 = wsv_ext2 = wsb + 1; 199 dw_ext0 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 200 dw_ext1 = dw_ext2 = dw0 - 1 - SQUISH_CONSTANT_4D; 201 } 202 } 203 204 //Contribution (0,0,0,0) 205 double attn0 = 2 - dx0 * dx0 - dy0 * dy0 - dz0 * dz0 - dw0 * dw0; 206 if (attn0 > 0) { 207 attn0 *= attn0; 208 value += attn0 * attn0 * extrapolate(xsb + 0, ysb + 0, zsb + 0, wsb + 0, dx0, dy0, dz0, dw0, perm); 209 } 210 211 //Contribution (1,0,0,0) 212 double dx1 = dx0 - 1 - SQUISH_CONSTANT_4D; 213 double dy1 = dy0 - 0 - SQUISH_CONSTANT_4D; 214 double dz1 = dz0 - 0 - SQUISH_CONSTANT_4D; 215 double dw1 = dw0 - 0 - SQUISH_CONSTANT_4D; 216 double attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1 - dw1 * dw1; 217 if (attn1 > 0) { 218 attn1 *= attn1; 219 value += attn1 * attn1 * extrapolate(xsb + 1, ysb + 0, zsb + 0, wsb + 0, dx1, dy1, dz1, dw1, perm); 220 } 221 222 //Contribution (0,1,0,0) 223 double dx2 = dx0 - 0 - SQUISH_CONSTANT_4D; 224 double dy2 = dy0 - 1 - SQUISH_CONSTANT_4D; 225 double dz2 = dz1; 226 double dw2 = dw1; 227 double attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2 - dw2 * dw2; 228 if (attn2 > 0) { 229 attn2 *= attn2; 230 value += attn2 * attn2 * extrapolate(xsb + 0, ysb + 1, zsb + 0, wsb + 0, dx2, dy2, dz2, dw2, perm); 231 } 232 233 //Contribution (0,0,1,0) 234 double dx3 = dx2; 235 double dy3 = dy1; 236 double dz3 = dz0 - 1 - SQUISH_CONSTANT_4D; 237 double dw3 = dw1; 238 double attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3 - dw3 * dw3; 239 if (attn3 > 0) { 240 attn3 *= attn3; 241 value += attn3 * attn3 * extrapolate(xsb + 0, ysb + 0, zsb + 1, wsb + 0, dx3, dy3, dz3, dw3, perm); 242 } 243 244 //Contribution (0,0,0,1) 245 double dx4 = dx2; 246 double dy4 = dy1; 247 double dz4 = dz1; 248 double dw4 = dw0 - 1 - SQUISH_CONSTANT_4D; 249 double attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4 - dw4 * dw4; 250 if (attn4 > 0) { 251 attn4 *= attn4; 252 value += attn4 * attn4 * extrapolate(xsb + 0, ysb + 0, zsb + 0, wsb + 1, dx4, dy4, dz4, dw4, perm); 253 } 254 } else if (inSum >= 3) { //We're inside the pentachoron (4-Simplex) at (1,1,1,1) 255 //Determine which two of (1,1,1,0), (1,1,0,1), (1,0,1,1), (0,1,1,1) are closest. 256 byte aPoint = 0x0E; 257 double aScore = xins; 258 byte bPoint = 0x0D; 259 double bScore = yins; 260 if (aScore <= bScore && zins < bScore) { 261 bScore = zins; 262 bPoint = 0x0B; 263 } else if (aScore > bScore && zins < aScore) { 264 aScore = zins; 265 aPoint = 0x0B; 266 } 267 if (aScore <= bScore && wins < bScore) { 268 bScore = wins; 269 bPoint = 0x07; 270 } else if (aScore > bScore && wins < aScore) { 271 aScore = wins; 272 aPoint = 0x07; 273 } 274 275 //Now we determine the three lattice points not part of the pentachoron that may contribute. 276 //This depends on the closest two pentachoron vertices, including (0,0,0,0) 277 double uins = 4 - inSum; 278 if (uins < aScore || uins < bScore) { //(1,1,1,1) is one of the closest two pentachoron 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 + 2; 283 xsv_ext1 = xsv_ext2 = xsb + 1; 284 dx_ext0 = dx0 - 2 - 4 * SQUISH_CONSTANT_4D; 285 dx_ext1 = dx_ext2 = dx0 - 1 - 4 * SQUISH_CONSTANT_4D; 286 } else { 287 xsv_ext0 = xsv_ext1 = xsv_ext2 = xsb; 288 dx_ext0 = dx_ext1 = dx_ext2 = dx0 - 4 * SQUISH_CONSTANT_4D; 289 } 290 291 if ((c & 0x02) != 0) { 292 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb + 1; 293 dy_ext0 = dy_ext1 = dy_ext2 = dy0 - 1 - 4 * SQUISH_CONSTANT_4D; 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 = ysv_ext2 = ysb; 303 dy_ext0 = dy_ext1 = dy_ext2 = dy0 - 4 * SQUISH_CONSTANT_4D; 304 } 305 306 if ((c & 0x04) != 0) { 307 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb + 1; 308 dz_ext0 = dz_ext1 = dz_ext2 = dz0 - 1 - 4 * SQUISH_CONSTANT_4D; 309 if ((c & 0x03) != 0x03) { 310 if ((c & 0x03) == 0) { 311 zsv_ext0 += 1; 312 dz_ext0 -= 1; 313 } else { 314 zsv_ext1 += 1; 315 dz_ext1 -= 1; 316 } 317 } else { 318 zsv_ext2 += 1; 319 dz_ext2 -= 1; 320 } 321 } else { 322 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb; 323 dz_ext0 = dz_ext1 = dz_ext2 = dz0 - 4 * SQUISH_CONSTANT_4D; 324 } 325 326 if ((c & 0x08) != 0) { 327 wsv_ext0 = wsv_ext1 = wsb + 1; 328 wsv_ext2 = wsb + 2; 329 dw_ext0 = dw_ext1 = dw0 - 1 - 4 * SQUISH_CONSTANT_4D; 330 dw_ext2 = dw0 - 2 - 4 * SQUISH_CONSTANT_4D; 331 } else { 332 wsv_ext0 = wsv_ext1 = wsv_ext2 = wsb; 333 dw_ext0 = dw_ext1 = dw_ext2 = dw0 - 4 * SQUISH_CONSTANT_4D; 334 } 335 } else { //(1,1,1,1) is not one of the closest two pentachoron vertices. 336 byte c = cast(byte)(aPoint & bPoint); //Our three extra vertices are determined by the closest two. 337 338 if ((c & 0x01) != 0) { 339 xsv_ext0 = xsv_ext2 = xsb + 1; 340 xsv_ext1 = xsb + 2; 341 dx_ext0 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 342 dx_ext1 = dx0 - 2 - 3 * SQUISH_CONSTANT_4D; 343 dx_ext2 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; 344 } else { 345 xsv_ext0 = xsv_ext1 = xsv_ext2 = xsb; 346 dx_ext0 = dx0 - 2 * SQUISH_CONSTANT_4D; 347 dx_ext1 = dx_ext2 = dx0 - 3 * SQUISH_CONSTANT_4D; 348 } 349 350 if ((c & 0x02) != 0) { 351 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb + 1; 352 dy_ext0 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 353 dy_ext1 = dy_ext2 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; 354 if ((c & 0x01) != 0) { 355 ysv_ext2 += 1; 356 dy_ext2 -= 1; 357 } else { 358 ysv_ext1 += 1; 359 dy_ext1 -= 1; 360 } 361 } else { 362 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb; 363 dy_ext0 = dy0 - 2 * SQUISH_CONSTANT_4D; 364 dy_ext1 = dy_ext2 = dy0 - 3 * SQUISH_CONSTANT_4D; 365 } 366 367 if ((c & 0x04) != 0) { 368 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb + 1; 369 dz_ext0 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 370 dz_ext1 = dz_ext2 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; 371 if ((c & 0x03) != 0) { 372 zsv_ext2 += 1; 373 dz_ext2 -= 1; 374 } else { 375 zsv_ext1 += 1; 376 dz_ext1 -= 1; 377 } 378 } else { 379 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb; 380 dz_ext0 = dz0 - 2 * SQUISH_CONSTANT_4D; 381 dz_ext1 = dz_ext2 = dz0 - 3 * SQUISH_CONSTANT_4D; 382 } 383 384 if ((c & 0x08) != 0) { 385 wsv_ext0 = wsv_ext1 = wsb + 1; 386 wsv_ext2 = wsb + 2; 387 dw_ext0 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 388 dw_ext1 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; 389 dw_ext2 = dw0 - 2 - 3 * SQUISH_CONSTANT_4D; 390 } else { 391 wsv_ext0 = wsv_ext1 = wsv_ext2 = wsb; 392 dw_ext0 = dw0 - 2 * SQUISH_CONSTANT_4D; 393 dw_ext1 = dw_ext2 = dw0 - 3 * SQUISH_CONSTANT_4D; 394 } 395 } 396 397 //Contribution (1,1,1,0) 398 double dx4 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; 399 double dy4 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; 400 double dz4 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; 401 double dw4 = dw0 - 3 * SQUISH_CONSTANT_4D; 402 double attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4 - dw4 * dw4; 403 if (attn4 > 0) { 404 attn4 *= attn4; 405 value += attn4 * attn4 * extrapolate(xsb + 1, ysb + 1, zsb + 1, wsb + 0, dx4, dy4, dz4, dw4, perm); 406 } 407 408 //Contribution (1,1,0,1) 409 double dx3 = dx4; 410 double dy3 = dy4; 411 double dz3 = dz0 - 3 * SQUISH_CONSTANT_4D; 412 double dw3 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; 413 double attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3 - dw3 * dw3; 414 if (attn3 > 0) { 415 attn3 *= attn3; 416 value += attn3 * attn3 * extrapolate(xsb + 1, ysb + 1, zsb + 0, wsb + 1, dx3, dy3, dz3, dw3, perm); 417 } 418 419 //Contribution (1,0,1,1) 420 double dx2 = dx4; 421 double dy2 = dy0 - 3 * SQUISH_CONSTANT_4D; 422 double dz2 = dz4; 423 double dw2 = dw3; 424 double attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2 - dw2 * dw2; 425 if (attn2 > 0) { 426 attn2 *= attn2; 427 value += attn2 * attn2 * extrapolate(xsb + 1, ysb + 0, zsb + 1, wsb + 1, dx2, dy2, dz2, dw2, perm); 428 } 429 430 //Contribution (0,1,1,1) 431 double dx1 = dx0 - 3 * SQUISH_CONSTANT_4D; 432 double dz1 = dz4; 433 double dy1 = dy4; 434 double dw1 = dw3; 435 double attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1 - dw1 * dw1; 436 if (attn1 > 0) { 437 attn1 *= attn1; 438 value += attn1 * attn1 * extrapolate(xsb + 0, ysb + 1, zsb + 1, wsb + 1, dx1, dy1, dz1, dw1, perm); 439 } 440 441 //Contribution (1,1,1,1) 442 dx0 = dx0 - 1 - 4 * SQUISH_CONSTANT_4D; 443 dy0 = dy0 - 1 - 4 * SQUISH_CONSTANT_4D; 444 dz0 = dz0 - 1 - 4 * SQUISH_CONSTANT_4D; 445 dw0 = dw0 - 1 - 4 * SQUISH_CONSTANT_4D; 446 double attn0 = 2 - dx0 * dx0 - dy0 * dy0 - dz0 * dz0 - dw0 * dw0; 447 if (attn0 > 0) { 448 attn0 *= attn0; 449 value += attn0 * attn0 * extrapolate(xsb + 1, ysb + 1, zsb + 1, wsb + 1, dx0, dy0, dz0, dw0, perm); 450 } 451 } else if (inSum <= 2) { //We're inside the first dispentachoron (Rectified 4-Simplex) 452 double aScore; 453 byte aPoint; 454 bool aIsBiggerSide = true; 455 double bScore; 456 byte bPoint; 457 bool bIsBiggerSide = true; 458 459 //Decide between (1,1,0,0) and (0,0,1,1) 460 if (xins + yins > zins + wins) { 461 aScore = xins + yins; 462 aPoint = 0x03; 463 } else { 464 aScore = zins + wins; 465 aPoint = 0x0C; 466 } 467 468 //Decide between (1,0,1,0) and (0,1,0,1) 469 if (xins + zins > yins + wins) { 470 bScore = xins + zins; 471 bPoint = 0x05; 472 } else { 473 bScore = yins + wins; 474 bPoint = 0x0A; 475 } 476 477 //Closer between (1,0,0,1) and (0,1,1,0) will replace the further of a and b, if closer. 478 if (xins + wins > yins + zins) { 479 double score = xins + wins; 480 if (aScore >= bScore && score > bScore) { 481 bScore = score; 482 bPoint = 0x09; 483 } else if (aScore < bScore && score > aScore) { 484 aScore = score; 485 aPoint = 0x09; 486 } 487 } else { 488 double score = yins + zins; 489 if (aScore >= bScore && score > bScore) { 490 bScore = score; 491 bPoint = 0x06; 492 } else if (aScore < bScore && score > aScore) { 493 aScore = score; 494 aPoint = 0x06; 495 } 496 } 497 498 //Decide if (1,0,0,0) is closer. 499 double p1 = 2 - inSum + xins; 500 if (aScore >= bScore && p1 > bScore) { 501 bScore = p1; 502 bPoint = 0x01; 503 bIsBiggerSide = false; 504 } else if (aScore < bScore && p1 > aScore) { 505 aScore = p1; 506 aPoint = 0x01; 507 aIsBiggerSide = false; 508 } 509 510 //Decide if (0,1,0,0) is closer. 511 double p2 = 2 - inSum + yins; 512 if (aScore >= bScore && p2 > bScore) { 513 bScore = p2; 514 bPoint = 0x02; 515 bIsBiggerSide = false; 516 } else if (aScore < bScore && p2 > aScore) { 517 aScore = p2; 518 aPoint = 0x02; 519 aIsBiggerSide = false; 520 } 521 522 //Decide if (0,0,1,0) is closer. 523 double p3 = 2 - inSum + zins; 524 if (aScore >= bScore && p3 > bScore) { 525 bScore = p3; 526 bPoint = 0x04; 527 bIsBiggerSide = false; 528 } else if (aScore < bScore && p3 > aScore) { 529 aScore = p3; 530 aPoint = 0x04; 531 aIsBiggerSide = false; 532 } 533 534 //Decide if (0,0,0,1) is closer. 535 double p4 = 2 - inSum + wins; 536 if (aScore >= bScore && p4 > bScore) { 537 bScore = p4; 538 bPoint = 0x08; 539 bIsBiggerSide = false; 540 } else if (aScore < bScore && p4 > aScore) { 541 aScore = p4; 542 aPoint = 0x08; 543 aIsBiggerSide = false; 544 } 545 546 //Where each of the two closest points are determines how the extra three vertices are calculated. 547 if (aIsBiggerSide == bIsBiggerSide) { 548 if (aIsBiggerSide) { //Both closest points on the bigger side 549 byte c1 = cast(byte)(aPoint | bPoint); 550 byte c2 = cast(byte)(aPoint & bPoint); 551 if ((c1 & 0x01) == 0) { 552 xsv_ext0 = xsb; 553 xsv_ext1 = xsb - 1; 554 dx_ext0 = dx0 - 3 * SQUISH_CONSTANT_4D; 555 dx_ext1 = dx0 + 1 - 2 * SQUISH_CONSTANT_4D; 556 } else { 557 xsv_ext0 = xsv_ext1 = xsb + 1; 558 dx_ext0 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; 559 dx_ext1 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 560 } 561 562 if ((c1 & 0x02) == 0) { 563 ysv_ext0 = ysb; 564 ysv_ext1 = ysb - 1; 565 dy_ext0 = dy0 - 3 * SQUISH_CONSTANT_4D; 566 dy_ext1 = dy0 + 1 - 2 * SQUISH_CONSTANT_4D; 567 } else { 568 ysv_ext0 = ysv_ext1 = ysb + 1; 569 dy_ext0 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; 570 dy_ext1 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 571 } 572 573 if ((c1 & 0x04) == 0) { 574 zsv_ext0 = zsb; 575 zsv_ext1 = zsb - 1; 576 dz_ext0 = dz0 - 3 * SQUISH_CONSTANT_4D; 577 dz_ext1 = dz0 + 1 - 2 * SQUISH_CONSTANT_4D; 578 } else { 579 zsv_ext0 = zsv_ext1 = zsb + 1; 580 dz_ext0 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; 581 dz_ext1 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 582 } 583 584 if ((c1 & 0x08) == 0) { 585 wsv_ext0 = wsb; 586 wsv_ext1 = wsb - 1; 587 dw_ext0 = dw0 - 3 * SQUISH_CONSTANT_4D; 588 dw_ext1 = dw0 + 1 - 2 * SQUISH_CONSTANT_4D; 589 } else { 590 wsv_ext0 = wsv_ext1 = wsb + 1; 591 dw_ext0 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; 592 dw_ext1 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 593 } 594 595 //One combination is a permutation of (0,0,0,2) based on c2 596 xsv_ext2 = xsb; 597 ysv_ext2 = ysb; 598 zsv_ext2 = zsb; 599 wsv_ext2 = wsb; 600 dx_ext2 = dx0 - 2 * SQUISH_CONSTANT_4D; 601 dy_ext2 = dy0 - 2 * SQUISH_CONSTANT_4D; 602 dz_ext2 = dz0 - 2 * SQUISH_CONSTANT_4D; 603 dw_ext2 = dw0 - 2 * SQUISH_CONSTANT_4D; 604 if ((c2 & 0x01) != 0) { 605 xsv_ext2 += 2; 606 dx_ext2 -= 2; 607 } else if ((c2 & 0x02) != 0) { 608 ysv_ext2 += 2; 609 dy_ext2 -= 2; 610 } else if ((c2 & 0x04) != 0) { 611 zsv_ext2 += 2; 612 dz_ext2 -= 2; 613 } else { 614 wsv_ext2 += 2; 615 dw_ext2 -= 2; 616 } 617 618 } else { //Both closest points on the smaller side 619 //One of the two extra points is (0,0,0,0) 620 xsv_ext2 = xsb; 621 ysv_ext2 = ysb; 622 zsv_ext2 = zsb; 623 wsv_ext2 = wsb; 624 dx_ext2 = dx0; 625 dy_ext2 = dy0; 626 dz_ext2 = dz0; 627 dw_ext2 = dw0; 628 629 //Other two points are based on the omitted axes. 630 byte c = cast(byte)(aPoint | bPoint); 631 632 if ((c & 0x01) == 0) { 633 xsv_ext0 = xsb - 1; 634 xsv_ext1 = xsb; 635 dx_ext0 = dx0 + 1 - SQUISH_CONSTANT_4D; 636 dx_ext1 = dx0 - SQUISH_CONSTANT_4D; 637 } else { 638 xsv_ext0 = xsv_ext1 = xsb + 1; 639 dx_ext0 = dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_4D; 640 } 641 642 if ((c & 0x02) == 0) { 643 ysv_ext0 = ysv_ext1 = ysb; 644 dy_ext0 = dy_ext1 = dy0 - SQUISH_CONSTANT_4D; 645 if ((c & 0x01) == 0x01) 646 { 647 ysv_ext0 -= 1; 648 dy_ext0 += 1; 649 } else { 650 ysv_ext1 -= 1; 651 dy_ext1 += 1; 652 } 653 } else { 654 ysv_ext0 = ysv_ext1 = ysb + 1; 655 dy_ext0 = dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_4D; 656 } 657 658 if ((c & 0x04) == 0) { 659 zsv_ext0 = zsv_ext1 = zsb; 660 dz_ext0 = dz_ext1 = dz0 - SQUISH_CONSTANT_4D; 661 if ((c & 0x03) == 0x03) 662 { 663 zsv_ext0 -= 1; 664 dz_ext0 += 1; 665 } else { 666 zsv_ext1 -= 1; 667 dz_ext1 += 1; 668 } 669 } else { 670 zsv_ext0 = zsv_ext1 = zsb + 1; 671 dz_ext0 = dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_4D; 672 } 673 674 if ((c & 0x08) == 0) 675 { 676 wsv_ext0 = wsb; 677 wsv_ext1 = wsb - 1; 678 dw_ext0 = dw0 - SQUISH_CONSTANT_4D; 679 dw_ext1 = dw0 + 1 - SQUISH_CONSTANT_4D; 680 } else { 681 wsv_ext0 = wsv_ext1 = wsb + 1; 682 dw_ext0 = dw_ext1 = dw0 - 1 - SQUISH_CONSTANT_4D; 683 } 684 685 } 686 } else { //One point on each "side" 687 byte c1, c2; 688 if (aIsBiggerSide) { 689 c1 = aPoint; 690 c2 = bPoint; 691 } else { 692 c1 = bPoint; 693 c2 = aPoint; 694 } 695 696 //Two contributions are the bigger-sided point with each 0 replaced with -1. 697 if ((c1 & 0x01) == 0) { 698 xsv_ext0 = xsb - 1; 699 xsv_ext1 = xsb; 700 dx_ext0 = dx0 + 1 - SQUISH_CONSTANT_4D; 701 dx_ext1 = dx0 - SQUISH_CONSTANT_4D; 702 } else { 703 xsv_ext0 = xsv_ext1 = xsb + 1; 704 dx_ext0 = dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_4D; 705 } 706 707 if ((c1 & 0x02) == 0) { 708 ysv_ext0 = ysv_ext1 = ysb; 709 dy_ext0 = dy_ext1 = dy0 - SQUISH_CONSTANT_4D; 710 if ((c1 & 0x01) == 0x01) { 711 ysv_ext0 -= 1; 712 dy_ext0 += 1; 713 } else { 714 ysv_ext1 -= 1; 715 dy_ext1 += 1; 716 } 717 } else { 718 ysv_ext0 = ysv_ext1 = ysb + 1; 719 dy_ext0 = dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_4D; 720 } 721 722 if ((c1 & 0x04) == 0) { 723 zsv_ext0 = zsv_ext1 = zsb; 724 dz_ext0 = dz_ext1 = dz0 - SQUISH_CONSTANT_4D; 725 if ((c1 & 0x03) == 0x03) { 726 zsv_ext0 -= 1; 727 dz_ext0 += 1; 728 } else { 729 zsv_ext1 -= 1; 730 dz_ext1 += 1; 731 } 732 } else { 733 zsv_ext0 = zsv_ext1 = zsb + 1; 734 dz_ext0 = dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_4D; 735 } 736 737 if ((c1 & 0x08) == 0) { 738 wsv_ext0 = wsb; 739 wsv_ext1 = wsb - 1; 740 dw_ext0 = dw0 - SQUISH_CONSTANT_4D; 741 dw_ext1 = dw0 + 1 - SQUISH_CONSTANT_4D; 742 } else { 743 wsv_ext0 = wsv_ext1 = wsb + 1; 744 dw_ext0 = dw_ext1 = dw0 - 1 - SQUISH_CONSTANT_4D; 745 } 746 747 //One contribution is a permutation of (0,0,0,2) based on the smaller-sided point 748 xsv_ext2 = xsb; 749 ysv_ext2 = ysb; 750 zsv_ext2 = zsb; 751 wsv_ext2 = wsb; 752 dx_ext2 = dx0 - 2 * SQUISH_CONSTANT_4D; 753 dy_ext2 = dy0 - 2 * SQUISH_CONSTANT_4D; 754 dz_ext2 = dz0 - 2 * SQUISH_CONSTANT_4D; 755 dw_ext2 = dw0 - 2 * SQUISH_CONSTANT_4D; 756 if ((c2 & 0x01) != 0) { 757 xsv_ext2 += 2; 758 dx_ext2 -= 2; 759 } else if ((c2 & 0x02) != 0) { 760 ysv_ext2 += 2; 761 dy_ext2 -= 2; 762 } else if ((c2 & 0x04) != 0) { 763 zsv_ext2 += 2; 764 dz_ext2 -= 2; 765 } else { 766 wsv_ext2 += 2; 767 dw_ext2 -= 2; 768 } 769 } 770 771 //Contribution (1,0,0,0) 772 double dx1 = dx0 - 1 - SQUISH_CONSTANT_4D; 773 double dy1 = dy0 - 0 - SQUISH_CONSTANT_4D; 774 double dz1 = dz0 - 0 - SQUISH_CONSTANT_4D; 775 double dw1 = dw0 - 0 - SQUISH_CONSTANT_4D; 776 double attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1 - dw1 * dw1; 777 if (attn1 > 0) { 778 attn1 *= attn1; 779 value += attn1 * attn1 * extrapolate(xsb + 1, ysb + 0, zsb + 0, wsb + 0, dx1, dy1, dz1, dw1, perm); 780 } 781 782 //Contribution (0,1,0,0) 783 double dx2 = dx0 - 0 - SQUISH_CONSTANT_4D; 784 double dy2 = dy0 - 1 - SQUISH_CONSTANT_4D; 785 double dz2 = dz1; 786 double dw2 = dw1; 787 double attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2 - dw2 * dw2; 788 if (attn2 > 0) { 789 attn2 *= attn2; 790 value += attn2 * attn2 * extrapolate(xsb + 0, ysb + 1, zsb + 0, wsb + 0, dx2, dy2, dz2, dw2, perm); 791 } 792 793 //Contribution (0,0,1,0) 794 double dx3 = dx2; 795 double dy3 = dy1; 796 double dz3 = dz0 - 1 - SQUISH_CONSTANT_4D; 797 double dw3 = dw1; 798 double attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3 - dw3 * dw3; 799 if (attn3 > 0) { 800 attn3 *= attn3; 801 value += attn3 * attn3 * extrapolate(xsb + 0, ysb + 0, zsb + 1, wsb + 0, dx3, dy3, dz3, dw3, perm); 802 } 803 804 //Contribution (0,0,0,1) 805 double dx4 = dx2; 806 double dy4 = dy1; 807 double dz4 = dz1; 808 double dw4 = dw0 - 1 - SQUISH_CONSTANT_4D; 809 double attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4 - dw4 * dw4; 810 if (attn4 > 0) { 811 attn4 *= attn4; 812 value += attn4 * attn4 * extrapolate(xsb + 0, ysb + 0, zsb + 0, wsb + 1, dx4, dy4, dz4, dw4, perm); 813 } 814 815 //Contribution (1,1,0,0) 816 double dx5 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 817 double dy5 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 818 double dz5 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; 819 double dw5 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; 820 double attn5 = 2 - dx5 * dx5 - dy5 * dy5 - dz5 * dz5 - dw5 * dw5; 821 if (attn5 > 0) { 822 attn5 *= attn5; 823 value += attn5 * attn5 * extrapolate(xsb + 1, ysb + 1, zsb + 0, wsb + 0, dx5, dy5, dz5, dw5, perm); 824 } 825 826 //Contribution (1,0,1,0) 827 double dx6 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 828 double dy6 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; 829 double dz6 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 830 double dw6 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; 831 double attn6 = 2 - dx6 * dx6 - dy6 * dy6 - dz6 * dz6 - dw6 * dw6; 832 if (attn6 > 0) { 833 attn6 *= attn6; 834 value += attn6 * attn6 * extrapolate(xsb + 1, ysb + 0, zsb + 1, wsb + 0, dx6, dy6, dz6, dw6, perm); 835 } 836 837 //Contribution (1,0,0,1) 838 double dx7 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 839 double dy7 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; 840 double dz7 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; 841 double dw7 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 842 double attn7 = 2 - dx7 * dx7 - dy7 * dy7 - dz7 * dz7 - dw7 * dw7; 843 if (attn7 > 0) { 844 attn7 *= attn7; 845 value += attn7 * attn7 * extrapolate(xsb + 1, ysb + 0, zsb + 0, wsb + 1, dx7, dy7, dz7, dw7, perm); 846 } 847 848 //Contribution (0,1,1,0) 849 double dx8 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; 850 double dy8 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 851 double dz8 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 852 double dw8 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; 853 double attn8 = 2 - dx8 * dx8 - dy8 * dy8 - dz8 * dz8 - dw8 * dw8; 854 if (attn8 > 0) { 855 attn8 *= attn8; 856 value += attn8 * attn8 * extrapolate(xsb + 0, ysb + 1, zsb + 1, wsb + 0, dx8, dy8, dz8, dw8, perm); 857 } 858 859 //Contribution (0,1,0,1) 860 double dx9 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; 861 double dy9 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 862 double dz9 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; 863 double dw9 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 864 double attn9 = 2 - dx9 * dx9 - dy9 * dy9 - dz9 * dz9 - dw9 * dw9; 865 if (attn9 > 0) { 866 attn9 *= attn9; 867 value += attn9 * attn9 * extrapolate(xsb + 0, ysb + 1, zsb + 0, wsb + 1, dx9, dy9, dz9, dw9, perm); 868 } 869 870 //Contribution (0,0,1,1) 871 double dx10 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; 872 double dy10 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; 873 double dz10 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 874 double dw10 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 875 double attn10 = 2 - dx10 * dx10 - dy10 * dy10 - dz10 * dz10 - dw10 * dw10; 876 if (attn10 > 0) { 877 attn10 *= attn10; 878 value += attn10 * attn10 * extrapolate(xsb + 0, ysb + 0, zsb + 1, wsb + 1, dx10, dy10, dz10, dw10, perm); 879 } 880 } else { //We're inside the second dispentachoron (Rectified 4-Simplex) 881 double aScore; 882 byte aPoint; 883 bool aIsBiggerSide = true; 884 double bScore; 885 byte bPoint; 886 bool bIsBiggerSide = true; 887 888 //Decide between (0,0,1,1) and (1,1,0,0) 889 if (xins + yins < zins + wins) { 890 aScore = xins + yins; 891 aPoint = 0x0C; 892 } else { 893 aScore = zins + wins; 894 aPoint = 0x03; 895 } 896 897 //Decide between (0,1,0,1) and (1,0,1,0) 898 if (xins + zins < yins + wins) { 899 bScore = xins + zins; 900 bPoint = 0x0A; 901 } else { 902 bScore = yins + wins; 903 bPoint = 0x05; 904 } 905 906 //Closer between (0,1,1,0) and (1,0,0,1) will replace the further of a and b, if closer. 907 if (xins + wins < yins + zins) { 908 double score = xins + wins; 909 if (aScore <= bScore && score < bScore) { 910 bScore = score; 911 bPoint = 0x06; 912 } else if (aScore > bScore && score < aScore) { 913 aScore = score; 914 aPoint = 0x06; 915 } 916 } else { 917 double score = yins + zins; 918 if (aScore <= bScore && score < bScore) { 919 bScore = score; 920 bPoint = 0x09; 921 } else if (aScore > bScore && score < aScore) { 922 aScore = score; 923 aPoint = 0x09; 924 } 925 } 926 927 //Decide if (0,1,1,1) is closer. 928 double p1 = 3 - inSum + xins; 929 if (aScore <= bScore && p1 < bScore) { 930 bScore = p1; 931 bPoint = 0x0E; 932 bIsBiggerSide = false; 933 } else if (aScore > bScore && p1 < aScore) { 934 aScore = p1; 935 aPoint = 0x0E; 936 aIsBiggerSide = false; 937 } 938 939 //Decide if (1,0,1,1) is closer. 940 double p2 = 3 - inSum + yins; 941 if (aScore <= bScore && p2 < bScore) { 942 bScore = p2; 943 bPoint = 0x0D; 944 bIsBiggerSide = false; 945 } else if (aScore > bScore && p2 < aScore) { 946 aScore = p2; 947 aPoint = 0x0D; 948 aIsBiggerSide = false; 949 } 950 951 //Decide if (1,1,0,1) is closer. 952 double p3 = 3 - inSum + zins; 953 if (aScore <= bScore && p3 < bScore) { 954 bScore = p3; 955 bPoint = 0x0B; 956 bIsBiggerSide = false; 957 } else if (aScore > bScore && p3 < aScore) { 958 aScore = p3; 959 aPoint = 0x0B; 960 aIsBiggerSide = false; 961 } 962 963 //Decide if (1,1,1,0) is closer. 964 double p4 = 3 - inSum + wins; 965 if (aScore <= bScore && p4 < bScore) { 966 bScore = p4; 967 bPoint = 0x07; 968 bIsBiggerSide = false; 969 } else if (aScore > bScore && p4 < aScore) { 970 aScore = p4; 971 aPoint = 0x07; 972 aIsBiggerSide = false; 973 } 974 975 //Where each of the two closest points are determines how the extra three vertices are calculated. 976 if (aIsBiggerSide == bIsBiggerSide) { 977 if (aIsBiggerSide) { //Both closest points on the bigger side 978 byte c1 = cast(byte)(aPoint & bPoint); 979 byte c2 = cast(byte)(aPoint | bPoint); 980 981 //Two contributions are permutations of (0,0,0,1) and (0,0,0,2) based on c1 982 xsv_ext0 = xsv_ext1 = xsb; 983 ysv_ext0 = ysv_ext1 = ysb; 984 zsv_ext0 = zsv_ext1 = zsb; 985 wsv_ext0 = wsv_ext1 = wsb; 986 dx_ext0 = dx0 - SQUISH_CONSTANT_4D; 987 dy_ext0 = dy0 - SQUISH_CONSTANT_4D; 988 dz_ext0 = dz0 - SQUISH_CONSTANT_4D; 989 dw_ext0 = dw0 - SQUISH_CONSTANT_4D; 990 dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_4D; 991 dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_4D; 992 dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_4D; 993 dw_ext1 = dw0 - 2 * SQUISH_CONSTANT_4D; 994 if ((c1 & 0x01) != 0) { 995 xsv_ext0 += 1; 996 dx_ext0 -= 1; 997 xsv_ext1 += 2; 998 dx_ext1 -= 2; 999 } else if ((c1 & 0x02) != 0) { 1000 ysv_ext0 += 1; 1001 dy_ext0 -= 1; 1002 ysv_ext1 += 2; 1003 dy_ext1 -= 2; 1004 } else if ((c1 & 0x04) != 0) { 1005 zsv_ext0 += 1; 1006 dz_ext0 -= 1; 1007 zsv_ext1 += 2; 1008 dz_ext1 -= 2; 1009 } else { 1010 wsv_ext0 += 1; 1011 dw_ext0 -= 1; 1012 wsv_ext1 += 2; 1013 dw_ext1 -= 2; 1014 } 1015 1016 //One contribution is a permutation of (1,1,1,-1) based on c2 1017 xsv_ext2 = xsb + 1; 1018 ysv_ext2 = ysb + 1; 1019 zsv_ext2 = zsb + 1; 1020 wsv_ext2 = wsb + 1; 1021 dx_ext2 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1022 dy_ext2 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1023 dz_ext2 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1024 dw_ext2 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1025 if ((c2 & 0x01) == 0) { 1026 xsv_ext2 -= 2; 1027 dx_ext2 += 2; 1028 } else if ((c2 & 0x02) == 0) { 1029 ysv_ext2 -= 2; 1030 dy_ext2 += 2; 1031 } else if ((c2 & 0x04) == 0) { 1032 zsv_ext2 -= 2; 1033 dz_ext2 += 2; 1034 } else { 1035 wsv_ext2 -= 2; 1036 dw_ext2 += 2; 1037 } 1038 } else { //Both closest points on the smaller side 1039 //One of the two extra points is (1,1,1,1) 1040 xsv_ext2 = xsb + 1; 1041 ysv_ext2 = ysb + 1; 1042 zsv_ext2 = zsb + 1; 1043 wsv_ext2 = wsb + 1; 1044 dx_ext2 = dx0 - 1 - 4 * SQUISH_CONSTANT_4D; 1045 dy_ext2 = dy0 - 1 - 4 * SQUISH_CONSTANT_4D; 1046 dz_ext2 = dz0 - 1 - 4 * SQUISH_CONSTANT_4D; 1047 dw_ext2 = dw0 - 1 - 4 * SQUISH_CONSTANT_4D; 1048 1049 //Other two points are based on the shared axes. 1050 byte c = cast(byte)(aPoint & bPoint); 1051 1052 if ((c & 0x01) != 0) { 1053 xsv_ext0 = xsb + 2; 1054 xsv_ext1 = xsb + 1; 1055 dx_ext0 = dx0 - 2 - 3 * SQUISH_CONSTANT_4D; 1056 dx_ext1 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; 1057 } else { 1058 xsv_ext0 = xsv_ext1 = xsb; 1059 dx_ext0 = dx_ext1 = dx0 - 3 * SQUISH_CONSTANT_4D; 1060 } 1061 1062 if ((c & 0x02) != 0) { 1063 ysv_ext0 = ysv_ext1 = ysb + 1; 1064 dy_ext0 = dy_ext1 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; 1065 if ((c & 0x01) == 0) 1066 { 1067 ysv_ext0 += 1; 1068 dy_ext0 -= 1; 1069 } else { 1070 ysv_ext1 += 1; 1071 dy_ext1 -= 1; 1072 } 1073 } else { 1074 ysv_ext0 = ysv_ext1 = ysb; 1075 dy_ext0 = dy_ext1 = dy0 - 3 * SQUISH_CONSTANT_4D; 1076 } 1077 1078 if ((c & 0x04) != 0) { 1079 zsv_ext0 = zsv_ext1 = zsb + 1; 1080 dz_ext0 = dz_ext1 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; 1081 if ((c & 0x03) == 0) 1082 { 1083 zsv_ext0 += 1; 1084 dz_ext0 -= 1; 1085 } else { 1086 zsv_ext1 += 1; 1087 dz_ext1 -= 1; 1088 } 1089 } else { 1090 zsv_ext0 = zsv_ext1 = zsb; 1091 dz_ext0 = dz_ext1 = dz0 - 3 * SQUISH_CONSTANT_4D; 1092 } 1093 1094 if ((c & 0x08) != 0) 1095 { 1096 wsv_ext0 = wsb + 1; 1097 wsv_ext1 = wsb + 2; 1098 dw_ext0 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; 1099 dw_ext1 = dw0 - 2 - 3 * SQUISH_CONSTANT_4D; 1100 } else { 1101 wsv_ext0 = wsv_ext1 = wsb; 1102 dw_ext0 = dw_ext1 = dw0 - 3 * SQUISH_CONSTANT_4D; 1103 } 1104 } 1105 } else { //One point on each "side" 1106 byte c1, c2; 1107 if (aIsBiggerSide) { 1108 c1 = aPoint; 1109 c2 = bPoint; 1110 } else { 1111 c1 = bPoint; 1112 c2 = aPoint; 1113 } 1114 1115 //Two contributions are the bigger-sided point with each 1 replaced with 2. 1116 if ((c1 & 0x01) != 0) { 1117 xsv_ext0 = xsb + 2; 1118 xsv_ext1 = xsb + 1; 1119 dx_ext0 = dx0 - 2 - 3 * SQUISH_CONSTANT_4D; 1120 dx_ext1 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; 1121 } else { 1122 xsv_ext0 = xsv_ext1 = xsb; 1123 dx_ext0 = dx_ext1 = dx0 - 3 * SQUISH_CONSTANT_4D; 1124 } 1125 1126 if ((c1 & 0x02) != 0) { 1127 ysv_ext0 = ysv_ext1 = ysb + 1; 1128 dy_ext0 = dy_ext1 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; 1129 if ((c1 & 0x01) == 0) { 1130 ysv_ext0 += 1; 1131 dy_ext0 -= 1; 1132 } else { 1133 ysv_ext1 += 1; 1134 dy_ext1 -= 1; 1135 } 1136 } else { 1137 ysv_ext0 = ysv_ext1 = ysb; 1138 dy_ext0 = dy_ext1 = dy0 - 3 * SQUISH_CONSTANT_4D; 1139 } 1140 1141 if ((c1 & 0x04) != 0) { 1142 zsv_ext0 = zsv_ext1 = zsb + 1; 1143 dz_ext0 = dz_ext1 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; 1144 if ((c1 & 0x03) == 0) { 1145 zsv_ext0 += 1; 1146 dz_ext0 -= 1; 1147 } else { 1148 zsv_ext1 += 1; 1149 dz_ext1 -= 1; 1150 } 1151 } else { 1152 zsv_ext0 = zsv_ext1 = zsb; 1153 dz_ext0 = dz_ext1 = dz0 - 3 * SQUISH_CONSTANT_4D; 1154 } 1155 1156 if ((c1 & 0x08) != 0) { 1157 wsv_ext0 = wsb + 1; 1158 wsv_ext1 = wsb + 2; 1159 dw_ext0 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; 1160 dw_ext1 = dw0 - 2 - 3 * SQUISH_CONSTANT_4D; 1161 } else { 1162 wsv_ext0 = wsv_ext1 = wsb; 1163 dw_ext0 = dw_ext1 = dw0 - 3 * SQUISH_CONSTANT_4D; 1164 } 1165 1166 //One contribution is a permutation of (1,1,1,-1) based on the smaller-sided point 1167 xsv_ext2 = xsb + 1; 1168 ysv_ext2 = ysb + 1; 1169 zsv_ext2 = zsb + 1; 1170 wsv_ext2 = wsb + 1; 1171 dx_ext2 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1172 dy_ext2 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1173 dz_ext2 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1174 dw_ext2 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1175 if ((c2 & 0x01) == 0) { 1176 xsv_ext2 -= 2; 1177 dx_ext2 += 2; 1178 } else if ((c2 & 0x02) == 0) { 1179 ysv_ext2 -= 2; 1180 dy_ext2 += 2; 1181 } else if ((c2 & 0x04) == 0) { 1182 zsv_ext2 -= 2; 1183 dz_ext2 += 2; 1184 } else { 1185 wsv_ext2 -= 2; 1186 dw_ext2 += 2; 1187 } 1188 } 1189 1190 //Contribution (1,1,1,0) 1191 double dx4 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; 1192 double dy4 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; 1193 double dz4 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; 1194 double dw4 = dw0 - 3 * SQUISH_CONSTANT_4D; 1195 double attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4 - dw4 * dw4; 1196 if (attn4 > 0) { 1197 attn4 *= attn4; 1198 value += attn4 * attn4 * extrapolate(xsb + 1, ysb + 1, zsb + 1, wsb + 0, dx4, dy4, dz4, dw4, perm); 1199 } 1200 1201 //Contribution (1,1,0,1) 1202 double dx3 = dx4; 1203 double dy3 = dy4; 1204 double dz3 = dz0 - 3 * SQUISH_CONSTANT_4D; 1205 double dw3 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; 1206 double attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3 - dw3 * dw3; 1207 if (attn3 > 0) { 1208 attn3 *= attn3; 1209 value += attn3 * attn3 * extrapolate(xsb + 1, ysb + 1, zsb + 0, wsb + 1, dx3, dy3, dz3, dw3, perm); 1210 } 1211 1212 //Contribution (1,0,1,1) 1213 double dx2 = dx4; 1214 double dy2 = dy0 - 3 * SQUISH_CONSTANT_4D; 1215 double dz2 = dz4; 1216 double dw2 = dw3; 1217 double attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2 - dw2 * dw2; 1218 if (attn2 > 0) { 1219 attn2 *= attn2; 1220 value += attn2 * attn2 * extrapolate(xsb + 1, ysb + 0, zsb + 1, wsb + 1, dx2, dy2, dz2, dw2, perm); 1221 } 1222 1223 //Contribution (0,1,1,1) 1224 double dx1 = dx0 - 3 * SQUISH_CONSTANT_4D; 1225 double dz1 = dz4; 1226 double dy1 = dy4; 1227 double dw1 = dw3; 1228 double attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1 - dw1 * dw1; 1229 if (attn1 > 0) { 1230 attn1 *= attn1; 1231 value += attn1 * attn1 * extrapolate(xsb + 0, ysb + 1, zsb + 1, wsb + 1, dx1, dy1, dz1, dw1, perm); 1232 } 1233 1234 //Contribution (1,1,0,0) 1235 double dx5 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1236 double dy5 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1237 double dz5 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; 1238 double dw5 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; 1239 double attn5 = 2 - dx5 * dx5 - dy5 * dy5 - dz5 * dz5 - dw5 * dw5; 1240 if (attn5 > 0) { 1241 attn5 *= attn5; 1242 value += attn5 * attn5 * extrapolate(xsb + 1, ysb + 1, zsb + 0, wsb + 0, dx5, dy5, dz5, dw5, perm); 1243 } 1244 1245 //Contribution (1,0,1,0) 1246 double dx6 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1247 double dy6 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; 1248 double dz6 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1249 double dw6 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; 1250 double attn6 = 2 - dx6 * dx6 - dy6 * dy6 - dz6 * dz6 - dw6 * dw6; 1251 if (attn6 > 0) { 1252 attn6 *= attn6; 1253 value += attn6 * attn6 * extrapolate(xsb + 1, ysb + 0, zsb + 1, wsb + 0, dx6, dy6, dz6, dw6, perm); 1254 } 1255 1256 //Contribution (1,0,0,1) 1257 double dx7 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1258 double dy7 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; 1259 double dz7 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; 1260 double dw7 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1261 double attn7 = 2 - dx7 * dx7 - dy7 * dy7 - dz7 * dz7 - dw7 * dw7; 1262 if (attn7 > 0) { 1263 attn7 *= attn7; 1264 value += attn7 * attn7 * extrapolate(xsb + 1, ysb + 0, zsb + 0, wsb + 1, dx7, dy7, dz7, dw7, perm); 1265 } 1266 1267 //Contribution (0,1,1,0) 1268 double dx8 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; 1269 double dy8 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1270 double dz8 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1271 double dw8 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; 1272 double attn8 = 2 - dx8 * dx8 - dy8 * dy8 - dz8 * dz8 - dw8 * dw8; 1273 if (attn8 > 0) { 1274 attn8 *= attn8; 1275 value += attn8 * attn8 * extrapolate(xsb + 0, ysb + 1, zsb + 1, wsb + 0, dx8, dy8, dz8, dw8, perm); 1276 } 1277 1278 //Contribution (0,1,0,1) 1279 double dx9 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; 1280 double dy9 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1281 double dz9 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; 1282 double dw9 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1283 double attn9 = 2 - dx9 * dx9 - dy9 * dy9 - dz9 * dz9 - dw9 * dw9; 1284 if (attn9 > 0) { 1285 attn9 *= attn9; 1286 value += attn9 * attn9 * extrapolate(xsb + 0, ysb + 1, zsb + 0, wsb + 1, dx9, dy9, dz9, dw9, perm); 1287 } 1288 1289 //Contribution (0,0,1,1) 1290 double dx10 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; 1291 double dy10 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; 1292 double dz10 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1293 double dw10 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1294 double attn10 = 2 - dx10 * dx10 - dy10 * dy10 - dz10 * dz10 - dw10 * dw10; 1295 if (attn10 > 0) { 1296 attn10 *= attn10; 1297 value += attn10 * attn10 * extrapolate(xsb + 0, ysb + 0, zsb + 1, wsb + 1, dx10, dy10, dz10, dw10, perm); 1298 } 1299 } 1300 1301 //First extra vertex 1302 double attn_ext0 = 2 - dx_ext0 * dx_ext0 - dy_ext0 * dy_ext0 - dz_ext0 * dz_ext0 - dw_ext0 * dw_ext0; 1303 if (attn_ext0 > 0) 1304 { 1305 attn_ext0 *= attn_ext0; 1306 value += attn_ext0 * attn_ext0 * extrapolate(xsv_ext0, ysv_ext0, zsv_ext0, wsv_ext0, dx_ext0, dy_ext0, dz_ext0, dw_ext0, perm); 1307 } 1308 1309 //Second extra vertex 1310 double attn_ext1 = 2 - dx_ext1 * dx_ext1 - dy_ext1 * dy_ext1 - dz_ext1 * dz_ext1 - dw_ext1 * dw_ext1; 1311 if (attn_ext1 > 0) 1312 { 1313 attn_ext1 *= attn_ext1; 1314 value += attn_ext1 * attn_ext1 * extrapolate(xsv_ext1, ysv_ext1, zsv_ext1, wsv_ext1, dx_ext1, dy_ext1, dz_ext1, dw_ext1, perm); 1315 } 1316 1317 //Third extra vertex 1318 double attn_ext2 = 2 - dx_ext2 * dx_ext2 - dy_ext2 * dy_ext2 - dz_ext2 * dz_ext2 - dw_ext2 * dw_ext2; 1319 if (attn_ext2 > 0) 1320 { 1321 attn_ext2 *= attn_ext2; 1322 value += attn_ext2 * attn_ext2 * extrapolate(xsv_ext2, ysv_ext2, zsv_ext2, wsv_ext2, dx_ext2, dy_ext2, dz_ext2, dw_ext2, perm); 1323 } 1324 1325 return value / NORM_CONSTANT_4D; 1326 } 1327 1328 /** 1329 extrapolates for 4D 1330 */ 1331 private @nogc @safe pure double extrapolate(int xsb, int ysb, int zsb, int wsb, double dx, double dy, double dz, double dw, const ref short[256] perm) 1332 { 1333 int index = perm[(perm[(perm[(perm[xsb & 0xFF] + ysb) & 0xFF] + zsb) & 0xFF] + wsb) & 0xFF] & 0xFC; 1334 return GRADIENTS.GRADIENTS_4D[index] * dx 1335 + GRADIENTS.GRADIENTS_4D[index + 1] * dy 1336 + GRADIENTS.GRADIENTS_4D[index + 2] * dz 1337 + GRADIENTS.GRADIENTS_4D[index + 3] * dw; 1338 }