Line | Branch | Exec | Source |
---|---|---|---|
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> | ||
19 | inline T pow4 (T x) | ||
20 | { | ||
21 | 10836 | x *= x; | |
22 | 10836 | return x*x; | |
23 | } | ||
24 | |||
25 | template <typename T> | ||
26 | inline T pow2 (T x) | ||
27 | { | ||
28 | 2408 | return x*x; | |
29 | } | ||
30 | |||
31 | template <typename T> | ||
32 | inline OSNoise::inttype fastFloori (T x) | ||
33 | { | ||
34 | 3612 | OSNoise::inttype ip = (OSNoise::inttype)x; | |
35 | |||
36 |
4/4✓ Branch 0 taken 444 times.
✓ Branch 1 taken 760 times.
✓ Branch 2 taken 760 times.
✓ Branch 3 taken 444 times.
|
1204 | if (x < 0.0) --ip; |
37 | |||
38 | return ip; | ||
39 | } | ||
40 | |||
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 | 14 | 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> | ||
63 | ✗ | 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 | 2408 | unsigned int index = mPermGradIndex[(mPerm[(mPerm[xsb & 0xFF] + ysb) & 0xFF] + zsb) & 0xFF]; | |
71 | 2408 | return sGradients[index] * dx + | |
72 | 2408 | sGradients[index + 1] * dy + | |
73 | 2408 | sGradients[index + 2] * dz; | |
74 | |||
75 | } | ||
76 | |||
77 | template <typename T> | ||
78 | ✗ | 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 | |||
92 | 14 | OSNoise::OSNoise(OSNoise::inttype seed) | |
93 | { | ||
94 | int source [256]; | ||
95 |
2/2✓ Branch 0 taken 3584 times.
✓ Branch 1 taken 14 times.
|
3598 | for (int i = 0; i < 256; ++i) { source[i] = i; } |
96 | LCG_STEP(seed); | ||
97 | LCG_STEP(seed); | ||
98 | LCG_STEP(seed); | ||
99 |
2/2✓ Branch 0 taken 3584 times.
✓ Branch 1 taken 14 times.
|
3598 | for (int i = 255; i >= 0; --i) { |
100 | LCG_STEP(seed); | ||
101 | 3584 | int r = (int)((seed + 31) % (i + 1)); | |
102 |
2/2✓ Branch 0 taken 1610 times.
✓ Branch 1 taken 1974 times.
|
3584 | if (r < 0) { r += (i + 1); } |
103 | 3584 | mPerm[i] = source[r]; | |
104 | 3584 | mPermGradIndex[i] = (int)((mPerm[i] % (72 / 3)) * 3); | |
105 | 3584 | source[r] = source[i]; | |
106 | } | ||
107 | 14 | } | |
108 | |||
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> | ||
119 | 1204 | 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 | 1204 | T stretchOffset = (x + y + z) * STRETCH_CONSTANT; | |
137 | 1204 | T xs = x + stretchOffset; | |
138 | 1204 | T ys = y + stretchOffset; | |
139 |
2/2✓ Branch 0 taken 380 times.
✓ Branch 1 taken 824 times.
|
1204 | 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 | 1204 | T xsbd = (T)xsb; | |
155 | 1204 | T ysbd = (T)ysb; | |
156 | 1204 | T zsbd = (T)zsb; | |
157 | #endif | ||
158 | |||
159 | // Skew out to get actual coordinates of rhombohedron origin. | ||
160 | 1204 | T squishOffset = (xsbd + ysbd + zsbd) * SQUISH_CONSTANT; | |
161 | 1204 | T xb = xsbd + squishOffset; | |
162 | 1204 | T yb = ysbd + squishOffset; | |
163 | 1204 | T zb = zsbd + squishOffset; | |
164 | |||
165 | // Positions relative to origin point. | ||
166 | 1204 | dx0 = x - xb; | |
167 | 1204 | dy0 = y - yb; | |
168 | 1204 | dz0 = z - zb; | |
169 | |||
170 | // Compute simplectic lattice coordinates relative to rhombohedral origin. | ||
171 | 1204 | xins = xs - xsbd; | |
172 | 1204 | yins = ys - ysbd; | |
173 | 1204 | 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 | 1204 | T inSum = xins + yins + zins; | |
184 | |||
185 |
4/4✓ Branch 0 taken 508 times.
✓ Branch 1 taken 696 times.
✓ Branch 2 taken 380 times.
✓ Branch 3 taken 128 times.
|
1204 | 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 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 316 times.
|
380 | if (p1 <= (T)1.0) { |
198 | 64 | aScore = (T)1.0 - p1; | |
199 | aPoint = 4; | ||
200 | aIsFurtherSide = false; | ||
201 | } else { | ||
202 | 316 | 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 | 380 | T p2 = xins + zins; | |
209 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 316 times.
|
380 | if (p2 <= (T)1.0) { |
210 | 64 | bScore = (T)1.0 - p2; | |
211 | bPoint = 2; | ||
212 | bIsFurtherSide = false; | ||
213 | } else { | ||
214 | 316 | 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 | 380 | T p3 = yins + zins; | |
222 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 380 times.
|
380 | 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 | 380 | T score = (T)1.0 - p3; | |
235 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 380 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
380 | if (aScore > bScore && bScore < score) { |
236 | bScore = score; | ||
237 | bPoint = 1; | ||
238 | bIsFurtherSide = false; | ||
239 |
2/4✓ Branch 0 taken 380 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 380 times.
|
380 | } 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 |
1/2✓ Branch 0 taken 380 times.
✗ Branch 1 not taken.
|
380 | if (aIsFurtherSide == bIsFurtherSide) { |
249 |
2/2✓ Branch 0 taken 316 times.
✓ Branch 1 taken 64 times.
|
380 | if (aIsFurtherSide) { |
250 | // Both closest points on (1,1,1) side. | ||
251 | |||
252 | // One of the two extra points is (1,1,1) | ||
253 | 316 | xsv_ext0 = xsb + 1; | |
254 | 316 | ysv_ext0 = ysb + 1; | |
255 | 316 | zsv_ext0 = zsb + 1; | |
256 | 316 | dx_ext0 = dx0 - (T)1.0 - (SQUISH_CONSTANT * (T)3.0); | |
257 | 316 | dy_ext0 = dy0 - (T)1.0 - (SQUISH_CONSTANT * (T)3.0); | |
258 | 316 | dz_ext0 = dz0 - (T)1.0 - (SQUISH_CONSTANT * (T)3.0); | |
259 | |||
260 | // Other extra point is based on the shared axis. | ||
261 | 316 | uint_fast8_t c = aPoint & bPoint; | |
262 |
1/2✓ Branch 0 taken 316 times.
✗ Branch 1 not taken.
|
316 | if (c & 0x01) { |
263 | 316 | xsv_ext1 = xsb + 2; | |
264 | ysv_ext1 = ysb; | ||
265 | zsv_ext1 = zsb; | ||
266 | 316 | dx_ext1 = dx0 - (T)2.0 - (SQUISH_CONSTANT * (T)2.0); | |
267 | 316 | dy_ext1 = dy0 - (SQUISH_CONSTANT * (T)2.0); | |
268 | 316 | 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 | 64 | uint_fast8_t c = aPoint | bPoint; | |
297 |
1/2✓ Branch 0 taken 64 times.
✗ Branch 1 not taken.
|
64 | if (!(c & 0x01)) { |
298 | 64 | xsv_ext1 = xsb - 1; | |
299 | 64 | ysv_ext1 = ysb + 1; | |
300 | 64 | zsv_ext1 = zsb + 1; | |
301 | 64 | dx_ext1 = dx0 + (T)1.0 - SQUISH_CONSTANT; | |
302 | 64 | dy_ext1 = dy0 - (T)1.0 - SQUISH_CONSTANT; | |
303 | 64 | 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 | 380 | contr_m[0] = contr_ext[0] = 0.0; | |
382 | |||
383 | // Contribution (0,0,1). | ||
384 | 380 | T dx1 = dx0 - (T)1.0 - SQUISH_CONSTANT; | |
385 | 380 | T dy1 = dy0 - SQUISH_CONSTANT; | |
386 | 380 | T dz1 = dz0 - SQUISH_CONSTANT; | |
387 | 380 | contr_m[1] = pow2(dx1) + pow2(dy1) + pow2(dz1); | |
388 | 380 | contr_ext[1] = extrapolate(xsb + 1, ysb, zsb, dx1, dy1, dz1); | |
389 | |||
390 | // Contribution (0,1,0). | ||
391 | 380 | T dx2 = dx0 - SQUISH_CONSTANT; | |
392 | 380 | T dy2 = dy0 - (T)1.0 - SQUISH_CONSTANT; | |
393 | T dz2 = dz1; | ||
394 | 380 | contr_m[2] = pow2(dx2) + pow2(dy2) + pow2(dz2); | |
395 | 380 | 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 | 380 | T dz3 = dz0 - (T)1.0 - SQUISH_CONSTANT; | |
401 | 380 | contr_m[3] = pow2(dx3) + pow2(dy3) + pow2(dz3); | |
402 | 380 | contr_ext[3] = extrapolate(xsb, ysb, zsb + 1, dx3, dy3, dz3); | |
403 | |||
404 | // Contribution (1,1,0). | ||
405 | 380 | T dx4 = dx0 - (T)1.0 - (SQUISH_CONSTANT * (T)2.0); | |
406 | 380 | T dy4 = dy0 - (T)1.0 - (SQUISH_CONSTANT * (T)2.0); | |
407 | 380 | T dz4 = dz0 - (SQUISH_CONSTANT * (T)2.0); | |
408 | 380 | contr_m[4] = pow2(dx4) + pow2(dy4) + pow2(dz4); | |
409 | 380 | contr_ext[4] = extrapolate(xsb + 1, ysb + 1, zsb, dx4, dy4, dz4); | |
410 | |||
411 | // Contribution (1,0,1). | ||
412 | T dx5 = dx4; | ||
413 | 380 | T dy5 = dy0 - (SQUISH_CONSTANT * (T)2.0); | |
414 | 380 | T dz5 = dz0 - (T)1.0 - (SQUISH_CONSTANT * (T)2.0); | |
415 | 380 | contr_m[5] = pow2(dx5) + pow2(dy5) + pow2(dz5); | |
416 | 380 | contr_ext[5] = extrapolate(xsb + 1, ysb, zsb + 1, dx5, dy5, dz5); | |
417 | |||
418 | // Contribution (0,1,1). | ||
419 | 380 | T dx6 = dx0 - (SQUISH_CONSTANT * (T)2.0); | |
420 | T dy6 = dy4; | ||
421 | T dz6 = dz5; | ||
422 | 380 | contr_m[6] = pow2(dx6) + pow2(dy6) + pow2(dz6); | |
423 | 380 | contr_ext[6] = extrapolate(xsb, ysb + 1, zsb + 1, dx6, dy6, dz6); | |
424 | |||
425 |
2/2✓ Branch 0 taken 696 times.
✓ Branch 1 taken 128 times.
|
824 | } 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 |
3/4✓ Branch 0 taken 316 times.
✓ Branch 1 taken 380 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 316 times.
|
696 | if (aScore < bScore && zins > aScore) { |
434 | aScore = zins; | ||
435 | aPoint = 4; | ||
436 |
3/4✓ Branch 0 taken 380 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 316 times.
✓ Branch 3 taken 64 times.
|
380 | } 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 | 696 | T wins = (T)1.0 - inSum; | |
444 |
3/4✓ Branch 0 taken 632 times.
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 632 times.
|
696 | 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 |
1/2✓ Branch 0 taken 64 times.
✗ Branch 1 not taken.
|
64 | uint_fast8_t c = ((bScore > aScore) ? bPoint : aPoint); |
449 | |||
450 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 64 times.
|
64 | 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 | 64 | xsv_ext0 = xsv_ext1 = xsb + 1; | |
457 | 64 | dx_ext0 = dx_ext1 = dx0 - (T)1.0; | |
458 | } | ||
459 | |||
460 |
1/2✓ Branch 0 taken 64 times.
✗ Branch 1 not taken.
|
64 | if (c != 2) { |
461 | ysv_ext0 = ysv_ext1 = ysb; | ||
462 | dy_ext0 = dy_ext1 = dy0; | ||
463 |
1/2✓ Branch 0 taken 64 times.
✗ Branch 1 not taken.
|
64 | if (c == 1) { |
464 | 64 | ysv_ext0 -= 1; | |
465 | 64 | 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 |
1/2✓ Branch 0 taken 64 times.
✗ Branch 1 not taken.
|
64 | if (c != 4) { |
476 | zsv_ext0 = zsb; | ||
477 | 64 | zsv_ext1 = zsb - 1; | |
478 | dz_ext0 = dz0; | ||
479 | 64 | 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 | 632 | uint_fast8_t c = (aPoint | bPoint); | |
489 | |||
490 |
2/2✓ Branch 0 taken 316 times.
✓ Branch 1 taken 316 times.
|
632 | if (c & 0x01) { |
491 | 316 | xsv_ext0 = xsv_ext1 = xsb + 1; | |
492 | 316 | dx_ext0 = dx0 - (T)1.0 - (SQUISH_CONSTANT * (T)2.0); | |
493 | 316 | dx_ext1 = dx0 - (T)1.0 - SQUISH_CONSTANT; | |
494 | } else { | ||
495 | xsv_ext0 = xsb; | ||
496 | 316 | xsv_ext1 = xsb - 1; | |
497 | 316 | dx_ext0 = dx0 - (SQUISH_CONSTANT * (T)2.0); | |
498 | 316 | dx_ext1 = dx0 + (T)1.0 - SQUISH_CONSTANT; | |
499 | } | ||
500 | |||
501 |
2/2✓ Branch 0 taken 316 times.
✓ Branch 1 taken 316 times.
|
632 | if (c & 0x02) { |
502 | 316 | ysv_ext0 = ysv_ext1 = ysb + 1; | |
503 | 316 | dy_ext0 = dy0 - (T)1.0 - (SQUISH_CONSTANT * (T)2.0); | |
504 | 316 | dy_ext1 = dy0 - (T)1.0 - SQUISH_CONSTANT; | |
505 | } else { | ||
506 | ysv_ext0 = ysb; | ||
507 | 316 | ysv_ext1 = ysb - 1; | |
508 | 316 | dy_ext0 = dy0 - (SQUISH_CONSTANT * (T)2.0); | |
509 | 316 | dy_ext1 = dy0 + (T)1.0 - SQUISH_CONSTANT; | |
510 | } | ||
511 | |||
512 |
1/2✓ Branch 0 taken 632 times.
✗ Branch 1 not taken.
|
632 | if (c & 0x04) { |
513 | 632 | zsv_ext0 = zsv_ext1 = zsb + 1; | |
514 | 632 | dz_ext0 = dz0 - (T)1.0 - (SQUISH_CONSTANT * (T)2.0); | |
515 | 632 | 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 | 696 | contr_m[0] = pow2(dx0) + pow2(dy0) + pow2(dz0); | |
527 | 696 | contr_ext[0] = extrapolate(xsb, ysb, zsb, dx0, dy0, dz0); | |
528 | } | ||
529 | |||
530 | // Contribution (0,0,1) | ||
531 | 696 | T dx1 = dx0 - (T)1.0 - SQUISH_CONSTANT; | |
532 | 696 | T dy1 = dy0 - SQUISH_CONSTANT; | |
533 | 696 | T dz1 = dz0 - SQUISH_CONSTANT; | |
534 | 696 | contr_m[1] = pow2(dx1) + pow2(dy1) + pow2(dz1); | |
535 | 696 | contr_ext[1] = extrapolate(xsb + 1, ysb, zsb, dx1, dy1, dz1); | |
536 | |||
537 | // Contribution (0,1,0) | ||
538 | 696 | T dx2 = dx0 - SQUISH_CONSTANT; | |
539 | 696 | T dy2 = dy0 - (T)1.0 - SQUISH_CONSTANT; | |
540 | T dz2 = dz1; | ||
541 | 696 | contr_m[2] = pow2(dx2) + pow2(dy2) + pow2(dz2); | |
542 | 696 | 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 | 696 | T dz3 = dz0 - (T)1.0 - SQUISH_CONSTANT; | |
548 | 696 | contr_m[3] = pow2(dx3) + pow2(dy3) + pow2(dz3); | |
549 | 696 | contr_ext[3] = extrapolate(xsb, ysb, zsb + 1, dx3, dy3, dz3); | |
550 | |||
551 | 696 | contr_m[4] = contr_m[5] = contr_m[6] = 0.0; | |
552 | 696 | 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 |
3/4✓ Branch 0 taken 128 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 64 times.
✓ Branch 3 taken 64 times.
|
128 | if (aScore <= bScore && zins < bScore) { |
564 | bScore = zins; | ||
565 | bPoint = 3; | ||
566 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
64 | } 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 | 128 | T wins = 3.0 - inSum; | |
574 |
3/4✓ Branch 0 taken 64 times.
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 64 times.
|
128 | 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 |
1/2✓ Branch 0 taken 64 times.
✗ Branch 1 not taken.
|
64 | uint_fast8_t c = ((bScore < aScore) ? bPoint : aPoint); |
579 | |||
580 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 64 times.
|
64 | 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 | 64 | dx_ext0 = dx_ext1 = dx0 - (SQUISH_CONSTANT * (T)3.0); | |
588 | } | ||
589 | |||
590 |
1/2✓ Branch 0 taken 64 times.
✗ Branch 1 not taken.
|
64 | if (c & 0x02) { |
591 | 64 | ysv_ext0 = ysv_ext1 = ysb + 1; | |
592 | 64 | dy_ext0 = dy_ext1 = dy0 - (T)1.0 - (SQUISH_CONSTANT * (T)3.0); | |
593 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 64 times.
|
64 | if (c & 0x01) { |
594 | ✗ | ysv_ext1 += 1; | |
595 | ✗ | dy_ext1 -= (T)1.0; | |
596 | } else { | ||
597 | 64 | ysv_ext0 += 1; | |
598 | 64 | 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 |
1/2✓ Branch 0 taken 64 times.
✗ Branch 1 not taken.
|
64 | if (c & 0x04) { |
606 | 64 | zsv_ext0 = zsb + 1; | |
607 | 64 | zsv_ext1 = zsb + 2; | |
608 | 64 | dz_ext0 = dz0 - (T)1.0 - (SQUISH_CONSTANT * (T)3.0); | |
609 | 64 | 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 | 64 | uint_fast8_t c = aPoint & bPoint; | |
619 | |||
620 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 64 times.
|
64 | 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 | 64 | dx_ext0 = dx0 - SQUISH_CONSTANT; | |
628 | 64 | dx_ext1 = dx0 - (SQUISH_CONSTANT * (T)2.0); | |
629 | } | ||
630 | |||
631 |
1/2✓ Branch 0 taken 64 times.
✗ Branch 1 not taken.
|
64 | if (c & 0x02) { |
632 | 64 | ysv_ext0 = ysb + 1; | |
633 | 64 | ysv_ext1 = ysb + 2; | |
634 | 64 | dy_ext0 = dy0 - (T)1.0 - SQUISH_CONSTANT; | |
635 | 64 | 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 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 64 times.
|
64 | 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 | 64 | dz_ext0 = dz0 - SQUISH_CONSTANT; | |
650 | 64 | dz_ext1 = dz0 - (SQUISH_CONSTANT * (T)2.0); | |
651 | } | ||
652 | } | ||
653 | |||
654 | // Contribution (1,1,0) | ||
655 | 128 | T dx3 = dx0 - (T)1.0 - (SQUISH_CONSTANT * (T)2.0); | |
656 | 128 | T dy3 = dy0 - (T)1.0 - (SQUISH_CONSTANT * (T)2.0); | |
657 | 128 | T dz3 = dz0 - (SQUISH_CONSTANT * (T)2.0); | |
658 | 128 | contr_m[3] = pow2(dx3) + pow2(dy3) + pow2(dz3); | |
659 | 128 | contr_ext[3] = extrapolate(xsb + 1, ysb + 1, zsb, dx3, dy3, dz3); | |
660 | |||
661 | // Contribution (1,0,1) | ||
662 | T dx2 = dx3; | ||
663 | 128 | T dy2 = dy0 - (SQUISH_CONSTANT * (T)2.0); | |
664 | 128 | T dz2 = dz0 - (T)1.0 - (SQUISH_CONSTANT * (T)2.0); | |
665 | 128 | contr_m[2] = pow2(dx2) + pow2(dy2) + pow2(dz2); | |
666 | 128 | contr_ext[2] = extrapolate(xsb + 1, ysb, zsb + 1, dx2, dy2, dz2); | |
667 | |||
668 | // Contribution (0,1,1) | ||
669 | { | ||
670 | 128 | T dx1 = dx0 - (SQUISH_CONSTANT * (T)2.0); | |
671 | T dy1 = dy3; | ||
672 | T dz1 = dz2; | ||
673 | 128 | contr_m[1] = pow2(dx1) + pow2(dy1) + pow2(dz1); | |
674 | 128 | contr_ext[1] = extrapolate(xsb, ysb + 1, zsb + 1, dx1, dy1, dz1); | |
675 | } | ||
676 | |||
677 | // Contribution (1,1,1) | ||
678 | { | ||
679 | 128 | dx0 = dx0 - (T)1.0 - (SQUISH_CONSTANT * (T)3.0); | |
680 | 128 | dy0 = dy0 - (T)1.0 - (SQUISH_CONSTANT * (T)3.0); | |
681 | 128 | dz0 = dz0 - (T)1.0 - (SQUISH_CONSTANT * (T)3.0); | |
682 | 128 | contr_m[0] = pow2(dx0) + pow2(dy0) + pow2(dz0); | |
683 | 128 | contr_ext[0] = extrapolate(xsb + 1, ysb + 1, zsb + 1, dx0, dy0, dz0); | |
684 | } | ||
685 | |||
686 | 128 | contr_m[4] = contr_m[5] = contr_m[6] = 0.0; | |
687 | 128 | contr_ext[4] = contr_ext[5] = contr_ext[6] = 0.0; | |
688 | |||
689 | } | ||
690 | |||
691 | // First extra vertex. | ||
692 | 1204 | contr_m[7] = pow2(dx_ext0) + pow2(dy_ext0) + pow2(dz_ext0); | |
693 | 1204 | contr_ext[7] = extrapolate(xsv_ext0, ysv_ext0, zsv_ext0, dx_ext0, dy_ext0, dz_ext0); | |
694 | |||
695 | // Second extra vertex. | ||
696 | 1204 | contr_m[8] = pow2(dx_ext1) + pow2(dy_ext1) + pow2(dz_ext1); | |
697 | 1204 | contr_ext[8] = extrapolate(xsv_ext1, ysv_ext1, zsv_ext1, dx_ext1, dy_ext1, dz_ext1); | |
698 | |||
699 | T value = 0.0; | ||
700 |
2/2✓ Branch 0 taken 10836 times.
✓ Branch 1 taken 1204 times.
|
12040 | for (int i=0; i<9; ++i) { |
701 |
2/2✓ Branch 0 taken 1268 times.
✓ Branch 1 taken 9568 times.
|
12104 | value += pow4(std::max((T)2.0 - contr_m[i], (T)0.0)) * contr_ext[i]; |
702 | } | ||
703 | |||
704 | 1204 | return (value * NORM_CONSTANT); | |
705 | } | ||
706 | |||
707 | template OPENVDB_AX_API 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 OPENVDB_AX_API 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 OPENVDB_AX_API double OSNoise::eval(const double x, const double y, const double z) const; | ||
714 | |||
715 | } // namespace OSN | ||
716 |