You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

open_simplex_noise.cpp 65KB

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