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 contour(float v)
324 ///////////////////////// [ New noise ] ////////////////////////////
327 float NoisePerlin2D(const NoiseParams *np, float x, float y, s32 seed)
337 for (size_t i = 0; i < np->octaves; i++) {
338 float noiseval = noise2d_gradient(x * f, y * f, seed + i,
339 np->flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED));
341 if (np->flags & NOISE_FLAG_ABSVALUE)
342 noiseval = std::fabs(noiseval);
349 return np->offset + a * np->scale;
353 float NoisePerlin3D(const NoiseParams *np, float x, float y, float z, s32 seed)
364 for (size_t i = 0; i < np->octaves; i++) {
365 float noiseval = noise3d_gradient(x * f, y * f, z * f, seed + i,
366 np->flags & NOISE_FLAG_EASED);
368 if (np->flags & NOISE_FLAG_ABSVALUE)
369 noiseval = std::fabs(noiseval);
376 return np->offset + a * np->scale;
380 Noise::Noise(const NoiseParams *np_, s32 seed, u32 sx, u32 sy, u32 sz)
394 delete[] gradient_buf;
395 delete[] persist_buf;
401 void Noise::allocBuffers()
410 this->noise_buf = NULL;
411 resizeNoiseBuf(sz > 1);
413 delete[] gradient_buf;
414 delete[] persist_buf;
418 size_t bufsize = sx * sy * sz;
419 this->persist_buf = NULL;
420 this->gradient_buf = new float[bufsize];
421 this->result = new float[bufsize];
422 } catch (std::bad_alloc &e) {
423 throw InvalidNoiseParamsException();
428 void Noise::setSize(u32 sx, u32 sy, u32 sz)
438 void Noise::setSpreadFactor(v3f spread)
440 this->np.spread = spread;
442 resizeNoiseBuf(sz > 1);
446 void Noise::setOctaves(int octaves)
448 this->np.octaves = octaves;
450 resizeNoiseBuf(sz > 1);
454 void Noise::resizeNoiseBuf(bool is3d)
456 // Maximum possible spread value factor
457 float ofactor = (np.lacunarity > 1.0) ?
458 pow(np.lacunarity, np.octaves - 1) :
461 // Noise lattice point count
462 // (int)(sz * spread * ofactor) is # of lattice points crossed due to length
463 float num_noise_points_x = sx * ofactor / np.spread.X;
464 float num_noise_points_y = sy * ofactor / np.spread.Y;
465 float num_noise_points_z = sz * ofactor / np.spread.Z;
467 // Protect against obviously invalid parameters
468 if (num_noise_points_x > 1000000000.f ||
469 num_noise_points_y > 1000000000.f ||
470 num_noise_points_z > 1000000000.f)
471 throw InvalidNoiseParamsException();
473 // Protect against an octave having a spread < 1, causing broken noise values
474 if (np.spread.X / ofactor < 1.0f ||
475 np.spread.Y / ofactor < 1.0f ||
476 np.spread.Z / ofactor < 1.0f) {
477 errorstream << "A noise parameter has too many octaves: "
478 << np.octaves << " octaves" << std::endl;
479 throw InvalidNoiseParamsException("A noise parameter has too many octaves");
482 // + 2 for the two initial endpoints
483 // + 1 for potentially crossing a boundary due to offset
484 size_t nlx = (size_t)std::ceil(num_noise_points_x) + 3;
485 size_t nly = (size_t)std::ceil(num_noise_points_y) + 3;
486 size_t nlz = is3d ? (size_t)std::ceil(num_noise_points_z) + 3 : 1;
490 noise_buf = new float[nlx * nly * nlz];
491 } catch (std::bad_alloc &e) {
492 throw InvalidNoiseParamsException();
498 * NB: This algorithm is not optimal in terms of space complexity. The entire
499 * integer lattice of noise points could be done as 2 lines instead, and for 3D,
500 * 2 lines + 2 planes.
501 * However, this would require the noise calls to be interposed with the
502 * interpolation loops, which may trash the icache, leading to lower overall
504 * Another optimization that could save half as many noise calls is to carry over
505 * values from the previous noise lattice as midpoints in the new lattice for the
508 #define idx(x, y) ((y) * nlx + (x))
509 void Noise::gradientMap2D(
511 float step_x, float step_y,
514 float v00, v01, v10, v11, u, v, orig_u;
515 u32 index, i, j, noisex, noisey;
519 bool eased = np.flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED);
520 Interp2dFxn interpolate = eased ?
521 biLinearInterpolation : biLinearInterpolationNoEase;
529 //calculate noise point lattice
530 nlx = (u32)(u + sx * step_x) + 2;
531 nly = (u32)(v + sy * step_y) + 2;
533 for (j = 0; j != nly; j++)
534 for (i = 0; i != nlx; i++)
535 noise_buf[index++] = noise2d(x0 + i, y0 + j, seed);
537 //calculate interpolations
540 for (j = 0; j != sy; j++) {
541 v00 = noise_buf[idx(0, noisey)];
542 v10 = noise_buf[idx(1, noisey)];
543 v01 = noise_buf[idx(0, noisey + 1)];
544 v11 = noise_buf[idx(1, noisey + 1)];
548 for (i = 0; i != sx; i++) {
549 gradient_buf[index++] = interpolate(v00, v10, v01, v11, u, v);
557 v10 = noise_buf[idx(noisex + 1, noisey)];
558 v11 = noise_buf[idx(noisex + 1, noisey + 1)];
572 #define idx(x, y, z) ((z) * nly * nlx + (y) * nlx + (x))
573 void Noise::gradientMap3D(
574 float x, float y, float z,
575 float step_x, float step_y, float step_z,
578 float v000, v010, v100, v110;
579 float v001, v011, v101, v111;
580 float u, v, w, orig_u, orig_v;
581 u32 index, i, j, k, noisex, noisey, noisez;
585 Interp3dFxn interpolate = (np.flags & NOISE_FLAG_EASED) ?
586 triLinearInterpolation : triLinearInterpolationNoEase;
597 //calculate noise point lattice
598 nlx = (u32)(u + sx * step_x) + 2;
599 nly = (u32)(v + sy * step_y) + 2;
600 nlz = (u32)(w + sz * step_z) + 2;
602 for (k = 0; k != nlz; k++)
603 for (j = 0; j != nly; j++)
604 for (i = 0; i != nlx; i++)
605 noise_buf[index++] = noise3d(x0 + i, y0 + j, z0 + k, seed);
607 //calculate interpolations
611 for (k = 0; k != sz; k++) {
614 for (j = 0; j != sy; j++) {
615 v000 = noise_buf[idx(0, noisey, noisez)];
616 v100 = noise_buf[idx(1, noisey, noisez)];
617 v010 = noise_buf[idx(0, noisey + 1, noisez)];
618 v110 = noise_buf[idx(1, noisey + 1, noisez)];
619 v001 = noise_buf[idx(0, noisey, noisez + 1)];
620 v101 = noise_buf[idx(1, noisey, noisez + 1)];
621 v011 = noise_buf[idx(0, noisey + 1, noisez + 1)];
622 v111 = noise_buf[idx(1, noisey + 1, noisez + 1)];
626 for (i = 0; i != sx; i++) {
627 gradient_buf[index++] = interpolate(
628 v000, v100, v010, v110,
629 v001, v101, v011, v111,
638 v100 = noise_buf[idx(noisex + 1, noisey, noisez)];
639 v110 = noise_buf[idx(noisex + 1, noisey + 1, noisez)];
642 v101 = noise_buf[idx(noisex + 1, noisey, noisez + 1)];
643 v111 = noise_buf[idx(noisex + 1, noisey + 1, noisez + 1)];
664 float *Noise::perlinMap2D(float x, float y, float *persistence_map)
666 float f = 1.0, g = 1.0;
667 size_t bufsize = sx * sy;
672 memset(result, 0, sizeof(float) * bufsize);
674 if (persistence_map) {
676 persist_buf = new float[bufsize];
677 for (size_t i = 0; i != bufsize; i++)
678 persist_buf[i] = 1.0;
681 for (size_t oct = 0; oct < np.octaves; oct++) {
682 gradientMap2D(x * f, y * f,
683 f / np.spread.X, f / np.spread.Y,
684 seed + np.seed + oct);
686 updateResults(g, persist_buf, persistence_map, bufsize);
692 if (std::fabs(np.offset - 0.f) > 0.00001 || std::fabs(np.scale - 1.f) > 0.00001) {
693 for (size_t i = 0; i != bufsize; i++)
694 result[i] = result[i] * np.scale + np.offset;
701 float *Noise::perlinMap3D(float x, float y, float z, float *persistence_map)
703 float f = 1.0, g = 1.0;
704 size_t bufsize = sx * sy * sz;
710 memset(result, 0, sizeof(float) * bufsize);
712 if (persistence_map) {
714 persist_buf = new float[bufsize];
715 for (size_t i = 0; i != bufsize; i++)
716 persist_buf[i] = 1.0;
719 for (size_t oct = 0; oct < np.octaves; oct++) {
720 gradientMap3D(x * f, y * f, z * f,
721 f / np.spread.X, f / np.spread.Y, f / np.spread.Z,
722 seed + np.seed + oct);
724 updateResults(g, persist_buf, persistence_map, bufsize);
730 if (std::fabs(np.offset - 0.f) > 0.00001 || std::fabs(np.scale - 1.f) > 0.00001) {
731 for (size_t i = 0; i != bufsize; i++)
732 result[i] = result[i] * np.scale + np.offset;
739 void Noise::updateResults(float g, float *gmap,
740 const float *persistence_map, size_t bufsize)
742 // This looks very ugly, but it is 50-70% faster than having
743 // conditional statements inside the loop
744 if (np.flags & NOISE_FLAG_ABSVALUE) {
745 if (persistence_map) {
746 for (size_t i = 0; i != bufsize; i++) {
747 result[i] += gmap[i] * std::fabs(gradient_buf[i]);
748 gmap[i] *= persistence_map[i];
751 for (size_t i = 0; i != bufsize; i++)
752 result[i] += g * std::fabs(gradient_buf[i]);
755 if (persistence_map) {
756 for (size_t i = 0; i != bufsize; i++) {
757 result[i] += gmap[i] * gradient_buf[i];
758 gmap[i] *= persistence_map[i];
761 for (size_t i = 0; i != bufsize; i++)
762 result[i] += g * gradient_buf[i];