3 * Copyright (C) 2010-2014 celeron55, Perttu Ahola <celeron55@gmail.com>
4 * Copyright (C) 2010-2014 kwolekr, Ryan Kwolek <kwolekr@minetest.net>
7 * Redistribution and use in source and binary forms, with or without modification, are
8 * permitted provided that the following conditions are met:
9 * 1. Redistributions of source code must retain the above copyright notice, this list of
10 * conditions and the following disclaimer.
11 * 2. Redistributions in binary form must reproduce the above copyright notice, this list
12 * of conditions and the following disclaimer in the documentation and/or other materials
13 * provided with the distribution.
15 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND ANY EXPRESS OR IMPLIED
16 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
17 * FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR
18 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
19 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
20 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
21 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
22 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
23 * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 #include <cstring> // memset
31 #include "util/numeric.h"
32 #include "util/string.h"
33 #include "exceptions.h"
35 #define NOISE_MAGIC_X 1619
36 #define NOISE_MAGIC_Y 31337
37 #define NOISE_MAGIC_Z 52591
38 #define NOISE_MAGIC_SEED 1013
40 typedef float (*Interp2dFxn)(
41 float v00, float v10, float v01, float v11,
44 typedef float (*Interp3dFxn)(
45 float v000, float v100, float v010, float v110,
46 float v001, float v101, float v011, float v111,
47 float x, float y, float z);
49 float cos_lookup[16] = {
50 1.0f, 0.9238f, 0.7071f, 0.3826f, .0f, -0.3826f, -0.7071f, -0.9238f,
51 1.0f, -0.9238f, -0.7071f, -0.3826f, .0f, 0.3826f, 0.7071f, 0.9238f
54 FlagDesc flagdesc_noiseparams[] = {
55 {"defaults", NOISE_FLAG_DEFAULTS},
56 {"eased", NOISE_FLAG_EASED},
57 {"absvalue", NOISE_FLAG_ABSVALUE},
58 {"pointbuffer", NOISE_FLAG_POINTBUFFER},
59 {"simplex", NOISE_FLAG_SIMPLEX},
63 ///////////////////////////////////////////////////////////////////////////////
65 PcgRandom::PcgRandom(u64 state, u64 seq)
70 void PcgRandom::seed(u64 state, u64 seq)
73 m_inc = (seq << 1u) | 1u;
82 u64 oldstate = m_state;
83 m_state = oldstate * 6364136223846793005ULL + m_inc;
85 u32 xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;
86 u32 rot = oldstate >> 59u;
87 return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));
91 u32 PcgRandom::range(u32 bound)
93 // If the bound is 0, we cover the whole RNG's range
98 This is an optimization of the expression:
99 0x100000000ull % bound
100 since 64-bit modulo operations typically much slower than 32.
102 u32 threshold = -bound % bound;
106 If the bound is not a multiple of the RNG's range, it may cause bias,
107 e.g. a RNG has a range from 0 to 3 and we take want a number 0 to 2.
108 Using rand() % 3, the number 0 would be twice as likely to appear.
109 With a very large RNG range, the effect becomes less prevalent but
112 This can be solved by modifying the range of the RNG to become a
113 multiple of bound by dropping values above the a threshold.
115 In our example, threshold == 4 % 3 == 1, so reject values < 1
116 (that is, 0), thus making the range == 3 with no bias.
118 This loop may look dangerous, but will always terminate due to the
119 RNG's property of uniformity.
121 while ((r = next()) < threshold)
128 s32 PcgRandom::range(s32 min, s32 max)
131 throw PrngException("Invalid range (max < min)");
133 // We have to cast to s64 because otherwise this could overflow,
134 // and signed overflow is undefined behavior.
135 u32 bound = (s64)max - (s64)min + 1;
136 return range(bound) + min;
140 void PcgRandom::bytes(void *out, size_t len)
142 u8 *outb = (u8 *)out;
147 if (bytes_left == 0) {
148 bytes_left = sizeof(u32);
160 s32 PcgRandom::randNormalDist(s32 min, s32 max, int num_trials)
163 for (int i = 0; i != num_trials; i++)
164 accum += range(min, max);
165 return myround((float)accum / num_trials);
168 ///////////////////////////////////////////////////////////////////////////////
170 float noise2d(int x, int y, s32 seed)
172 unsigned int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y
173 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
175 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
176 return 1.f - (float)(int)n / 0x40000000;
180 float noise3d(int x, int y, int z, s32 seed)
182 unsigned int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y + NOISE_MAGIC_Z * z
183 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
185 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
186 return 1.f - (float)(int)n / 0x40000000;
190 inline float dotProduct(float vx, float vy, float wx, float wy)
192 return vx * wx + vy * wy;
196 inline float linearInterpolation(float v0, float v1, float t)
198 return v0 + (v1 - v0) * t;
202 inline float biLinearInterpolation(
203 float v00, float v10,
204 float v01, float v11,
207 float tx = easeCurve(x);
208 float ty = easeCurve(y);
209 float u = linearInterpolation(v00, v10, tx);
210 float v = linearInterpolation(v01, v11, tx);
211 return linearInterpolation(u, v, ty);
215 inline float biLinearInterpolationNoEase(
216 float v00, float v10,
217 float v01, float v11,
220 float u = linearInterpolation(v00, v10, x);
221 float v = linearInterpolation(v01, v11, x);
222 return linearInterpolation(u, v, y);
226 float triLinearInterpolation(
227 float v000, float v100, float v010, float v110,
228 float v001, float v101, float v011, float v111,
229 float x, float y, float z)
231 float tx = easeCurve(x);
232 float ty = easeCurve(y);
233 float tz = easeCurve(z);
234 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, tx, ty);
235 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, tx, ty);
236 return linearInterpolation(u, v, tz);
239 float triLinearInterpolationNoEase(
240 float v000, float v100, float v010, float v110,
241 float v001, float v101, float v011, float v111,
242 float x, float y, float z)
244 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, x, y);
245 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, x, y);
246 return linearInterpolation(u, v, z);
249 float noise2d_gradient(float x, float y, s32 seed, bool eased)
251 // Calculate the integer coordinates
254 // Calculate the remaining part of the coordinates
255 float xl = x - (float)x0;
256 float yl = y - (float)y0;
257 // Get values for corners of square
258 float v00 = noise2d(x0, y0, seed);
259 float v10 = noise2d(x0+1, y0, seed);
260 float v01 = noise2d(x0, y0+1, seed);
261 float v11 = noise2d(x0+1, y0+1, seed);
264 return biLinearInterpolation(v00, v10, v01, v11, xl, yl);
266 return biLinearInterpolationNoEase(v00, v10, v01, v11, xl, yl);
270 float noise3d_gradient(float x, float y, float z, s32 seed, bool eased)
272 // Calculate the integer coordinates
276 // Calculate the remaining part of the coordinates
277 float xl = x - (float)x0;
278 float yl = y - (float)y0;
279 float zl = z - (float)z0;
280 // Get values for corners of cube
281 float v000 = noise3d(x0, y0, z0, seed);
282 float v100 = noise3d(x0 + 1, y0, z0, seed);
283 float v010 = noise3d(x0, y0 + 1, z0, seed);
284 float v110 = noise3d(x0 + 1, y0 + 1, z0, seed);
285 float v001 = noise3d(x0, y0, z0 + 1, seed);
286 float v101 = noise3d(x0 + 1, y0, z0 + 1, seed);
287 float v011 = noise3d(x0, y0 + 1, z0 + 1, seed);
288 float v111 = noise3d(x0 + 1, y0 + 1, z0 + 1, seed);
291 return triLinearInterpolation(
292 v000, v100, v010, v110,
293 v001, v101, v011, v111,
297 return triLinearInterpolationNoEase(
298 v000, v100, v010, v110,
299 v001, v101, v011, v111,
304 float noise2d_perlin(float x, float y, s32 seed,
305 int octaves, float persistence, bool eased)
310 for (int i = 0; i < octaves; i++)
312 a += g * noise2d_gradient(x * f, y * f, seed + i, eased);
320 float noise2d_perlin_abs(float x, float y, s32 seed,
321 int octaves, float persistence, bool eased)
326 for (int i = 0; i < octaves; i++) {
327 a += g * std::fabs(noise2d_gradient(x * f, y * f, seed + i, eased));
335 float noise3d_perlin(float x, float y, float z, s32 seed,
336 int octaves, float persistence, bool eased)
341 for (int i = 0; i < octaves; i++) {
342 a += g * noise3d_gradient(x * f, y * f, z * f, seed + i, eased);
350 float noise3d_perlin_abs(float x, float y, float z, s32 seed,
351 int octaves, float persistence, bool eased)
356 for (int i = 0; i < octaves; i++) {
357 a += g * std::fabs(noise3d_gradient(x * f, y * f, z * f, seed + i, eased));
365 float contour(float v)
374 ///////////////////////// [ New noise ] ////////////////////////////
377 float NoisePerlin2D(NoiseParams *np, float x, float y, s32 seed)
387 for (size_t i = 0; i < np->octaves; i++) {
388 float noiseval = noise2d_gradient(x * f, y * f, seed + i,
389 np->flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED));
391 if (np->flags & NOISE_FLAG_ABSVALUE)
392 noiseval = std::fabs(noiseval);
399 return np->offset + a * np->scale;
403 float NoisePerlin3D(NoiseParams *np, float x, float y, float z, s32 seed)
414 for (size_t i = 0; i < np->octaves; i++) {
415 float noiseval = noise3d_gradient(x * f, y * f, z * f, seed + i,
416 np->flags & NOISE_FLAG_EASED);
418 if (np->flags & NOISE_FLAG_ABSVALUE)
419 noiseval = std::fabs(noiseval);
426 return np->offset + a * np->scale;
430 Noise::Noise(NoiseParams *np_, s32 seed, u32 sx, u32 sy, u32 sz)
432 memcpy(&np, np_, sizeof(np));
444 delete[] gradient_buf;
445 delete[] persist_buf;
451 void Noise::allocBuffers()
460 this->noise_buf = NULL;
461 resizeNoiseBuf(sz > 1);
463 delete[] gradient_buf;
464 delete[] persist_buf;
468 size_t bufsize = sx * sy * sz;
469 this->persist_buf = NULL;
470 this->gradient_buf = new float[bufsize];
471 this->result = new float[bufsize];
472 } catch (std::bad_alloc &e) {
473 throw InvalidNoiseParamsException();
478 void Noise::setSize(u32 sx, u32 sy, u32 sz)
488 void Noise::setSpreadFactor(v3f spread)
490 this->np.spread = spread;
492 resizeNoiseBuf(sz > 1);
496 void Noise::setOctaves(int octaves)
498 this->np.octaves = octaves;
500 resizeNoiseBuf(sz > 1);
504 void Noise::resizeNoiseBuf(bool is3d)
506 //maximum possible spread value factor
507 float ofactor = (np.lacunarity > 1.0) ?
508 pow(np.lacunarity, np.octaves - 1) :
511 // noise lattice point count
512 // (int)(sz * spread * ofactor) is # of lattice points crossed due to length
513 float num_noise_points_x = sx * ofactor / np.spread.X;
514 float num_noise_points_y = sy * ofactor / np.spread.Y;
515 float num_noise_points_z = sz * ofactor / np.spread.Z;
517 // protect against obviously invalid parameters
518 if (num_noise_points_x > 1000000000.f ||
519 num_noise_points_y > 1000000000.f ||
520 num_noise_points_z > 1000000000.f)
521 throw InvalidNoiseParamsException();
523 // + 2 for the two initial endpoints
524 // + 1 for potentially crossing a boundary due to offset
525 size_t nlx = (size_t)std::ceil(num_noise_points_x) + 3;
526 size_t nly = (size_t)std::ceil(num_noise_points_y) + 3;
527 size_t nlz = is3d ? (size_t)std::ceil(num_noise_points_z) + 3 : 1;
531 noise_buf = new float[nlx * nly * nlz];
532 } catch (std::bad_alloc &e) {
533 throw InvalidNoiseParamsException();
539 * NB: This algorithm is not optimal in terms of space complexity. The entire
540 * integer lattice of noise points could be done as 2 lines instead, and for 3D,
541 * 2 lines + 2 planes.
542 * However, this would require the noise calls to be interposed with the
543 * interpolation loops, which may trash the icache, leading to lower overall
545 * Another optimization that could save half as many noise calls is to carry over
546 * values from the previous noise lattice as midpoints in the new lattice for the
549 #define idx(x, y) ((y) * nlx + (x))
550 void Noise::gradientMap2D(
552 float step_x, float step_y,
555 float v00, v01, v10, v11, u, v, orig_u;
556 u32 index, i, j, noisex, noisey;
560 bool eased = np.flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED);
561 Interp2dFxn interpolate = eased ?
562 biLinearInterpolation : biLinearInterpolationNoEase;
570 //calculate noise point lattice
571 nlx = (u32)(u + sx * step_x) + 2;
572 nly = (u32)(v + sy * step_y) + 2;
574 for (j = 0; j != nly; j++)
575 for (i = 0; i != nlx; i++)
576 noise_buf[index++] = noise2d(x0 + i, y0 + j, seed);
578 //calculate interpolations
581 for (j = 0; j != sy; j++) {
582 v00 = noise_buf[idx(0, noisey)];
583 v10 = noise_buf[idx(1, noisey)];
584 v01 = noise_buf[idx(0, noisey + 1)];
585 v11 = noise_buf[idx(1, noisey + 1)];
589 for (i = 0; i != sx; i++) {
590 gradient_buf[index++] = interpolate(v00, v10, v01, v11, u, v);
598 v10 = noise_buf[idx(noisex + 1, noisey)];
599 v11 = noise_buf[idx(noisex + 1, noisey + 1)];
613 #define idx(x, y, z) ((z) * nly * nlx + (y) * nlx + (x))
614 void Noise::gradientMap3D(
615 float x, float y, float z,
616 float step_x, float step_y, float step_z,
619 float v000, v010, v100, v110;
620 float v001, v011, v101, v111;
621 float u, v, w, orig_u, orig_v;
622 u32 index, i, j, k, noisex, noisey, noisez;
626 Interp3dFxn interpolate = (np.flags & NOISE_FLAG_EASED) ?
627 triLinearInterpolation : triLinearInterpolationNoEase;
638 //calculate noise point lattice
639 nlx = (u32)(u + sx * step_x) + 2;
640 nly = (u32)(v + sy * step_y) + 2;
641 nlz = (u32)(w + sz * step_z) + 2;
643 for (k = 0; k != nlz; k++)
644 for (j = 0; j != nly; j++)
645 for (i = 0; i != nlx; i++)
646 noise_buf[index++] = noise3d(x0 + i, y0 + j, z0 + k, seed);
648 //calculate interpolations
652 for (k = 0; k != sz; k++) {
655 for (j = 0; j != sy; j++) {
656 v000 = noise_buf[idx(0, noisey, noisez)];
657 v100 = noise_buf[idx(1, noisey, noisez)];
658 v010 = noise_buf[idx(0, noisey + 1, noisez)];
659 v110 = noise_buf[idx(1, noisey + 1, noisez)];
660 v001 = noise_buf[idx(0, noisey, noisez + 1)];
661 v101 = noise_buf[idx(1, noisey, noisez + 1)];
662 v011 = noise_buf[idx(0, noisey + 1, noisez + 1)];
663 v111 = noise_buf[idx(1, noisey + 1, noisez + 1)];
667 for (i = 0; i != sx; i++) {
668 gradient_buf[index++] = interpolate(
669 v000, v100, v010, v110,
670 v001, v101, v011, v111,
679 v100 = noise_buf[idx(noisex + 1, noisey, noisez)];
680 v110 = noise_buf[idx(noisex + 1, noisey + 1, noisez)];
683 v101 = noise_buf[idx(noisex + 1, noisey, noisez + 1)];
684 v111 = noise_buf[idx(noisex + 1, noisey + 1, noisez + 1)];
705 float *Noise::perlinMap2D(float x, float y, float *persistence_map)
707 float f = 1.0, g = 1.0;
708 size_t bufsize = sx * sy;
713 memset(result, 0, sizeof(float) * bufsize);
715 if (persistence_map) {
717 persist_buf = new float[bufsize];
718 for (size_t i = 0; i != bufsize; i++)
719 persist_buf[i] = 1.0;
722 for (size_t oct = 0; oct < np.octaves; oct++) {
723 gradientMap2D(x * f, y * f,
724 f / np.spread.X, f / np.spread.Y,
725 seed + np.seed + oct);
727 updateResults(g, persist_buf, persistence_map, bufsize);
733 if (std::fabs(np.offset - 0.f) > 0.00001 || std::fabs(np.scale - 1.f) > 0.00001) {
734 for (size_t i = 0; i != bufsize; i++)
735 result[i] = result[i] * np.scale + np.offset;
742 float *Noise::perlinMap3D(float x, float y, float z, float *persistence_map)
744 float f = 1.0, g = 1.0;
745 size_t bufsize = sx * sy * sz;
751 memset(result, 0, sizeof(float) * bufsize);
753 if (persistence_map) {
755 persist_buf = new float[bufsize];
756 for (size_t i = 0; i != bufsize; i++)
757 persist_buf[i] = 1.0;
760 for (size_t oct = 0; oct < np.octaves; oct++) {
761 gradientMap3D(x * f, y * f, z * f,
762 f / np.spread.X, f / np.spread.Y, f / np.spread.Z,
763 seed + np.seed + oct);
765 updateResults(g, persist_buf, persistence_map, bufsize);
771 if (std::fabs(np.offset - 0.f) > 0.00001 || std::fabs(np.scale - 1.f) > 0.00001) {
772 for (size_t i = 0; i != bufsize; i++)
773 result[i] = result[i] * np.scale + np.offset;
780 void Noise::updateResults(float g, float *gmap,
781 const float *persistence_map, size_t bufsize)
783 // This looks very ugly, but it is 50-70% faster than having
784 // conditional statements inside the loop
785 if (np.flags & NOISE_FLAG_ABSVALUE) {
786 if (persistence_map) {
787 for (size_t i = 0; i != bufsize; i++) {
788 result[i] += gmap[i] * std::fabs(gradient_buf[i]);
789 gmap[i] *= persistence_map[i];
792 for (size_t i = 0; i != bufsize; i++)
793 result[i] += g * std::fabs(gradient_buf[i]);
796 if (persistence_map) {
797 for (size_t i = 0; i != bufsize; i++) {
798 result[i] += gmap[i] * gradient_buf[i];
799 gmap[i] *= persistence_map[i];
802 for (size_t i = 0; i != bufsize; i++)
803 result[i] += g * gradient_buf[i];