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.osimplex3d; 5 6 import dosimplex.util; 7 8 /// "Config" enum for the noise generator 9 private enum : double { 10 STRETCH_CONSTANT_3D = -1.0 / 6, ///(1/Math.sqrt(3+1)-1)/3; 11 SQUISH_CONSTANT_3D = 1.0 / 3, ///(Math.sqrt(3+1)-1)/3; 12 NORM_CONSTANT_3D = 103., 13 } 14 15 /** 16 3D OpenSimplexNoise implementation. 17 */ 18 @nogc @safe pure double osNoise3D(double x, double y, double z, const ref short[256] perm, const ref short[256] permGradIndex3D) 19 { 20 21 //Place input coordinates on simplectic honeycomb. 22 double stretchOffset = (x + y + z) * STRETCH_CONSTANT_3D; 23 double xs = x + stretchOffset; 24 double ys = y + stretchOffset; 25 double zs = z + stretchOffset; 26 27 //Floor to get simplectic honeycomb coordinates of rhombohedron (stretched cube) super-cell origin. 28 int xsb = fastFloor(xs); 29 int ysb = fastFloor(ys); 30 int zsb = fastFloor(zs); 31 32 //Skew out to get actual coordinates of rhombohedron origin. We'll need these later. 33 double squishOffset = (xsb + ysb + zsb) * SQUISH_CONSTANT_3D; 34 double xb = xsb + squishOffset; 35 double yb = ysb + squishOffset; 36 double zb = zsb + squishOffset; 37 38 //Compute simplectic honeycomb coordinates relative to rhombohedral origin. 39 double xins = xs - xsb; 40 double yins = ys - ysb; 41 double zins = zs - zsb; 42 43 //Sum those together to get a value that determines which region we're in. 44 double inSum = xins + yins + zins; 45 46 //Positions relative to origin point. 47 double dx0 = x - xb; 48 double dy0 = y - yb; 49 double dz0 = z - zb; 50 51 //We'll be defining these inside the next block and using them afterwards. 52 double dx_ext0, dy_ext0, dz_ext0; 53 double dx_ext1, dy_ext1, dz_ext1; 54 int xsv_ext0, ysv_ext0, zsv_ext0; 55 int xsv_ext1, ysv_ext1, zsv_ext1; 56 57 double value = 0; 58 if (inSum <= 1) { //We're inside the tetrahedron (3-Simplex) at (0,0,0) 59 60 //Determine which two of (0,0,1), (0,1,0), (1,0,0) are closest. 61 byte aPoint = 0x01; 62 double aScore = xins; 63 byte bPoint = 0x02; 64 double bScore = yins; 65 if (aScore >= bScore && zins > bScore) { 66 bScore = zins; 67 bPoint = 0x04; 68 } else if (aScore < bScore && zins > aScore) { 69 aScore = zins; 70 aPoint = 0x04; 71 } 72 73 //Now we determine the two lattice points not part of the tetrahedron that may contribute. 74 //This depends on the closest two tetrahedral vertices, including (0,0,0) 75 double wins = 1 - inSum; 76 if (wins > aScore || wins > bScore) { //(0,0,0) is one of the closest two tetrahedral vertices. 77 byte c = (bScore > aScore ? bPoint : aPoint); //Our other closest vertex is the closest out of a and b. 78 79 if ((c & 0x01) == 0) { 80 xsv_ext0 = xsb - 1; 81 xsv_ext1 = xsb; 82 dx_ext0 = dx0 + 1; 83 dx_ext1 = dx0; 84 } else { 85 xsv_ext0 = xsv_ext1 = xsb + 1; 86 dx_ext0 = dx_ext1 = dx0 - 1; 87 } 88 89 if ((c & 0x02) == 0) { 90 ysv_ext0 = ysv_ext1 = ysb; 91 dy_ext0 = dy_ext1 = dy0; 92 if ((c & 0x01) == 0) { 93 ysv_ext1 -= 1; 94 dy_ext1 += 1; 95 } else { 96 ysv_ext0 -= 1; 97 dy_ext0 += 1; 98 } 99 } else { 100 ysv_ext0 = ysv_ext1 = ysb + 1; 101 dy_ext0 = dy_ext1 = dy0 - 1; 102 } 103 104 if ((c & 0x04) == 0) { 105 zsv_ext0 = zsb; 106 zsv_ext1 = zsb - 1; 107 dz_ext0 = dz0; 108 dz_ext1 = dz0 + 1; 109 } else { 110 zsv_ext0 = zsv_ext1 = zsb + 1; 111 dz_ext0 = dz_ext1 = dz0 - 1; 112 } 113 } else { //(0,0,0) is not one of the closest two tetrahedral vertices. 114 byte c = cast(byte)(aPoint | bPoint); //Our two extra vertices are determined by the closest two. 115 116 if ((c & 0x01) == 0) { 117 xsv_ext0 = xsb; 118 xsv_ext1 = xsb - 1; 119 dx_ext0 = dx0 - 2 * SQUISH_CONSTANT_3D; 120 dx_ext1 = dx0 + 1 - SQUISH_CONSTANT_3D; 121 } else { 122 xsv_ext0 = xsv_ext1 = xsb + 1; 123 dx_ext0 = dx0 - 1 - 2 * SQUISH_CONSTANT_3D; 124 dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_3D; 125 } 126 127 if ((c & 0x02) == 0) { 128 ysv_ext0 = ysb; 129 ysv_ext1 = ysb - 1; 130 dy_ext0 = dy0 - 2 * SQUISH_CONSTANT_3D; 131 dy_ext1 = dy0 + 1 - SQUISH_CONSTANT_3D; 132 } else { 133 ysv_ext0 = ysv_ext1 = ysb + 1; 134 dy_ext0 = dy0 - 1 - 2 * SQUISH_CONSTANT_3D; 135 dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_3D; 136 } 137 138 if ((c & 0x04) == 0) { 139 zsv_ext0 = zsb; 140 zsv_ext1 = zsb - 1; 141 dz_ext0 = dz0 - 2 * SQUISH_CONSTANT_3D; 142 dz_ext1 = dz0 + 1 - SQUISH_CONSTANT_3D; 143 } else { 144 zsv_ext0 = zsv_ext1 = zsb + 1; 145 dz_ext0 = dz0 - 1 - 2 * SQUISH_CONSTANT_3D; 146 dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_3D; 147 } 148 } 149 150 //Contribution (0,0,0) 151 double attn0 = 2 - dx0 * dx0 - dy0 * dy0 - dz0 * dz0; 152 if (attn0 > 0) { 153 attn0 *= attn0; 154 value += attn0 * attn0 * extrapolate(xsb + 0, ysb + 0, zsb + 0, dx0, dy0, dz0, perm, permGradIndex3D); 155 } 156 157 //Contribution (1,0,0) 158 double dx1 = dx0 - 1 - SQUISH_CONSTANT_3D; 159 double dy1 = dy0 - 0 - SQUISH_CONSTANT_3D; 160 double dz1 = dz0 - 0 - SQUISH_CONSTANT_3D; 161 double attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1; 162 if (attn1 > 0) { 163 attn1 *= attn1; 164 value += attn1 * attn1 * extrapolate(xsb + 1, ysb + 0, zsb + 0, dx1, dy1, dz1, perm, permGradIndex3D); 165 } 166 167 //Contribution (0,1,0) 168 double dx2 = dx0 - 0 - SQUISH_CONSTANT_3D; 169 double dy2 = dy0 - 1 - SQUISH_CONSTANT_3D; 170 double dz2 = dz1; 171 double attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2; 172 if (attn2 > 0) { 173 attn2 *= attn2; 174 value += attn2 * attn2 * extrapolate(xsb + 0, ysb + 1, zsb + 0, dx2, dy2, dz2, perm, permGradIndex3D); 175 } 176 177 //Contribution (0,0,1) 178 double dx3 = dx2; 179 double dy3 = dy1; 180 double dz3 = dz0 - 1 - SQUISH_CONSTANT_3D; 181 double attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3; 182 if (attn3 > 0) { 183 attn3 *= attn3; 184 value += attn3 * attn3 * extrapolate(xsb + 0, ysb + 0, zsb + 1, dx3, dy3, dz3, perm, permGradIndex3D); 185 } 186 } else if (inSum >= 2) { //We're inside the tetrahedron (3-Simplex) at (1,1,1) 187 188 //Determine which two tetrahedral vertices are the closest, out of (1,1,0), (1,0,1), (0,1,1) but not (1,1,1). 189 byte aPoint = 0x06; 190 double aScore = xins; 191 byte bPoint = 0x05; 192 double bScore = yins; 193 if (aScore <= bScore && zins < bScore) { 194 bScore = zins; 195 bPoint = 0x03; 196 } else if (aScore > bScore && zins < aScore) { 197 aScore = zins; 198 aPoint = 0x03; 199 } 200 201 //Now we determine the two lattice points not part of the tetrahedron that may contribute. 202 //This depends on the closest two tetrahedral vertices, including (1,1,1) 203 double wins = 3 - inSum; 204 if (wins < aScore || wins < bScore) { //(1,1,1) is one of the closest two tetrahedral vertices. 205 byte c = (bScore < aScore ? bPoint : aPoint); //Our other closest vertex is the closest out of a and b. 206 207 if ((c & 0x01) != 0) { 208 xsv_ext0 = xsb + 2; 209 xsv_ext1 = xsb + 1; 210 dx_ext0 = dx0 - 2 - 3 * SQUISH_CONSTANT_3D; 211 dx_ext1 = dx0 - 1 - 3 * SQUISH_CONSTANT_3D; 212 } else { 213 xsv_ext0 = xsv_ext1 = xsb; 214 dx_ext0 = dx_ext1 = dx0 - 3 * SQUISH_CONSTANT_3D; 215 } 216 217 if ((c & 0x02) != 0) { 218 ysv_ext0 = ysv_ext1 = ysb + 1; 219 dy_ext0 = dy_ext1 = dy0 - 1 - 3 * SQUISH_CONSTANT_3D; 220 if ((c & 0x01) != 0) { 221 ysv_ext1 += 1; 222 dy_ext1 -= 1; 223 } else { 224 ysv_ext0 += 1; 225 dy_ext0 -= 1; 226 } 227 } else { 228 ysv_ext0 = ysv_ext1 = ysb; 229 dy_ext0 = dy_ext1 = dy0 - 3 * SQUISH_CONSTANT_3D; 230 } 231 232 if ((c & 0x04) != 0) { 233 zsv_ext0 = zsb + 1; 234 zsv_ext1 = zsb + 2; 235 dz_ext0 = dz0 - 1 - 3 * SQUISH_CONSTANT_3D; 236 dz_ext1 = dz0 - 2 - 3 * SQUISH_CONSTANT_3D; 237 } else { 238 zsv_ext0 = zsv_ext1 = zsb; 239 dz_ext0 = dz_ext1 = dz0 - 3 * SQUISH_CONSTANT_3D; 240 } 241 } else { //(1,1,1) is not one of the closest two tetrahedral vertices. 242 byte c = cast(byte)(aPoint & bPoint); //Our two extra vertices are determined by the closest two. 243 244 if ((c & 0x01) != 0) { 245 xsv_ext0 = xsb + 1; 246 xsv_ext1 = xsb + 2; 247 dx_ext0 = dx0 - 1 - SQUISH_CONSTANT_3D; 248 dx_ext1 = dx0 - 2 - 2 * SQUISH_CONSTANT_3D; 249 } else { 250 xsv_ext0 = xsv_ext1 = xsb; 251 dx_ext0 = dx0 - SQUISH_CONSTANT_3D; 252 dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D; 253 } 254 255 if ((c & 0x02) != 0) { 256 ysv_ext0 = ysb + 1; 257 ysv_ext1 = ysb + 2; 258 dy_ext0 = dy0 - 1 - SQUISH_CONSTANT_3D; 259 dy_ext1 = dy0 - 2 - 2 * SQUISH_CONSTANT_3D; 260 } else { 261 ysv_ext0 = ysv_ext1 = ysb; 262 dy_ext0 = dy0 - SQUISH_CONSTANT_3D; 263 dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D; 264 } 265 266 if ((c & 0x04) != 0) { 267 zsv_ext0 = zsb + 1; 268 zsv_ext1 = zsb + 2; 269 dz_ext0 = dz0 - 1 - SQUISH_CONSTANT_3D; 270 dz_ext1 = dz0 - 2 - 2 * SQUISH_CONSTANT_3D; 271 } else { 272 zsv_ext0 = zsv_ext1 = zsb; 273 dz_ext0 = dz0 - SQUISH_CONSTANT_3D; 274 dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D; 275 } 276 } 277 278 //Contribution (1,1,0) 279 double dx3 = dx0 - 1 - 2 * SQUISH_CONSTANT_3D; 280 double dy3 = dy0 - 1 - 2 * SQUISH_CONSTANT_3D; 281 double dz3 = dz0 - 0 - 2 * SQUISH_CONSTANT_3D; 282 double attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3; 283 if (attn3 > 0) { 284 attn3 *= attn3; 285 value += attn3 * attn3 * extrapolate(xsb + 1, ysb + 1, zsb + 0, dx3, dy3, dz3, perm, permGradIndex3D); 286 } 287 288 //Contribution (1,0,1) 289 double dx2 = dx3; 290 double dy2 = dy0 - 0 - 2 * SQUISH_CONSTANT_3D; 291 double dz2 = dz0 - 1 - 2 * SQUISH_CONSTANT_3D; 292 double attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2; 293 if (attn2 > 0) { 294 attn2 *= attn2; 295 value += attn2 * attn2 * extrapolate(xsb + 1, ysb + 0, zsb + 1, dx2, dy2, dz2, perm, permGradIndex3D); 296 } 297 298 //Contribution (0,1,1) 299 double dx1 = dx0 - 0 - 2 * SQUISH_CONSTANT_3D; 300 double dy1 = dy3; 301 double dz1 = dz2; 302 double attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1; 303 if (attn1 > 0) { 304 attn1 *= attn1; 305 value += attn1 * attn1 * extrapolate(xsb + 0, ysb + 1, zsb + 1, dx1, dy1, dz1, perm, permGradIndex3D); 306 } 307 308 //Contribution (1,1,1) 309 dx0 = dx0 - 1 - 3 * SQUISH_CONSTANT_3D; 310 dy0 = dy0 - 1 - 3 * SQUISH_CONSTANT_3D; 311 dz0 = dz0 - 1 - 3 * SQUISH_CONSTANT_3D; 312 double attn0 = 2 - dx0 * dx0 - dy0 * dy0 - dz0 * dz0; 313 if (attn0 > 0) { 314 attn0 *= attn0; 315 value += attn0 * attn0 * extrapolate(xsb + 1, ysb + 1, zsb + 1, dx0, dy0, dz0, perm, permGradIndex3D); 316 } 317 } else { //We're inside the octahedron (Rectified 3-Simplex) in between. 318 double aScore; 319 byte aPoint; 320 bool aIsFurtherSide; 321 double bScore; 322 byte bPoint; 323 bool bIsFurtherSide; 324 325 //Decide between point (0,0,1) and (1,1,0) as closest 326 double p1 = xins + yins; 327 if (p1 > 1) { 328 aScore = p1 - 1; 329 aPoint = 0x03; 330 aIsFurtherSide = true; 331 } else { 332 aScore = 1 - p1; 333 aPoint = 0x04; 334 aIsFurtherSide = false; 335 } 336 337 //Decide between point (0,1,0) and (1,0,1) as closest 338 double p2 = xins + zins; 339 if (p2 > 1) { 340 bScore = p2 - 1; 341 bPoint = 0x05; 342 bIsFurtherSide = true; 343 } else { 344 bScore = 1 - p2; 345 bPoint = 0x02; 346 bIsFurtherSide = false; 347 } 348 349 //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. 350 double p3 = yins + zins; 351 if (p3 > 1) { 352 double score = p3 - 1; 353 if (aScore <= bScore && aScore < score) { 354 aScore = score; 355 aPoint = 0x06; 356 aIsFurtherSide = true; 357 } else if (aScore > bScore && bScore < score) { 358 bScore = score; 359 bPoint = 0x06; 360 bIsFurtherSide = true; 361 } 362 } else { 363 double score = 1 - p3; 364 if (aScore <= bScore && aScore < score) { 365 aScore = score; 366 aPoint = 0x01; 367 aIsFurtherSide = false; 368 } else if (aScore > bScore && bScore < score) { 369 bScore = score; 370 bPoint = 0x01; 371 bIsFurtherSide = false; 372 } 373 } 374 375 //Where each of the two closest points are determines how the extra two vertices are calculated. 376 if (aIsFurtherSide == bIsFurtherSide) { 377 if (aIsFurtherSide) { //Both closest points on (1,1,1) side 378 379 //One of the two extra points is (1,1,1) 380 dx_ext0 = dx0 - 1 - 3 * SQUISH_CONSTANT_3D; 381 dy_ext0 = dy0 - 1 - 3 * SQUISH_CONSTANT_3D; 382 dz_ext0 = dz0 - 1 - 3 * SQUISH_CONSTANT_3D; 383 xsv_ext0 = xsb + 1; 384 ysv_ext0 = ysb + 1; 385 zsv_ext0 = zsb + 1; 386 387 //Other extra point is based on the shared axis. 388 byte c = cast(byte)(aPoint & bPoint); 389 if ((c & 0x01) != 0) { 390 dx_ext1 = dx0 - 2 - 2 * SQUISH_CONSTANT_3D; 391 dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D; 392 dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D; 393 xsv_ext1 = xsb + 2; 394 ysv_ext1 = ysb; 395 zsv_ext1 = zsb; 396 } else if ((c & 0x02) != 0) { 397 dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D; 398 dy_ext1 = dy0 - 2 - 2 * SQUISH_CONSTANT_3D; 399 dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D; 400 xsv_ext1 = xsb; 401 ysv_ext1 = ysb + 2; 402 zsv_ext1 = zsb; 403 } else { 404 dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D; 405 dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D; 406 dz_ext1 = dz0 - 2 - 2 * SQUISH_CONSTANT_3D; 407 xsv_ext1 = xsb; 408 ysv_ext1 = ysb; 409 zsv_ext1 = zsb + 2; 410 } 411 } else {//Both closest points on (0,0,0) side 412 413 //One of the two extra points is (0,0,0) 414 dx_ext0 = dx0; 415 dy_ext0 = dy0; 416 dz_ext0 = dz0; 417 xsv_ext0 = xsb; 418 ysv_ext0 = ysb; 419 zsv_ext0 = zsb; 420 421 //Other extra point is based on the omitted axis. 422 byte c = cast(byte)(aPoint | bPoint); 423 if ((c & 0x01) == 0) { 424 dx_ext1 = dx0 + 1 - SQUISH_CONSTANT_3D; 425 dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_3D; 426 dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_3D; 427 xsv_ext1 = xsb - 1; 428 ysv_ext1 = ysb + 1; 429 zsv_ext1 = zsb + 1; 430 } else if ((c & 0x02) == 0) { 431 dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_3D; 432 dy_ext1 = dy0 + 1 - SQUISH_CONSTANT_3D; 433 dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_3D; 434 xsv_ext1 = xsb + 1; 435 ysv_ext1 = ysb - 1; 436 zsv_ext1 = zsb + 1; 437 } else { 438 dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_3D; 439 dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_3D; 440 dz_ext1 = dz0 + 1 - SQUISH_CONSTANT_3D; 441 xsv_ext1 = xsb + 1; 442 ysv_ext1 = ysb + 1; 443 zsv_ext1 = zsb - 1; 444 } 445 } 446 } else { //One point on (0,0,0) side, one point on (1,1,1) side 447 byte c1, c2; 448 if (aIsFurtherSide) { 449 c1 = aPoint; 450 c2 = bPoint; 451 } else { 452 c1 = bPoint; 453 c2 = aPoint; 454 } 455 456 //One contribution is a permutation of (1,1,-1) 457 if ((c1 & 0x01) == 0) { 458 dx_ext0 = dx0 + 1 - SQUISH_CONSTANT_3D; 459 dy_ext0 = dy0 - 1 - SQUISH_CONSTANT_3D; 460 dz_ext0 = dz0 - 1 - SQUISH_CONSTANT_3D; 461 xsv_ext0 = xsb - 1; 462 ysv_ext0 = ysb + 1; 463 zsv_ext0 = zsb + 1; 464 } else if ((c1 & 0x02) == 0) { 465 dx_ext0 = dx0 - 1 - SQUISH_CONSTANT_3D; 466 dy_ext0 = dy0 + 1 - SQUISH_CONSTANT_3D; 467 dz_ext0 = dz0 - 1 - SQUISH_CONSTANT_3D; 468 xsv_ext0 = xsb + 1; 469 ysv_ext0 = ysb - 1; 470 zsv_ext0 = zsb + 1; 471 } else { 472 dx_ext0 = dx0 - 1 - SQUISH_CONSTANT_3D; 473 dy_ext0 = dy0 - 1 - SQUISH_CONSTANT_3D; 474 dz_ext0 = dz0 + 1 - SQUISH_CONSTANT_3D; 475 xsv_ext0 = xsb + 1; 476 ysv_ext0 = ysb + 1; 477 zsv_ext0 = zsb - 1; 478 } 479 480 //One contribution is a permutation of (0,0,2) 481 dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D; 482 dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D; 483 dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D; 484 xsv_ext1 = xsb; 485 ysv_ext1 = ysb; 486 zsv_ext1 = zsb; 487 if ((c2 & 0x01) != 0) { 488 dx_ext1 -= 2; 489 xsv_ext1 += 2; 490 } else if ((c2 & 0x02) != 0) { 491 dy_ext1 -= 2; 492 ysv_ext1 += 2; 493 } else { 494 dz_ext1 -= 2; 495 zsv_ext1 += 2; 496 } 497 } 498 499 //Contribution (1,0,0) 500 double dx1 = dx0 - 1 - SQUISH_CONSTANT_3D; 501 double dy1 = dy0 - 0 - SQUISH_CONSTANT_3D; 502 double dz1 = dz0 - 0 - SQUISH_CONSTANT_3D; 503 double attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1; 504 if (attn1 > 0) { 505 attn1 *= attn1; 506 value += attn1 * attn1 * extrapolate(xsb + 1, ysb + 0, zsb + 0, dx1, dy1, dz1, perm, permGradIndex3D); 507 } 508 509 //Contribution (0,1,0) 510 double dx2 = dx0 - 0 - SQUISH_CONSTANT_3D; 511 double dy2 = dy0 - 1 - SQUISH_CONSTANT_3D; 512 double dz2 = dz1; 513 double attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2; 514 if (attn2 > 0) { 515 attn2 *= attn2; 516 value += attn2 * attn2 * extrapolate(xsb + 0, ysb + 1, zsb + 0, dx2, dy2, dz2, perm, permGradIndex3D); 517 } 518 519 //Contribution (0,0,1) 520 double dx3 = dx2; 521 double dy3 = dy1; 522 double dz3 = dz0 - 1 - SQUISH_CONSTANT_3D; 523 double attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3; 524 if (attn3 > 0) { 525 attn3 *= attn3; 526 value += attn3 * attn3 * extrapolate(xsb + 0, ysb + 0, zsb + 1, dx3, dy3, dz3, perm, permGradIndex3D); 527 } 528 529 //Contribution (1,1,0) 530 double dx4 = dx0 - 1 - 2 * SQUISH_CONSTANT_3D; 531 double dy4 = dy0 - 1 - 2 * SQUISH_CONSTANT_3D; 532 double dz4 = dz0 - 0 - 2 * SQUISH_CONSTANT_3D; 533 double attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4; 534 if (attn4 > 0) { 535 attn4 *= attn4; 536 value += attn4 * attn4 * extrapolate(xsb + 1, ysb + 1, zsb + 0, dx4, dy4, dz4, perm, permGradIndex3D); 537 } 538 539 //Contribution (1,0,1) 540 double dx5 = dx4; 541 double dy5 = dy0 - 0 - 2 * SQUISH_CONSTANT_3D; 542 double dz5 = dz0 - 1 - 2 * SQUISH_CONSTANT_3D; 543 double attn5 = 2 - dx5 * dx5 - dy5 * dy5 - dz5 * dz5; 544 if (attn5 > 0) { 545 attn5 *= attn5; 546 value += attn5 * attn5 * extrapolate(xsb + 1, ysb + 0, zsb + 1, dx5, dy5, dz5, perm, permGradIndex3D); 547 } 548 549 //Contribution (0,1,1) 550 double dx6 = dx0 - 0 - 2 * SQUISH_CONSTANT_3D; 551 double dy6 = dy4; 552 double dz6 = dz5; 553 double attn6 = 2 - dx6 * dx6 - dy6 * dy6 - dz6 * dz6; 554 if (attn6 > 0) { 555 attn6 *= attn6; 556 value += attn6 * attn6 * extrapolate(xsb + 0, ysb + 1, zsb + 1, dx6, dy6, dz6, perm, permGradIndex3D); 557 } 558 } 559 560 //First extra vertex 561 double attn_ext0 = 2 - dx_ext0 * dx_ext0 - dy_ext0 * dy_ext0 - dz_ext0 * dz_ext0; 562 if (attn_ext0 > 0) 563 { 564 attn_ext0 *= attn_ext0; 565 value += attn_ext0 * attn_ext0 * extrapolate(xsv_ext0, ysv_ext0, zsv_ext0, dx_ext0, dy_ext0, dz_ext0, perm, permGradIndex3D); 566 } 567 568 //Second extra vertex 569 double attn_ext1 = 2 - dx_ext1 * dx_ext1 - dy_ext1 * dy_ext1 - dz_ext1 * dz_ext1; 570 if (attn_ext1 > 0) 571 { 572 attn_ext1 *= attn_ext1; 573 value += attn_ext1 * attn_ext1 * extrapolate(xsv_ext1, ysv_ext1, zsv_ext1, dx_ext1, dy_ext1, dz_ext1, perm, permGradIndex3D); 574 } 575 576 return value / NORM_CONSTANT_3D; 577 } 578 579 /** 580 extrapolates for 3D 581 */ 582 @nogc @safe pure double extrapolate(int xsb, int ysb, int zsb, double dx, double dy, double dz, const ref short[256] perm, const ref short[256] permGradIndex3D) 583 { 584 int index = permGradIndex3D[(perm[(perm[xsb & 0xFF] + ysb) & 0xFF] + zsb) & 0xFF]; 585 return GRADIENTS.GRADIENTS_3D[index] * dx + GRADIENTS.GRADIENTS_3D[index + 1] * dy + GRADIENTS.GRADIENTS_3D[index + 2] * dz; 586 }