open-simplex-noise.c 65 KB

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