Guest User

Untitled

a guest
Nov 17th, 2017
57
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 73.24 KB | None | 0 0
  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.  
  17. public 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 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 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 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 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. }
Add Comment
Please, Sign In to add comment