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 <string.h> // 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.0, 0.9238, 0.7071, 0.3826, 0, -0.3826, -0.7071, -0.9238,
51 1.0, -0.9238, -0.7071, -0.3826, 0, 0.3826, 0.7071, 0.9238
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
97 If the bound is not a multiple of the RNG's range, it may cause bias,
98 e.g. a RNG has a range from 0 to 3 and we take want a number 0 to 2.
99 Using rand() % 3, the number 0 would be twice as likely to appear.
100 With a very large RNG range, the effect becomes less prevalent but
101 still present. This can be solved by modifying the range of the RNG
102 to become a multiple of bound by dropping values above the a threshhold.
103 In our example, threshhold == 4 - 3 = 1 % 3 == 1, so reject 0, thus
104 making the range 3 with no bias.
106 This loop looks dangerous, but will always terminate due to the
107 RNG's property of uniformity.
109 u32 threshhold = -bound % bound;
112 while ((r = next()) < threshhold)
119 s32 PcgRandom::range(s32 min, s32 max)
122 throw PrngException("Invalid range (max < min)");
124 u32 bound = max - min + 1;
125 return range(bound) + min;
129 void PcgRandom::bytes(void *out, size_t len)
131 u8 *outb = (u8 *)out;
136 if (bytes_left == 0) {
137 bytes_left = sizeof(u32);
149 s32 PcgRandom::randNormalDist(s32 min, s32 max, int num_trials)
152 for (int i = 0; i != num_trials; i++)
153 accum += range(min, max);
154 return myround((float)accum / num_trials);
157 ///////////////////////////////////////////////////////////////////////////////
159 float noise2d(int x, int y, int seed)
161 int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y
162 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
164 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
165 return 1.f - (float)n / 0x40000000;
169 float noise3d(int x, int y, int z, int seed)
171 int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y + NOISE_MAGIC_Z * z
172 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
174 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
175 return 1.f - (float)n / 0x40000000;
179 inline float dotProduct(float vx, float vy, float wx, float wy)
181 return vx * wx + vy * wy;
185 inline float linearInterpolation(float v0, float v1, float t)
187 return v0 + (v1 - v0) * t;
191 inline float biLinearInterpolation(
192 float v00, float v10,
193 float v01, float v11,
196 float tx = easeCurve(x);
197 float ty = easeCurve(y);
198 float u = linearInterpolation(v00, v10, tx);
199 float v = linearInterpolation(v01, v11, tx);
200 return linearInterpolation(u, v, ty);
204 inline float biLinearInterpolationNoEase(
205 float v00, float v10,
206 float v01, float v11,
209 float u = linearInterpolation(v00, v10, x);
210 float v = linearInterpolation(v01, v11, x);
211 return linearInterpolation(u, v, y);
215 float triLinearInterpolation(
216 float v000, float v100, float v010, float v110,
217 float v001, float v101, float v011, float v111,
218 float x, float y, float z)
220 float tx = easeCurve(x);
221 float ty = easeCurve(y);
222 float tz = easeCurve(z);
223 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, tx, ty);
224 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, tx, ty);
225 return linearInterpolation(u, v, tz);
228 float triLinearInterpolationNoEase(
229 float v000, float v100, float v010, float v110,
230 float v001, float v101, float v011, float v111,
231 float x, float y, float z)
233 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, x, y);
234 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, x, y);
235 return linearInterpolation(u, v, z);
238 float noise2d_gradient(float x, float y, int seed, bool eased)
240 // Calculate the integer coordinates
243 // Calculate the remaining part of the coordinates
244 float xl = x - (float)x0;
245 float yl = y - (float)y0;
246 // Get values for corners of square
247 float v00 = noise2d(x0, y0, seed);
248 float v10 = noise2d(x0+1, y0, seed);
249 float v01 = noise2d(x0, y0+1, seed);
250 float v11 = noise2d(x0+1, y0+1, seed);
253 return biLinearInterpolation(v00, v10, v01, v11, xl, yl);
255 return biLinearInterpolationNoEase(v00, v10, v01, v11, xl, yl);
259 float noise3d_gradient(float x, float y, float z, int seed, bool eased)
261 // Calculate the integer coordinates
265 // Calculate the remaining part of the coordinates
266 float xl = x - (float)x0;
267 float yl = y - (float)y0;
268 float zl = z - (float)z0;
269 // Get values for corners of cube
270 float v000 = noise3d(x0, y0, z0, seed);
271 float v100 = noise3d(x0 + 1, y0, z0, seed);
272 float v010 = noise3d(x0, y0 + 1, z0, seed);
273 float v110 = noise3d(x0 + 1, y0 + 1, z0, seed);
274 float v001 = noise3d(x0, y0, z0 + 1, seed);
275 float v101 = noise3d(x0 + 1, y0, z0 + 1, seed);
276 float v011 = noise3d(x0, y0 + 1, z0 + 1, seed);
277 float v111 = noise3d(x0 + 1, y0 + 1, z0 + 1, seed);
280 return triLinearInterpolation(
281 v000, v100, v010, v110,
282 v001, v101, v011, v111,
285 return triLinearInterpolationNoEase(
286 v000, v100, v010, v110,
287 v001, v101, v011, v111,
293 float noise2d_perlin(float x, float y, int seed,
294 int octaves, float persistence, bool eased)
299 for (int i = 0; i < octaves; i++)
301 a += g * noise2d_gradient(x * f, y * f, seed + i, eased);
309 float noise2d_perlin_abs(float x, float y, int seed,
310 int octaves, float persistence, bool eased)
315 for (int i = 0; i < octaves; i++) {
316 a += g * fabs(noise2d_gradient(x * f, y * f, seed + i, eased));
324 float noise3d_perlin(float x, float y, float z, int seed,
325 int octaves, float persistence, bool eased)
330 for (int i = 0; i < octaves; i++) {
331 a += g * noise3d_gradient(x * f, y * f, z * f, seed + i, eased);
339 float noise3d_perlin_abs(float x, float y, float z, int seed,
340 int octaves, float persistence, bool eased)
345 for (int i = 0; i < octaves; i++) {
346 a += g * fabs(noise3d_gradient(x * f, y * f, z * f, seed + i, eased));
354 float contour(float v)
363 ///////////////////////// [ New noise ] ////////////////////////////
366 float NoisePerlin2D(NoiseParams *np, float x, float y, int seed)
376 for (size_t i = 0; i < np->octaves; i++) {
377 float noiseval = noise2d_gradient(x * f, y * f, seed + i,
378 np->flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED));
380 if (np->flags & NOISE_FLAG_ABSVALUE)
381 noiseval = fabs(noiseval);
388 return np->offset + a * np->scale;
392 float NoisePerlin3D(NoiseParams *np, float x, float y, float z, int seed)
403 for (size_t i = 0; i < np->octaves; i++) {
404 float noiseval = noise3d_gradient(x * f, y * f, z * f, seed + i,
405 np->flags & NOISE_FLAG_EASED);
407 if (np->flags & NOISE_FLAG_ABSVALUE)
408 noiseval = fabs(noiseval);
415 return np->offset + a * np->scale;
419 Noise::Noise(NoiseParams *np_, int seed, u32 sx, u32 sy, u32 sz)
421 memcpy(&np, np_, sizeof(np));
427 this->persist_buf = NULL;
428 this->gradient_buf = NULL;
437 delete[] gradient_buf;
438 delete[] persist_buf;
444 void Noise::allocBuffers()
453 this->noise_buf = NULL;
454 resizeNoiseBuf(sz > 1);
456 delete[] gradient_buf;
457 delete[] persist_buf;
461 size_t bufsize = sx * sy * sz;
462 this->persist_buf = NULL;
463 this->gradient_buf = new float[bufsize];
464 this->result = new float[bufsize];
465 } catch (std::bad_alloc &e) {
466 throw InvalidNoiseParamsException();
471 void Noise::setSize(u32 sx, u32 sy, u32 sz)
481 void Noise::setSpreadFactor(v3f spread)
483 this->np.spread = spread;
485 resizeNoiseBuf(sz > 1);
489 void Noise::setOctaves(int octaves)
491 this->np.octaves = octaves;
493 resizeNoiseBuf(sz > 1);
497 void Noise::resizeNoiseBuf(bool is3d)
499 //maximum possible spread value factor
500 float ofactor = (np.lacunarity > 1.0) ?
501 pow(np.lacunarity, np.octaves - 1) :
504 // noise lattice point count
505 // (int)(sz * spread * ofactor) is # of lattice points crossed due to length
506 float num_noise_points_x = sx * ofactor / np.spread.X;
507 float num_noise_points_y = sy * ofactor / np.spread.Y;
508 float num_noise_points_z = sz * ofactor / np.spread.Z;
510 // protect against obviously invalid parameters
511 if (num_noise_points_x > 1000000000.f ||
512 num_noise_points_y > 1000000000.f ||
513 num_noise_points_z > 1000000000.f)
514 throw InvalidNoiseParamsException();
516 // + 2 for the two initial endpoints
517 // + 1 for potentially crossing a boundary due to offset
518 size_t nlx = (size_t)ceil(num_noise_points_x) + 3;
519 size_t nly = (size_t)ceil(num_noise_points_y) + 3;
520 size_t nlz = is3d ? (size_t)ceil(num_noise_points_z) + 3 : 1;
524 noise_buf = new float[nlx * nly * nlz];
525 } catch (std::bad_alloc &e) {
526 throw InvalidNoiseParamsException();
532 * NB: This algorithm is not optimal in terms of space complexity. The entire
533 * integer lattice of noise points could be done as 2 lines instead, and for 3D,
534 * 2 lines + 2 planes.
535 * However, this would require the noise calls to be interposed with the
536 * interpolation loops, which may trash the icache, leading to lower overall
538 * Another optimization that could save half as many noise calls is to carry over
539 * values from the previous noise lattice as midpoints in the new lattice for the
542 #define idx(x, y) ((y) * nlx + (x))
543 void Noise::gradientMap2D(
545 float step_x, float step_y,
548 float v00, v01, v10, v11, u, v, orig_u;
549 u32 index, i, j, noisex, noisey;
553 bool eased = np.flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED);
554 Interp2dFxn interpolate = eased ?
555 biLinearInterpolation : biLinearInterpolationNoEase;
563 //calculate noise point lattice
564 nlx = (u32)(u + sx * step_x) + 2;
565 nly = (u32)(v + sy * step_y) + 2;
567 for (j = 0; j != nly; j++)
568 for (i = 0; i != nlx; i++)
569 noise_buf[index++] = noise2d(x0 + i, y0 + j, seed);
571 //calculate interpolations
574 for (j = 0; j != sy; j++) {
575 v00 = noise_buf[idx(0, noisey)];
576 v10 = noise_buf[idx(1, noisey)];
577 v01 = noise_buf[idx(0, noisey + 1)];
578 v11 = noise_buf[idx(1, noisey + 1)];
582 for (i = 0; i != sx; i++) {
583 gradient_buf[index++] = interpolate(v00, v10, v01, v11, u, v);
591 v10 = noise_buf[idx(noisex + 1, noisey)];
592 v11 = noise_buf[idx(noisex + 1, noisey + 1)];
606 #define idx(x, y, z) ((z) * nly * nlx + (y) * nlx + (x))
607 void Noise::gradientMap3D(
608 float x, float y, float z,
609 float step_x, float step_y, float step_z,
612 float v000, v010, v100, v110;
613 float v001, v011, v101, v111;
614 float u, v, w, orig_u, orig_v;
615 u32 index, i, j, k, noisex, noisey, noisez;
619 Interp3dFxn interpolate = (np.flags & NOISE_FLAG_EASED) ?
620 triLinearInterpolation : triLinearInterpolationNoEase;
631 //calculate noise point lattice
632 nlx = (u32)(u + sx * step_x) + 2;
633 nly = (u32)(v + sy * step_y) + 2;
634 nlz = (u32)(w + sz * step_z) + 2;
636 for (k = 0; k != nlz; k++)
637 for (j = 0; j != nly; j++)
638 for (i = 0; i != nlx; i++)
639 noise_buf[index++] = noise3d(x0 + i, y0 + j, z0 + k, seed);
641 //calculate interpolations
645 for (k = 0; k != sz; k++) {
648 for (j = 0; j != sy; j++) {
649 v000 = noise_buf[idx(0, noisey, noisez)];
650 v100 = noise_buf[idx(1, noisey, noisez)];
651 v010 = noise_buf[idx(0, noisey + 1, noisez)];
652 v110 = noise_buf[idx(1, noisey + 1, noisez)];
653 v001 = noise_buf[idx(0, noisey, noisez + 1)];
654 v101 = noise_buf[idx(1, noisey, noisez + 1)];
655 v011 = noise_buf[idx(0, noisey + 1, noisez + 1)];
656 v111 = noise_buf[idx(1, noisey + 1, noisez + 1)];
660 for (i = 0; i != sx; i++) {
661 gradient_buf[index++] = interpolate(
662 v000, v100, v010, v110,
663 v001, v101, v011, v111,
672 v100 = noise_buf[idx(noisex + 1, noisey, noisez)];
673 v110 = noise_buf[idx(noisex + 1, noisey + 1, noisez)];
676 v101 = noise_buf[idx(noisex + 1, noisey, noisez + 1)];
677 v111 = noise_buf[idx(noisex + 1, noisey + 1, noisez + 1)];
698 float *Noise::perlinMap2D(float x, float y, float *persistence_map)
700 float f = 1.0, g = 1.0;
701 size_t bufsize = sx * sy;
706 memset(result, 0, sizeof(float) * bufsize);
708 if (persistence_map) {
710 persist_buf = new float[bufsize];
711 for (size_t i = 0; i != bufsize; i++)
712 persist_buf[i] = 1.0;
715 for (size_t oct = 0; oct < np.octaves; oct++) {
716 gradientMap2D(x * f, y * f,
717 f / np.spread.X, f / np.spread.Y,
718 seed + np.seed + oct);
720 updateResults(g, persist_buf, persistence_map, bufsize);
726 if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
727 for (size_t i = 0; i != bufsize; i++)
728 result[i] = result[i] * np.scale + np.offset;
735 float *Noise::perlinMap3D(float x, float y, float z, float *persistence_map)
737 float f = 1.0, g = 1.0;
738 size_t bufsize = sx * sy * sz;
744 memset(result, 0, sizeof(float) * bufsize);
746 if (persistence_map) {
748 persist_buf = new float[bufsize];
749 for (size_t i = 0; i != bufsize; i++)
750 persist_buf[i] = 1.0;
753 for (size_t oct = 0; oct < np.octaves; oct++) {
754 gradientMap3D(x * f, y * f, z * f,
755 f / np.spread.X, f / np.spread.Y, f / np.spread.Z,
756 seed + np.seed + oct);
758 updateResults(g, persist_buf, persistence_map, bufsize);
764 if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
765 for (size_t i = 0; i != bufsize; i++)
766 result[i] = result[i] * np.scale + np.offset;
773 void Noise::updateResults(float g, float *gmap,
774 float *persistence_map, size_t bufsize)
776 // This looks very ugly, but it is 50-70% faster than having
777 // conditional statements inside the loop
778 if (np.flags & NOISE_FLAG_ABSVALUE) {
779 if (persistence_map) {
780 for (size_t i = 0; i != bufsize; i++) {
781 result[i] += gmap[i] * fabs(gradient_buf[i]);
782 gmap[i] *= persistence_map[i];
785 for (size_t i = 0; i != bufsize; i++)
786 result[i] += g * fabs(gradient_buf[i]);
789 if (persistence_map) {
790 for (size_t i = 0; i != bufsize; i++) {
791 result[i] += gmap[i] * gradient_buf[i];
792 gmap[i] *= persistence_map[i];
795 for (size_t i = 0; i != bufsize; i++)
796 result[i] += g * gradient_buf[i];