File indexing completed on 2024-09-22 04:08:01
0001 /* 0002 * OpenSimplex (Simplectic) Noise in C. 0003 * Ported by Stephen M. Cameron from Kurt Spencer's java implementation 0004 * 0005 * v1.1 (October 5, 2014) 0006 * - Added 2D and 4D implementations. 0007 * - Proper gradient sets for all dimensions, from a 0008 * dimensionally-generalizable scheme with an actual 0009 * rhyme and reason behind it. 0010 * - Removed default permutation array in favor of 0011 * default seed. 0012 * - Changed seed-based constructor to be independent 0013 * of any particular randomization library, so results 0014 * will be the same when ported to other languages. 0015 */ 0016 #include <math.h> 0017 #include <stdlib.h> 0018 #include <stdint.h> 0019 #include <string.h> 0020 #include <errno.h> 0021 0022 #include "open-simplex-noise.h" 0023 0024 #define STRETCH_CONSTANT_2D (-0.211324865405187) /* (1 / sqrt(2 + 1) - 1 ) / 2; */ 0025 #define SQUISH_CONSTANT_2D (0.366025403784439) /* (sqrt(2 + 1) -1) / 2; */ 0026 #define STRETCH_CONSTANT_3D (-1.0 / 6.0) /* (1 / sqrt(3 + 1) - 1) / 3; */ 0027 #define SQUISH_CONSTANT_3D (1.0 / 3.0) /* (sqrt(3+1)-1)/3; */ 0028 #define STRETCH_CONSTANT_4D (-0.138196601125011) /* (1 / sqrt(4 + 1) - 1) / 4; */ 0029 #define SQUISH_CONSTANT_4D (0.309016994374947) /* (sqrt(4 + 1) - 1) / 4; */ 0030 0031 #define NORM_CONSTANT_2D (47.0) 0032 #define NORM_CONSTANT_3D (103.0) 0033 #define NORM_CONSTANT_4D (30.0) 0034 0035 #define DEFAULT_SEED (0LL) 0036 0037 struct osn_context { 0038 int16_t *perm; 0039 int16_t *permGradIndex3D; 0040 }; 0041 0042 #define ARRAYSIZE(x) (sizeof((x)) / sizeof((x)[0])) 0043 0044 /* 0045 * Gradients for 2D. They approximate the directions to the 0046 * vertices of an octagon from the center. 0047 */ 0048 static const int8_t gradients2D[] = { 0049 5, 2, 2, 5, 0050 -5, 2, -2, 5, 0051 5, -2, 2, -5, 0052 -5, -2, -2, -5, 0053 }; 0054 0055 /* 0056 * Gradients for 3D. They approximate the directions to the 0057 * vertices of a rhombicuboctahedron from the center, skewed so 0058 * that the triangular and square facets can be inscribed inside 0059 * circles of the same radius. 0060 */ 0061 static const signed char gradients3D[] = { 0062 -11, 4, 4, -4, 11, 4, -4, 4, 11, 0063 11, 4, 4, 4, 11, 4, 4, 4, 11, 0064 -11, -4, 4, -4, -11, 4, -4, -4, 11, 0065 11, -4, 4, 4, -11, 4, 4, -4, 11, 0066 -11, 4, -4, -4, 11, -4, -4, 4, -11, 0067 11, 4, -4, 4, 11, -4, 4, 4, -11, 0068 -11, -4, -4, -4, -11, -4, -4, -4, -11, 0069 11, -4, -4, 4, -11, -4, 4, -4, -11, 0070 }; 0071 0072 /* 0073 * Gradients for 4D. They approximate the directions to the 0074 * vertices of a disprismatotesseractihexadecachoron from the center, 0075 * skewed so that the tetrahedral and cubic facets can be inscribed inside 0076 * spheres of the same radius. 0077 */ 0078 static const signed char gradients4D[] = { 0079 3, 1, 1, 1, 1, 3, 1, 1, 1, 1, 3, 1, 1, 1, 1, 3, 0080 -3, 1, 1, 1, -1, 3, 1, 1, -1, 1, 3, 1, -1, 1, 1, 3, 0081 3, -1, 1, 1, 1, -3, 1, 1, 1, -1, 3, 1, 1, -1, 1, 3, 0082 -3, -1, 1, 1, -1, -3, 1, 1, -1, -1, 3, 1, -1, -1, 1, 3, 0083 3, 1, -1, 1, 1, 3, -1, 1, 1, 1, -3, 1, 1, 1, -1, 3, 0084 -3, 1, -1, 1, -1, 3, -1, 1, -1, 1, -3, 1, -1, 1, -1, 3, 0085 3, -1, -1, 1, 1, -3, -1, 1, 1, -1, -3, 1, 1, -1, -1, 3, 0086 -3, -1, -1, 1, -1, -3, -1, 1, -1, -1, -3, 1, -1, -1, -1, 3, 0087 3, 1, 1, -1, 1, 3, 1, -1, 1, 1, 3, -1, 1, 1, 1, -3, 0088 -3, 1, 1, -1, -1, 3, 1, -1, -1, 1, 3, -1, -1, 1, 1, -3, 0089 3, -1, 1, -1, 1, -3, 1, -1, 1, -1, 3, -1, 1, -1, 1, -3, 0090 -3, -1, 1, -1, -1, -3, 1, -1, -1, -1, 3, -1, -1, -1, 1, -3, 0091 3, 1, -1, -1, 1, 3, -1, -1, 1, 1, -3, -1, 1, 1, -1, -3, 0092 -3, 1, -1, -1, -1, 3, -1, -1, -1, 1, -3, -1, -1, 1, -1, -3, 0093 3, -1, -1, -1, 1, -3, -1, -1, 1, -1, -3, -1, 1, -1, -1, -3, 0094 -3, -1, -1, -1, -1, -3, -1, -1, -1, -1, -3, -1, -1, -1, -1, -3, 0095 }; 0096 0097 static double extrapolate2(struct osn_context *ctx, int xsb, int ysb, double dx, double dy) 0098 { 0099 int16_t *perm = ctx->perm; 0100 int index = perm[(perm[xsb & 0xFF] + ysb) & 0xFF] & 0x0E; 0101 return gradients2D[index] * dx 0102 + gradients2D[index + 1] * dy; 0103 } 0104 0105 static double extrapolate3(struct osn_context *ctx, int xsb, int ysb, int zsb, double dx, double dy, double dz) 0106 { 0107 int16_t *perm = ctx->perm; 0108 int16_t *permGradIndex3D = ctx->permGradIndex3D; 0109 int index = permGradIndex3D[(perm[(perm[xsb & 0xFF] + ysb) & 0xFF] + zsb) & 0xFF]; 0110 return gradients3D[index] * dx 0111 + gradients3D[index + 1] * dy 0112 + gradients3D[index + 2] * dz; 0113 } 0114 0115 static double extrapolate4(struct osn_context *ctx, int xsb, int ysb, int zsb, int wsb, double dx, double dy, double dz, double dw) 0116 { 0117 int16_t *perm = ctx->perm; 0118 int index = perm[(perm[(perm[(perm[xsb & 0xFF] + ysb) & 0xFF] + zsb) & 0xFF] + wsb) & 0xFF] & 0xFC; 0119 return gradients4D[index] * dx 0120 + gradients4D[index + 1] * dy 0121 + gradients4D[index + 2] * dz 0122 + gradients4D[index + 3] * dw; 0123 } 0124 0125 static INLINE int fastFloor(double x) { 0126 int xi = (int) x; 0127 return x < xi ? xi - 1 : xi; 0128 } 0129 0130 static int allocate_perm(struct osn_context *ctx, int nperm, int ngrad) 0131 { 0132 if (ctx->perm) 0133 free(ctx->perm); 0134 if (ctx->permGradIndex3D) 0135 free(ctx->permGradIndex3D); 0136 ctx->perm = (int16_t *) malloc(sizeof(*ctx->perm) * nperm); 0137 if (!ctx->perm) 0138 return -ENOMEM; 0139 ctx->permGradIndex3D = (int16_t *) malloc(sizeof(*ctx->permGradIndex3D) * ngrad); 0140 if (!ctx->permGradIndex3D) { 0141 free(ctx->perm); 0142 return -ENOMEM; 0143 } 0144 return 0; 0145 } 0146 0147 int open_simplex_noise_init_perm(struct osn_context *ctx, int16_t p[], int nelements) 0148 { 0149 int i, rc; 0150 0151 rc = allocate_perm(ctx, nelements, 256); 0152 if (rc) 0153 return rc; 0154 memcpy(ctx->perm, p, sizeof(*ctx->perm) * nelements); 0155 0156 for (i = 0; i < 256; i++) { 0157 /* Since 3D has 24 gradients, simple bitmask won't work, so precompute modulo array. */ 0158 ctx->permGradIndex3D[i] = (int16_t)((ctx->perm[i] % (ARRAYSIZE(gradients3D) / 3)) * 3); 0159 } 0160 return 0; 0161 } 0162 0163 /* 0164 * Initializes using a permutation array generated from a 64-bit seed. 0165 * Generates a proper permutation (i.e. doesn't merely perform N successive pair 0166 * swaps on a base array). Uses a simple 64-bit LCG. 0167 */ 0168 int open_simplex_noise(int64_t seed, struct osn_context **ctx) 0169 { 0170 int rc; 0171 int16_t source[256]; 0172 int i; 0173 int16_t *perm; 0174 int16_t *permGradIndex3D; 0175 int r; 0176 0177 *ctx = (struct osn_context *) malloc(sizeof(**ctx)); 0178 if (!(*ctx)) 0179 return -ENOMEM; 0180 (*ctx)->perm = NULL; 0181 (*ctx)->permGradIndex3D = NULL; 0182 0183 rc = allocate_perm(*ctx, 256, 256); 0184 if (rc) { 0185 free(*ctx); 0186 return rc; 0187 } 0188 0189 perm = (*ctx)->perm; 0190 permGradIndex3D = (*ctx)->permGradIndex3D; 0191 0192 for (i = 0; i < 256; i++) 0193 source[i] = (int16_t) i; 0194 seed = seed * 6364136223846793005LL + 1442695040888963407LL; 0195 seed = seed * 6364136223846793005LL + 1442695040888963407LL; 0196 seed = seed * 6364136223846793005LL + 1442695040888963407LL; 0197 for (i = 255; i >= 0; i--) { 0198 seed = seed * 6364136223846793005LL + 1442695040888963407LL; 0199 r = (int)((seed + 31) % (i + 1)); 0200 if (r < 0) 0201 r += (i + 1); 0202 perm[i] = source[r]; 0203 permGradIndex3D[i] = (short)((perm[i] % (ARRAYSIZE(gradients3D) / 3)) * 3); 0204 source[r] = source[i]; 0205 } 0206 return 0; 0207 } 0208 0209 void open_simplex_noise_free(struct osn_context *ctx) 0210 { 0211 if (!ctx) 0212 return; 0213 if (ctx->perm) { 0214 free(ctx->perm); 0215 ctx->perm = NULL; 0216 } 0217 if (ctx->permGradIndex3D) { 0218 free(ctx->permGradIndex3D); 0219 ctx->permGradIndex3D = NULL; 0220 } 0221 free(ctx); 0222 } 0223 0224 /* 2D OpenSimplex (Simplectic) Noise. */ 0225 double open_simplex_noise2(struct osn_context *ctx, double x, double y) 0226 { 0227 0228 /* Place input coordinates onto grid. */ 0229 double stretchOffset = (x + y) * STRETCH_CONSTANT_2D; 0230 double xs = x + stretchOffset; 0231 double ys = y + stretchOffset; 0232 0233 /* Floor to get grid coordinates of rhombus (stretched square) super-cell origin. */ 0234 int xsb = fastFloor(xs); 0235 int ysb = fastFloor(ys); 0236 0237 /* Skew out to get actual coordinates of rhombus origin. We'll need these later. */ 0238 double squishOffset = (xsb + ysb) * SQUISH_CONSTANT_2D; 0239 double xb = xsb + squishOffset; 0240 double yb = ysb + squishOffset; 0241 0242 /* Compute grid coordinates relative to rhombus origin. */ 0243 double xins = xs - xsb; 0244 double yins = ys - ysb; 0245 0246 /* Sum those together to get a value that determines which region we're in. */ 0247 double inSum = xins + yins; 0248 0249 /* Positions relative to origin point. */ 0250 double dx0 = x - xb; 0251 double dy0 = y - yb; 0252 0253 /* We'll be defining these inside the next block and using them afterwards. */ 0254 double dx_ext, dy_ext; 0255 int xsv_ext, ysv_ext; 0256 0257 double dx1; 0258 double dy1; 0259 double attn1; 0260 double dx2; 0261 double dy2; 0262 double attn2; 0263 double zins; 0264 double attn0; 0265 double attn_ext; 0266 0267 double value = 0; 0268 0269 /* Contribution (1,0) */ 0270 dx1 = dx0 - 1 - SQUISH_CONSTANT_2D; 0271 dy1 = dy0 - 0 - SQUISH_CONSTANT_2D; 0272 attn1 = 2 - dx1 * dx1 - dy1 * dy1; 0273 if (attn1 > 0) { 0274 attn1 *= attn1; 0275 value += attn1 * attn1 * extrapolate2(ctx, xsb + 1, ysb + 0, dx1, dy1); 0276 } 0277 0278 /* Contribution (0,1) */ 0279 dx2 = dx0 - 0 - SQUISH_CONSTANT_2D; 0280 dy2 = dy0 - 1 - SQUISH_CONSTANT_2D; 0281 attn2 = 2 - dx2 * dx2 - dy2 * dy2; 0282 if (attn2 > 0) { 0283 attn2 *= attn2; 0284 value += attn2 * attn2 * extrapolate2(ctx, xsb + 0, ysb + 1, dx2, dy2); 0285 } 0286 0287 if (inSum <= 1) { /* We're inside the triangle (2-Simplex) at (0,0) */ 0288 zins = 1 - inSum; 0289 if (zins > xins || zins > yins) { /* (0,0) is one of the closest two triangular vertices */ 0290 if (xins > yins) { 0291 xsv_ext = xsb + 1; 0292 ysv_ext = ysb - 1; 0293 dx_ext = dx0 - 1; 0294 dy_ext = dy0 + 1; 0295 } else { 0296 xsv_ext = xsb - 1; 0297 ysv_ext = ysb + 1; 0298 dx_ext = dx0 + 1; 0299 dy_ext = dy0 - 1; 0300 } 0301 } else { /* (1,0) and (0,1) are the closest two vertices. */ 0302 xsv_ext = xsb + 1; 0303 ysv_ext = ysb + 1; 0304 dx_ext = dx0 - 1 - 2 * SQUISH_CONSTANT_2D; 0305 dy_ext = dy0 - 1 - 2 * SQUISH_CONSTANT_2D; 0306 } 0307 } else { /* We're inside the triangle (2-Simplex) at (1,1) */ 0308 zins = 2 - inSum; 0309 if (zins < xins || zins < yins) { /* (0,0) is one of the closest two triangular vertices */ 0310 if (xins > yins) { 0311 xsv_ext = xsb + 2; 0312 ysv_ext = ysb + 0; 0313 dx_ext = dx0 - 2 - 2 * SQUISH_CONSTANT_2D; 0314 dy_ext = dy0 + 0 - 2 * SQUISH_CONSTANT_2D; 0315 } else { 0316 xsv_ext = xsb + 0; 0317 ysv_ext = ysb + 2; 0318 dx_ext = dx0 + 0 - 2 * SQUISH_CONSTANT_2D; 0319 dy_ext = dy0 - 2 - 2 * SQUISH_CONSTANT_2D; 0320 } 0321 } else { /* (1,0) and (0,1) are the closest two vertices. */ 0322 dx_ext = dx0; 0323 dy_ext = dy0; 0324 xsv_ext = xsb; 0325 ysv_ext = ysb; 0326 } 0327 xsb += 1; 0328 ysb += 1; 0329 dx0 = dx0 - 1 - 2 * SQUISH_CONSTANT_2D; 0330 dy0 = dy0 - 1 - 2 * SQUISH_CONSTANT_2D; 0331 } 0332 0333 /* Contribution (0,0) or (1,1) */ 0334 attn0 = 2 - dx0 * dx0 - dy0 * dy0; 0335 if (attn0 > 0) { 0336 attn0 *= attn0; 0337 value += attn0 * attn0 * extrapolate2(ctx, xsb, ysb, dx0, dy0); 0338 } 0339 0340 /* Extra Vertex */ 0341 attn_ext = 2 - dx_ext * dx_ext - dy_ext * dy_ext; 0342 if (attn_ext > 0) { 0343 attn_ext *= attn_ext; 0344 value += attn_ext * attn_ext * extrapolate2(ctx, xsv_ext, ysv_ext, dx_ext, dy_ext); 0345 } 0346 0347 return value / NORM_CONSTANT_2D; 0348 } 0349 0350 /* 0351 * 3D OpenSimplex (Simplectic) Noise 0352 */ 0353 double open_simplex_noise3(struct osn_context *ctx, double x, double y, double z) 0354 { 0355 0356 /* Place input coordinates on simplectic honeycomb. */ 0357 double stretchOffset = (x + y + z) * STRETCH_CONSTANT_3D; 0358 double xs = x + stretchOffset; 0359 double ys = y + stretchOffset; 0360 double zs = z + stretchOffset; 0361 0362 /* Floor to get simplectic honeycomb coordinates of rhombohedron (stretched cube) super-cell origin. */ 0363 int xsb = fastFloor(xs); 0364 int ysb = fastFloor(ys); 0365 int zsb = fastFloor(zs); 0366 0367 /* Skew out to get actual coordinates of rhombohedron origin. We'll need these later. */ 0368 double squishOffset = (xsb + ysb + zsb) * SQUISH_CONSTANT_3D; 0369 double xb = xsb + squishOffset; 0370 double yb = ysb + squishOffset; 0371 double zb = zsb + squishOffset; 0372 0373 /* Compute simplectic honeycomb coordinates relative to rhombohedral origin. */ 0374 double xins = xs - xsb; 0375 double yins = ys - ysb; 0376 double zins = zs - zsb; 0377 0378 /* Sum those together to get a value that determines which region we're in. */ 0379 double inSum = xins + yins + zins; 0380 0381 /* Positions relative to origin point. */ 0382 double dx0 = x - xb; 0383 double dy0 = y - yb; 0384 double dz0 = z - zb; 0385 0386 /* We'll be defining these inside the next block and using them afterwards. */ 0387 double dx_ext0, dy_ext0, dz_ext0; 0388 double dx_ext1, dy_ext1, dz_ext1; 0389 int xsv_ext0, ysv_ext0, zsv_ext0; 0390 int xsv_ext1, ysv_ext1, zsv_ext1; 0391 0392 double wins; 0393 int8_t c, c1, c2; 0394 int8_t aPoint, bPoint; 0395 double aScore, bScore; 0396 int aIsFurtherSide; 0397 int bIsFurtherSide; 0398 double p1, p2, p3; 0399 double score; 0400 double attn0, attn1, attn2, attn3, attn4, attn5, attn6; 0401 double dx1, dy1, dz1; 0402 double dx2, dy2, dz2; 0403 double dx3, dy3, dz3; 0404 double dx4, dy4, dz4; 0405 double dx5, dy5, dz5; 0406 double dx6, dy6, dz6; 0407 double attn_ext0, attn_ext1; 0408 0409 double value = 0; 0410 if (inSum <= 1) { /* We're inside the tetrahedron (3-Simplex) at (0,0,0) */ 0411 0412 /* Determine which two of (0,0,1), (0,1,0), (1,0,0) are closest. */ 0413 aPoint = 0x01; 0414 aScore = xins; 0415 bPoint = 0x02; 0416 bScore = yins; 0417 if (aScore >= bScore && zins > bScore) { 0418 bScore = zins; 0419 bPoint = 0x04; 0420 } else if (aScore < bScore && zins > aScore) { 0421 aScore = zins; 0422 aPoint = 0x04; 0423 } 0424 0425 /* Now we determine the two lattice points not part of the tetrahedron that may contribute. 0426 This depends on the closest two tetrahedral vertices, including (0,0,0) */ 0427 wins = 1 - inSum; 0428 if (wins > aScore || wins > bScore) { /* (0,0,0) is one of the closest two tetrahedral vertices. */ 0429 c = (bScore > aScore ? bPoint : aPoint); /* Our other closest vertex is the closest out of a and b. */ 0430 0431 if ((c & 0x01) == 0) { 0432 xsv_ext0 = xsb - 1; 0433 xsv_ext1 = xsb; 0434 dx_ext0 = dx0 + 1; 0435 dx_ext1 = dx0; 0436 } else { 0437 xsv_ext0 = xsv_ext1 = xsb + 1; 0438 dx_ext0 = dx_ext1 = dx0 - 1; 0439 } 0440 0441 if ((c & 0x02) == 0) { 0442 ysv_ext0 = ysv_ext1 = ysb; 0443 dy_ext0 = dy_ext1 = dy0; 0444 if ((c & 0x01) == 0) { 0445 ysv_ext1 -= 1; 0446 dy_ext1 += 1; 0447 } else { 0448 ysv_ext0 -= 1; 0449 dy_ext0 += 1; 0450 } 0451 } else { 0452 ysv_ext0 = ysv_ext1 = ysb + 1; 0453 dy_ext0 = dy_ext1 = dy0 - 1; 0454 } 0455 0456 if ((c & 0x04) == 0) { 0457 zsv_ext0 = zsb; 0458 zsv_ext1 = zsb - 1; 0459 dz_ext0 = dz0; 0460 dz_ext1 = dz0 + 1; 0461 } else { 0462 zsv_ext0 = zsv_ext1 = zsb + 1; 0463 dz_ext0 = dz_ext1 = dz0 - 1; 0464 } 0465 } else { /* (0,0,0) is not one of the closest two tetrahedral vertices. */ 0466 c = (int8_t)(aPoint | bPoint); /* Our two extra vertices are determined by the closest two. */ 0467 0468 if ((c & 0x01) == 0) { 0469 xsv_ext0 = xsb; 0470 xsv_ext1 = xsb - 1; 0471 dx_ext0 = dx0 - 2 * SQUISH_CONSTANT_3D; 0472 dx_ext1 = dx0 + 1 - SQUISH_CONSTANT_3D; 0473 } else { 0474 xsv_ext0 = xsv_ext1 = xsb + 1; 0475 dx_ext0 = dx0 - 1 - 2 * SQUISH_CONSTANT_3D; 0476 dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_3D; 0477 } 0478 0479 if ((c & 0x02) == 0) { 0480 ysv_ext0 = ysb; 0481 ysv_ext1 = ysb - 1; 0482 dy_ext0 = dy0 - 2 * SQUISH_CONSTANT_3D; 0483 dy_ext1 = dy0 + 1 - SQUISH_CONSTANT_3D; 0484 } else { 0485 ysv_ext0 = ysv_ext1 = ysb + 1; 0486 dy_ext0 = dy0 - 1 - 2 * SQUISH_CONSTANT_3D; 0487 dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_3D; 0488 } 0489 0490 if ((c & 0x04) == 0) { 0491 zsv_ext0 = zsb; 0492 zsv_ext1 = zsb - 1; 0493 dz_ext0 = dz0 - 2 * SQUISH_CONSTANT_3D; 0494 dz_ext1 = dz0 + 1 - SQUISH_CONSTANT_3D; 0495 } else { 0496 zsv_ext0 = zsv_ext1 = zsb + 1; 0497 dz_ext0 = dz0 - 1 - 2 * SQUISH_CONSTANT_3D; 0498 dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_3D; 0499 } 0500 } 0501 0502 /* Contribution (0,0,0) */ 0503 attn0 = 2 - dx0 * dx0 - dy0 * dy0 - dz0 * dz0; 0504 if (attn0 > 0) { 0505 attn0 *= attn0; 0506 value += attn0 * attn0 * extrapolate3(ctx, xsb + 0, ysb + 0, zsb + 0, dx0, dy0, dz0); 0507 } 0508 0509 /* Contribution (1,0,0) */ 0510 dx1 = dx0 - 1 - SQUISH_CONSTANT_3D; 0511 dy1 = dy0 - 0 - SQUISH_CONSTANT_3D; 0512 dz1 = dz0 - 0 - SQUISH_CONSTANT_3D; 0513 attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1; 0514 if (attn1 > 0) { 0515 attn1 *= attn1; 0516 value += attn1 * attn1 * extrapolate3(ctx, xsb + 1, ysb + 0, zsb + 0, dx1, dy1, dz1); 0517 } 0518 0519 /* Contribution (0,1,0) */ 0520 dx2 = dx0 - 0 - SQUISH_CONSTANT_3D; 0521 dy2 = dy0 - 1 - SQUISH_CONSTANT_3D; 0522 dz2 = dz1; 0523 attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2; 0524 if (attn2 > 0) { 0525 attn2 *= attn2; 0526 value += attn2 * attn2 * extrapolate3(ctx, xsb + 0, ysb + 1, zsb + 0, dx2, dy2, dz2); 0527 } 0528 0529 /* Contribution (0,0,1) */ 0530 dx3 = dx2; 0531 dy3 = dy1; 0532 dz3 = dz0 - 1 - SQUISH_CONSTANT_3D; 0533 attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3; 0534 if (attn3 > 0) { 0535 attn3 *= attn3; 0536 value += attn3 * attn3 * extrapolate3(ctx, xsb + 0, ysb + 0, zsb + 1, dx3, dy3, dz3); 0537 } 0538 } else if (inSum >= 2) { /* We're inside the tetrahedron (3-Simplex) at (1,1,1) */ 0539 0540 /* Determine which two tetrahedral vertices are the closest, out of (1,1,0), (1,0,1), (0,1,1) but not (1,1,1). */ 0541 aPoint = 0x06; 0542 aScore = xins; 0543 bPoint = 0x05; 0544 bScore = yins; 0545 if (aScore <= bScore && zins < bScore) { 0546 bScore = zins; 0547 bPoint = 0x03; 0548 } else if (aScore > bScore && zins < aScore) { 0549 aScore = zins; 0550 aPoint = 0x03; 0551 } 0552 0553 /* Now we determine the two lattice points not part of the tetrahedron that may contribute. 0554 This depends on the closest two tetrahedral vertices, including (1,1,1) */ 0555 wins = 3 - inSum; 0556 if (wins < aScore || wins < bScore) { /* (1,1,1) is one of the closest two tetrahedral vertices. */ 0557 c = (bScore < aScore ? bPoint : aPoint); /* Our other closest vertex is the closest out of a and b. */ 0558 0559 if ((c & 0x01) != 0) { 0560 xsv_ext0 = xsb + 2; 0561 xsv_ext1 = xsb + 1; 0562 dx_ext0 = dx0 - 2 - 3 * SQUISH_CONSTANT_3D; 0563 dx_ext1 = dx0 - 1 - 3 * SQUISH_CONSTANT_3D; 0564 } else { 0565 xsv_ext0 = xsv_ext1 = xsb; 0566 dx_ext0 = dx_ext1 = dx0 - 3 * SQUISH_CONSTANT_3D; 0567 } 0568 0569 if ((c & 0x02) != 0) { 0570 ysv_ext0 = ysv_ext1 = ysb + 1; 0571 dy_ext0 = dy_ext1 = dy0 - 1 - 3 * SQUISH_CONSTANT_3D; 0572 if ((c & 0x01) != 0) { 0573 ysv_ext1 += 1; 0574 dy_ext1 -= 1; 0575 } else { 0576 ysv_ext0 += 1; 0577 dy_ext0 -= 1; 0578 } 0579 } else { 0580 ysv_ext0 = ysv_ext1 = ysb; 0581 dy_ext0 = dy_ext1 = dy0 - 3 * SQUISH_CONSTANT_3D; 0582 } 0583 0584 if ((c & 0x04) != 0) { 0585 zsv_ext0 = zsb + 1; 0586 zsv_ext1 = zsb + 2; 0587 dz_ext0 = dz0 - 1 - 3 * SQUISH_CONSTANT_3D; 0588 dz_ext1 = dz0 - 2 - 3 * SQUISH_CONSTANT_3D; 0589 } else { 0590 zsv_ext0 = zsv_ext1 = zsb; 0591 dz_ext0 = dz_ext1 = dz0 - 3 * SQUISH_CONSTANT_3D; 0592 } 0593 } else { /* (1,1,1) is not one of the closest two tetrahedral vertices. */ 0594 c = (int8_t)(aPoint & bPoint); /* Our two extra vertices are determined by the closest two. */ 0595 0596 if ((c & 0x01) != 0) { 0597 xsv_ext0 = xsb + 1; 0598 xsv_ext1 = xsb + 2; 0599 dx_ext0 = dx0 - 1 - SQUISH_CONSTANT_3D; 0600 dx_ext1 = dx0 - 2 - 2 * SQUISH_CONSTANT_3D; 0601 } else { 0602 xsv_ext0 = xsv_ext1 = xsb; 0603 dx_ext0 = dx0 - SQUISH_CONSTANT_3D; 0604 dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D; 0605 } 0606 0607 if ((c & 0x02) != 0) { 0608 ysv_ext0 = ysb + 1; 0609 ysv_ext1 = ysb + 2; 0610 dy_ext0 = dy0 - 1 - SQUISH_CONSTANT_3D; 0611 dy_ext1 = dy0 - 2 - 2 * SQUISH_CONSTANT_3D; 0612 } else { 0613 ysv_ext0 = ysv_ext1 = ysb; 0614 dy_ext0 = dy0 - SQUISH_CONSTANT_3D; 0615 dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D; 0616 } 0617 0618 if ((c & 0x04) != 0) { 0619 zsv_ext0 = zsb + 1; 0620 zsv_ext1 = zsb + 2; 0621 dz_ext0 = dz0 - 1 - SQUISH_CONSTANT_3D; 0622 dz_ext1 = dz0 - 2 - 2 * SQUISH_CONSTANT_3D; 0623 } else { 0624 zsv_ext0 = zsv_ext1 = zsb; 0625 dz_ext0 = dz0 - SQUISH_CONSTANT_3D; 0626 dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D; 0627 } 0628 } 0629 0630 /* Contribution (1,1,0) */ 0631 dx3 = dx0 - 1 - 2 * SQUISH_CONSTANT_3D; 0632 dy3 = dy0 - 1 - 2 * SQUISH_CONSTANT_3D; 0633 dz3 = dz0 - 0 - 2 * SQUISH_CONSTANT_3D; 0634 attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3; 0635 if (attn3 > 0) { 0636 attn3 *= attn3; 0637 value += attn3 * attn3 * extrapolate3(ctx, xsb + 1, ysb + 1, zsb + 0, dx3, dy3, dz3); 0638 } 0639 0640 /* Contribution (1,0,1) */ 0641 dx2 = dx3; 0642 dy2 = dy0 - 0 - 2 * SQUISH_CONSTANT_3D; 0643 dz2 = dz0 - 1 - 2 * SQUISH_CONSTANT_3D; 0644 attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2; 0645 if (attn2 > 0) { 0646 attn2 *= attn2; 0647 value += attn2 * attn2 * extrapolate3(ctx, xsb + 1, ysb + 0, zsb + 1, dx2, dy2, dz2); 0648 } 0649 0650 /* Contribution (0,1,1) */ 0651 dx1 = dx0 - 0 - 2 * SQUISH_CONSTANT_3D; 0652 dy1 = dy3; 0653 dz1 = dz2; 0654 attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1; 0655 if (attn1 > 0) { 0656 attn1 *= attn1; 0657 value += attn1 * attn1 * extrapolate3(ctx, xsb + 0, ysb + 1, zsb + 1, dx1, dy1, dz1); 0658 } 0659 0660 /* Contribution (1,1,1) */ 0661 dx0 = dx0 - 1 - 3 * SQUISH_CONSTANT_3D; 0662 dy0 = dy0 - 1 - 3 * SQUISH_CONSTANT_3D; 0663 dz0 = dz0 - 1 - 3 * SQUISH_CONSTANT_3D; 0664 attn0 = 2 - dx0 * dx0 - dy0 * dy0 - dz0 * dz0; 0665 if (attn0 > 0) { 0666 attn0 *= attn0; 0667 value += attn0 * attn0 * extrapolate3(ctx, xsb + 1, ysb + 1, zsb + 1, dx0, dy0, dz0); 0668 } 0669 } else { /* We're inside the octahedron (Rectified 3-Simplex) in between. 0670 Decide between point (0,0,1) and (1,1,0) as closest */ 0671 p1 = xins + yins; 0672 if (p1 > 1) { 0673 aScore = p1 - 1; 0674 aPoint = 0x03; 0675 aIsFurtherSide = 1; 0676 } else { 0677 aScore = 1 - p1; 0678 aPoint = 0x04; 0679 aIsFurtherSide = 0; 0680 } 0681 0682 /* Decide between point (0,1,0) and (1,0,1) as closest */ 0683 p2 = xins + zins; 0684 if (p2 > 1) { 0685 bScore = p2 - 1; 0686 bPoint = 0x05; 0687 bIsFurtherSide = 1; 0688 } else { 0689 bScore = 1 - p2; 0690 bPoint = 0x02; 0691 bIsFurtherSide = 0; 0692 } 0693 0694 /* 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. */ 0695 p3 = yins + zins; 0696 if (p3 > 1) { 0697 score = p3 - 1; 0698 if (aScore <= bScore && aScore < score) { 0699 aScore = score; 0700 aPoint = 0x06; 0701 aIsFurtherSide = 1; 0702 } else if (aScore > bScore && bScore < score) { 0703 bScore = score; 0704 bPoint = 0x06; 0705 bIsFurtherSide = 1; 0706 } 0707 } else { 0708 score = 1 - p3; 0709 if (aScore <= bScore && aScore < score) { 0710 aScore = score; 0711 aPoint = 0x01; 0712 aIsFurtherSide = 0; 0713 } else if (aScore > bScore && bScore < score) { 0714 bScore = score; 0715 bPoint = 0x01; 0716 bIsFurtherSide = 0; 0717 } 0718 } 0719 0720 /* Where each of the two closest points are determines how the extra two vertices are calculated. */ 0721 if (aIsFurtherSide == bIsFurtherSide) { 0722 if (aIsFurtherSide) { /* Both closest points on (1,1,1) side */ 0723 0724 /* One of the two extra points is (1,1,1) */ 0725 dx_ext0 = dx0 - 1 - 3 * SQUISH_CONSTANT_3D; 0726 dy_ext0 = dy0 - 1 - 3 * SQUISH_CONSTANT_3D; 0727 dz_ext0 = dz0 - 1 - 3 * SQUISH_CONSTANT_3D; 0728 xsv_ext0 = xsb + 1; 0729 ysv_ext0 = ysb + 1; 0730 zsv_ext0 = zsb + 1; 0731 0732 /* Other extra point is based on the shared axis. */ 0733 c = (int8_t)(aPoint & bPoint); 0734 if ((c & 0x01) != 0) { 0735 dx_ext1 = dx0 - 2 - 2 * SQUISH_CONSTANT_3D; 0736 dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D; 0737 dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D; 0738 xsv_ext1 = xsb + 2; 0739 ysv_ext1 = ysb; 0740 zsv_ext1 = zsb; 0741 } else if ((c & 0x02) != 0) { 0742 dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D; 0743 dy_ext1 = dy0 - 2 - 2 * SQUISH_CONSTANT_3D; 0744 dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D; 0745 xsv_ext1 = xsb; 0746 ysv_ext1 = ysb + 2; 0747 zsv_ext1 = zsb; 0748 } else { 0749 dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D; 0750 dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D; 0751 dz_ext1 = dz0 - 2 - 2 * SQUISH_CONSTANT_3D; 0752 xsv_ext1 = xsb; 0753 ysv_ext1 = ysb; 0754 zsv_ext1 = zsb + 2; 0755 } 0756 } else { /* Both closest points on (0,0,0) side */ 0757 0758 /* One of the two extra points is (0,0,0) */ 0759 dx_ext0 = dx0; 0760 dy_ext0 = dy0; 0761 dz_ext0 = dz0; 0762 xsv_ext0 = xsb; 0763 ysv_ext0 = ysb; 0764 zsv_ext0 = zsb; 0765 0766 /* Other extra point is based on the omitted axis. */ 0767 c = (int8_t)(aPoint | bPoint); 0768 if ((c & 0x01) == 0) { 0769 dx_ext1 = dx0 + 1 - SQUISH_CONSTANT_3D; 0770 dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_3D; 0771 dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_3D; 0772 xsv_ext1 = xsb - 1; 0773 ysv_ext1 = ysb + 1; 0774 zsv_ext1 = zsb + 1; 0775 } else if ((c & 0x02) == 0) { 0776 dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_3D; 0777 dy_ext1 = dy0 + 1 - SQUISH_CONSTANT_3D; 0778 dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_3D; 0779 xsv_ext1 = xsb + 1; 0780 ysv_ext1 = ysb - 1; 0781 zsv_ext1 = zsb + 1; 0782 } else { 0783 dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_3D; 0784 dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_3D; 0785 dz_ext1 = dz0 + 1 - SQUISH_CONSTANT_3D; 0786 xsv_ext1 = xsb + 1; 0787 ysv_ext1 = ysb + 1; 0788 zsv_ext1 = zsb - 1; 0789 } 0790 } 0791 } else { /* One point on (0,0,0) side, one point on (1,1,1) side */ 0792 if (aIsFurtherSide) { 0793 c1 = aPoint; 0794 c2 = bPoint; 0795 } else { 0796 c1 = bPoint; 0797 c2 = aPoint; 0798 } 0799 0800 /* One contribution is a permutation of (1,1,-1) */ 0801 if ((c1 & 0x01) == 0) { 0802 dx_ext0 = dx0 + 1 - SQUISH_CONSTANT_3D; 0803 dy_ext0 = dy0 - 1 - SQUISH_CONSTANT_3D; 0804 dz_ext0 = dz0 - 1 - SQUISH_CONSTANT_3D; 0805 xsv_ext0 = xsb - 1; 0806 ysv_ext0 = ysb + 1; 0807 zsv_ext0 = zsb + 1; 0808 } else if ((c1 & 0x02) == 0) { 0809 dx_ext0 = dx0 - 1 - SQUISH_CONSTANT_3D; 0810 dy_ext0 = dy0 + 1 - SQUISH_CONSTANT_3D; 0811 dz_ext0 = dz0 - 1 - SQUISH_CONSTANT_3D; 0812 xsv_ext0 = xsb + 1; 0813 ysv_ext0 = ysb - 1; 0814 zsv_ext0 = zsb + 1; 0815 } else { 0816 dx_ext0 = dx0 - 1 - SQUISH_CONSTANT_3D; 0817 dy_ext0 = dy0 - 1 - SQUISH_CONSTANT_3D; 0818 dz_ext0 = dz0 + 1 - SQUISH_CONSTANT_3D; 0819 xsv_ext0 = xsb + 1; 0820 ysv_ext0 = ysb + 1; 0821 zsv_ext0 = zsb - 1; 0822 } 0823 0824 /* One contribution is a permutation of (0,0,2) */ 0825 dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D; 0826 dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D; 0827 dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D; 0828 xsv_ext1 = xsb; 0829 ysv_ext1 = ysb; 0830 zsv_ext1 = zsb; 0831 if ((c2 & 0x01) != 0) { 0832 dx_ext1 -= 2; 0833 xsv_ext1 += 2; 0834 } else if ((c2 & 0x02) != 0) { 0835 dy_ext1 -= 2; 0836 ysv_ext1 += 2; 0837 } else { 0838 dz_ext1 -= 2; 0839 zsv_ext1 += 2; 0840 } 0841 } 0842 0843 /* Contribution (1,0,0) */ 0844 dx1 = dx0 - 1 - SQUISH_CONSTANT_3D; 0845 dy1 = dy0 - 0 - SQUISH_CONSTANT_3D; 0846 dz1 = dz0 - 0 - SQUISH_CONSTANT_3D; 0847 attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1; 0848 if (attn1 > 0) { 0849 attn1 *= attn1; 0850 value += attn1 * attn1 * extrapolate3(ctx, xsb + 1, ysb + 0, zsb + 0, dx1, dy1, dz1); 0851 } 0852 0853 /* Contribution (0,1,0) */ 0854 dx2 = dx0 - 0 - SQUISH_CONSTANT_3D; 0855 dy2 = dy0 - 1 - SQUISH_CONSTANT_3D; 0856 dz2 = dz1; 0857 attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2; 0858 if (attn2 > 0) { 0859 attn2 *= attn2; 0860 value += attn2 * attn2 * extrapolate3(ctx, xsb + 0, ysb + 1, zsb + 0, dx2, dy2, dz2); 0861 } 0862 0863 /* Contribution (0,0,1) */ 0864 dx3 = dx2; 0865 dy3 = dy1; 0866 dz3 = dz0 - 1 - SQUISH_CONSTANT_3D; 0867 attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3; 0868 if (attn3 > 0) { 0869 attn3 *= attn3; 0870 value += attn3 * attn3 * extrapolate3(ctx, xsb + 0, ysb + 0, zsb + 1, dx3, dy3, dz3); 0871 } 0872 0873 /* Contribution (1,1,0) */ 0874 dx4 = dx0 - 1 - 2 * SQUISH_CONSTANT_3D; 0875 dy4 = dy0 - 1 - 2 * SQUISH_CONSTANT_3D; 0876 dz4 = dz0 - 0 - 2 * SQUISH_CONSTANT_3D; 0877 attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4; 0878 if (attn4 > 0) { 0879 attn4 *= attn4; 0880 value += attn4 * attn4 * extrapolate3(ctx, xsb + 1, ysb + 1, zsb + 0, dx4, dy4, dz4); 0881 } 0882 0883 /* Contribution (1,0,1) */ 0884 dx5 = dx4; 0885 dy5 = dy0 - 0 - 2 * SQUISH_CONSTANT_3D; 0886 dz5 = dz0 - 1 - 2 * SQUISH_CONSTANT_3D; 0887 attn5 = 2 - dx5 * dx5 - dy5 * dy5 - dz5 * dz5; 0888 if (attn5 > 0) { 0889 attn5 *= attn5; 0890 value += attn5 * attn5 * extrapolate3(ctx, xsb + 1, ysb + 0, zsb + 1, dx5, dy5, dz5); 0891 } 0892 0893 /* Contribution (0,1,1) */ 0894 dx6 = dx0 - 0 - 2 * SQUISH_CONSTANT_3D; 0895 dy6 = dy4; 0896 dz6 = dz5; 0897 attn6 = 2 - dx6 * dx6 - dy6 * dy6 - dz6 * dz6; 0898 if (attn6 > 0) { 0899 attn6 *= attn6; 0900 value += attn6 * attn6 * extrapolate3(ctx, xsb + 0, ysb + 1, zsb + 1, dx6, dy6, dz6); 0901 } 0902 } 0903 0904 /* First extra vertex */ 0905 attn_ext0 = 2 - dx_ext0 * dx_ext0 - dy_ext0 * dy_ext0 - dz_ext0 * dz_ext0; 0906 if (attn_ext0 > 0) 0907 { 0908 attn_ext0 *= attn_ext0; 0909 value += attn_ext0 * attn_ext0 * extrapolate3(ctx, xsv_ext0, ysv_ext0, zsv_ext0, dx_ext0, dy_ext0, dz_ext0); 0910 } 0911 0912 /* Second extra vertex */ 0913 attn_ext1 = 2 - dx_ext1 * dx_ext1 - dy_ext1 * dy_ext1 - dz_ext1 * dz_ext1; 0914 if (attn_ext1 > 0) 0915 { 0916 attn_ext1 *= attn_ext1; 0917 value += attn_ext1 * attn_ext1 * extrapolate3(ctx, xsv_ext1, ysv_ext1, zsv_ext1, dx_ext1, dy_ext1, dz_ext1); 0918 } 0919 0920 return value / NORM_CONSTANT_3D; 0921 } 0922 0923 /* 0924 * 4D OpenSimplex (Simplectic) Noise. 0925 */ 0926 double open_simplex_noise4(struct osn_context *ctx, double x, double y, double z, double w) 0927 { 0928 double uins; 0929 double dx1, dy1, dz1, dw1; 0930 double dx2, dy2, dz2, dw2; 0931 double dx3, dy3, dz3, dw3; 0932 double dx4, dy4, dz4, dw4; 0933 double dx5, dy5, dz5, dw5; 0934 double dx6, dy6, dz6, dw6; 0935 double dx7, dy7, dz7, dw7; 0936 double dx8, dy8, dz8, dw8; 0937 double dx9, dy9, dz9, dw9; 0938 double dx10, dy10, dz10, dw10; 0939 double attn0, attn1, attn2, attn3, attn4; 0940 double attn5, attn6, attn7, attn8, attn9, attn10; 0941 double attn_ext0, attn_ext1, attn_ext2; 0942 int8_t c, c1, c2; 0943 int8_t aPoint, bPoint; 0944 double aScore, bScore; 0945 int aIsBiggerSide; 0946 int bIsBiggerSide; 0947 double p1, p2, p3, p4; 0948 double score; 0949 0950 /* Place input coordinates on simplectic honeycomb. */ 0951 double stretchOffset = (x + y + z + w) * STRETCH_CONSTANT_4D; 0952 double xs = x + stretchOffset; 0953 double ys = y + stretchOffset; 0954 double zs = z + stretchOffset; 0955 double ws = w + stretchOffset; 0956 0957 /* Floor to get simplectic honeycomb coordinates of rhombo-hypercube super-cell origin. */ 0958 int xsb = fastFloor(xs); 0959 int ysb = fastFloor(ys); 0960 int zsb = fastFloor(zs); 0961 int wsb = fastFloor(ws); 0962 0963 /* Skew out to get actual coordinates of stretched rhombo-hypercube origin. We'll need these later. */ 0964 double squishOffset = (xsb + ysb + zsb + wsb) * SQUISH_CONSTANT_4D; 0965 double xb = xsb + squishOffset; 0966 double yb = ysb + squishOffset; 0967 double zb = zsb + squishOffset; 0968 double wb = wsb + squishOffset; 0969 0970 /* Compute simplectic honeycomb coordinates relative to rhombo-hypercube origin. */ 0971 double xins = xs - xsb; 0972 double yins = ys - ysb; 0973 double zins = zs - zsb; 0974 double wins = ws - wsb; 0975 0976 /* Sum those together to get a value that determines which region we're in. */ 0977 double inSum = xins + yins + zins + wins; 0978 0979 /* Positions relative to origin point. */ 0980 double dx0 = x - xb; 0981 double dy0 = y - yb; 0982 double dz0 = z - zb; 0983 double dw0 = w - wb; 0984 0985 /* We'll be defining these inside the next block and using them afterwards. */ 0986 double dx_ext0, dy_ext0, dz_ext0, dw_ext0; 0987 double dx_ext1, dy_ext1, dz_ext1, dw_ext1; 0988 double dx_ext2, dy_ext2, dz_ext2, dw_ext2; 0989 int xsv_ext0, ysv_ext0, zsv_ext0, wsv_ext0; 0990 int xsv_ext1, ysv_ext1, zsv_ext1, wsv_ext1; 0991 int xsv_ext2, ysv_ext2, zsv_ext2, wsv_ext2; 0992 0993 double value = 0; 0994 if (inSum <= 1) { /* We're inside the pentachoron (4-Simplex) at (0,0,0,0) */ 0995 0996 /* Determine which two of (0,0,0,1), (0,0,1,0), (0,1,0,0), (1,0,0,0) are closest. */ 0997 aPoint = 0x01; 0998 aScore = xins; 0999 bPoint = 0x02; 1000 bScore = yins; 1001 if (aScore >= bScore && zins > bScore) { 1002 bScore = zins; 1003 bPoint = 0x04; 1004 } else if (aScore < bScore && zins > aScore) { 1005 aScore = zins; 1006 aPoint = 0x04; 1007 } 1008 if (aScore >= bScore && wins > bScore) { 1009 bScore = wins; 1010 bPoint = 0x08; 1011 } else if (aScore < bScore && wins > aScore) { 1012 aScore = wins; 1013 aPoint = 0x08; 1014 } 1015 1016 /* Now we determine the three lattice points not part of the pentachoron that may contribute. 1017 This depends on the closest two pentachoron vertices, including (0,0,0,0) */ 1018 uins = 1 - inSum; 1019 if (uins > aScore || uins > bScore) { /* (0,0,0,0) is one of the closest two pentachoron vertices. */ 1020 c = (bScore > aScore ? bPoint : aPoint); /* Our other closest vertex is the closest out of a and b. */ 1021 if ((c & 0x01) == 0) { 1022 xsv_ext0 = xsb - 1; 1023 xsv_ext1 = xsv_ext2 = xsb; 1024 dx_ext0 = dx0 + 1; 1025 dx_ext1 = dx_ext2 = dx0; 1026 } else { 1027 xsv_ext0 = xsv_ext1 = xsv_ext2 = xsb + 1; 1028 dx_ext0 = dx_ext1 = dx_ext2 = dx0 - 1; 1029 } 1030 1031 if ((c & 0x02) == 0) { 1032 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb; 1033 dy_ext0 = dy_ext1 = dy_ext2 = dy0; 1034 if ((c & 0x01) == 0x01) { 1035 ysv_ext0 -= 1; 1036 dy_ext0 += 1; 1037 } else { 1038 ysv_ext1 -= 1; 1039 dy_ext1 += 1; 1040 } 1041 } else { 1042 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb + 1; 1043 dy_ext0 = dy_ext1 = dy_ext2 = dy0 - 1; 1044 } 1045 1046 if ((c & 0x04) == 0) { 1047 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb; 1048 dz_ext0 = dz_ext1 = dz_ext2 = dz0; 1049 if ((c & 0x03) != 0) { 1050 if ((c & 0x03) == 0x03) { 1051 zsv_ext0 -= 1; 1052 dz_ext0 += 1; 1053 } else { 1054 zsv_ext1 -= 1; 1055 dz_ext1 += 1; 1056 } 1057 } else { 1058 zsv_ext2 -= 1; 1059 dz_ext2 += 1; 1060 } 1061 } else { 1062 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb + 1; 1063 dz_ext0 = dz_ext1 = dz_ext2 = dz0 - 1; 1064 } 1065 1066 if ((c & 0x08) == 0) { 1067 wsv_ext0 = wsv_ext1 = wsb; 1068 wsv_ext2 = wsb - 1; 1069 dw_ext0 = dw_ext1 = dw0; 1070 dw_ext2 = dw0 + 1; 1071 } else { 1072 wsv_ext0 = wsv_ext1 = wsv_ext2 = wsb + 1; 1073 dw_ext0 = dw_ext1 = dw_ext2 = dw0 - 1; 1074 } 1075 } else { /* (0,0,0,0) is not one of the closest two pentachoron vertices. */ 1076 c = (int8_t)(aPoint | bPoint); /* Our three extra vertices are determined by the closest two. */ 1077 1078 if ((c & 0x01) == 0) { 1079 xsv_ext0 = xsv_ext2 = xsb; 1080 xsv_ext1 = xsb - 1; 1081 dx_ext0 = dx0 - 2 * SQUISH_CONSTANT_4D; 1082 dx_ext1 = dx0 + 1 - SQUISH_CONSTANT_4D; 1083 dx_ext2 = dx0 - SQUISH_CONSTANT_4D; 1084 } else { 1085 xsv_ext0 = xsv_ext1 = xsv_ext2 = xsb + 1; 1086 dx_ext0 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1087 dx_ext1 = dx_ext2 = dx0 - 1 - SQUISH_CONSTANT_4D; 1088 } 1089 1090 if ((c & 0x02) == 0) { 1091 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb; 1092 dy_ext0 = dy0 - 2 * SQUISH_CONSTANT_4D; 1093 dy_ext1 = dy_ext2 = dy0 - SQUISH_CONSTANT_4D; 1094 if ((c & 0x01) == 0x01) { 1095 ysv_ext1 -= 1; 1096 dy_ext1 += 1; 1097 } else { 1098 ysv_ext2 -= 1; 1099 dy_ext2 += 1; 1100 } 1101 } else { 1102 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb + 1; 1103 dy_ext0 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1104 dy_ext1 = dy_ext2 = dy0 - 1 - SQUISH_CONSTANT_4D; 1105 } 1106 1107 if ((c & 0x04) == 0) { 1108 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb; 1109 dz_ext0 = dz0 - 2 * SQUISH_CONSTANT_4D; 1110 dz_ext1 = dz_ext2 = dz0 - SQUISH_CONSTANT_4D; 1111 if ((c & 0x03) == 0x03) { 1112 zsv_ext1 -= 1; 1113 dz_ext1 += 1; 1114 } else { 1115 zsv_ext2 -= 1; 1116 dz_ext2 += 1; 1117 } 1118 } else { 1119 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb + 1; 1120 dz_ext0 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1121 dz_ext1 = dz_ext2 = dz0 - 1 - SQUISH_CONSTANT_4D; 1122 } 1123 1124 if ((c & 0x08) == 0) { 1125 wsv_ext0 = wsv_ext1 = wsb; 1126 wsv_ext2 = wsb - 1; 1127 dw_ext0 = dw0 - 2 * SQUISH_CONSTANT_4D; 1128 dw_ext1 = dw0 - SQUISH_CONSTANT_4D; 1129 dw_ext2 = dw0 + 1 - SQUISH_CONSTANT_4D; 1130 } else { 1131 wsv_ext0 = wsv_ext1 = wsv_ext2 = wsb + 1; 1132 dw_ext0 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1133 dw_ext1 = dw_ext2 = dw0 - 1 - SQUISH_CONSTANT_4D; 1134 } 1135 } 1136 1137 /* Contribution (0,0,0,0) */ 1138 attn0 = 2 - dx0 * dx0 - dy0 * dy0 - dz0 * dz0 - dw0 * dw0; 1139 if (attn0 > 0) { 1140 attn0 *= attn0; 1141 value += attn0 * attn0 * extrapolate4(ctx, xsb + 0, ysb + 0, zsb + 0, wsb + 0, dx0, dy0, dz0, dw0); 1142 } 1143 1144 /* Contribution (1,0,0,0) */ 1145 dx1 = dx0 - 1 - SQUISH_CONSTANT_4D; 1146 dy1 = dy0 - 0 - SQUISH_CONSTANT_4D; 1147 dz1 = dz0 - 0 - SQUISH_CONSTANT_4D; 1148 dw1 = dw0 - 0 - SQUISH_CONSTANT_4D; 1149 attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1 - dw1 * dw1; 1150 if (attn1 > 0) { 1151 attn1 *= attn1; 1152 value += attn1 * attn1 * extrapolate4(ctx, xsb + 1, ysb + 0, zsb + 0, wsb + 0, dx1, dy1, dz1, dw1); 1153 } 1154 1155 /* Contribution (0,1,0,0) */ 1156 dx2 = dx0 - 0 - SQUISH_CONSTANT_4D; 1157 dy2 = dy0 - 1 - SQUISH_CONSTANT_4D; 1158 dz2 = dz1; 1159 dw2 = dw1; 1160 attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2 - dw2 * dw2; 1161 if (attn2 > 0) { 1162 attn2 *= attn2; 1163 value += attn2 * attn2 * extrapolate4(ctx, xsb + 0, ysb + 1, zsb + 0, wsb + 0, dx2, dy2, dz2, dw2); 1164 } 1165 1166 /* Contribution (0,0,1,0) */ 1167 dx3 = dx2; 1168 dy3 = dy1; 1169 dz3 = dz0 - 1 - SQUISH_CONSTANT_4D; 1170 dw3 = dw1; 1171 attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3 - dw3 * dw3; 1172 if (attn3 > 0) { 1173 attn3 *= attn3; 1174 value += attn3 * attn3 * extrapolate4(ctx, xsb + 0, ysb + 0, zsb + 1, wsb + 0, dx3, dy3, dz3, dw3); 1175 } 1176 1177 /* Contribution (0,0,0,1) */ 1178 dx4 = dx2; 1179 dy4 = dy1; 1180 dz4 = dz1; 1181 dw4 = dw0 - 1 - SQUISH_CONSTANT_4D; 1182 attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4 - dw4 * dw4; 1183 if (attn4 > 0) { 1184 attn4 *= attn4; 1185 value += attn4 * attn4 * extrapolate4(ctx, xsb + 0, ysb + 0, zsb + 0, wsb + 1, dx4, dy4, dz4, dw4); 1186 } 1187 } else if (inSum >= 3) { /* We're inside the pentachoron (4-Simplex) at (1,1,1,1) 1188 Determine which two of (1,1,1,0), (1,1,0,1), (1,0,1,1), (0,1,1,1) are closest. */ 1189 aPoint = 0x0E; 1190 aScore = xins; 1191 bPoint = 0x0D; 1192 bScore = yins; 1193 if (aScore <= bScore && zins < bScore) { 1194 bScore = zins; 1195 bPoint = 0x0B; 1196 } else if (aScore > bScore && zins < aScore) { 1197 aScore = zins; 1198 aPoint = 0x0B; 1199 } 1200 if (aScore <= bScore && wins < bScore) { 1201 bScore = wins; 1202 bPoint = 0x07; 1203 } else if (aScore > bScore && wins < aScore) { 1204 aScore = wins; 1205 aPoint = 0x07; 1206 } 1207 1208 /* Now we determine the three lattice points not part of the pentachoron that may contribute. 1209 This depends on the closest two pentachoron vertices, including (0,0,0,0) */ 1210 uins = 4 - inSum; 1211 if (uins < aScore || uins < bScore) { /* (1,1,1,1) is one of the closest two pentachoron vertices. */ 1212 c = (bScore < aScore ? bPoint : aPoint); /* Our other closest vertex is the closest out of a and b. */ 1213 1214 if ((c & 0x01) != 0) { 1215 xsv_ext0 = xsb + 2; 1216 xsv_ext1 = xsv_ext2 = xsb + 1; 1217 dx_ext0 = dx0 - 2 - 4 * SQUISH_CONSTANT_4D; 1218 dx_ext1 = dx_ext2 = dx0 - 1 - 4 * SQUISH_CONSTANT_4D; 1219 } else { 1220 xsv_ext0 = xsv_ext1 = xsv_ext2 = xsb; 1221 dx_ext0 = dx_ext1 = dx_ext2 = dx0 - 4 * SQUISH_CONSTANT_4D; 1222 } 1223 1224 if ((c & 0x02) != 0) { 1225 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb + 1; 1226 dy_ext0 = dy_ext1 = dy_ext2 = dy0 - 1 - 4 * SQUISH_CONSTANT_4D; 1227 if ((c & 0x01) != 0) { 1228 ysv_ext1 += 1; 1229 dy_ext1 -= 1; 1230 } else { 1231 ysv_ext0 += 1; 1232 dy_ext0 -= 1; 1233 } 1234 } else { 1235 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb; 1236 dy_ext0 = dy_ext1 = dy_ext2 = dy0 - 4 * SQUISH_CONSTANT_4D; 1237 } 1238 1239 if ((c & 0x04) != 0) { 1240 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb + 1; 1241 dz_ext0 = dz_ext1 = dz_ext2 = dz0 - 1 - 4 * SQUISH_CONSTANT_4D; 1242 if ((c & 0x03) != 0x03) { 1243 if ((c & 0x03) == 0) { 1244 zsv_ext0 += 1; 1245 dz_ext0 -= 1; 1246 } else { 1247 zsv_ext1 += 1; 1248 dz_ext1 -= 1; 1249 } 1250 } else { 1251 zsv_ext2 += 1; 1252 dz_ext2 -= 1; 1253 } 1254 } else { 1255 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb; 1256 dz_ext0 = dz_ext1 = dz_ext2 = dz0 - 4 * SQUISH_CONSTANT_4D; 1257 } 1258 1259 if ((c & 0x08) != 0) { 1260 wsv_ext0 = wsv_ext1 = wsb + 1; 1261 wsv_ext2 = wsb + 2; 1262 dw_ext0 = dw_ext1 = dw0 - 1 - 4 * SQUISH_CONSTANT_4D; 1263 dw_ext2 = dw0 - 2 - 4 * SQUISH_CONSTANT_4D; 1264 } else { 1265 wsv_ext0 = wsv_ext1 = wsv_ext2 = wsb; 1266 dw_ext0 = dw_ext1 = dw_ext2 = dw0 - 4 * SQUISH_CONSTANT_4D; 1267 } 1268 } else { /* (1,1,1,1) is not one of the closest two pentachoron vertices. */ 1269 c = (int8_t)(aPoint & bPoint); /* Our three extra vertices are determined by the closest two. */ 1270 1271 if ((c & 0x01) != 0) { 1272 xsv_ext0 = xsv_ext2 = xsb + 1; 1273 xsv_ext1 = xsb + 2; 1274 dx_ext0 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1275 dx_ext1 = dx0 - 2 - 3 * SQUISH_CONSTANT_4D; 1276 dx_ext2 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; 1277 } else { 1278 xsv_ext0 = xsv_ext1 = xsv_ext2 = xsb; 1279 dx_ext0 = dx0 - 2 * SQUISH_CONSTANT_4D; 1280 dx_ext1 = dx_ext2 = dx0 - 3 * SQUISH_CONSTANT_4D; 1281 } 1282 1283 if ((c & 0x02) != 0) { 1284 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb + 1; 1285 dy_ext0 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1286 dy_ext1 = dy_ext2 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; 1287 if ((c & 0x01) != 0) { 1288 ysv_ext2 += 1; 1289 dy_ext2 -= 1; 1290 } else { 1291 ysv_ext1 += 1; 1292 dy_ext1 -= 1; 1293 } 1294 } else { 1295 ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb; 1296 dy_ext0 = dy0 - 2 * SQUISH_CONSTANT_4D; 1297 dy_ext1 = dy_ext2 = dy0 - 3 * SQUISH_CONSTANT_4D; 1298 } 1299 1300 if ((c & 0x04) != 0) { 1301 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb + 1; 1302 dz_ext0 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1303 dz_ext1 = dz_ext2 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; 1304 if ((c & 0x03) != 0) { 1305 zsv_ext2 += 1; 1306 dz_ext2 -= 1; 1307 } else { 1308 zsv_ext1 += 1; 1309 dz_ext1 -= 1; 1310 } 1311 } else { 1312 zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb; 1313 dz_ext0 = dz0 - 2 * SQUISH_CONSTANT_4D; 1314 dz_ext1 = dz_ext2 = dz0 - 3 * SQUISH_CONSTANT_4D; 1315 } 1316 1317 if ((c & 0x08) != 0) { 1318 wsv_ext0 = wsv_ext1 = wsb + 1; 1319 wsv_ext2 = wsb + 2; 1320 dw_ext0 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1321 dw_ext1 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; 1322 dw_ext2 = dw0 - 2 - 3 * SQUISH_CONSTANT_4D; 1323 } else { 1324 wsv_ext0 = wsv_ext1 = wsv_ext2 = wsb; 1325 dw_ext0 = dw0 - 2 * SQUISH_CONSTANT_4D; 1326 dw_ext1 = dw_ext2 = dw0 - 3 * SQUISH_CONSTANT_4D; 1327 } 1328 } 1329 1330 /* Contribution (1,1,1,0) */ 1331 dx4 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; 1332 dy4 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; 1333 dz4 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; 1334 dw4 = dw0 - 3 * SQUISH_CONSTANT_4D; 1335 attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4 - dw4 * dw4; 1336 if (attn4 > 0) { 1337 attn4 *= attn4; 1338 value += attn4 * attn4 * extrapolate4(ctx, xsb + 1, ysb + 1, zsb + 1, wsb + 0, dx4, dy4, dz4, dw4); 1339 } 1340 1341 /* Contribution (1,1,0,1) */ 1342 dx3 = dx4; 1343 dy3 = dy4; 1344 dz3 = dz0 - 3 * SQUISH_CONSTANT_4D; 1345 dw3 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; 1346 attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3 - dw3 * dw3; 1347 if (attn3 > 0) { 1348 attn3 *= attn3; 1349 value += attn3 * attn3 * extrapolate4(ctx, xsb + 1, ysb + 1, zsb + 0, wsb + 1, dx3, dy3, dz3, dw3); 1350 } 1351 1352 /* Contribution (1,0,1,1) */ 1353 dx2 = dx4; 1354 dy2 = dy0 - 3 * SQUISH_CONSTANT_4D; 1355 dz2 = dz4; 1356 dw2 = dw3; 1357 attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2 - dw2 * dw2; 1358 if (attn2 > 0) { 1359 attn2 *= attn2; 1360 value += attn2 * attn2 * extrapolate4(ctx, xsb + 1, ysb + 0, zsb + 1, wsb + 1, dx2, dy2, dz2, dw2); 1361 } 1362 1363 /* Contribution (0,1,1,1) */ 1364 dx1 = dx0 - 3 * SQUISH_CONSTANT_4D; 1365 dz1 = dz4; 1366 dy1 = dy4; 1367 dw1 = dw3; 1368 attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1 - dw1 * dw1; 1369 if (attn1 > 0) { 1370 attn1 *= attn1; 1371 value += attn1 * attn1 * extrapolate4(ctx, xsb + 0, ysb + 1, zsb + 1, wsb + 1, dx1, dy1, dz1, dw1); 1372 } 1373 1374 /* Contribution (1,1,1,1) */ 1375 dx0 = dx0 - 1 - 4 * SQUISH_CONSTANT_4D; 1376 dy0 = dy0 - 1 - 4 * SQUISH_CONSTANT_4D; 1377 dz0 = dz0 - 1 - 4 * SQUISH_CONSTANT_4D; 1378 dw0 = dw0 - 1 - 4 * SQUISH_CONSTANT_4D; 1379 attn0 = 2 - dx0 * dx0 - dy0 * dy0 - dz0 * dz0 - dw0 * dw0; 1380 if (attn0 > 0) { 1381 attn0 *= attn0; 1382 value += attn0 * attn0 * extrapolate4(ctx, xsb + 1, ysb + 1, zsb + 1, wsb + 1, dx0, dy0, dz0, dw0); 1383 } 1384 } else if (inSum <= 2) { /* We're inside the first dispentachoron (Rectified 4-Simplex) */ 1385 aIsBiggerSide = 1; 1386 bIsBiggerSide = 1; 1387 1388 /* Decide between (1,1,0,0) and (0,0,1,1) */ 1389 if (xins + yins > zins + wins) { 1390 aScore = xins + yins; 1391 aPoint = 0x03; 1392 } else { 1393 aScore = zins + wins; 1394 aPoint = 0x0C; 1395 } 1396 1397 /* Decide between (1,0,1,0) and (0,1,0,1) */ 1398 if (xins + zins > yins + wins) { 1399 bScore = xins + zins; 1400 bPoint = 0x05; 1401 } else { 1402 bScore = yins + wins; 1403 bPoint = 0x0A; 1404 } 1405 1406 /* Closer between (1,0,0,1) and (0,1,1,0) will replace the further of a and b, if closer. */ 1407 if (xins + wins > yins + zins) { 1408 score = xins + wins; 1409 if (aScore >= bScore && score > bScore) { 1410 bScore = score; 1411 bPoint = 0x09; 1412 } else if (aScore < bScore && score > aScore) { 1413 aScore = score; 1414 aPoint = 0x09; 1415 } 1416 } else { 1417 score = yins + zins; 1418 if (aScore >= bScore && score > bScore) { 1419 bScore = score; 1420 bPoint = 0x06; 1421 } else if (aScore < bScore && score > aScore) { 1422 aScore = score; 1423 aPoint = 0x06; 1424 } 1425 } 1426 1427 /* Decide if (1,0,0,0) is closer. */ 1428 p1 = 2 - inSum + xins; 1429 if (aScore >= bScore && p1 > bScore) { 1430 bScore = p1; 1431 bPoint = 0x01; 1432 bIsBiggerSide = 0; 1433 } else if (aScore < bScore && p1 > aScore) { 1434 aScore = p1; 1435 aPoint = 0x01; 1436 aIsBiggerSide = 0; 1437 } 1438 1439 /* Decide if (0,1,0,0) is closer. */ 1440 p2 = 2 - inSum + yins; 1441 if (aScore >= bScore && p2 > bScore) { 1442 bScore = p2; 1443 bPoint = 0x02; 1444 bIsBiggerSide = 0; 1445 } else if (aScore < bScore && p2 > aScore) { 1446 aScore = p2; 1447 aPoint = 0x02; 1448 aIsBiggerSide = 0; 1449 } 1450 1451 /* Decide if (0,0,1,0) is closer. */ 1452 p3 = 2 - inSum + zins; 1453 if (aScore >= bScore && p3 > bScore) { 1454 bScore = p3; 1455 bPoint = 0x04; 1456 bIsBiggerSide = 0; 1457 } else if (aScore < bScore && p3 > aScore) { 1458 aScore = p3; 1459 aPoint = 0x04; 1460 aIsBiggerSide = 0; 1461 } 1462 1463 /* Decide if (0,0,0,1) is closer. */ 1464 p4 = 2 - inSum + wins; 1465 if (aScore >= bScore && p4 > bScore) { 1466 bScore = p4; 1467 bPoint = 0x08; 1468 bIsBiggerSide = 0; 1469 } else if (aScore < bScore && p4 > aScore) { 1470 aScore = p4; 1471 aPoint = 0x08; 1472 aIsBiggerSide = 0; 1473 } 1474 1475 /* Where each of the two closest points are determines how the extra three vertices are calculated. */ 1476 if (aIsBiggerSide == bIsBiggerSide) { 1477 if (aIsBiggerSide) { /* Both closest points on the bigger side */ 1478 c1 = (int8_t)(aPoint | bPoint); 1479 c2 = (int8_t)(aPoint & bPoint); 1480 if ((c1 & 0x01) == 0) { 1481 xsv_ext0 = xsb; 1482 xsv_ext1 = xsb - 1; 1483 dx_ext0 = dx0 - 3 * SQUISH_CONSTANT_4D; 1484 dx_ext1 = dx0 + 1 - 2 * SQUISH_CONSTANT_4D; 1485 } else { 1486 xsv_ext0 = xsv_ext1 = xsb + 1; 1487 dx_ext0 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; 1488 dx_ext1 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1489 } 1490 1491 if ((c1 & 0x02) == 0) { 1492 ysv_ext0 = ysb; 1493 ysv_ext1 = ysb - 1; 1494 dy_ext0 = dy0 - 3 * SQUISH_CONSTANT_4D; 1495 dy_ext1 = dy0 + 1 - 2 * SQUISH_CONSTANT_4D; 1496 } else { 1497 ysv_ext0 = ysv_ext1 = ysb + 1; 1498 dy_ext0 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; 1499 dy_ext1 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1500 } 1501 1502 if ((c1 & 0x04) == 0) { 1503 zsv_ext0 = zsb; 1504 zsv_ext1 = zsb - 1; 1505 dz_ext0 = dz0 - 3 * SQUISH_CONSTANT_4D; 1506 dz_ext1 = dz0 + 1 - 2 * SQUISH_CONSTANT_4D; 1507 } else { 1508 zsv_ext0 = zsv_ext1 = zsb + 1; 1509 dz_ext0 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; 1510 dz_ext1 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1511 } 1512 1513 if ((c1 & 0x08) == 0) { 1514 wsv_ext0 = wsb; 1515 wsv_ext1 = wsb - 1; 1516 dw_ext0 = dw0 - 3 * SQUISH_CONSTANT_4D; 1517 dw_ext1 = dw0 + 1 - 2 * SQUISH_CONSTANT_4D; 1518 } else { 1519 wsv_ext0 = wsv_ext1 = wsb + 1; 1520 dw_ext0 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; 1521 dw_ext1 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1522 } 1523 1524 /* One combination is a permutation of (0,0,0,2) based on c2 */ 1525 xsv_ext2 = xsb; 1526 ysv_ext2 = ysb; 1527 zsv_ext2 = zsb; 1528 wsv_ext2 = wsb; 1529 dx_ext2 = dx0 - 2 * SQUISH_CONSTANT_4D; 1530 dy_ext2 = dy0 - 2 * SQUISH_CONSTANT_4D; 1531 dz_ext2 = dz0 - 2 * SQUISH_CONSTANT_4D; 1532 dw_ext2 = dw0 - 2 * SQUISH_CONSTANT_4D; 1533 if ((c2 & 0x01) != 0) { 1534 xsv_ext2 += 2; 1535 dx_ext2 -= 2; 1536 } else if ((c2 & 0x02) != 0) { 1537 ysv_ext2 += 2; 1538 dy_ext2 -= 2; 1539 } else if ((c2 & 0x04) != 0) { 1540 zsv_ext2 += 2; 1541 dz_ext2 -= 2; 1542 } else { 1543 wsv_ext2 += 2; 1544 dw_ext2 -= 2; 1545 } 1546 1547 } else { /* Both closest points on the smaller side */ 1548 /* One of the two extra points is (0,0,0,0) */ 1549 xsv_ext2 = xsb; 1550 ysv_ext2 = ysb; 1551 zsv_ext2 = zsb; 1552 wsv_ext2 = wsb; 1553 dx_ext2 = dx0; 1554 dy_ext2 = dy0; 1555 dz_ext2 = dz0; 1556 dw_ext2 = dw0; 1557 1558 /* Other two points are based on the omitted axes. */ 1559 c = (int8_t)(aPoint | bPoint); 1560 1561 if ((c & 0x01) == 0) { 1562 xsv_ext0 = xsb - 1; 1563 xsv_ext1 = xsb; 1564 dx_ext0 = dx0 + 1 - SQUISH_CONSTANT_4D; 1565 dx_ext1 = dx0 - SQUISH_CONSTANT_4D; 1566 } else { 1567 xsv_ext0 = xsv_ext1 = xsb + 1; 1568 dx_ext0 = dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_4D; 1569 } 1570 1571 if ((c & 0x02) == 0) { 1572 ysv_ext0 = ysv_ext1 = ysb; 1573 dy_ext0 = dy_ext1 = dy0 - SQUISH_CONSTANT_4D; 1574 if ((c & 0x01) == 0x01) 1575 { 1576 ysv_ext0 -= 1; 1577 dy_ext0 += 1; 1578 } else { 1579 ysv_ext1 -= 1; 1580 dy_ext1 += 1; 1581 } 1582 } else { 1583 ysv_ext0 = ysv_ext1 = ysb + 1; 1584 dy_ext0 = dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_4D; 1585 } 1586 1587 if ((c & 0x04) == 0) { 1588 zsv_ext0 = zsv_ext1 = zsb; 1589 dz_ext0 = dz_ext1 = dz0 - SQUISH_CONSTANT_4D; 1590 if ((c & 0x03) == 0x03) 1591 { 1592 zsv_ext0 -= 1; 1593 dz_ext0 += 1; 1594 } else { 1595 zsv_ext1 -= 1; 1596 dz_ext1 += 1; 1597 } 1598 } else { 1599 zsv_ext0 = zsv_ext1 = zsb + 1; 1600 dz_ext0 = dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_4D; 1601 } 1602 1603 if ((c & 0x08) == 0) 1604 { 1605 wsv_ext0 = wsb; 1606 wsv_ext1 = wsb - 1; 1607 dw_ext0 = dw0 - SQUISH_CONSTANT_4D; 1608 dw_ext1 = dw0 + 1 - SQUISH_CONSTANT_4D; 1609 } else { 1610 wsv_ext0 = wsv_ext1 = wsb + 1; 1611 dw_ext0 = dw_ext1 = dw0 - 1 - SQUISH_CONSTANT_4D; 1612 } 1613 1614 } 1615 } else { /* One point on each "side" */ 1616 if (aIsBiggerSide) { 1617 c1 = aPoint; 1618 c2 = bPoint; 1619 } else { 1620 c1 = bPoint; 1621 c2 = aPoint; 1622 } 1623 1624 /* Two contributions are the bigger-sided point with each 0 replaced with -1. */ 1625 if ((c1 & 0x01) == 0) { 1626 xsv_ext0 = xsb - 1; 1627 xsv_ext1 = xsb; 1628 dx_ext0 = dx0 + 1 - SQUISH_CONSTANT_4D; 1629 dx_ext1 = dx0 - SQUISH_CONSTANT_4D; 1630 } else { 1631 xsv_ext0 = xsv_ext1 = xsb + 1; 1632 dx_ext0 = dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_4D; 1633 } 1634 1635 if ((c1 & 0x02) == 0) { 1636 ysv_ext0 = ysv_ext1 = ysb; 1637 dy_ext0 = dy_ext1 = dy0 - SQUISH_CONSTANT_4D; 1638 if ((c1 & 0x01) == 0x01) { 1639 ysv_ext0 -= 1; 1640 dy_ext0 += 1; 1641 } else { 1642 ysv_ext1 -= 1; 1643 dy_ext1 += 1; 1644 } 1645 } else { 1646 ysv_ext0 = ysv_ext1 = ysb + 1; 1647 dy_ext0 = dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_4D; 1648 } 1649 1650 if ((c1 & 0x04) == 0) { 1651 zsv_ext0 = zsv_ext1 = zsb; 1652 dz_ext0 = dz_ext1 = dz0 - SQUISH_CONSTANT_4D; 1653 if ((c1 & 0x03) == 0x03) { 1654 zsv_ext0 -= 1; 1655 dz_ext0 += 1; 1656 } else { 1657 zsv_ext1 -= 1; 1658 dz_ext1 += 1; 1659 } 1660 } else { 1661 zsv_ext0 = zsv_ext1 = zsb + 1; 1662 dz_ext0 = dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_4D; 1663 } 1664 1665 if ((c1 & 0x08) == 0) { 1666 wsv_ext0 = wsb; 1667 wsv_ext1 = wsb - 1; 1668 dw_ext0 = dw0 - SQUISH_CONSTANT_4D; 1669 dw_ext1 = dw0 + 1 - SQUISH_CONSTANT_4D; 1670 } else { 1671 wsv_ext0 = wsv_ext1 = wsb + 1; 1672 dw_ext0 = dw_ext1 = dw0 - 1 - SQUISH_CONSTANT_4D; 1673 } 1674 1675 /* One contribution is a permutation of (0,0,0,2) based on the smaller-sided point */ 1676 xsv_ext2 = xsb; 1677 ysv_ext2 = ysb; 1678 zsv_ext2 = zsb; 1679 wsv_ext2 = wsb; 1680 dx_ext2 = dx0 - 2 * SQUISH_CONSTANT_4D; 1681 dy_ext2 = dy0 - 2 * SQUISH_CONSTANT_4D; 1682 dz_ext2 = dz0 - 2 * SQUISH_CONSTANT_4D; 1683 dw_ext2 = dw0 - 2 * SQUISH_CONSTANT_4D; 1684 if ((c2 & 0x01) != 0) { 1685 xsv_ext2 += 2; 1686 dx_ext2 -= 2; 1687 } else if ((c2 & 0x02) != 0) { 1688 ysv_ext2 += 2; 1689 dy_ext2 -= 2; 1690 } else if ((c2 & 0x04) != 0) { 1691 zsv_ext2 += 2; 1692 dz_ext2 -= 2; 1693 } else { 1694 wsv_ext2 += 2; 1695 dw_ext2 -= 2; 1696 } 1697 } 1698 1699 /* Contribution (1,0,0,0) */ 1700 dx1 = dx0 - 1 - SQUISH_CONSTANT_4D; 1701 dy1 = dy0 - 0 - SQUISH_CONSTANT_4D; 1702 dz1 = dz0 - 0 - SQUISH_CONSTANT_4D; 1703 dw1 = dw0 - 0 - SQUISH_CONSTANT_4D; 1704 attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1 - dw1 * dw1; 1705 if (attn1 > 0) { 1706 attn1 *= attn1; 1707 value += attn1 * attn1 * extrapolate4(ctx, xsb + 1, ysb + 0, zsb + 0, wsb + 0, dx1, dy1, dz1, dw1); 1708 } 1709 1710 /* Contribution (0,1,0,0) */ 1711 dx2 = dx0 - 0 - SQUISH_CONSTANT_4D; 1712 dy2 = dy0 - 1 - SQUISH_CONSTANT_4D; 1713 dz2 = dz1; 1714 dw2 = dw1; 1715 attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2 - dw2 * dw2; 1716 if (attn2 > 0) { 1717 attn2 *= attn2; 1718 value += attn2 * attn2 * extrapolate4(ctx, xsb + 0, ysb + 1, zsb + 0, wsb + 0, dx2, dy2, dz2, dw2); 1719 } 1720 1721 /* Contribution (0,0,1,0) */ 1722 dx3 = dx2; 1723 dy3 = dy1; 1724 dz3 = dz0 - 1 - SQUISH_CONSTANT_4D; 1725 dw3 = dw1; 1726 attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3 - dw3 * dw3; 1727 if (attn3 > 0) { 1728 attn3 *= attn3; 1729 value += attn3 * attn3 * extrapolate4(ctx, xsb + 0, ysb + 0, zsb + 1, wsb + 0, dx3, dy3, dz3, dw3); 1730 } 1731 1732 /* Contribution (0,0,0,1) */ 1733 dx4 = dx2; 1734 dy4 = dy1; 1735 dz4 = dz1; 1736 dw4 = dw0 - 1 - SQUISH_CONSTANT_4D; 1737 attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4 - dw4 * dw4; 1738 if (attn4 > 0) { 1739 attn4 *= attn4; 1740 value += attn4 * attn4 * extrapolate4(ctx, xsb + 0, ysb + 0, zsb + 0, wsb + 1, dx4, dy4, dz4, dw4); 1741 } 1742 1743 /* Contribution (1,1,0,0) */ 1744 dx5 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1745 dy5 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1746 dz5 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; 1747 dw5 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; 1748 attn5 = 2 - dx5 * dx5 - dy5 * dy5 - dz5 * dz5 - dw5 * dw5; 1749 if (attn5 > 0) { 1750 attn5 *= attn5; 1751 value += attn5 * attn5 * extrapolate4(ctx, xsb + 1, ysb + 1, zsb + 0, wsb + 0, dx5, dy5, dz5, dw5); 1752 } 1753 1754 /* Contribution (1,0,1,0) */ 1755 dx6 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1756 dy6 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; 1757 dz6 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1758 dw6 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; 1759 attn6 = 2 - dx6 * dx6 - dy6 * dy6 - dz6 * dz6 - dw6 * dw6; 1760 if (attn6 > 0) { 1761 attn6 *= attn6; 1762 value += attn6 * attn6 * extrapolate4(ctx, xsb + 1, ysb + 0, zsb + 1, wsb + 0, dx6, dy6, dz6, dw6); 1763 } 1764 1765 /* Contribution (1,0,0,1) */ 1766 dx7 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1767 dy7 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; 1768 dz7 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; 1769 dw7 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1770 attn7 = 2 - dx7 * dx7 - dy7 * dy7 - dz7 * dz7 - dw7 * dw7; 1771 if (attn7 > 0) { 1772 attn7 *= attn7; 1773 value += attn7 * attn7 * extrapolate4(ctx, xsb + 1, ysb + 0, zsb + 0, wsb + 1, dx7, dy7, dz7, dw7); 1774 } 1775 1776 /* Contribution (0,1,1,0) */ 1777 dx8 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; 1778 dy8 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1779 dz8 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1780 dw8 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; 1781 attn8 = 2 - dx8 * dx8 - dy8 * dy8 - dz8 * dz8 - dw8 * dw8; 1782 if (attn8 > 0) { 1783 attn8 *= attn8; 1784 value += attn8 * attn8 * extrapolate4(ctx, xsb + 0, ysb + 1, zsb + 1, wsb + 0, dx8, dy8, dz8, dw8); 1785 } 1786 1787 /* Contribution (0,1,0,1) */ 1788 dx9 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; 1789 dy9 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1790 dz9 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; 1791 dw9 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1792 attn9 = 2 - dx9 * dx9 - dy9 * dy9 - dz9 * dz9 - dw9 * dw9; 1793 if (attn9 > 0) { 1794 attn9 *= attn9; 1795 value += attn9 * attn9 * extrapolate4(ctx, xsb + 0, ysb + 1, zsb + 0, wsb + 1, dx9, dy9, dz9, dw9); 1796 } 1797 1798 /* Contribution (0,0,1,1) */ 1799 dx10 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; 1800 dy10 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; 1801 dz10 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1802 dw10 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1803 attn10 = 2 - dx10 * dx10 - dy10 * dy10 - dz10 * dz10 - dw10 * dw10; 1804 if (attn10 > 0) { 1805 attn10 *= attn10; 1806 value += attn10 * attn10 * extrapolate4(ctx, xsb + 0, ysb + 0, zsb + 1, wsb + 1, dx10, dy10, dz10, dw10); 1807 } 1808 } else { /* We're inside the second dispentachoron (Rectified 4-Simplex) */ 1809 aIsBiggerSide = 1; 1810 bIsBiggerSide = 1; 1811 1812 /* Decide between (0,0,1,1) and (1,1,0,0) */ 1813 if (xins + yins < zins + wins) { 1814 aScore = xins + yins; 1815 aPoint = 0x0C; 1816 } else { 1817 aScore = zins + wins; 1818 aPoint = 0x03; 1819 } 1820 1821 /* Decide between (0,1,0,1) and (1,0,1,0) */ 1822 if (xins + zins < yins + wins) { 1823 bScore = xins + zins; 1824 bPoint = 0x0A; 1825 } else { 1826 bScore = yins + wins; 1827 bPoint = 0x05; 1828 } 1829 1830 /* Closer between (0,1,1,0) and (1,0,0,1) will replace the further of a and b, if closer. */ 1831 if (xins + wins < yins + zins) { 1832 score = xins + wins; 1833 if (aScore <= bScore && score < bScore) { 1834 bScore = score; 1835 bPoint = 0x06; 1836 } else if (aScore > bScore && score < aScore) { 1837 aScore = score; 1838 aPoint = 0x06; 1839 } 1840 } else { 1841 score = yins + zins; 1842 if (aScore <= bScore && score < bScore) { 1843 bScore = score; 1844 bPoint = 0x09; 1845 } else if (aScore > bScore && score < aScore) { 1846 aScore = score; 1847 aPoint = 0x09; 1848 } 1849 } 1850 1851 /* Decide if (0,1,1,1) is closer. */ 1852 p1 = 3 - inSum + xins; 1853 if (aScore <= bScore && p1 < bScore) { 1854 bScore = p1; 1855 bPoint = 0x0E; 1856 bIsBiggerSide = 0; 1857 } else if (aScore > bScore && p1 < aScore) { 1858 aScore = p1; 1859 aPoint = 0x0E; 1860 aIsBiggerSide = 0; 1861 } 1862 1863 /* Decide if (1,0,1,1) is closer. */ 1864 p2 = 3 - inSum + yins; 1865 if (aScore <= bScore && p2 < bScore) { 1866 bScore = p2; 1867 bPoint = 0x0D; 1868 bIsBiggerSide = 0; 1869 } else if (aScore > bScore && p2 < aScore) { 1870 aScore = p2; 1871 aPoint = 0x0D; 1872 aIsBiggerSide = 0; 1873 } 1874 1875 /* Decide if (1,1,0,1) is closer. */ 1876 p3 = 3 - inSum + zins; 1877 if (aScore <= bScore && p3 < bScore) { 1878 bScore = p3; 1879 bPoint = 0x0B; 1880 bIsBiggerSide = 0; 1881 } else if (aScore > bScore && p3 < aScore) { 1882 aScore = p3; 1883 aPoint = 0x0B; 1884 aIsBiggerSide = 0; 1885 } 1886 1887 /* Decide if (1,1,1,0) is closer. */ 1888 p4 = 3 - inSum + wins; 1889 if (aScore <= bScore && p4 < bScore) { 1890 bScore = p4; 1891 bPoint = 0x07; 1892 bIsBiggerSide = 0; 1893 } else if (aScore > bScore && p4 < aScore) { 1894 aScore = p4; 1895 aPoint = 0x07; 1896 aIsBiggerSide = 0; 1897 } 1898 1899 /* Where each of the two closest points are determines how the extra three vertices are calculated. */ 1900 if (aIsBiggerSide == bIsBiggerSide) { 1901 if (aIsBiggerSide) { /* Both closest points on the bigger side */ 1902 c1 = (int8_t)(aPoint & bPoint); 1903 c2 = (int8_t)(aPoint | bPoint); 1904 1905 /* Two contributions are permutations of (0,0,0,1) and (0,0,0,2) based on c1 */ 1906 xsv_ext0 = xsv_ext1 = xsb; 1907 ysv_ext0 = ysv_ext1 = ysb; 1908 zsv_ext0 = zsv_ext1 = zsb; 1909 wsv_ext0 = wsv_ext1 = wsb; 1910 dx_ext0 = dx0 - SQUISH_CONSTANT_4D; 1911 dy_ext0 = dy0 - SQUISH_CONSTANT_4D; 1912 dz_ext0 = dz0 - SQUISH_CONSTANT_4D; 1913 dw_ext0 = dw0 - SQUISH_CONSTANT_4D; 1914 dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_4D; 1915 dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_4D; 1916 dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_4D; 1917 dw_ext1 = dw0 - 2 * SQUISH_CONSTANT_4D; 1918 if ((c1 & 0x01) != 0) { 1919 xsv_ext0 += 1; 1920 dx_ext0 -= 1; 1921 xsv_ext1 += 2; 1922 dx_ext1 -= 2; 1923 } else if ((c1 & 0x02) != 0) { 1924 ysv_ext0 += 1; 1925 dy_ext0 -= 1; 1926 ysv_ext1 += 2; 1927 dy_ext1 -= 2; 1928 } else if ((c1 & 0x04) != 0) { 1929 zsv_ext0 += 1; 1930 dz_ext0 -= 1; 1931 zsv_ext1 += 2; 1932 dz_ext1 -= 2; 1933 } else { 1934 wsv_ext0 += 1; 1935 dw_ext0 -= 1; 1936 wsv_ext1 += 2; 1937 dw_ext1 -= 2; 1938 } 1939 1940 /* One contribution is a permutation of (1,1,1,-1) based on c2 */ 1941 xsv_ext2 = xsb + 1; 1942 ysv_ext2 = ysb + 1; 1943 zsv_ext2 = zsb + 1; 1944 wsv_ext2 = wsb + 1; 1945 dx_ext2 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 1946 dy_ext2 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 1947 dz_ext2 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 1948 dw_ext2 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 1949 if ((c2 & 0x01) == 0) { 1950 xsv_ext2 -= 2; 1951 dx_ext2 += 2; 1952 } else if ((c2 & 0x02) == 0) { 1953 ysv_ext2 -= 2; 1954 dy_ext2 += 2; 1955 } else if ((c2 & 0x04) == 0) { 1956 zsv_ext2 -= 2; 1957 dz_ext2 += 2; 1958 } else { 1959 wsv_ext2 -= 2; 1960 dw_ext2 += 2; 1961 } 1962 } else { /* Both closest points on the smaller side */ 1963 /* One of the two extra points is (1,1,1,1) */ 1964 xsv_ext2 = xsb + 1; 1965 ysv_ext2 = ysb + 1; 1966 zsv_ext2 = zsb + 1; 1967 wsv_ext2 = wsb + 1; 1968 dx_ext2 = dx0 - 1 - 4 * SQUISH_CONSTANT_4D; 1969 dy_ext2 = dy0 - 1 - 4 * SQUISH_CONSTANT_4D; 1970 dz_ext2 = dz0 - 1 - 4 * SQUISH_CONSTANT_4D; 1971 dw_ext2 = dw0 - 1 - 4 * SQUISH_CONSTANT_4D; 1972 1973 /* Other two points are based on the shared axes. */ 1974 c = (int8_t)(aPoint & bPoint); 1975 1976 if ((c & 0x01) != 0) { 1977 xsv_ext0 = xsb + 2; 1978 xsv_ext1 = xsb + 1; 1979 dx_ext0 = dx0 - 2 - 3 * SQUISH_CONSTANT_4D; 1980 dx_ext1 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; 1981 } else { 1982 xsv_ext0 = xsv_ext1 = xsb; 1983 dx_ext0 = dx_ext1 = dx0 - 3 * SQUISH_CONSTANT_4D; 1984 } 1985 1986 if ((c & 0x02) != 0) { 1987 ysv_ext0 = ysv_ext1 = ysb + 1; 1988 dy_ext0 = dy_ext1 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; 1989 if ((c & 0x01) == 0) 1990 { 1991 ysv_ext0 += 1; 1992 dy_ext0 -= 1; 1993 } else { 1994 ysv_ext1 += 1; 1995 dy_ext1 -= 1; 1996 } 1997 } else { 1998 ysv_ext0 = ysv_ext1 = ysb; 1999 dy_ext0 = dy_ext1 = dy0 - 3 * SQUISH_CONSTANT_4D; 2000 } 2001 2002 if ((c & 0x04) != 0) { 2003 zsv_ext0 = zsv_ext1 = zsb + 1; 2004 dz_ext0 = dz_ext1 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; 2005 if ((c & 0x03) == 0) 2006 { 2007 zsv_ext0 += 1; 2008 dz_ext0 -= 1; 2009 } else { 2010 zsv_ext1 += 1; 2011 dz_ext1 -= 1; 2012 } 2013 } else { 2014 zsv_ext0 = zsv_ext1 = zsb; 2015 dz_ext0 = dz_ext1 = dz0 - 3 * SQUISH_CONSTANT_4D; 2016 } 2017 2018 if ((c & 0x08) != 0) 2019 { 2020 wsv_ext0 = wsb + 1; 2021 wsv_ext1 = wsb + 2; 2022 dw_ext0 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; 2023 dw_ext1 = dw0 - 2 - 3 * SQUISH_CONSTANT_4D; 2024 } else { 2025 wsv_ext0 = wsv_ext1 = wsb; 2026 dw_ext0 = dw_ext1 = dw0 - 3 * SQUISH_CONSTANT_4D; 2027 } 2028 } 2029 } else { /* One point on each "side" */ 2030 if (aIsBiggerSide) { 2031 c1 = aPoint; 2032 c2 = bPoint; 2033 } else { 2034 c1 = bPoint; 2035 c2 = aPoint; 2036 } 2037 2038 /* Two contributions are the bigger-sided point with each 1 replaced with 2. */ 2039 if ((c1 & 0x01) != 0) { 2040 xsv_ext0 = xsb + 2; 2041 xsv_ext1 = xsb + 1; 2042 dx_ext0 = dx0 - 2 - 3 * SQUISH_CONSTANT_4D; 2043 dx_ext1 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; 2044 } else { 2045 xsv_ext0 = xsv_ext1 = xsb; 2046 dx_ext0 = dx_ext1 = dx0 - 3 * SQUISH_CONSTANT_4D; 2047 } 2048 2049 if ((c1 & 0x02) != 0) { 2050 ysv_ext0 = ysv_ext1 = ysb + 1; 2051 dy_ext0 = dy_ext1 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; 2052 if ((c1 & 0x01) == 0) { 2053 ysv_ext0 += 1; 2054 dy_ext0 -= 1; 2055 } else { 2056 ysv_ext1 += 1; 2057 dy_ext1 -= 1; 2058 } 2059 } else { 2060 ysv_ext0 = ysv_ext1 = ysb; 2061 dy_ext0 = dy_ext1 = dy0 - 3 * SQUISH_CONSTANT_4D; 2062 } 2063 2064 if ((c1 & 0x04) != 0) { 2065 zsv_ext0 = zsv_ext1 = zsb + 1; 2066 dz_ext0 = dz_ext1 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; 2067 if ((c1 & 0x03) == 0) { 2068 zsv_ext0 += 1; 2069 dz_ext0 -= 1; 2070 } else { 2071 zsv_ext1 += 1; 2072 dz_ext1 -= 1; 2073 } 2074 } else { 2075 zsv_ext0 = zsv_ext1 = zsb; 2076 dz_ext0 = dz_ext1 = dz0 - 3 * SQUISH_CONSTANT_4D; 2077 } 2078 2079 if ((c1 & 0x08) != 0) { 2080 wsv_ext0 = wsb + 1; 2081 wsv_ext1 = wsb + 2; 2082 dw_ext0 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; 2083 dw_ext1 = dw0 - 2 - 3 * SQUISH_CONSTANT_4D; 2084 } else { 2085 wsv_ext0 = wsv_ext1 = wsb; 2086 dw_ext0 = dw_ext1 = dw0 - 3 * SQUISH_CONSTANT_4D; 2087 } 2088 2089 /* One contribution is a permutation of (1,1,1,-1) based on the smaller-sided point */ 2090 xsv_ext2 = xsb + 1; 2091 ysv_ext2 = ysb + 1; 2092 zsv_ext2 = zsb + 1; 2093 wsv_ext2 = wsb + 1; 2094 dx_ext2 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 2095 dy_ext2 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 2096 dz_ext2 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 2097 dw_ext2 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 2098 if ((c2 & 0x01) == 0) { 2099 xsv_ext2 -= 2; 2100 dx_ext2 += 2; 2101 } else if ((c2 & 0x02) == 0) { 2102 ysv_ext2 -= 2; 2103 dy_ext2 += 2; 2104 } else if ((c2 & 0x04) == 0) { 2105 zsv_ext2 -= 2; 2106 dz_ext2 += 2; 2107 } else { 2108 wsv_ext2 -= 2; 2109 dw_ext2 += 2; 2110 } 2111 } 2112 2113 /* Contribution (1,1,1,0) */ 2114 dx4 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; 2115 dy4 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; 2116 dz4 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; 2117 dw4 = dw0 - 3 * SQUISH_CONSTANT_4D; 2118 attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4 - dw4 * dw4; 2119 if (attn4 > 0) { 2120 attn4 *= attn4; 2121 value += attn4 * attn4 * extrapolate4(ctx, xsb + 1, ysb + 1, zsb + 1, wsb + 0, dx4, dy4, dz4, dw4); 2122 } 2123 2124 /* Contribution (1,1,0,1) */ 2125 dx3 = dx4; 2126 dy3 = dy4; 2127 dz3 = dz0 - 3 * SQUISH_CONSTANT_4D; 2128 dw3 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; 2129 attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3 - dw3 * dw3; 2130 if (attn3 > 0) { 2131 attn3 *= attn3; 2132 value += attn3 * attn3 * extrapolate4(ctx, xsb + 1, ysb + 1, zsb + 0, wsb + 1, dx3, dy3, dz3, dw3); 2133 } 2134 2135 /* Contribution (1,0,1,1) */ 2136 dx2 = dx4; 2137 dy2 = dy0 - 3 * SQUISH_CONSTANT_4D; 2138 dz2 = dz4; 2139 dw2 = dw3; 2140 attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2 - dw2 * dw2; 2141 if (attn2 > 0) { 2142 attn2 *= attn2; 2143 value += attn2 * attn2 * extrapolate4(ctx, xsb + 1, ysb + 0, zsb + 1, wsb + 1, dx2, dy2, dz2, dw2); 2144 } 2145 2146 /* Contribution (0,1,1,1) */ 2147 dx1 = dx0 - 3 * SQUISH_CONSTANT_4D; 2148 dz1 = dz4; 2149 dy1 = dy4; 2150 dw1 = dw3; 2151 attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1 - dw1 * dw1; 2152 if (attn1 > 0) { 2153 attn1 *= attn1; 2154 value += attn1 * attn1 * extrapolate4(ctx, xsb + 0, ysb + 1, zsb + 1, wsb + 1, dx1, dy1, dz1, dw1); 2155 } 2156 2157 /* Contribution (1,1,0,0) */ 2158 dx5 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 2159 dy5 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 2160 dz5 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; 2161 dw5 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; 2162 attn5 = 2 - dx5 * dx5 - dy5 * dy5 - dz5 * dz5 - dw5 * dw5; 2163 if (attn5 > 0) { 2164 attn5 *= attn5; 2165 value += attn5 * attn5 * extrapolate4(ctx, xsb + 1, ysb + 1, zsb + 0, wsb + 0, dx5, dy5, dz5, dw5); 2166 } 2167 2168 /* Contribution (1,0,1,0) */ 2169 dx6 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 2170 dy6 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; 2171 dz6 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 2172 dw6 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; 2173 attn6 = 2 - dx6 * dx6 - dy6 * dy6 - dz6 * dz6 - dw6 * dw6; 2174 if (attn6 > 0) { 2175 attn6 *= attn6; 2176 value += attn6 * attn6 * extrapolate4(ctx, xsb + 1, ysb + 0, zsb + 1, wsb + 0, dx6, dy6, dz6, dw6); 2177 } 2178 2179 /* Contribution (1,0,0,1) */ 2180 dx7 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; 2181 dy7 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; 2182 dz7 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; 2183 dw7 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 2184 attn7 = 2 - dx7 * dx7 - dy7 * dy7 - dz7 * dz7 - dw7 * dw7; 2185 if (attn7 > 0) { 2186 attn7 *= attn7; 2187 value += attn7 * attn7 * extrapolate4(ctx, xsb + 1, ysb + 0, zsb + 0, wsb + 1, dx7, dy7, dz7, dw7); 2188 } 2189 2190 /* Contribution (0,1,1,0) */ 2191 dx8 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; 2192 dy8 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 2193 dz8 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 2194 dw8 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; 2195 attn8 = 2 - dx8 * dx8 - dy8 * dy8 - dz8 * dz8 - dw8 * dw8; 2196 if (attn8 > 0) { 2197 attn8 *= attn8; 2198 value += attn8 * attn8 * extrapolate4(ctx, xsb + 0, ysb + 1, zsb + 1, wsb + 0, dx8, dy8, dz8, dw8); 2199 } 2200 2201 /* Contribution (0,1,0,1) */ 2202 dx9 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; 2203 dy9 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; 2204 dz9 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; 2205 dw9 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 2206 attn9 = 2 - dx9 * dx9 - dy9 * dy9 - dz9 * dz9 - dw9 * dw9; 2207 if (attn9 > 0) { 2208 attn9 *= attn9; 2209 value += attn9 * attn9 * extrapolate4(ctx, xsb + 0, ysb + 1, zsb + 0, wsb + 1, dx9, dy9, dz9, dw9); 2210 } 2211 2212 /* Contribution (0,0,1,1) */ 2213 dx10 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; 2214 dy10 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; 2215 dz10 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; 2216 dw10 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; 2217 attn10 = 2 - dx10 * dx10 - dy10 * dy10 - dz10 * dz10 - dw10 * dw10; 2218 if (attn10 > 0) { 2219 attn10 *= attn10; 2220 value += attn10 * attn10 * extrapolate4(ctx, xsb + 0, ysb + 0, zsb + 1, wsb + 1, dx10, dy10, dz10, dw10); 2221 } 2222 } 2223 2224 /* First extra vertex */ 2225 attn_ext0 = 2 - dx_ext0 * dx_ext0 - dy_ext0 * dy_ext0 - dz_ext0 * dz_ext0 - dw_ext0 * dw_ext0; 2226 if (attn_ext0 > 0) 2227 { 2228 attn_ext0 *= attn_ext0; 2229 value += attn_ext0 * attn_ext0 * extrapolate4(ctx, xsv_ext0, ysv_ext0, zsv_ext0, wsv_ext0, dx_ext0, dy_ext0, dz_ext0, dw_ext0); 2230 } 2231 2232 /* Second extra vertex */ 2233 attn_ext1 = 2 - dx_ext1 * dx_ext1 - dy_ext1 * dy_ext1 - dz_ext1 * dz_ext1 - dw_ext1 * dw_ext1; 2234 if (attn_ext1 > 0) 2235 { 2236 attn_ext1 *= attn_ext1; 2237 value += attn_ext1 * attn_ext1 * extrapolate4(ctx, xsv_ext1, ysv_ext1, zsv_ext1, wsv_ext1, dx_ext1, dy_ext1, dz_ext1, dw_ext1); 2238 } 2239 2240 /* Third extra vertex */ 2241 attn_ext2 = 2 - dx_ext2 * dx_ext2 - dy_ext2 * dy_ext2 - dz_ext2 * dz_ext2 - dw_ext2 * dw_ext2; 2242 if (attn_ext2 > 0) 2243 { 2244 attn_ext2 *= attn_ext2; 2245 value += attn_ext2 * attn_ext2 * extrapolate4(ctx, xsv_ext2, ysv_ext2, zsv_ext2, wsv_ext2, dx_ext2, dy_ext2, dz_ext2, dw_ext2); 2246 } 2247 2248 return value / NORM_CONSTANT_4D; 2249 } 2250