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 FlagDesc flagdesc_noiseparams[] = {
50 {"defaults", NOISE_FLAG_DEFAULTS},
51 {"eased", NOISE_FLAG_EASED},
52 {"absvalue", NOISE_FLAG_ABSVALUE},
53 {"pointbuffer", NOISE_FLAG_POINTBUFFER},
54 {"simplex", NOISE_FLAG_SIMPLEX},
58 ///////////////////////////////////////////////////////////////////////////////
60 PcgRandom::PcgRandom(u64 state, u64 seq)
65 void PcgRandom::seed(u64 state, u64 seq)
68 m_inc = (seq << 1u) | 1u;
77 u64 oldstate = m_state;
78 m_state = oldstate * 6364136223846793005ULL + m_inc;
80 u32 xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;
81 u32 rot = oldstate >> 59u;
82 return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));
86 u32 PcgRandom::range(u32 bound)
88 // If the bound is 0, we cover the whole RNG's range
93 This is an optimization of the expression:
94 0x100000000ull % bound
95 since 64-bit modulo operations typically much slower than 32.
97 u32 threshold = -bound % bound;
101 If the bound is not a multiple of the RNG's range, it may cause bias,
102 e.g. a RNG has a range from 0 to 3 and we take want a number 0 to 2.
103 Using rand() % 3, the number 0 would be twice as likely to appear.
104 With a very large RNG range, the effect becomes less prevalent but
107 This can be solved by modifying the range of the RNG to become a
108 multiple of bound by dropping values above the a threshold.
110 In our example, threshold == 4 % 3 == 1, so reject values < 1
111 (that is, 0), thus making the range == 3 with no bias.
113 This loop may look dangerous, but will always terminate due to the
114 RNG's property of uniformity.
116 while ((r = next()) < threshold)
123 s32 PcgRandom::range(s32 min, s32 max)
126 throw PrngException("Invalid range (max < min)");
128 // We have to cast to s64 because otherwise this could overflow,
129 // and signed overflow is undefined behavior.
130 u32 bound = (s64)max - (s64)min + 1;
131 return range(bound) + min;
135 void PcgRandom::bytes(void *out, size_t len)
137 u8 *outb = (u8 *)out;
142 if (bytes_left == 0) {
143 bytes_left = sizeof(u32);
155 s32 PcgRandom::randNormalDist(s32 min, s32 max, int num_trials)
158 for (int i = 0; i != num_trials; i++)
159 accum += range(min, max);
160 return myround((float)accum / num_trials);
163 ///////////////////////////////////////////////////////////////////////////////
165 float noise2d(int x, int y, s32 seed)
167 unsigned int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y
168 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
170 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
171 return 1.f - (float)(int)n / 0x40000000;
175 float noise3d(int x, int y, int z, s32 seed)
177 unsigned int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y + NOISE_MAGIC_Z * z
178 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
180 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
181 return 1.f - (float)(int)n / 0x40000000;
185 inline float dotProduct(float vx, float vy, float wx, float wy)
187 return vx * wx + vy * wy;
191 inline float linearInterpolation(float v0, float v1, float t)
193 return v0 + (v1 - v0) * t;
197 inline float biLinearInterpolation(
198 float v00, float v10,
199 float v01, float v11,
202 float tx = easeCurve(x);
203 float ty = easeCurve(y);
204 float u = linearInterpolation(v00, v10, tx);
205 float v = linearInterpolation(v01, v11, tx);
206 return linearInterpolation(u, v, ty);
210 inline float biLinearInterpolationNoEase(
211 float v00, float v10,
212 float v01, float v11,
215 float u = linearInterpolation(v00, v10, x);
216 float v = linearInterpolation(v01, v11, x);
217 return linearInterpolation(u, v, y);
221 float triLinearInterpolation(
222 float v000, float v100, float v010, float v110,
223 float v001, float v101, float v011, float v111,
224 float x, float y, float z)
226 float tx = easeCurve(x);
227 float ty = easeCurve(y);
228 float tz = easeCurve(z);
229 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, tx, ty);
230 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, tx, ty);
231 return linearInterpolation(u, v, tz);
234 float triLinearInterpolationNoEase(
235 float v000, float v100, float v010, float v110,
236 float v001, float v101, float v011, float v111,
237 float x, float y, float z)
239 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, x, y);
240 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, x, y);
241 return linearInterpolation(u, v, z);
244 float noise2d_gradient(float x, float y, s32 seed, bool eased)
246 // Calculate the integer coordinates
249 // Calculate the remaining part of the coordinates
250 float xl = x - (float)x0;
251 float yl = y - (float)y0;
252 // Get values for corners of square
253 float v00 = noise2d(x0, y0, seed);
254 float v10 = noise2d(x0+1, y0, seed);
255 float v01 = noise2d(x0, y0+1, seed);
256 float v11 = noise2d(x0+1, y0+1, seed);
259 return biLinearInterpolation(v00, v10, v01, v11, xl, yl);
261 return biLinearInterpolationNoEase(v00, v10, v01, v11, xl, yl);
265 float noise3d_gradient(float x, float y, float z, s32 seed, bool eased)
267 // Calculate the integer coordinates
271 // Calculate the remaining part of the coordinates
272 float xl = x - (float)x0;
273 float yl = y - (float)y0;
274 float zl = z - (float)z0;
275 // Get values for corners of cube
276 float v000 = noise3d(x0, y0, z0, seed);
277 float v100 = noise3d(x0 + 1, y0, z0, seed);
278 float v010 = noise3d(x0, y0 + 1, z0, seed);
279 float v110 = noise3d(x0 + 1, y0 + 1, z0, seed);
280 float v001 = noise3d(x0, y0, z0 + 1, seed);
281 float v101 = noise3d(x0 + 1, y0, z0 + 1, seed);
282 float v011 = noise3d(x0, y0 + 1, z0 + 1, seed);
283 float v111 = noise3d(x0 + 1, y0 + 1, z0 + 1, seed);
286 return triLinearInterpolation(
287 v000, v100, v010, v110,
288 v001, v101, v011, v111,
292 return triLinearInterpolationNoEase(
293 v000, v100, v010, v110,
294 v001, v101, v011, v111,
299 float noise2d_perlin(float x, float y, s32 seed,
300 int octaves, float persistence, bool eased)
305 for (int i = 0; i < octaves; i++)
307 a += g * noise2d_gradient(x * f, y * f, seed + i, eased);
315 float noise2d_perlin_abs(float x, float y, s32 seed,
316 int octaves, float persistence, bool eased)
321 for (int i = 0; i < octaves; i++) {
322 a += g * std::fabs(noise2d_gradient(x * f, y * f, seed + i, eased));
330 float noise3d_perlin(float x, float y, float z, s32 seed,
331 int octaves, float persistence, bool eased)
336 for (int i = 0; i < octaves; i++) {
337 a += g * noise3d_gradient(x * f, y * f, z * f, seed + i, eased);
345 float noise3d_perlin_abs(float x, float y, float z, s32 seed,
346 int octaves, float persistence, bool eased)
351 for (int i = 0; i < octaves; i++) {
352 a += g * std::fabs(noise3d_gradient(x * f, y * f, z * f, seed + i, eased));
360 float contour(float v)
369 ///////////////////////// [ New noise ] ////////////////////////////
372 float NoisePerlin2D(NoiseParams *np, float x, float y, s32 seed)
382 for (size_t i = 0; i < np->octaves; i++) {
383 float noiseval = noise2d_gradient(x * f, y * f, seed + i,
384 np->flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED));
386 if (np->flags & NOISE_FLAG_ABSVALUE)
387 noiseval = std::fabs(noiseval);
394 return np->offset + a * np->scale;
398 float NoisePerlin3D(NoiseParams *np, float x, float y, float z, s32 seed)
409 for (size_t i = 0; i < np->octaves; i++) {
410 float noiseval = noise3d_gradient(x * f, y * f, z * f, seed + i,
411 np->flags & NOISE_FLAG_EASED);
413 if (np->flags & NOISE_FLAG_ABSVALUE)
414 noiseval = std::fabs(noiseval);
421 return np->offset + a * np->scale;
425 Noise::Noise(NoiseParams *np_, s32 seed, u32 sx, u32 sy, u32 sz)
427 memcpy(&np, np_, sizeof(np));
439 delete[] gradient_buf;
440 delete[] persist_buf;
446 void Noise::allocBuffers()
455 this->noise_buf = NULL;
456 resizeNoiseBuf(sz > 1);
458 delete[] gradient_buf;
459 delete[] persist_buf;
463 size_t bufsize = sx * sy * sz;
464 this->persist_buf = NULL;
465 this->gradient_buf = new float[bufsize];
466 this->result = new float[bufsize];
467 } catch (std::bad_alloc &e) {
468 throw InvalidNoiseParamsException();
473 void Noise::setSize(u32 sx, u32 sy, u32 sz)
483 void Noise::setSpreadFactor(v3f spread)
485 this->np.spread = spread;
487 resizeNoiseBuf(sz > 1);
491 void Noise::setOctaves(int octaves)
493 this->np.octaves = octaves;
495 resizeNoiseBuf(sz > 1);
499 void Noise::resizeNoiseBuf(bool is3d)
501 // Maximum possible spread value factor
502 float ofactor = (np.lacunarity > 1.0) ?
503 pow(np.lacunarity, np.octaves - 1) :
506 // Noise lattice point count
507 // (int)(sz * spread * ofactor) is # of lattice points crossed due to length
508 float num_noise_points_x = sx * ofactor / np.spread.X;
509 float num_noise_points_y = sy * ofactor / np.spread.Y;
510 float num_noise_points_z = sz * ofactor / np.spread.Z;
512 // Protect against obviously invalid parameters
513 if (num_noise_points_x > 1000000000.f ||
514 num_noise_points_y > 1000000000.f ||
515 num_noise_points_z > 1000000000.f)
516 throw InvalidNoiseParamsException();
518 // Protect against an octave having a spread < 1, causing broken noise values
519 if (np.spread.X / ofactor < 1.0f ||
520 np.spread.Y / ofactor < 1.0f ||
521 np.spread.Z / ofactor < 1.0f) {
522 errorstream << "A noise parameter has too many octaves: "
523 << np.octaves << " octaves" << std::endl;
524 throw InvalidNoiseParamsException("A noise parameter has too many octaves");
527 // + 2 for the two initial endpoints
528 // + 1 for potentially crossing a boundary due to offset
529 size_t nlx = (size_t)std::ceil(num_noise_points_x) + 3;
530 size_t nly = (size_t)std::ceil(num_noise_points_y) + 3;
531 size_t nlz = is3d ? (size_t)std::ceil(num_noise_points_z) + 3 : 1;
535 noise_buf = new float[nlx * nly * nlz];
536 } catch (std::bad_alloc &e) {
537 throw InvalidNoiseParamsException();
543 * NB: This algorithm is not optimal in terms of space complexity. The entire
544 * integer lattice of noise points could be done as 2 lines instead, and for 3D,
545 * 2 lines + 2 planes.
546 * However, this would require the noise calls to be interposed with the
547 * interpolation loops, which may trash the icache, leading to lower overall
549 * Another optimization that could save half as many noise calls is to carry over
550 * values from the previous noise lattice as midpoints in the new lattice for the
553 #define idx(x, y) ((y) * nlx + (x))
554 void Noise::gradientMap2D(
556 float step_x, float step_y,
559 float v00, v01, v10, v11, u, v, orig_u;
560 u32 index, i, j, noisex, noisey;
564 bool eased = np.flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED);
565 Interp2dFxn interpolate = eased ?
566 biLinearInterpolation : biLinearInterpolationNoEase;
574 //calculate noise point lattice
575 nlx = (u32)(u + sx * step_x) + 2;
576 nly = (u32)(v + sy * step_y) + 2;
578 for (j = 0; j != nly; j++)
579 for (i = 0; i != nlx; i++)
580 noise_buf[index++] = noise2d(x0 + i, y0 + j, seed);
582 //calculate interpolations
585 for (j = 0; j != sy; j++) {
586 v00 = noise_buf[idx(0, noisey)];
587 v10 = noise_buf[idx(1, noisey)];
588 v01 = noise_buf[idx(0, noisey + 1)];
589 v11 = noise_buf[idx(1, noisey + 1)];
593 for (i = 0; i != sx; i++) {
594 gradient_buf[index++] = interpolate(v00, v10, v01, v11, u, v);
602 v10 = noise_buf[idx(noisex + 1, noisey)];
603 v11 = noise_buf[idx(noisex + 1, noisey + 1)];
617 #define idx(x, y, z) ((z) * nly * nlx + (y) * nlx + (x))
618 void Noise::gradientMap3D(
619 float x, float y, float z,
620 float step_x, float step_y, float step_z,
623 float v000, v010, v100, v110;
624 float v001, v011, v101, v111;
625 float u, v, w, orig_u, orig_v;
626 u32 index, i, j, k, noisex, noisey, noisez;
630 Interp3dFxn interpolate = (np.flags & NOISE_FLAG_EASED) ?
631 triLinearInterpolation : triLinearInterpolationNoEase;
642 //calculate noise point lattice
643 nlx = (u32)(u + sx * step_x) + 2;
644 nly = (u32)(v + sy * step_y) + 2;
645 nlz = (u32)(w + sz * step_z) + 2;
647 for (k = 0; k != nlz; k++)
648 for (j = 0; j != nly; j++)
649 for (i = 0; i != nlx; i++)
650 noise_buf[index++] = noise3d(x0 + i, y0 + j, z0 + k, seed);
652 //calculate interpolations
656 for (k = 0; k != sz; k++) {
659 for (j = 0; j != sy; j++) {
660 v000 = noise_buf[idx(0, noisey, noisez)];
661 v100 = noise_buf[idx(1, noisey, noisez)];
662 v010 = noise_buf[idx(0, noisey + 1, noisez)];
663 v110 = noise_buf[idx(1, noisey + 1, noisez)];
664 v001 = noise_buf[idx(0, noisey, noisez + 1)];
665 v101 = noise_buf[idx(1, noisey, noisez + 1)];
666 v011 = noise_buf[idx(0, noisey + 1, noisez + 1)];
667 v111 = noise_buf[idx(1, noisey + 1, noisez + 1)];
671 for (i = 0; i != sx; i++) {
672 gradient_buf[index++] = interpolate(
673 v000, v100, v010, v110,
674 v001, v101, v011, v111,
683 v100 = noise_buf[idx(noisex + 1, noisey, noisez)];
684 v110 = noise_buf[idx(noisex + 1, noisey + 1, noisez)];
687 v101 = noise_buf[idx(noisex + 1, noisey, noisez + 1)];
688 v111 = noise_buf[idx(noisex + 1, noisey + 1, noisez + 1)];
709 float *Noise::perlinMap2D(float x, float y, float *persistence_map)
711 float f = 1.0, g = 1.0;
712 size_t bufsize = sx * sy;
717 memset(result, 0, sizeof(float) * bufsize);
719 if (persistence_map) {
721 persist_buf = new float[bufsize];
722 for (size_t i = 0; i != bufsize; i++)
723 persist_buf[i] = 1.0;
726 for (size_t oct = 0; oct < np.octaves; oct++) {
727 gradientMap2D(x * f, y * f,
728 f / np.spread.X, f / np.spread.Y,
729 seed + np.seed + oct);
731 updateResults(g, persist_buf, persistence_map, bufsize);
737 if (std::fabs(np.offset - 0.f) > 0.00001 || std::fabs(np.scale - 1.f) > 0.00001) {
738 for (size_t i = 0; i != bufsize; i++)
739 result[i] = result[i] * np.scale + np.offset;
746 float *Noise::perlinMap3D(float x, float y, float z, float *persistence_map)
748 float f = 1.0, g = 1.0;
749 size_t bufsize = sx * sy * sz;
755 memset(result, 0, sizeof(float) * bufsize);
757 if (persistence_map) {
759 persist_buf = new float[bufsize];
760 for (size_t i = 0; i != bufsize; i++)
761 persist_buf[i] = 1.0;
764 for (size_t oct = 0; oct < np.octaves; oct++) {
765 gradientMap3D(x * f, y * f, z * f,
766 f / np.spread.X, f / np.spread.Y, f / np.spread.Z,
767 seed + np.seed + oct);
769 updateResults(g, persist_buf, persistence_map, bufsize);
775 if (std::fabs(np.offset - 0.f) > 0.00001 || std::fabs(np.scale - 1.f) > 0.00001) {
776 for (size_t i = 0; i != bufsize; i++)
777 result[i] = result[i] * np.scale + np.offset;
784 void Noise::updateResults(float g, float *gmap,
785 const float *persistence_map, size_t bufsize)
787 // This looks very ugly, but it is 50-70% faster than having
788 // conditional statements inside the loop
789 if (np.flags & NOISE_FLAG_ABSVALUE) {
790 if (persistence_map) {
791 for (size_t i = 0; i != bufsize; i++) {
792 result[i] += gmap[i] * std::fabs(gradient_buf[i]);
793 gmap[i] *= persistence_map[i];
796 for (size_t i = 0; i != bufsize; i++)
797 result[i] += g * std::fabs(gradient_buf[i]);
800 if (persistence_map) {
801 for (size_t i = 0; i != bufsize; i++) {
802 result[i] += gmap[i] * gradient_buf[i];
803 gmap[i] *= persistence_map[i];
806 for (size_t i = 0; i != bufsize; i++)
807 result[i] += g * gradient_buf[i];