open-simplex-noise.c 65 KB

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