1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3
4 /// @file math/OpenSimplexNoise.cc
5
6 #include "OpenSimplexNoise.h"
7
8 #include <algorithm>
9 #include <cmath>
10 #include <type_traits>
11
12 // see OpenSimplexNoise.h for details about the origin on this code
13
14 namespace OSN {
15
16 namespace {
17
18 template <typename T>
pow4(T x)19 inline T pow4 (T x)
20 {
21 x *= x;
22 return x*x;
23 }
24
25 template <typename T>
pow2(T x)26 inline T pow2 (T x)
27 {
28 return x*x;
29 }
30
31 template <typename T>
fastFloori(T x)32 inline OSNoise::inttype fastFloori (T x)
33 {
34 OSNoise::inttype ip = (OSNoise::inttype)x;
35
36 if (x < 0.0) --ip;
37
38 return ip;
39 }
40
LCG_STEP(int64_t & x)41 inline void LCG_STEP (int64_t & x)
42 {
43 // Magic constants are attributed to Donald Knuth's MMIX implementation.
44 static const int64_t MULTIPLIER = 6364136223846793005LL;
45 static const int64_t INCREMENT = 1442695040888963407LL;
46 x = ((x * MULTIPLIER) + INCREMENT);
47 }
48
49 } // anonymous namespace
50
51 // Array of gradient values for 3D. They approximate the directions to the
52 // vertices of a rhombicuboctahedron from its center, skewed so that the
53 // triangular and square facets can be inscribed in circles of the same radius.
54 // New gradient set 2014-10-06.
55 const int OSNoise::sGradients [] = {
56 -11, 4, 4, -4, 11, 4, -4, 4, 11, 11, 4, 4, 4, 11, 4, 4, 4, 11,
57 -11,-4, 4, -4,-11, 4, -4,-4, 11, 11,-4, 4, 4,-11, 4, 4,-4, 11,
58 -11, 4,-4, -4, 11,-4, -4, 4,-11, 11, 4,-4, 4, 11,-4, 4, 4,-11,
59 -11,-4,-4, -4,-11,-4, -4,-4,-11, 11,-4,-4, 4,-11,-4, 4,-4,-11
60 };
61
62 template <typename T>
extrapolate(const OSNoise::inttype xsb,const OSNoise::inttype ysb,const OSNoise::inttype zsb,const T dx,const T dy,const T dz) const63 inline T OSNoise::extrapolate(const OSNoise::inttype xsb,
64 const OSNoise::inttype ysb,
65 const OSNoise::inttype zsb,
66 const T dx,
67 const T dy,
68 const T dz) const
69 {
70 unsigned int index = mPermGradIndex[(mPerm[(mPerm[xsb & 0xFF] + ysb) & 0xFF] + zsb) & 0xFF];
71 return sGradients[index] * dx +
72 sGradients[index + 1] * dy +
73 sGradients[index + 2] * dz;
74
75 }
76
77 template <typename T>
extrapolate(const OSNoise::inttype xsb,const OSNoise::inttype ysb,const OSNoise::inttype zsb,const T dx,const T dy,const T dz,T (& de)[3]) const78 inline T OSNoise::extrapolate(const OSNoise::inttype xsb,
79 const OSNoise::inttype ysb,
80 const OSNoise::inttype zsb,
81 const T dx,
82 const T dy,
83 const T dz,
84 T (&de) [3]) const
85 {
86 unsigned int index = mPermGradIndex[(mPerm[(mPerm[xsb & 0xFF] + ysb) & 0xFF] + zsb) & 0xFF];
87 return (de[0] = sGradients[index]) * dx +
88 (de[1] = sGradients[index + 1]) * dy +
89 (de[2] = sGradients[index + 2]) * dz;
90 }
91
OSNoise(OSNoise::inttype seed)92 OSNoise::OSNoise(OSNoise::inttype seed)
93 {
94 int source [256];
95 for (int i = 0; i < 256; ++i) { source[i] = i; }
96 LCG_STEP(seed);
97 LCG_STEP(seed);
98 LCG_STEP(seed);
99 for (int i = 255; i >= 0; --i) {
100 LCG_STEP(seed);
101 int r = (int)((seed + 31) % (i + 1));
102 if (r < 0) { r += (i + 1); }
103 mPerm[i] = source[r];
104 mPermGradIndex[i] = (int)((mPerm[i] % (72 / 3)) * 3);
105 source[r] = source[i];
106 }
107 }
108
OSNoise(const int * p)109 OSNoise::OSNoise(const int * p)
110 {
111 // Copy the supplied permutation array into this instance.
112 for (int i = 0; i < 256; ++i) {
113 mPerm[i] = p[i];
114 mPermGradIndex[i] = (int)((mPerm[i] % (72 / 3)) * 3);
115 }
116 }
117
118 template <typename T>
eval(const T x,const T y,const T z) const119 T OSNoise::eval(const T x, const T y, const T z) const
120 {
121 static_assert(std::is_floating_point<T>::value, "OpenSimplexNoise can only be used with floating-point types");
122
123 static const T STRETCH_CONSTANT = (T)(-1.0 / 6.0); // (1 / sqrt(3 + 1) - 1) / 3
124 static const T SQUISH_CONSTANT = (T)(1.0 / 3.0); // (sqrt(3 + 1) - 1) / 3
125 static const T NORM_CONSTANT = (T)(1.0 / 103.0);
126
127 OSNoise::inttype xsb, ysb, zsb;
128 T dx0, dy0, dz0;
129 T xins, yins, zins;
130
131 // Parameters for the individual contributions
132 T contr_m [9], contr_ext [9];
133
134 {
135 // Place input coordinates on simplectic lattice.
136 T stretchOffset = (x + y + z) * STRETCH_CONSTANT;
137 T xs = x + stretchOffset;
138 T ys = y + stretchOffset;
139 T zs = z + stretchOffset;
140
141 // Floor to get simplectic lattice coordinates of rhombohedron
142 // (stretched cube) super-cell.
143 #ifdef __FAST_MATH__
144 T xsbd = std::floor(xs);
145 T ysbd = std::floor(ys);
146 T zsbd = std::floor(zs);
147 xsb = (OSNoise::inttype)xsbd;
148 ysb = (OSNoise::inttype)ysbd;
149 zsb = (OSNoise::inttype)zsbd;
150 #else
151 xsb = fastFloori(xs);
152 ysb = fastFloori(ys);
153 zsb = fastFloori(zs);
154 T xsbd = (T)xsb;
155 T ysbd = (T)ysb;
156 T zsbd = (T)zsb;
157 #endif
158
159 // Skew out to get actual coordinates of rhombohedron origin.
160 T squishOffset = (xsbd + ysbd + zsbd) * SQUISH_CONSTANT;
161 T xb = xsbd + squishOffset;
162 T yb = ysbd + squishOffset;
163 T zb = zsbd + squishOffset;
164
165 // Positions relative to origin point.
166 dx0 = x - xb;
167 dy0 = y - yb;
168 dz0 = z - zb;
169
170 // Compute simplectic lattice coordinates relative to rhombohedral origin.
171 xins = xs - xsbd;
172 yins = ys - ysbd;
173 zins = zs - zsbd;
174 }
175
176 // These are given values inside the next block, and used afterwards.
177 OSNoise::inttype xsv_ext0, ysv_ext0, zsv_ext0;
178 OSNoise::inttype xsv_ext1, ysv_ext1, zsv_ext1;
179 T dx_ext0, dy_ext0, dz_ext0;
180 T dx_ext1, dy_ext1, dz_ext1;
181
182 // Sum together to get a value that determines which cell we are in.
183 T inSum = xins + yins + zins;
184
185 if (inSum > (T)1.0 && inSum < (T)2.0) {
186 // The point is inside the octahedron (rectified 3-Simplex) inbetween.
187
188 T aScore;
189 uint_fast8_t aPoint;
190 bool aIsFurtherSide;
191 T bScore;
192 uint_fast8_t bPoint;
193 bool bIsFurtherSide;
194
195 // Decide between point (1,0,0) and (0,1,1) as closest.
196 T p1 = xins + yins;
197 if (p1 <= (T)1.0) {
198 aScore = (T)1.0 - p1;
199 aPoint = 4;
200 aIsFurtherSide = false;
201 } else {
202 aScore = p1 - (T)1.0;
203 aPoint = 3;
204 aIsFurtherSide = true;
205 }
206
207 // Decide between point (0,1,0) and (1,0,1) as closest.
208 T p2 = xins + zins;
209 if (p2 <= (T)1.0) {
210 bScore = (T)1.0 - p2;
211 bPoint = 2;
212 bIsFurtherSide = false;
213 } else {
214 bScore = p2 - (T)1.0;
215 bPoint = 5;
216 bIsFurtherSide = true;
217 }
218
219 // The closest out of the two (0,0,1) and (1,1,0) will replace the
220 // furthest out of the two decided above if closer.
221 T p3 = yins + zins;
222 if (p3 > (T)1.0) {
223 T score = p3 - (T)1.0;
224 if (aScore > bScore && bScore < score) {
225 bScore = score;
226 bPoint = 6;
227 bIsFurtherSide = true;
228 } else if (aScore <= bScore && aScore < score) {
229 aScore = score;
230 aPoint = 6;
231 aIsFurtherSide = true;
232 }
233 } else {
234 T score = (T)1.0 - p3;
235 if (aScore > bScore && bScore < score) {
236 bScore = score;
237 bPoint = 1;
238 bIsFurtherSide = false;
239 } else if (aScore <= bScore && aScore < score) {
240 aScore = score;
241 aPoint = 1;
242 aIsFurtherSide = false;
243 }
244 }
245
246 // Where each of the two closest points are determines how the
247 // extra two vertices are calculated.
248 if (aIsFurtherSide == bIsFurtherSide) {
249 if (aIsFurtherSide) {
250 // Both closest points on (1,1,1) side.
251
252 // One of the two extra points is (1,1,1)
253 xsv_ext0 = xsb + 1;
254 ysv_ext0 = ysb + 1;
255 zsv_ext0 = zsb + 1;
256 dx_ext0 = dx0 - (T)1.0 - (SQUISH_CONSTANT * (T)3.0);
257 dy_ext0 = dy0 - (T)1.0 - (SQUISH_CONSTANT * (T)3.0);
258 dz_ext0 = dz0 - (T)1.0 - (SQUISH_CONSTANT * (T)3.0);
259
260 // Other extra point is based on the shared axis.
261 uint_fast8_t c = aPoint & bPoint;
262 if (c & 0x01) {
263 xsv_ext1 = xsb + 2;
264 ysv_ext1 = ysb;
265 zsv_ext1 = zsb;
266 dx_ext1 = dx0 - (T)2.0 - (SQUISH_CONSTANT * (T)2.0);
267 dy_ext1 = dy0 - (SQUISH_CONSTANT * (T)2.0);
268 dz_ext1 = dz0 - (SQUISH_CONSTANT * (T)2.0);
269 } else if (c & 0x02) {
270 xsv_ext1 = xsb;
271 ysv_ext1 = ysb + 2;
272 zsv_ext1 = zsb;
273 dx_ext1 = dx0 - (SQUISH_CONSTANT * (T)2.0);
274 dy_ext1 = dy0 - (T)2.0 - (SQUISH_CONSTANT * (T)2.0);
275 dz_ext1 = dz0 - (SQUISH_CONSTANT * (T)2.0);
276 } else {
277 xsv_ext1 = xsb;
278 ysv_ext1 = ysb;
279 zsv_ext1 = zsb + 2;
280 dx_ext1 = dx0 - (SQUISH_CONSTANT * (T)2.0);
281 dy_ext1 = dy0 - (SQUISH_CONSTANT * (T)2.0);
282 dz_ext1 = dz0 - (T)2.0 - (SQUISH_CONSTANT * (T)2.0);
283 }
284 } else {
285 // Both closest points are on the (0,0,0) side.
286
287 // One of the two extra points is (0,0,0).
288 xsv_ext0 = xsb;
289 ysv_ext0 = ysb;
290 zsv_ext0 = zsb;
291 dx_ext0 = dx0;
292 dy_ext0 = dy0;
293 dz_ext0 = dz0;
294
295 // The other extra point is based on the omitted axis.
296 uint_fast8_t c = aPoint | bPoint;
297 if (!(c & 0x01)) {
298 xsv_ext1 = xsb - 1;
299 ysv_ext1 = ysb + 1;
300 zsv_ext1 = zsb + 1;
301 dx_ext1 = dx0 + (T)1.0 - SQUISH_CONSTANT;
302 dy_ext1 = dy0 - (T)1.0 - SQUISH_CONSTANT;
303 dz_ext1 = dz0 - (T)1.0 - SQUISH_CONSTANT;
304 } else if (!(c & 0x02)) {
305 xsv_ext1 = xsb + 1;
306 ysv_ext1 = ysb - 1;
307 zsv_ext1 = zsb + 1;
308 dx_ext1 = dx0 - (T)1.0 - SQUISH_CONSTANT;
309 dy_ext1 = dy0 + (T)1.0 - SQUISH_CONSTANT;
310 dz_ext1 = dz0 - (T)1.0 - SQUISH_CONSTANT;
311 } else {
312 xsv_ext1 = xsb + 1;
313 ysv_ext1 = ysb + 1;
314 zsv_ext1 = zsb - 1;
315 dx_ext1 = dx0 - (T)1.0 - SQUISH_CONSTANT;
316 dy_ext1 = dy0 - (T)1.0 - SQUISH_CONSTANT;
317 dz_ext1 = dz0 + (T)1.0 - SQUISH_CONSTANT;
318 }
319 }
320 } else {
321 // One point is on the (0,0,0) side, one point is on the (1,1,1) side.
322
323 uint_fast8_t c1, c2;
324 if (aIsFurtherSide) {
325 c1 = aPoint;
326 c2 = bPoint;
327 } else {
328 c1 = bPoint;
329 c2 = aPoint;
330 }
331
332 // One contribution is a permutation of (1,1,-1).
333 if (!(c1 & 0x01)) {
334 xsv_ext0 = xsb - 1;
335 ysv_ext0 = ysb + 1;
336 zsv_ext0 = zsb + 1;
337 dx_ext0 = dx0 + (T)1.0 - SQUISH_CONSTANT;
338 dy_ext0 = dy0 - (T)1.0 - SQUISH_CONSTANT;
339 dz_ext0 = dz0 - (T)1.0 - SQUISH_CONSTANT;
340 } else if (!(c1 & 0x02)) {
341 xsv_ext0 = xsb + 1;
342 ysv_ext0 = ysb - 1;
343 zsv_ext0 = zsb + 1;
344 dx_ext0 = dx0 - (T)1.0 - SQUISH_CONSTANT;
345 dy_ext0 = dy0 + (T)1.0 - SQUISH_CONSTANT;
346 dz_ext0 = dz0 - (T)1.0 - SQUISH_CONSTANT;
347 } else {
348 xsv_ext0 = xsb + 1;
349 ysv_ext0 = ysb + 1;
350 zsv_ext0 = zsb - 1;
351 dx_ext0 = dx0 - (T)1.0 - SQUISH_CONSTANT;
352 dy_ext0 = dy0 - (T)1.0 - SQUISH_CONSTANT;
353 dz_ext0 = dz0 + (T)1.0 - SQUISH_CONSTANT;
354 }
355
356 // One contribution is a permutation of (0,0,2).
357 if (c2 & 0x01) {
358 xsv_ext1 = xsb + 2;
359 ysv_ext1 = ysb;
360 zsv_ext1 = zsb;
361 dx_ext1 = dx0 - (T)2.0 - (SQUISH_CONSTANT * (T)2.0);
362 dy_ext1 = dy0 - (SQUISH_CONSTANT * (T)2.0);
363 dz_ext1 = dz0 - (SQUISH_CONSTANT * (T)2.0);
364 } else if (c2 & 0x02) {
365 xsv_ext1 = xsb;
366 ysv_ext1 = ysb + 2;
367 zsv_ext1 = zsb;
368 dx_ext1 = dx0 - (SQUISH_CONSTANT * (T)2.0);
369 dy_ext1 = dy0 - (T)2.0 - (SQUISH_CONSTANT * (T)2.0);
370 dz_ext1 = dz0 - (SQUISH_CONSTANT * (T)2.0);
371 } else {
372 xsv_ext1 = xsb;
373 ysv_ext1 = ysb;
374 zsv_ext1 = zsb + 2;
375 dx_ext1 = dx0 - (SQUISH_CONSTANT * (T)2.0);
376 dy_ext1 = dy0 - (SQUISH_CONSTANT * (T)2.0);
377 dz_ext1 = dz0 - (T)2.0 - (SQUISH_CONSTANT * (T)2.0);
378 }
379 }
380
381 contr_m[0] = contr_ext[0] = 0.0;
382
383 // Contribution (0,0,1).
384 T dx1 = dx0 - (T)1.0 - SQUISH_CONSTANT;
385 T dy1 = dy0 - SQUISH_CONSTANT;
386 T dz1 = dz0 - SQUISH_CONSTANT;
387 contr_m[1] = pow2(dx1) + pow2(dy1) + pow2(dz1);
388 contr_ext[1] = extrapolate(xsb + 1, ysb, zsb, dx1, dy1, dz1);
389
390 // Contribution (0,1,0).
391 T dx2 = dx0 - SQUISH_CONSTANT;
392 T dy2 = dy0 - (T)1.0 - SQUISH_CONSTANT;
393 T dz2 = dz1;
394 contr_m[2] = pow2(dx2) + pow2(dy2) + pow2(dz2);
395 contr_ext[2] = extrapolate(xsb, ysb + 1, zsb, dx2, dy2, dz2);
396
397 // Contribution (1,0,0).
398 T dx3 = dx2;
399 T dy3 = dy1;
400 T dz3 = dz0 - (T)1.0 - SQUISH_CONSTANT;
401 contr_m[3] = pow2(dx3) + pow2(dy3) + pow2(dz3);
402 contr_ext[3] = extrapolate(xsb, ysb, zsb + 1, dx3, dy3, dz3);
403
404 // Contribution (1,1,0).
405 T dx4 = dx0 - (T)1.0 - (SQUISH_CONSTANT * (T)2.0);
406 T dy4 = dy0 - (T)1.0 - (SQUISH_CONSTANT * (T)2.0);
407 T dz4 = dz0 - (SQUISH_CONSTANT * (T)2.0);
408 contr_m[4] = pow2(dx4) + pow2(dy4) + pow2(dz4);
409 contr_ext[4] = extrapolate(xsb + 1, ysb + 1, zsb, dx4, dy4, dz4);
410
411 // Contribution (1,0,1).
412 T dx5 = dx4;
413 T dy5 = dy0 - (SQUISH_CONSTANT * (T)2.0);
414 T dz5 = dz0 - (T)1.0 - (SQUISH_CONSTANT * (T)2.0);
415 contr_m[5] = pow2(dx5) + pow2(dy5) + pow2(dz5);
416 contr_ext[5] = extrapolate(xsb + 1, ysb, zsb + 1, dx5, dy5, dz5);
417
418 // Contribution (0,1,1).
419 T dx6 = dx0 - (SQUISH_CONSTANT * (T)2.0);
420 T dy6 = dy4;
421 T dz6 = dz5;
422 contr_m[6] = pow2(dx6) + pow2(dy6) + pow2(dz6);
423 contr_ext[6] = extrapolate(xsb, ysb + 1, zsb + 1, dx6, dy6, dz6);
424
425 } else if (inSum <= (T)1.0) {
426 // The point is inside the tetrahedron (3-Simplex) at (0,0,0)
427
428 // Determine which of (0,0,1), (0,1,0), (1,0,0) are closest.
429 uint_fast8_t aPoint = 1;
430 T aScore = xins;
431 uint_fast8_t bPoint = 2;
432 T bScore = yins;
433 if (aScore < bScore && zins > aScore) {
434 aScore = zins;
435 aPoint = 4;
436 } else if (aScore >= bScore && zins > bScore) {
437 bScore = zins;
438 bPoint = 4;
439 }
440
441 // Determine the two lattice points not part of the tetrahedron that may contribute.
442 // This depends on the closest two tetrahedral vertices, including (0,0,0).
443 T wins = (T)1.0 - inSum;
444 if (wins > aScore || wins > bScore) {
445 // (0,0,0) is one of the closest two tetrahedral vertices.
446
447 // The other closest vertex is the closer of a and b.
448 uint_fast8_t c = ((bScore > aScore) ? bPoint : aPoint);
449
450 if (c != 1) {
451 xsv_ext0 = xsb - 1;
452 xsv_ext1 = xsb;
453 dx_ext0 = dx0 + (T)1.0;
454 dx_ext1 = dx0;
455 } else {
456 xsv_ext0 = xsv_ext1 = xsb + 1;
457 dx_ext0 = dx_ext1 = dx0 - (T)1.0;
458 }
459
460 if (c != 2) {
461 ysv_ext0 = ysv_ext1 = ysb;
462 dy_ext0 = dy_ext1 = dy0;
463 if (c == 1) {
464 ysv_ext0 -= 1;
465 dy_ext0 += (T)1.0;
466 } else {
467 ysv_ext1 -= 1;
468 dy_ext1 += (T)1.0;
469 }
470 } else {
471 ysv_ext0 = ysv_ext1 = ysb + 1;
472 dy_ext0 = dy_ext1 = dy0 - (T)1.0;
473 }
474
475 if (c != 4) {
476 zsv_ext0 = zsb;
477 zsv_ext1 = zsb - 1;
478 dz_ext0 = dz0;
479 dz_ext1 = dz0 + (T)1.0;
480 } else {
481 zsv_ext0 = zsv_ext1 = zsb + 1;
482 dz_ext0 = dz_ext1 = dz0 - (T)1.0;
483 }
484 } else {
485 // (0,0,0) is not one of the closest two tetrahedral vertices.
486
487 // The two extra vertices are determined by the closest two.
488 uint_fast8_t c = (aPoint | bPoint);
489
490 if (c & 0x01) {
491 xsv_ext0 = xsv_ext1 = xsb + 1;
492 dx_ext0 = dx0 - (T)1.0 - (SQUISH_CONSTANT * (T)2.0);
493 dx_ext1 = dx0 - (T)1.0 - SQUISH_CONSTANT;
494 } else {
495 xsv_ext0 = xsb;
496 xsv_ext1 = xsb - 1;
497 dx_ext0 = dx0 - (SQUISH_CONSTANT * (T)2.0);
498 dx_ext1 = dx0 + (T)1.0 - SQUISH_CONSTANT;
499 }
500
501 if (c & 0x02) {
502 ysv_ext0 = ysv_ext1 = ysb + 1;
503 dy_ext0 = dy0 - (T)1.0 - (SQUISH_CONSTANT * (T)2.0);
504 dy_ext1 = dy0 - (T)1.0 - SQUISH_CONSTANT;
505 } else {
506 ysv_ext0 = ysb;
507 ysv_ext1 = ysb - 1;
508 dy_ext0 = dy0 - (SQUISH_CONSTANT * (T)2.0);
509 dy_ext1 = dy0 + (T)1.0 - SQUISH_CONSTANT;
510 }
511
512 if (c & 0x04) {
513 zsv_ext0 = zsv_ext1 = zsb + 1;
514 dz_ext0 = dz0 - (T)1.0 - (SQUISH_CONSTANT * (T)2.0);
515 dz_ext1 = dz0 - (T)1.0 - SQUISH_CONSTANT;
516 } else {
517 zsv_ext0 = zsb;
518 zsv_ext1 = zsb - 1;
519 dz_ext0 = dz0 - (SQUISH_CONSTANT * (T)2.0);
520 dz_ext1 = dz0 + (T)1.0 - SQUISH_CONSTANT;
521 }
522 }
523
524 // Contribution (0,0,0)
525 {
526 contr_m[0] = pow2(dx0) + pow2(dy0) + pow2(dz0);
527 contr_ext[0] = extrapolate(xsb, ysb, zsb, dx0, dy0, dz0);
528 }
529
530 // Contribution (0,0,1)
531 T dx1 = dx0 - (T)1.0 - SQUISH_CONSTANT;
532 T dy1 = dy0 - SQUISH_CONSTANT;
533 T dz1 = dz0 - SQUISH_CONSTANT;
534 contr_m[1] = pow2(dx1) + pow2(dy1) + pow2(dz1);
535 contr_ext[1] = extrapolate(xsb + 1, ysb, zsb, dx1, dy1, dz1);
536
537 // Contribution (0,1,0)
538 T dx2 = dx0 - SQUISH_CONSTANT;
539 T dy2 = dy0 - (T)1.0 - SQUISH_CONSTANT;
540 T dz2 = dz1;
541 contr_m[2] = pow2(dx2) + pow2(dy2) + pow2(dz2);
542 contr_ext[2] = extrapolate(xsb, ysb + 1, zsb, dx2, dy2, dz2);
543
544 // Contribution (1,0,0)
545 T dx3 = dx2;
546 T dy3 = dy1;
547 T dz3 = dz0 - (T)1.0 - SQUISH_CONSTANT;
548 contr_m[3] = pow2(dx3) + pow2(dy3) + pow2(dz3);
549 contr_ext[3] = extrapolate(xsb, ysb, zsb + 1, dx3, dy3, dz3);
550
551 contr_m[4] = contr_m[5] = contr_m[6] = 0.0;
552 contr_ext[4] = contr_ext[5] = contr_ext[6] = 0.0;
553
554 } else {
555 // The point is inside the tetrahedron (3-Simplex) at (1,1,1)
556
557 // Determine which two tetrahedral vertices are the closest
558 // out of (1,1,0), (1,0,1), and (0,1,1), but not (1,1,1).
559 uint_fast8_t aPoint = 6;
560 T aScore = xins;
561 uint_fast8_t bPoint = 5;
562 T bScore = yins;
563 if (aScore <= bScore && zins < bScore) {
564 bScore = zins;
565 bPoint = 3;
566 } else if (aScore > bScore && zins < aScore) {
567 aScore = zins;
568 aPoint = 3;
569 }
570
571 // Determine the two lattice points not part of the tetrahedron that may contribute.
572 // This depends on the closest two tetrahedral vertices, including (1,1,1).
573 T wins = 3.0 - inSum;
574 if (wins < aScore || wins < bScore) {
575 // (1,1,1) is one of the closest two tetrahedral vertices.
576
577 // The other closest vertex is the closest of a and b.
578 uint_fast8_t c = ((bScore < aScore) ? bPoint : aPoint);
579
580 if (c & 0x01) {
581 xsv_ext0 = xsb + 2;
582 xsv_ext1 = xsb + 1;
583 dx_ext0 = dx0 - (T)2.0 - (SQUISH_CONSTANT * (T)3.0);
584 dx_ext1 = dx0 - (T)1.0 - (SQUISH_CONSTANT * (T)3.0);
585 } else {
586 xsv_ext0 = xsv_ext1 = xsb;
587 dx_ext0 = dx_ext1 = dx0 - (SQUISH_CONSTANT * (T)3.0);
588 }
589
590 if (c & 0x02) {
591 ysv_ext0 = ysv_ext1 = ysb + 1;
592 dy_ext0 = dy_ext1 = dy0 - (T)1.0 - (SQUISH_CONSTANT * (T)3.0);
593 if (c & 0x01) {
594 ysv_ext1 += 1;
595 dy_ext1 -= (T)1.0;
596 } else {
597 ysv_ext0 += 1;
598 dy_ext0 -= (T)1.0;
599 }
600 } else {
601 ysv_ext0 = ysv_ext1 = ysb;
602 dy_ext0 = dy_ext1 = dy0 - (SQUISH_CONSTANT * (T)3.0);
603 }
604
605 if (c & 0x04) {
606 zsv_ext0 = zsb + 1;
607 zsv_ext1 = zsb + 2;
608 dz_ext0 = dz0 - (T)1.0 - (SQUISH_CONSTANT * (T)3.0);
609 dz_ext1 = dz0 - (T)2.0 - (SQUISH_CONSTANT * (T)3.0);
610 } else {
611 zsv_ext0 = zsv_ext1 = zsb;
612 dz_ext0 = dz_ext1 = dz0 - (SQUISH_CONSTANT * (T)3.0);
613 }
614 } else {
615 // (1,1,1) is not one of the closest two tetrahedral vertices.
616
617 // The two extra vertices are determined by the closest two.
618 uint_fast8_t c = aPoint & bPoint;
619
620 if (c & 0x01) {
621 xsv_ext0 = xsb + 1;
622 xsv_ext1 = xsb + 2;
623 dx_ext0 = dx0 - (T)1.0 - SQUISH_CONSTANT;
624 dx_ext1 = dx0 - (T)2.0 - (SQUISH_CONSTANT * (T)2.0);
625 } else {
626 xsv_ext0 = xsv_ext1 = xsb;
627 dx_ext0 = dx0 - SQUISH_CONSTANT;
628 dx_ext1 = dx0 - (SQUISH_CONSTANT * (T)2.0);
629 }
630
631 if (c & 0x02) {
632 ysv_ext0 = ysb + 1;
633 ysv_ext1 = ysb + 2;
634 dy_ext0 = dy0 - (T)1.0 - SQUISH_CONSTANT;
635 dy_ext1 = dy0 - (T)2.0 - (SQUISH_CONSTANT * (T)2.0);
636 } else {
637 ysv_ext0 = ysv_ext1 = ysb;
638 dy_ext0 = dy0 - SQUISH_CONSTANT;
639 dy_ext1 = dy0 - (SQUISH_CONSTANT * (T)2.0);
640 }
641
642 if (c & 0x04) {
643 zsv_ext0 = zsb + 1;
644 zsv_ext1 = zsb + 2;
645 dz_ext0 = dz0 - (T)1.0 - SQUISH_CONSTANT;
646 dz_ext1 = dz0 - (T)2.0 - (SQUISH_CONSTANT * (T)2.0);
647 } else {
648 zsv_ext0 = zsv_ext1 = zsb;
649 dz_ext0 = dz0 - SQUISH_CONSTANT;
650 dz_ext1 = dz0 - (SQUISH_CONSTANT * (T)2.0);
651 }
652 }
653
654 // Contribution (1,1,0)
655 T dx3 = dx0 - (T)1.0 - (SQUISH_CONSTANT * (T)2.0);
656 T dy3 = dy0 - (T)1.0 - (SQUISH_CONSTANT * (T)2.0);
657 T dz3 = dz0 - (SQUISH_CONSTANT * (T)2.0);
658 contr_m[3] = pow2(dx3) + pow2(dy3) + pow2(dz3);
659 contr_ext[3] = extrapolate(xsb + 1, ysb + 1, zsb, dx3, dy3, dz3);
660
661 // Contribution (1,0,1)
662 T dx2 = dx3;
663 T dy2 = dy0 - (SQUISH_CONSTANT * (T)2.0);
664 T dz2 = dz0 - (T)1.0 - (SQUISH_CONSTANT * (T)2.0);
665 contr_m[2] = pow2(dx2) + pow2(dy2) + pow2(dz2);
666 contr_ext[2] = extrapolate(xsb + 1, ysb, zsb + 1, dx2, dy2, dz2);
667
668 // Contribution (0,1,1)
669 {
670 T dx1 = dx0 - (SQUISH_CONSTANT * (T)2.0);
671 T dy1 = dy3;
672 T dz1 = dz2;
673 contr_m[1] = pow2(dx1) + pow2(dy1) + pow2(dz1);
674 contr_ext[1] = extrapolate(xsb, ysb + 1, zsb + 1, dx1, dy1, dz1);
675 }
676
677 // Contribution (1,1,1)
678 {
679 dx0 = dx0 - (T)1.0 - (SQUISH_CONSTANT * (T)3.0);
680 dy0 = dy0 - (T)1.0 - (SQUISH_CONSTANT * (T)3.0);
681 dz0 = dz0 - (T)1.0 - (SQUISH_CONSTANT * (T)3.0);
682 contr_m[0] = pow2(dx0) + pow2(dy0) + pow2(dz0);
683 contr_ext[0] = extrapolate(xsb + 1, ysb + 1, zsb + 1, dx0, dy0, dz0);
684 }
685
686 contr_m[4] = contr_m[5] = contr_m[6] = 0.0;
687 contr_ext[4] = contr_ext[5] = contr_ext[6] = 0.0;
688
689 }
690
691 // First extra vertex.
692 contr_m[7] = pow2(dx_ext0) + pow2(dy_ext0) + pow2(dz_ext0);
693 contr_ext[7] = extrapolate(xsv_ext0, ysv_ext0, zsv_ext0, dx_ext0, dy_ext0, dz_ext0);
694
695 // Second extra vertex.
696 contr_m[8] = pow2(dx_ext1) + pow2(dy_ext1) + pow2(dz_ext1);
697 contr_ext[8] = extrapolate(xsv_ext1, ysv_ext1, zsv_ext1, dx_ext1, dy_ext1, dz_ext1);
698
699 T value = 0.0;
700 for (int i=0; i<9; ++i) {
701 value += pow4(std::max((T)2.0 - contr_m[i], (T)0.0)) * contr_ext[i];
702 }
703
704 return (value * NORM_CONSTANT);
705 }
706
707 template double OSNoise::extrapolate(const OSNoise::inttype xsb, const OSNoise::inttype ysb, const OSNoise::inttype zsb,
708 const double dx, const double dy, const double dz) const;
709 template double OSNoise::extrapolate(const OSNoise::inttype xsb, const OSNoise::inttype ysb, const OSNoise::inttype zsb,
710 const double dx, const double dy, const double dz,
711 double (&de) [3]) const;
712
713 template double OSNoise::eval(const double x, const double y, const double z) const;
714
715 } // namespace OSN
716