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