Krita Source Code Documentation
Loading...
Searching...
No Matches
open-simplex-noise.c
Go to the documentation of this file.
1/*
2 * OpenSimplex (Simplectic) Noise in C.
3 * Ported by Stephen M. Cameron from Kurt Spencer's java implementation
4 *
5 * v1.1 (October 5, 2014)
6 * - Added 2D and 4D implementations.
7 * - Proper gradient sets for all dimensions, from a
8 * dimensionally-generalizable scheme with an actual
9 * rhyme and reason behind it.
10 * - Removed default permutation array in favor of
11 * default seed.
12 * - Changed seed-based constructor to be independent
13 * of any particular randomization library, so results
14 * will be the same when ported to other languages.
15 */
16#include <math.h>
17#include <stdlib.h>
18#include <stdint.h>
19#include <string.h>
20#include <errno.h>
21
22#include "open-simplex-noise.h"
23
24#define STRETCH_CONSTANT_2D (-0.211324865405187) /* (1 / sqrt(2 + 1) - 1 ) / 2; */
25#define SQUISH_CONSTANT_2D (0.366025403784439) /* (sqrt(2 + 1) -1) / 2; */
26#define STRETCH_CONSTANT_3D (-1.0 / 6.0) /* (1 / sqrt(3 + 1) - 1) / 3; */
27#define SQUISH_CONSTANT_3D (1.0 / 3.0) /* (sqrt(3+1)-1)/3; */
28#define STRETCH_CONSTANT_4D (-0.138196601125011) /* (1 / sqrt(4 + 1) - 1) / 4; */
29#define SQUISH_CONSTANT_4D (0.309016994374947) /* (sqrt(4 + 1) - 1) / 4; */
30
31#define NORM_CONSTANT_2D (47.0)
32#define NORM_CONSTANT_3D (103.0)
33#define NORM_CONSTANT_4D (30.0)
34
35#define DEFAULT_SEED (0LL)
36
37struct osn_context {
38 int16_t *perm;
39 int16_t *permGradIndex3D;
40};
41
42#define ARRAYSIZE(x) (sizeof((x)) / sizeof((x)[0]))
43
44/*
45 * Gradients for 2D. They approximate the directions to the
46 * vertices of an octagon from the center.
47 */
48static const int8_t gradients2D[] = {
49 5, 2, 2, 5,
50 -5, 2, -2, 5,
51 5, -2, 2, -5,
52 -5, -2, -2, -5,
53};
54
55/*
56 * Gradients for 3D. They approximate the directions to the
57 * vertices of a rhombicuboctahedron from the center, skewed so
58 * that the triangular and square facets can be inscribed inside
59 * circles of the same radius.
60 */
61static const signed char gradients3D[] = {
62 -11, 4, 4, -4, 11, 4, -4, 4, 11,
63 11, 4, 4, 4, 11, 4, 4, 4, 11,
64 -11, -4, 4, -4, -11, 4, -4, -4, 11,
65 11, -4, 4, 4, -11, 4, 4, -4, 11,
66 -11, 4, -4, -4, 11, -4, -4, 4, -11,
67 11, 4, -4, 4, 11, -4, 4, 4, -11,
68 -11, -4, -4, -4, -11, -4, -4, -4, -11,
69 11, -4, -4, 4, -11, -4, 4, -4, -11,
70};
71
72/*
73 * Gradients for 4D. They approximate the directions to the
74 * vertices of a disprismatotesseractihexadecachoron from the center,
75 * skewed so that the tetrahedral and cubic facets can be inscribed inside
76 * spheres of the same radius.
77 */
78static const signed char gradients4D[] = {
79 3, 1, 1, 1, 1, 3, 1, 1, 1, 1, 3, 1, 1, 1, 1, 3,
80 -3, 1, 1, 1, -1, 3, 1, 1, -1, 1, 3, 1, -1, 1, 1, 3,
81 3, -1, 1, 1, 1, -3, 1, 1, 1, -1, 3, 1, 1, -1, 1, 3,
82 -3, -1, 1, 1, -1, -3, 1, 1, -1, -1, 3, 1, -1, -1, 1, 3,
83 3, 1, -1, 1, 1, 3, -1, 1, 1, 1, -3, 1, 1, 1, -1, 3,
84 -3, 1, -1, 1, -1, 3, -1, 1, -1, 1, -3, 1, -1, 1, -1, 3,
85 3, -1, -1, 1, 1, -3, -1, 1, 1, -1, -3, 1, 1, -1, -1, 3,
86 -3, -1, -1, 1, -1, -3, -1, 1, -1, -1, -3, 1, -1, -1, -1, 3,
87 3, 1, 1, -1, 1, 3, 1, -1, 1, 1, 3, -1, 1, 1, 1, -3,
88 -3, 1, 1, -1, -1, 3, 1, -1, -1, 1, 3, -1, -1, 1, 1, -3,
89 3, -1, 1, -1, 1, -3, 1, -1, 1, -1, 3, -1, 1, -1, 1, -3,
90 -3, -1, 1, -1, -1, -3, 1, -1, -1, -1, 3, -1, -1, -1, 1, -3,
91 3, 1, -1, -1, 1, 3, -1, -1, 1, 1, -3, -1, 1, 1, -1, -3,
92 -3, 1, -1, -1, -1, 3, -1, -1, -1, 1, -3, -1, -1, 1, -1, -3,
93 3, -1, -1, -1, 1, -3, -1, -1, 1, -1, -3, -1, 1, -1, -1, -3,
94 -3, -1, -1, -1, -1, -3, -1, -1, -1, -1, -3, -1, -1, -1, -1, -3,
95};
96
97static double extrapolate2(struct osn_context *ctx, int xsb, int ysb, double dx, double dy)
98{
99 int16_t *perm = ctx->perm;
100 int index = perm[(perm[xsb & 0xFF] + ysb) & 0xFF] & 0x0E;
101 return gradients2D[index] * dx
102 + gradients2D[index + 1] * dy;
103}
104
105static double extrapolate3(struct osn_context *ctx, int xsb, int ysb, int zsb, double dx, double dy, double dz)
106{
107 int16_t *perm = ctx->perm;
108 int16_t *permGradIndex3D = ctx->permGradIndex3D;
109 int index = permGradIndex3D[(perm[(perm[xsb & 0xFF] + ysb) & 0xFF] + zsb) & 0xFF];
110 return gradients3D[index] * dx
111 + gradients3D[index + 1] * dy
112 + gradients3D[index + 2] * dz;
113}
114
115static double extrapolate4(struct osn_context *ctx, int xsb, int ysb, int zsb, int wsb, double dx, double dy, double dz, double dw)
116{
117 int16_t *perm = ctx->perm;
118 int index = perm[(perm[(perm[(perm[xsb & 0xFF] + ysb) & 0xFF] + zsb) & 0xFF] + wsb) & 0xFF] & 0xFC;
119 return gradients4D[index] * dx
120 + gradients4D[index + 1] * dy
121 + gradients4D[index + 2] * dz
122 + gradients4D[index + 3] * dw;
123}
124
125static INLINE int fastFloor(double x) {
126 int xi = (int) x;
127 return x < xi ? xi - 1 : xi;
128}
129
130static int allocate_perm(struct osn_context *ctx, int nperm, int ngrad)
131{
132 if (ctx->perm)
133 free(ctx->perm);
134 if (ctx->permGradIndex3D)
135 free(ctx->permGradIndex3D);
136 ctx->perm = (int16_t *) malloc(sizeof(*ctx->perm) * nperm);
137 if (!ctx->perm)
138 return -ENOMEM;
139 ctx->permGradIndex3D = (int16_t *) malloc(sizeof(*ctx->permGradIndex3D) * ngrad);
140 if (!ctx->permGradIndex3D) {
141 free(ctx->perm);
142 return -ENOMEM;
143 }
144 return 0;
145}
146
147int open_simplex_noise_init_perm(struct osn_context *ctx, int16_t p[], int nelements)
148{
149 int i, rc;
150
151 rc = allocate_perm(ctx, nelements, 256);
152 if (rc)
153 return rc;
154 memcpy(ctx->perm, p, sizeof(*ctx->perm) * nelements);
155
156 for (i = 0; i < 256; i++) {
157 /* Since 3D has 24 gradients, simple bitmask won't work, so precompute modulo array. */
158 ctx->permGradIndex3D[i] = (int16_t)((ctx->perm[i] % (ARRAYSIZE(gradients3D) / 3)) * 3);
159 }
160 return 0;
161}
162
163/*
164 * Initializes using a permutation array generated from a 64-bit seed.
165 * Generates a proper permutation (i.e. doesn't merely perform N successive pair
166 * swaps on a base array). Uses a simple 64-bit LCG.
167 */
168int open_simplex_noise(int64_t seed, struct osn_context **ctx)
169{
170 int rc;
171 int16_t source[256];
172 int i;
173 int16_t *perm;
174 int16_t *permGradIndex3D;
175 int r;
176
177 *ctx = (struct osn_context *) malloc(sizeof(**ctx));
178 if (!(*ctx))
179 return -ENOMEM;
180 (*ctx)->perm = NULL;
181 (*ctx)->permGradIndex3D = NULL;
182
183 rc = allocate_perm(*ctx, 256, 256);
184 if (rc) {
185 free(*ctx);
186 return rc;
187 }
188
189 perm = (*ctx)->perm;
190 permGradIndex3D = (*ctx)->permGradIndex3D;
191
192 for (i = 0; i < 256; i++)
193 source[i] = (int16_t) i;
194 seed = seed * 6364136223846793005LL + 1442695040888963407LL;
195 seed = seed * 6364136223846793005LL + 1442695040888963407LL;
196 seed = seed * 6364136223846793005LL + 1442695040888963407LL;
197 for (i = 255; i >= 0; i--) {
198 seed = seed * 6364136223846793005LL + 1442695040888963407LL;
199 r = (int)((seed + 31) % (i + 1));
200 if (r < 0)
201 r += (i + 1);
202 perm[i] = source[r];
203 permGradIndex3D[i] = (short)((perm[i] % (ARRAYSIZE(gradients3D) / 3)) * 3);
204 source[r] = source[i];
205 }
206 return 0;
207}
208
209void open_simplex_noise_free(struct osn_context *ctx)
210{
211 if (!ctx)
212 return;
213 if (ctx->perm) {
214 free(ctx->perm);
215 ctx->perm = NULL;
216 }
217 if (ctx->permGradIndex3D) {
218 free(ctx->permGradIndex3D);
219 ctx->permGradIndex3D = NULL;
220 }
221 free(ctx);
222}
223
224/* 2D OpenSimplex (Simplectic) Noise. */
225double open_simplex_noise2(struct osn_context *ctx, double x, double y)
226{
227
228 /* Place input coordinates onto grid. */
229 double stretchOffset = (x + y) * STRETCH_CONSTANT_2D;
230 double xs = x + stretchOffset;
231 double ys = y + stretchOffset;
232
233 /* Floor to get grid coordinates of rhombus (stretched square) super-cell origin. */
234 int xsb = fastFloor(xs);
235 int ysb = fastFloor(ys);
236
237 /* Skew out to get actual coordinates of rhombus origin. We'll need these later. */
238 double squishOffset = (xsb + ysb) * SQUISH_CONSTANT_2D;
239 double xb = xsb + squishOffset;
240 double yb = ysb + squishOffset;
241
242 /* Compute grid coordinates relative to rhombus origin. */
243 double xins = xs - xsb;
244 double yins = ys - ysb;
245
246 /* Sum those together to get a value that determines which region we're in. */
247 double inSum = xins + yins;
248
249 /* Positions relative to origin point. */
250 double dx0 = x - xb;
251 double dy0 = y - yb;
252
253 /* We'll be defining these inside the next block and using them afterwards. */
254 double dx_ext, dy_ext;
255 int xsv_ext, ysv_ext;
256
257 double dx1;
258 double dy1;
259 double attn1;
260 double dx2;
261 double dy2;
262 double attn2;
263 double zins;
264 double attn0;
265 double attn_ext;
266
267 double value = 0;
268
269 /* Contribution (1,0) */
270 dx1 = dx0 - 1 - SQUISH_CONSTANT_2D;
271 dy1 = dy0 - 0 - SQUISH_CONSTANT_2D;
272 attn1 = 2 - dx1 * dx1 - dy1 * dy1;
273 if (attn1 > 0) {
274 attn1 *= attn1;
275 value += attn1 * attn1 * extrapolate2(ctx, xsb + 1, ysb + 0, dx1, dy1);
276 }
277
278 /* Contribution (0,1) */
279 dx2 = dx0 - 0 - SQUISH_CONSTANT_2D;
280 dy2 = dy0 - 1 - SQUISH_CONSTANT_2D;
281 attn2 = 2 - dx2 * dx2 - dy2 * dy2;
282 if (attn2 > 0) {
283 attn2 *= attn2;
284 value += attn2 * attn2 * extrapolate2(ctx, xsb + 0, ysb + 1, dx2, dy2);
285 }
286
287 if (inSum <= 1) { /* We're inside the triangle (2-Simplex) at (0,0) */
288 zins = 1 - inSum;
289 if (zins > xins || zins > yins) { /* (0,0) is one of the closest two triangular vertices */
290 if (xins > yins) {
291 xsv_ext = xsb + 1;
292 ysv_ext = ysb - 1;
293 dx_ext = dx0 - 1;
294 dy_ext = dy0 + 1;
295 } else {
296 xsv_ext = xsb - 1;
297 ysv_ext = ysb + 1;
298 dx_ext = dx0 + 1;
299 dy_ext = dy0 - 1;
300 }
301 } else { /* (1,0) and (0,1) are the closest two vertices. */
302 xsv_ext = xsb + 1;
303 ysv_ext = ysb + 1;
304 dx_ext = dx0 - 1 - 2 * SQUISH_CONSTANT_2D;
305 dy_ext = dy0 - 1 - 2 * SQUISH_CONSTANT_2D;
306 }
307 } else { /* We're inside the triangle (2-Simplex) at (1,1) */
308 zins = 2 - inSum;
309 if (zins < xins || zins < yins) { /* (0,0) is one of the closest two triangular vertices */
310 if (xins > yins) {
311 xsv_ext = xsb + 2;
312 ysv_ext = ysb + 0;
313 dx_ext = dx0 - 2 - 2 * SQUISH_CONSTANT_2D;
314 dy_ext = dy0 + 0 - 2 * SQUISH_CONSTANT_2D;
315 } else {
316 xsv_ext = xsb + 0;
317 ysv_ext = ysb + 2;
318 dx_ext = dx0 + 0 - 2 * SQUISH_CONSTANT_2D;
319 dy_ext = dy0 - 2 - 2 * SQUISH_CONSTANT_2D;
320 }
321 } else { /* (1,0) and (0,1) are the closest two vertices. */
322 dx_ext = dx0;
323 dy_ext = dy0;
324 xsv_ext = xsb;
325 ysv_ext = ysb;
326 }
327 xsb += 1;
328 ysb += 1;
329 dx0 = dx0 - 1 - 2 * SQUISH_CONSTANT_2D;
330 dy0 = dy0 - 1 - 2 * SQUISH_CONSTANT_2D;
331 }
332
333 /* Contribution (0,0) or (1,1) */
334 attn0 = 2 - dx0 * dx0 - dy0 * dy0;
335 if (attn0 > 0) {
336 attn0 *= attn0;
337 value += attn0 * attn0 * extrapolate2(ctx, xsb, ysb, dx0, dy0);
338 }
339
340 /* Extra Vertex */
341 attn_ext = 2 - dx_ext * dx_ext - dy_ext * dy_ext;
342 if (attn_ext > 0) {
343 attn_ext *= attn_ext;
344 value += attn_ext * attn_ext * extrapolate2(ctx, xsv_ext, ysv_ext, dx_ext, dy_ext);
345 }
346
347 return value / NORM_CONSTANT_2D;
348}
349
350/*
351 * 3D OpenSimplex (Simplectic) Noise
352 */
353double open_simplex_noise3(struct osn_context *ctx, double x, double y, double z)
354{
355
356 /* Place input coordinates on simplectic honeycomb. */
357 double stretchOffset = (x + y + z) * STRETCH_CONSTANT_3D;
358 double xs = x + stretchOffset;
359 double ys = y + stretchOffset;
360 double zs = z + stretchOffset;
361
362 /* Floor to get simplectic honeycomb coordinates of rhombohedron (stretched cube) super-cell origin. */
363 int xsb = fastFloor(xs);
364 int ysb = fastFloor(ys);
365 int zsb = fastFloor(zs);
366
367 /* Skew out to get actual coordinates of rhombohedron origin. We'll need these later. */
368 double squishOffset = (xsb + ysb + zsb) * SQUISH_CONSTANT_3D;
369 double xb = xsb + squishOffset;
370 double yb = ysb + squishOffset;
371 double zb = zsb + squishOffset;
372
373 /* Compute simplectic honeycomb coordinates relative to rhombohedral origin. */
374 double xins = xs - xsb;
375 double yins = ys - ysb;
376 double zins = zs - zsb;
377
378 /* Sum those together to get a value that determines which region we're in. */
379 double inSum = xins + yins + zins;
380
381 /* Positions relative to origin point. */
382 double dx0 = x - xb;
383 double dy0 = y - yb;
384 double dz0 = z - zb;
385
386 /* We'll be defining these inside the next block and using them afterwards. */
387 double dx_ext0, dy_ext0, dz_ext0;
388 double dx_ext1, dy_ext1, dz_ext1;
389 int xsv_ext0, ysv_ext0, zsv_ext0;
390 int xsv_ext1, ysv_ext1, zsv_ext1;
391
392 double wins;
393 int8_t c, c1, c2;
394 int8_t aPoint, bPoint;
395 double aScore, bScore;
396 int aIsFurtherSide;
397 int bIsFurtherSide;
398 double p1, p2, p3;
399 double score;
400 double attn0, attn1, attn2, attn3, attn4, attn5, attn6;
401 double dx1, dy1, dz1;
402 double dx2, dy2, dz2;
403 double dx3, dy3, dz3;
404 double dx4, dy4, dz4;
405 double dx5, dy5, dz5;
406 double dx6, dy6, dz6;
407 double attn_ext0, attn_ext1;
408
409 double value = 0;
410 if (inSum <= 1) { /* We're inside the tetrahedron (3-Simplex) at (0,0,0) */
411
412 /* Determine which two of (0,0,1), (0,1,0), (1,0,0) are closest. */
413 aPoint = 0x01;
414 aScore = xins;
415 bPoint = 0x02;
416 bScore = yins;
417 if (aScore >= bScore && zins > bScore) {
418 bScore = zins;
419 bPoint = 0x04;
420 } else if (aScore < bScore && zins > aScore) {
421 aScore = zins;
422 aPoint = 0x04;
423 }
424
425 /* Now we determine the two lattice points not part of the tetrahedron that may contribute.
426 This depends on the closest two tetrahedral vertices, including (0,0,0) */
427 wins = 1 - inSum;
428 if (wins > aScore || wins > bScore) { /* (0,0,0) is one of the closest two tetrahedral vertices. */
429 c = (bScore > aScore ? bPoint : aPoint); /* Our other closest vertex is the closest out of a and b. */
430
431 if ((c & 0x01) == 0) {
432 xsv_ext0 = xsb - 1;
433 xsv_ext1 = xsb;
434 dx_ext0 = dx0 + 1;
435 dx_ext1 = dx0;
436 } else {
437 xsv_ext0 = xsv_ext1 = xsb + 1;
438 dx_ext0 = dx_ext1 = dx0 - 1;
439 }
440
441 if ((c & 0x02) == 0) {
442 ysv_ext0 = ysv_ext1 = ysb;
443 dy_ext0 = dy_ext1 = dy0;
444 if ((c & 0x01) == 0) {
445 ysv_ext1 -= 1;
446 dy_ext1 += 1;
447 } else {
448 ysv_ext0 -= 1;
449 dy_ext0 += 1;
450 }
451 } else {
452 ysv_ext0 = ysv_ext1 = ysb + 1;
453 dy_ext0 = dy_ext1 = dy0 - 1;
454 }
455
456 if ((c & 0x04) == 0) {
457 zsv_ext0 = zsb;
458 zsv_ext1 = zsb - 1;
459 dz_ext0 = dz0;
460 dz_ext1 = dz0 + 1;
461 } else {
462 zsv_ext0 = zsv_ext1 = zsb + 1;
463 dz_ext0 = dz_ext1 = dz0 - 1;
464 }
465 } else { /* (0,0,0) is not one of the closest two tetrahedral vertices. */
466 c = (int8_t)(aPoint | bPoint); /* Our two extra vertices are determined by the closest two. */
467
468 if ((c & 0x01) == 0) {
469 xsv_ext0 = xsb;
470 xsv_ext1 = xsb - 1;
471 dx_ext0 = dx0 - 2 * SQUISH_CONSTANT_3D;
472 dx_ext1 = dx0 + 1 - SQUISH_CONSTANT_3D;
473 } else {
474 xsv_ext0 = xsv_ext1 = xsb + 1;
475 dx_ext0 = dx0 - 1 - 2 * SQUISH_CONSTANT_3D;
476 dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_3D;
477 }
478
479 if ((c & 0x02) == 0) {
480 ysv_ext0 = ysb;
481 ysv_ext1 = ysb - 1;
482 dy_ext0 = dy0 - 2 * SQUISH_CONSTANT_3D;
483 dy_ext1 = dy0 + 1 - SQUISH_CONSTANT_3D;
484 } else {
485 ysv_ext0 = ysv_ext1 = ysb + 1;
486 dy_ext0 = dy0 - 1 - 2 * SQUISH_CONSTANT_3D;
487 dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_3D;
488 }
489
490 if ((c & 0x04) == 0) {
491 zsv_ext0 = zsb;
492 zsv_ext1 = zsb - 1;
493 dz_ext0 = dz0 - 2 * SQUISH_CONSTANT_3D;
494 dz_ext1 = dz0 + 1 - SQUISH_CONSTANT_3D;
495 } else {
496 zsv_ext0 = zsv_ext1 = zsb + 1;
497 dz_ext0 = dz0 - 1 - 2 * SQUISH_CONSTANT_3D;
498 dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_3D;
499 }
500 }
501
502 /* Contribution (0,0,0) */
503 attn0 = 2 - dx0 * dx0 - dy0 * dy0 - dz0 * dz0;
504 if (attn0 > 0) {
505 attn0 *= attn0;
506 value += attn0 * attn0 * extrapolate3(ctx, xsb + 0, ysb + 0, zsb + 0, dx0, dy0, dz0);
507 }
508
509 /* Contribution (1,0,0) */
510 dx1 = dx0 - 1 - SQUISH_CONSTANT_3D;
511 dy1 = dy0 - 0 - SQUISH_CONSTANT_3D;
512 dz1 = dz0 - 0 - SQUISH_CONSTANT_3D;
513 attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1;
514 if (attn1 > 0) {
515 attn1 *= attn1;
516 value += attn1 * attn1 * extrapolate3(ctx, xsb + 1, ysb + 0, zsb + 0, dx1, dy1, dz1);
517 }
518
519 /* Contribution (0,1,0) */
520 dx2 = dx0 - 0 - SQUISH_CONSTANT_3D;
521 dy2 = dy0 - 1 - SQUISH_CONSTANT_3D;
522 dz2 = dz1;
523 attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2;
524 if (attn2 > 0) {
525 attn2 *= attn2;
526 value += attn2 * attn2 * extrapolate3(ctx, xsb + 0, ysb + 1, zsb + 0, dx2, dy2, dz2);
527 }
528
529 /* Contribution (0,0,1) */
530 dx3 = dx2;
531 dy3 = dy1;
532 dz3 = dz0 - 1 - SQUISH_CONSTANT_3D;
533 attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3;
534 if (attn3 > 0) {
535 attn3 *= attn3;
536 value += attn3 * attn3 * extrapolate3(ctx, xsb + 0, ysb + 0, zsb + 1, dx3, dy3, dz3);
537 }
538 } else if (inSum >= 2) { /* We're inside the tetrahedron (3-Simplex) at (1,1,1) */
539
540 /* Determine which two tetrahedral vertices are the closest, out of (1,1,0), (1,0,1), (0,1,1) but not (1,1,1). */
541 aPoint = 0x06;
542 aScore = xins;
543 bPoint = 0x05;
544 bScore = yins;
545 if (aScore <= bScore && zins < bScore) {
546 bScore = zins;
547 bPoint = 0x03;
548 } else if (aScore > bScore && zins < aScore) {
549 aScore = zins;
550 aPoint = 0x03;
551 }
552
553 /* Now we determine the two lattice points not part of the tetrahedron that may contribute.
554 This depends on the closest two tetrahedral vertices, including (1,1,1) */
555 wins = 3 - inSum;
556 if (wins < aScore || wins < bScore) { /* (1,1,1) is one of the closest two tetrahedral vertices. */
557 c = (bScore < aScore ? bPoint : aPoint); /* Our other closest vertex is the closest out of a and b. */
558
559 if ((c & 0x01) != 0) {
560 xsv_ext0 = xsb + 2;
561 xsv_ext1 = xsb + 1;
562 dx_ext0 = dx0 - 2 - 3 * SQUISH_CONSTANT_3D;
563 dx_ext1 = dx0 - 1 - 3 * SQUISH_CONSTANT_3D;
564 } else {
565 xsv_ext0 = xsv_ext1 = xsb;
566 dx_ext0 = dx_ext1 = dx0 - 3 * SQUISH_CONSTANT_3D;
567 }
568
569 if ((c & 0x02) != 0) {
570 ysv_ext0 = ysv_ext1 = ysb + 1;
571 dy_ext0 = dy_ext1 = dy0 - 1 - 3 * SQUISH_CONSTANT_3D;
572 if ((c & 0x01) != 0) {
573 ysv_ext1 += 1;
574 dy_ext1 -= 1;
575 } else {
576 ysv_ext0 += 1;
577 dy_ext0 -= 1;
578 }
579 } else {
580 ysv_ext0 = ysv_ext1 = ysb;
581 dy_ext0 = dy_ext1 = dy0 - 3 * SQUISH_CONSTANT_3D;
582 }
583
584 if ((c & 0x04) != 0) {
585 zsv_ext0 = zsb + 1;
586 zsv_ext1 = zsb + 2;
587 dz_ext0 = dz0 - 1 - 3 * SQUISH_CONSTANT_3D;
588 dz_ext1 = dz0 - 2 - 3 * SQUISH_CONSTANT_3D;
589 } else {
590 zsv_ext0 = zsv_ext1 = zsb;
591 dz_ext0 = dz_ext1 = dz0 - 3 * SQUISH_CONSTANT_3D;
592 }
593 } else { /* (1,1,1) is not one of the closest two tetrahedral vertices. */
594 c = (int8_t)(aPoint & bPoint); /* Our two extra vertices are determined by the closest two. */
595
596 if ((c & 0x01) != 0) {
597 xsv_ext0 = xsb + 1;
598 xsv_ext1 = xsb + 2;
599 dx_ext0 = dx0 - 1 - SQUISH_CONSTANT_3D;
600 dx_ext1 = dx0 - 2 - 2 * SQUISH_CONSTANT_3D;
601 } else {
602 xsv_ext0 = xsv_ext1 = xsb;
603 dx_ext0 = dx0 - SQUISH_CONSTANT_3D;
604 dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D;
605 }
606
607 if ((c & 0x02) != 0) {
608 ysv_ext0 = ysb + 1;
609 ysv_ext1 = ysb + 2;
610 dy_ext0 = dy0 - 1 - SQUISH_CONSTANT_3D;
611 dy_ext1 = dy0 - 2 - 2 * SQUISH_CONSTANT_3D;
612 } else {
613 ysv_ext0 = ysv_ext1 = ysb;
614 dy_ext0 = dy0 - SQUISH_CONSTANT_3D;
615 dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D;
616 }
617
618 if ((c & 0x04) != 0) {
619 zsv_ext0 = zsb + 1;
620 zsv_ext1 = zsb + 2;
621 dz_ext0 = dz0 - 1 - SQUISH_CONSTANT_3D;
622 dz_ext1 = dz0 - 2 - 2 * SQUISH_CONSTANT_3D;
623 } else {
624 zsv_ext0 = zsv_ext1 = zsb;
625 dz_ext0 = dz0 - SQUISH_CONSTANT_3D;
626 dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D;
627 }
628 }
629
630 /* Contribution (1,1,0) */
631 dx3 = dx0 - 1 - 2 * SQUISH_CONSTANT_3D;
632 dy3 = dy0 - 1 - 2 * SQUISH_CONSTANT_3D;
633 dz3 = dz0 - 0 - 2 * SQUISH_CONSTANT_3D;
634 attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3;
635 if (attn3 > 0) {
636 attn3 *= attn3;
637 value += attn3 * attn3 * extrapolate3(ctx, xsb + 1, ysb + 1, zsb + 0, dx3, dy3, dz3);
638 }
639
640 /* Contribution (1,0,1) */
641 dx2 = dx3;
642 dy2 = dy0 - 0 - 2 * SQUISH_CONSTANT_3D;
643 dz2 = dz0 - 1 - 2 * SQUISH_CONSTANT_3D;
644 attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2;
645 if (attn2 > 0) {
646 attn2 *= attn2;
647 value += attn2 * attn2 * extrapolate3(ctx, xsb + 1, ysb + 0, zsb + 1, dx2, dy2, dz2);
648 }
649
650 /* Contribution (0,1,1) */
651 dx1 = dx0 - 0 - 2 * SQUISH_CONSTANT_3D;
652 dy1 = dy3;
653 dz1 = dz2;
654 attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1;
655 if (attn1 > 0) {
656 attn1 *= attn1;
657 value += attn1 * attn1 * extrapolate3(ctx, xsb + 0, ysb + 1, zsb + 1, dx1, dy1, dz1);
658 }
659
660 /* Contribution (1,1,1) */
661 dx0 = dx0 - 1 - 3 * SQUISH_CONSTANT_3D;
662 dy0 = dy0 - 1 - 3 * SQUISH_CONSTANT_3D;
663 dz0 = dz0 - 1 - 3 * SQUISH_CONSTANT_3D;
664 attn0 = 2 - dx0 * dx0 - dy0 * dy0 - dz0 * dz0;
665 if (attn0 > 0) {
666 attn0 *= attn0;
667 value += attn0 * attn0 * extrapolate3(ctx, xsb + 1, ysb + 1, zsb + 1, dx0, dy0, dz0);
668 }
669 } else { /* We're inside the octahedron (Rectified 3-Simplex) in between.
670 Decide between point (0,0,1) and (1,1,0) as closest */
671 p1 = xins + yins;
672 if (p1 > 1) {
673 aScore = p1 - 1;
674 aPoint = 0x03;
675 aIsFurtherSide = 1;
676 } else {
677 aScore = 1 - p1;
678 aPoint = 0x04;
679 aIsFurtherSide = 0;
680 }
681
682 /* Decide between point (0,1,0) and (1,0,1) as closest */
683 p2 = xins + zins;
684 if (p2 > 1) {
685 bScore = p2 - 1;
686 bPoint = 0x05;
687 bIsFurtherSide = 1;
688 } else {
689 bScore = 1 - p2;
690 bPoint = 0x02;
691 bIsFurtherSide = 0;
692 }
693
694 /* 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. */
695 p3 = yins + zins;
696 if (p3 > 1) {
697 score = p3 - 1;
698 if (aScore <= bScore && aScore < score) {
699 aScore = score;
700 aPoint = 0x06;
701 aIsFurtherSide = 1;
702 } else if (aScore > bScore && bScore < score) {
703 bScore = score;
704 bPoint = 0x06;
705 bIsFurtherSide = 1;
706 }
707 } else {
708 score = 1 - p3;
709 if (aScore <= bScore && aScore < score) {
710 aScore = score;
711 aPoint = 0x01;
712 aIsFurtherSide = 0;
713 } else if (aScore > bScore && bScore < score) {
714 bScore = score;
715 bPoint = 0x01;
716 bIsFurtherSide = 0;
717 }
718 }
719
720 /* Where each of the two closest points are determines how the extra two vertices are calculated. */
721 if (aIsFurtherSide == bIsFurtherSide) {
722 if (aIsFurtherSide) { /* Both closest points on (1,1,1) side */
723
724 /* One of the two extra points is (1,1,1) */
725 dx_ext0 = dx0 - 1 - 3 * SQUISH_CONSTANT_3D;
726 dy_ext0 = dy0 - 1 - 3 * SQUISH_CONSTANT_3D;
727 dz_ext0 = dz0 - 1 - 3 * SQUISH_CONSTANT_3D;
728 xsv_ext0 = xsb + 1;
729 ysv_ext0 = ysb + 1;
730 zsv_ext0 = zsb + 1;
731
732 /* Other extra point is based on the shared axis. */
733 c = (int8_t)(aPoint & bPoint);
734 if ((c & 0x01) != 0) {
735 dx_ext1 = dx0 - 2 - 2 * SQUISH_CONSTANT_3D;
736 dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D;
737 dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D;
738 xsv_ext1 = xsb + 2;
739 ysv_ext1 = ysb;
740 zsv_ext1 = zsb;
741 } else if ((c & 0x02) != 0) {
742 dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D;
743 dy_ext1 = dy0 - 2 - 2 * SQUISH_CONSTANT_3D;
744 dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D;
745 xsv_ext1 = xsb;
746 ysv_ext1 = ysb + 2;
747 zsv_ext1 = zsb;
748 } else {
749 dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D;
750 dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D;
751 dz_ext1 = dz0 - 2 - 2 * SQUISH_CONSTANT_3D;
752 xsv_ext1 = xsb;
753 ysv_ext1 = ysb;
754 zsv_ext1 = zsb + 2;
755 }
756 } else { /* Both closest points on (0,0,0) side */
757
758 /* One of the two extra points is (0,0,0) */
759 dx_ext0 = dx0;
760 dy_ext0 = dy0;
761 dz_ext0 = dz0;
762 xsv_ext0 = xsb;
763 ysv_ext0 = ysb;
764 zsv_ext0 = zsb;
765
766 /* Other extra point is based on the omitted axis. */
767 c = (int8_t)(aPoint | bPoint);
768 if ((c & 0x01) == 0) {
769 dx_ext1 = dx0 + 1 - SQUISH_CONSTANT_3D;
770 dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_3D;
771 dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_3D;
772 xsv_ext1 = xsb - 1;
773 ysv_ext1 = ysb + 1;
774 zsv_ext1 = zsb + 1;
775 } else if ((c & 0x02) == 0) {
776 dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_3D;
777 dy_ext1 = dy0 + 1 - SQUISH_CONSTANT_3D;
778 dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_3D;
779 xsv_ext1 = xsb + 1;
780 ysv_ext1 = ysb - 1;
781 zsv_ext1 = zsb + 1;
782 } else {
783 dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_3D;
784 dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_3D;
785 dz_ext1 = dz0 + 1 - SQUISH_CONSTANT_3D;
786 xsv_ext1 = xsb + 1;
787 ysv_ext1 = ysb + 1;
788 zsv_ext1 = zsb - 1;
789 }
790 }
791 } else { /* One point on (0,0,0) side, one point on (1,1,1) side */
792 if (aIsFurtherSide) {
793 c1 = aPoint;
794 c2 = bPoint;
795 } else {
796 c1 = bPoint;
797 c2 = aPoint;
798 }
799
800 /* One contribution is a permutation of (1,1,-1) */
801 if ((c1 & 0x01) == 0) {
802 dx_ext0 = dx0 + 1 - SQUISH_CONSTANT_3D;
803 dy_ext0 = dy0 - 1 - SQUISH_CONSTANT_3D;
804 dz_ext0 = dz0 - 1 - SQUISH_CONSTANT_3D;
805 xsv_ext0 = xsb - 1;
806 ysv_ext0 = ysb + 1;
807 zsv_ext0 = zsb + 1;
808 } else if ((c1 & 0x02) == 0) {
809 dx_ext0 = dx0 - 1 - SQUISH_CONSTANT_3D;
810 dy_ext0 = dy0 + 1 - SQUISH_CONSTANT_3D;
811 dz_ext0 = dz0 - 1 - SQUISH_CONSTANT_3D;
812 xsv_ext0 = xsb + 1;
813 ysv_ext0 = ysb - 1;
814 zsv_ext0 = zsb + 1;
815 } else {
816 dx_ext0 = dx0 - 1 - SQUISH_CONSTANT_3D;
817 dy_ext0 = dy0 - 1 - SQUISH_CONSTANT_3D;
818 dz_ext0 = dz0 + 1 - SQUISH_CONSTANT_3D;
819 xsv_ext0 = xsb + 1;
820 ysv_ext0 = ysb + 1;
821 zsv_ext0 = zsb - 1;
822 }
823
824 /* One contribution is a permutation of (0,0,2) */
825 dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D;
826 dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D;
827 dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D;
828 xsv_ext1 = xsb;
829 ysv_ext1 = ysb;
830 zsv_ext1 = zsb;
831 if ((c2 & 0x01) != 0) {
832 dx_ext1 -= 2;
833 xsv_ext1 += 2;
834 } else if ((c2 & 0x02) != 0) {
835 dy_ext1 -= 2;
836 ysv_ext1 += 2;
837 } else {
838 dz_ext1 -= 2;
839 zsv_ext1 += 2;
840 }
841 }
842
843 /* Contribution (1,0,0) */
844 dx1 = dx0 - 1 - SQUISH_CONSTANT_3D;
845 dy1 = dy0 - 0 - SQUISH_CONSTANT_3D;
846 dz1 = dz0 - 0 - SQUISH_CONSTANT_3D;
847 attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1;
848 if (attn1 > 0) {
849 attn1 *= attn1;
850 value += attn1 * attn1 * extrapolate3(ctx, xsb + 1, ysb + 0, zsb + 0, dx1, dy1, dz1);
851 }
852
853 /* Contribution (0,1,0) */
854 dx2 = dx0 - 0 - SQUISH_CONSTANT_3D;
855 dy2 = dy0 - 1 - SQUISH_CONSTANT_3D;
856 dz2 = dz1;
857 attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2;
858 if (attn2 > 0) {
859 attn2 *= attn2;
860 value += attn2 * attn2 * extrapolate3(ctx, xsb + 0, ysb + 1, zsb + 0, dx2, dy2, dz2);
861 }
862
863 /* Contribution (0,0,1) */
864 dx3 = dx2;
865 dy3 = dy1;
866 dz3 = dz0 - 1 - SQUISH_CONSTANT_3D;
867 attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3;
868 if (attn3 > 0) {
869 attn3 *= attn3;
870 value += attn3 * attn3 * extrapolate3(ctx, xsb + 0, ysb + 0, zsb + 1, dx3, dy3, dz3);
871 }
872
873 /* Contribution (1,1,0) */
874 dx4 = dx0 - 1 - 2 * SQUISH_CONSTANT_3D;
875 dy4 = dy0 - 1 - 2 * SQUISH_CONSTANT_3D;
876 dz4 = dz0 - 0 - 2 * SQUISH_CONSTANT_3D;
877 attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4;
878 if (attn4 > 0) {
879 attn4 *= attn4;
880 value += attn4 * attn4 * extrapolate3(ctx, xsb + 1, ysb + 1, zsb + 0, dx4, dy4, dz4);
881 }
882
883 /* Contribution (1,0,1) */
884 dx5 = dx4;
885 dy5 = dy0 - 0 - 2 * SQUISH_CONSTANT_3D;
886 dz5 = dz0 - 1 - 2 * SQUISH_CONSTANT_3D;
887 attn5 = 2 - dx5 * dx5 - dy5 * dy5 - dz5 * dz5;
888 if (attn5 > 0) {
889 attn5 *= attn5;
890 value += attn5 * attn5 * extrapolate3(ctx, xsb + 1, ysb + 0, zsb + 1, dx5, dy5, dz5);
891 }
892
893 /* Contribution (0,1,1) */
894 dx6 = dx0 - 0 - 2 * SQUISH_CONSTANT_3D;
895 dy6 = dy4;
896 dz6 = dz5;
897 attn6 = 2 - dx6 * dx6 - dy6 * dy6 - dz6 * dz6;
898 if (attn6 > 0) {
899 attn6 *= attn6;
900 value += attn6 * attn6 * extrapolate3(ctx, xsb + 0, ysb + 1, zsb + 1, dx6, dy6, dz6);
901 }
902 }
903
904 /* First extra vertex */
905 attn_ext0 = 2 - dx_ext0 * dx_ext0 - dy_ext0 * dy_ext0 - dz_ext0 * dz_ext0;
906 if (attn_ext0 > 0)
907 {
908 attn_ext0 *= attn_ext0;
909 value += attn_ext0 * attn_ext0 * extrapolate3(ctx, xsv_ext0, ysv_ext0, zsv_ext0, dx_ext0, dy_ext0, dz_ext0);
910 }
911
912 /* Second extra vertex */
913 attn_ext1 = 2 - dx_ext1 * dx_ext1 - dy_ext1 * dy_ext1 - dz_ext1 * dz_ext1;
914 if (attn_ext1 > 0)
915 {
916 attn_ext1 *= attn_ext1;
917 value += attn_ext1 * attn_ext1 * extrapolate3(ctx, xsv_ext1, ysv_ext1, zsv_ext1, dx_ext1, dy_ext1, dz_ext1);
918 }
919
920 return value / NORM_CONSTANT_3D;
921}
922
923/*
924 * 4D OpenSimplex (Simplectic) Noise.
925 */
926double open_simplex_noise4(struct osn_context *ctx, double x, double y, double z, double w)
927{
928 double uins;
929 double dx1, dy1, dz1, dw1;
930 double dx2, dy2, dz2, dw2;
931 double dx3, dy3, dz3, dw3;
932 double dx4, dy4, dz4, dw4;
933 double dx5, dy5, dz5, dw5;
934 double dx6, dy6, dz6, dw6;
935 double dx7, dy7, dz7, dw7;
936 double dx8, dy8, dz8, dw8;
937 double dx9, dy9, dz9, dw9;
938 double dx10, dy10, dz10, dw10;
939 double attn0, attn1, attn2, attn3, attn4;
940 double attn5, attn6, attn7, attn8, attn9, attn10;
941 double attn_ext0, attn_ext1, attn_ext2;
942 int8_t c, c1, c2;
943 int8_t aPoint, bPoint;
944 double aScore, bScore;
945 int aIsBiggerSide;
946 int bIsBiggerSide;
947 double p1, p2, p3, p4;
948 double score;
949
950 /* Place input coordinates on simplectic honeycomb. */
951 double stretchOffset = (x + y + z + w) * STRETCH_CONSTANT_4D;
952 double xs = x + stretchOffset;
953 double ys = y + stretchOffset;
954 double zs = z + stretchOffset;
955 double ws = w + stretchOffset;
956
957 /* Floor to get simplectic honeycomb coordinates of rhombo-hypercube super-cell origin. */
958 int xsb = fastFloor(xs);
959 int ysb = fastFloor(ys);
960 int zsb = fastFloor(zs);
961 int wsb = fastFloor(ws);
962
963 /* Skew out to get actual coordinates of stretched rhombo-hypercube origin. We'll need these later. */
964 double squishOffset = (xsb + ysb + zsb + wsb) * SQUISH_CONSTANT_4D;
965 double xb = xsb + squishOffset;
966 double yb = ysb + squishOffset;
967 double zb = zsb + squishOffset;
968 double wb = wsb + squishOffset;
969
970 /* Compute simplectic honeycomb coordinates relative to rhombo-hypercube origin. */
971 double xins = xs - xsb;
972 double yins = ys - ysb;
973 double zins = zs - zsb;
974 double wins = ws - wsb;
975
976 /* Sum those together to get a value that determines which region we're in. */
977 double inSum = xins + yins + zins + wins;
978
979 /* Positions relative to origin point. */
980 double dx0 = x - xb;
981 double dy0 = y - yb;
982 double dz0 = z - zb;
983 double dw0 = w - wb;
984
985 /* We'll be defining these inside the next block and using them afterwards. */
986 double dx_ext0, dy_ext0, dz_ext0, dw_ext0;
987 double dx_ext1, dy_ext1, dz_ext1, dw_ext1;
988 double dx_ext2, dy_ext2, dz_ext2, dw_ext2;
989 int xsv_ext0, ysv_ext0, zsv_ext0, wsv_ext0;
990 int xsv_ext1, ysv_ext1, zsv_ext1, wsv_ext1;
991 int xsv_ext2, ysv_ext2, zsv_ext2, wsv_ext2;
992
993 double value = 0;
994 if (inSum <= 1) { /* We're inside the pentachoron (4-Simplex) at (0,0,0,0) */
995
996 /* Determine which two of (0,0,0,1), (0,0,1,0), (0,1,0,0), (1,0,0,0) are closest. */
997 aPoint = 0x01;
998 aScore = xins;
999 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
value
float value(const T *src, size_t ch)
Definition JPEGXLImport.cpp:93
p
const Params2D p
Definition KisBezierUtils.cpp:703
p2
QPointF p2
Definition KisBezierUtils.cpp:499
p3
QPointF p3
Definition KisBezierUtils.cpp:499
p1
QPointF p1
Definition KisBezierUtils.cpp:499
source
KisMagneticGraph::vertex_descriptor source(typename KisMagneticGraph::edge_descriptor e, KisMagneticGraph g)
Definition KisMagneticGraph.h:251
gradients3D
static const signed char gradients3D[]
Definition open-simplex-noise.c:61
NORM_CONSTANT_4D
#define NORM_CONSTANT_4D
Definition open-simplex-noise.c:33
gradients2D
static const int8_t gradients2D[]
Definition open-simplex-noise.c:48
open_simplex_noise
int open_simplex_noise(int64_t seed, struct osn_context **ctx)
Definition open-simplex-noise.c:168
open_simplex_noise_free
void open_simplex_noise_free(struct osn_context *ctx)
Definition open-simplex-noise.c:209
SQUISH_CONSTANT_4D
#define SQUISH_CONSTANT_4D
Definition open-simplex-noise.c:29
open_simplex_noise2
double open_simplex_noise2(struct osn_context *ctx, double x, double y)
Definition open-simplex-noise.c:225
extrapolate2
static double extrapolate2(struct osn_context *ctx, int xsb, int ysb, double dx, double dy)
Definition open-simplex-noise.c:97
SQUISH_CONSTANT_2D
#define SQUISH_CONSTANT_2D
Definition open-simplex-noise.c:25
extrapolate3
static double extrapolate3(struct osn_context *ctx, int xsb, int ysb, int zsb, double dx, double dy, double dz)
Definition open-simplex-noise.c:105
extrapolate4
static double extrapolate4(struct osn_context *ctx, int xsb, int ysb, int zsb, int wsb, double dx, double dy, double dz, double dw)
Definition open-simplex-noise.c:115
NORM_CONSTANT_3D
#define NORM_CONSTANT_3D
Definition open-simplex-noise.c:32
SQUISH_CONSTANT_3D
#define SQUISH_CONSTANT_3D
Definition open-simplex-noise.c:27
open_simplex_noise_init_perm
int open_simplex_noise_init_perm(struct osn_context *ctx, int16_t p[], int nelements)
Definition open-simplex-noise.c:147
STRETCH_CONSTANT_3D
#define STRETCH_CONSTANT_3D
Definition open-simplex-noise.c:26
open_simplex_noise3
double open_simplex_noise3(struct osn_context *ctx, double x, double y, double z)
Definition open-simplex-noise.c:353
ARRAYSIZE
#define ARRAYSIZE(x)
Definition open-simplex-noise.c:42
fastFloor
static INLINE int fastFloor(double x)
Definition open-simplex-noise.c:125
STRETCH_CONSTANT_2D
#define STRETCH_CONSTANT_2D
Definition open-simplex-noise.c:24
allocate_perm
static int allocate_perm(struct osn_context *ctx, int nperm, int ngrad)
Definition open-simplex-noise.c:130
gradients4D
static const signed char gradients4D[]
Definition open-simplex-noise.c:78
NORM_CONSTANT_2D
#define NORM_CONSTANT_2D
Definition open-simplex-noise.c:31
STRETCH_CONSTANT_4D
#define STRETCH_CONSTANT_4D
Definition open-simplex-noise.c:28
open_simplex_noise4
double open_simplex_noise4(struct osn_context *ctx, double x, double y, double z, double w)
Definition open-simplex-noise.c:926
INLINE
#define INLINE
Definition open-simplex-noise.h:33
osn_context
Definition open-simplex-noise.c:37
osn_context::permGradIndex3D
int16_t * permGradIndex3D
Definition open-simplex-noise.c:39
osn_context::perm
int16_t * perm
Definition open-simplex-noise.c:38