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
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 u32 bound = max - min + 1;
134 return range(bound) + min;
138 void PcgRandom::bytes(void *out, size_t len)
140 u8 *outb = (u8 *)out;
145 if (bytes_left == 0) {
146 bytes_left = sizeof(u32);
158 s32 PcgRandom::randNormalDist(s32 min, s32 max, int num_trials)
161 for (int i = 0; i != num_trials; i++)
162 accum += range(min, max);
163 return myround((float)accum / num_trials);
166 ///////////////////////////////////////////////////////////////////////////////
168 float noise2d(int x, int y, s32 seed)
170 unsigned int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y
171 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
173 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
174 return 1.f - (float)(int)n / 0x40000000;
178 float noise3d(int x, int y, int z, s32 seed)
180 unsigned int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y + NOISE_MAGIC_Z * z
181 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
183 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
184 return 1.f - (float)(int)n / 0x40000000;
188 inline float dotProduct(float vx, float vy, float wx, float wy)
190 return vx * wx + vy * wy;
194 inline float linearInterpolation(float v0, float v1, float t)
196 return v0 + (v1 - v0) * t;
200 inline float biLinearInterpolation(
201 float v00, float v10,
202 float v01, float v11,
205 float tx = easeCurve(x);
206 float ty = easeCurve(y);
207 float u = linearInterpolation(v00, v10, tx);
208 float v = linearInterpolation(v01, v11, tx);
209 return linearInterpolation(u, v, ty);
213 inline float biLinearInterpolationNoEase(
214 float v00, float v10,
215 float v01, float v11,
218 float u = linearInterpolation(v00, v10, x);
219 float v = linearInterpolation(v01, v11, x);
220 return linearInterpolation(u, v, y);
224 float triLinearInterpolation(
225 float v000, float v100, float v010, float v110,
226 float v001, float v101, float v011, float v111,
227 float x, float y, float z)
229 float tx = easeCurve(x);
230 float ty = easeCurve(y);
231 float tz = easeCurve(z);
232 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, tx, ty);
233 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, tx, ty);
234 return linearInterpolation(u, v, tz);
237 float triLinearInterpolationNoEase(
238 float v000, float v100, float v010, float v110,
239 float v001, float v101, float v011, float v111,
240 float x, float y, float z)
242 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, x, y);
243 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, x, y);
244 return linearInterpolation(u, v, z);
247 float noise2d_gradient(float x, float y, s32 seed, bool eased)
249 // Calculate the integer coordinates
252 // Calculate the remaining part of the coordinates
253 float xl = x - (float)x0;
254 float yl = y - (float)y0;
255 // Get values for corners of square
256 float v00 = noise2d(x0, y0, seed);
257 float v10 = noise2d(x0+1, y0, seed);
258 float v01 = noise2d(x0, y0+1, seed);
259 float v11 = noise2d(x0+1, y0+1, seed);
262 return biLinearInterpolation(v00, v10, v01, v11, xl, yl);
264 return biLinearInterpolationNoEase(v00, v10, v01, v11, xl, yl);
268 float noise3d_gradient(float x, float y, float z, s32 seed, bool eased)
270 // Calculate the integer coordinates
274 // Calculate the remaining part of the coordinates
275 float xl = x - (float)x0;
276 float yl = y - (float)y0;
277 float zl = z - (float)z0;
278 // Get values for corners of cube
279 float v000 = noise3d(x0, y0, z0, seed);
280 float v100 = noise3d(x0 + 1, y0, z0, seed);
281 float v010 = noise3d(x0, y0 + 1, z0, seed);
282 float v110 = noise3d(x0 + 1, y0 + 1, z0, seed);
283 float v001 = noise3d(x0, y0, z0 + 1, seed);
284 float v101 = noise3d(x0 + 1, y0, z0 + 1, seed);
285 float v011 = noise3d(x0, y0 + 1, z0 + 1, seed);
286 float v111 = noise3d(x0 + 1, y0 + 1, z0 + 1, seed);
289 return triLinearInterpolation(
290 v000, v100, v010, v110,
291 v001, v101, v011, v111,
294 return triLinearInterpolationNoEase(
295 v000, v100, v010, v110,
296 v001, v101, v011, v111,
302 float noise2d_perlin(float x, float y, s32 seed,
303 int octaves, float persistence, bool eased)
308 for (int i = 0; i < octaves; i++)
310 a += g * noise2d_gradient(x * f, y * f, seed + i, eased);
318 float noise2d_perlin_abs(float x, float y, s32 seed,
319 int octaves, float persistence, bool eased)
324 for (int i = 0; i < octaves; i++) {
325 a += g * fabs(noise2d_gradient(x * f, y * f, seed + i, eased));
333 float noise3d_perlin(float x, float y, float z, s32 seed,
334 int octaves, float persistence, bool eased)
339 for (int i = 0; i < octaves; i++) {
340 a += g * noise3d_gradient(x * f, y * f, z * f, seed + i, eased);
348 float noise3d_perlin_abs(float x, float y, float z, s32 seed,
349 int octaves, float persistence, bool eased)
354 for (int i = 0; i < octaves; i++) {
355 a += g * fabs(noise3d_gradient(x * f, y * f, z * f, seed + i, eased));
363 float contour(float v)
372 ///////////////////////// [ New noise ] ////////////////////////////
375 float NoisePerlin2D(NoiseParams *np, float x, float y, s32 seed)
385 for (size_t i = 0; i < np->octaves; i++) {
386 float noiseval = noise2d_gradient(x * f, y * f, seed + i,
387 np->flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED));
389 if (np->flags & NOISE_FLAG_ABSVALUE)
390 noiseval = fabs(noiseval);
397 return np->offset + a * np->scale;
401 float NoisePerlin3D(NoiseParams *np, float x, float y, float z, s32 seed)
412 for (size_t i = 0; i < np->octaves; i++) {
413 float noiseval = noise3d_gradient(x * f, y * f, z * f, seed + i,
414 np->flags & NOISE_FLAG_EASED);
416 if (np->flags & NOISE_FLAG_ABSVALUE)
417 noiseval = fabs(noiseval);
424 return np->offset + a * np->scale;
428 Noise::Noise(NoiseParams *np_, s32 seed, u32 sx, u32 sy, u32 sz)
430 memcpy(&np, np_, sizeof(np));
436 this->persist_buf = NULL;
437 this->gradient_buf = NULL;
446 delete[] gradient_buf;
447 delete[] persist_buf;
453 void Noise::allocBuffers()
462 this->noise_buf = NULL;
463 resizeNoiseBuf(sz > 1);
465 delete[] gradient_buf;
466 delete[] persist_buf;
470 size_t bufsize = sx * sy * sz;
471 this->persist_buf = NULL;
472 this->gradient_buf = new float[bufsize];
473 this->result = new float[bufsize];
474 } catch (std::bad_alloc &e) {
475 throw InvalidNoiseParamsException();
480 void Noise::setSize(u32 sx, u32 sy, u32 sz)
490 void Noise::setSpreadFactor(v3f spread)
492 this->np.spread = spread;
494 resizeNoiseBuf(sz > 1);
498 void Noise::setOctaves(int octaves)
500 this->np.octaves = octaves;
502 resizeNoiseBuf(sz > 1);
506 void Noise::resizeNoiseBuf(bool is3d)
508 //maximum possible spread value factor
509 float ofactor = (np.lacunarity > 1.0) ?
510 pow(np.lacunarity, np.octaves - 1) :
513 // noise lattice point count
514 // (int)(sz * spread * ofactor) is # of lattice points crossed due to length
515 float num_noise_points_x = sx * ofactor / np.spread.X;
516 float num_noise_points_y = sy * ofactor / np.spread.Y;
517 float num_noise_points_z = sz * ofactor / np.spread.Z;
519 // protect against obviously invalid parameters
520 if (num_noise_points_x > 1000000000.f ||
521 num_noise_points_y > 1000000000.f ||
522 num_noise_points_z > 1000000000.f)
523 throw InvalidNoiseParamsException();
525 // + 2 for the two initial endpoints
526 // + 1 for potentially crossing a boundary due to offset
527 size_t nlx = (size_t)ceil(num_noise_points_x) + 3;
528 size_t nly = (size_t)ceil(num_noise_points_y) + 3;
529 size_t nlz = is3d ? (size_t)ceil(num_noise_points_z) + 3 : 1;
533 noise_buf = new float[nlx * nly * nlz];
534 } catch (std::bad_alloc &e) {
535 throw InvalidNoiseParamsException();
541 * NB: This algorithm is not optimal in terms of space complexity. The entire
542 * integer lattice of noise points could be done as 2 lines instead, and for 3D,
543 * 2 lines + 2 planes.
544 * However, this would require the noise calls to be interposed with the
545 * interpolation loops, which may trash the icache, leading to lower overall
547 * Another optimization that could save half as many noise calls is to carry over
548 * values from the previous noise lattice as midpoints in the new lattice for the
551 #define idx(x, y) ((y) * nlx + (x))
552 void Noise::gradientMap2D(
554 float step_x, float step_y,
557 float v00, v01, v10, v11, u, v, orig_u;
558 u32 index, i, j, noisex, noisey;
562 bool eased = np.flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED);
563 Interp2dFxn interpolate = eased ?
564 biLinearInterpolation : biLinearInterpolationNoEase;
572 //calculate noise point lattice
573 nlx = (u32)(u + sx * step_x) + 2;
574 nly = (u32)(v + sy * step_y) + 2;
576 for (j = 0; j != nly; j++)
577 for (i = 0; i != nlx; i++)
578 noise_buf[index++] = noise2d(x0 + i, y0 + j, seed);
580 //calculate interpolations
583 for (j = 0; j != sy; j++) {
584 v00 = noise_buf[idx(0, noisey)];
585 v10 = noise_buf[idx(1, noisey)];
586 v01 = noise_buf[idx(0, noisey + 1)];
587 v11 = noise_buf[idx(1, noisey + 1)];
591 for (i = 0; i != sx; i++) {
592 gradient_buf[index++] = interpolate(v00, v10, v01, v11, u, v);
600 v10 = noise_buf[idx(noisex + 1, noisey)];
601 v11 = noise_buf[idx(noisex + 1, noisey + 1)];
615 #define idx(x, y, z) ((z) * nly * nlx + (y) * nlx + (x))
616 void Noise::gradientMap3D(
617 float x, float y, float z,
618 float step_x, float step_y, float step_z,
621 float v000, v010, v100, v110;
622 float v001, v011, v101, v111;
623 float u, v, w, orig_u, orig_v;
624 u32 index, i, j, k, noisex, noisey, noisez;
628 Interp3dFxn interpolate = (np.flags & NOISE_FLAG_EASED) ?
629 triLinearInterpolation : triLinearInterpolationNoEase;
640 //calculate noise point lattice
641 nlx = (u32)(u + sx * step_x) + 2;
642 nly = (u32)(v + sy * step_y) + 2;
643 nlz = (u32)(w + sz * step_z) + 2;
645 for (k = 0; k != nlz; k++)
646 for (j = 0; j != nly; j++)
647 for (i = 0; i != nlx; i++)
648 noise_buf[index++] = noise3d(x0 + i, y0 + j, z0 + k, seed);
650 //calculate interpolations
654 for (k = 0; k != sz; k++) {
657 for (j = 0; j != sy; j++) {
658 v000 = noise_buf[idx(0, noisey, noisez)];
659 v100 = noise_buf[idx(1, noisey, noisez)];
660 v010 = noise_buf[idx(0, noisey + 1, noisez)];
661 v110 = noise_buf[idx(1, noisey + 1, noisez)];
662 v001 = noise_buf[idx(0, noisey, noisez + 1)];
663 v101 = noise_buf[idx(1, noisey, noisez + 1)];
664 v011 = noise_buf[idx(0, noisey + 1, noisez + 1)];
665 v111 = noise_buf[idx(1, noisey + 1, noisez + 1)];
669 for (i = 0; i != sx; i++) {
670 gradient_buf[index++] = interpolate(
671 v000, v100, v010, v110,
672 v001, v101, v011, v111,
681 v100 = noise_buf[idx(noisex + 1, noisey, noisez)];
682 v110 = noise_buf[idx(noisex + 1, noisey + 1, noisez)];
685 v101 = noise_buf[idx(noisex + 1, noisey, noisez + 1)];
686 v111 = noise_buf[idx(noisex + 1, noisey + 1, noisez + 1)];
707 float *Noise::perlinMap2D(float x, float y, float *persistence_map)
709 float f = 1.0, g = 1.0;
710 size_t bufsize = sx * sy;
715 memset(result, 0, sizeof(float) * bufsize);
717 if (persistence_map) {
719 persist_buf = new float[bufsize];
720 for (size_t i = 0; i != bufsize; i++)
721 persist_buf[i] = 1.0;
724 for (size_t oct = 0; oct < np.octaves; oct++) {
725 gradientMap2D(x * f, y * f,
726 f / np.spread.X, f / np.spread.Y,
727 seed + np.seed + oct);
729 updateResults(g, persist_buf, persistence_map, bufsize);
735 if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
736 for (size_t i = 0; i != bufsize; i++)
737 result[i] = result[i] * np.scale + np.offset;
744 float *Noise::perlinMap3D(float x, float y, float z, float *persistence_map)
746 float f = 1.0, g = 1.0;
747 size_t bufsize = sx * sy * sz;
753 memset(result, 0, sizeof(float) * bufsize);
755 if (persistence_map) {
757 persist_buf = new float[bufsize];
758 for (size_t i = 0; i != bufsize; i++)
759 persist_buf[i] = 1.0;
762 for (size_t oct = 0; oct < np.octaves; oct++) {
763 gradientMap3D(x * f, y * f, z * f,
764 f / np.spread.X, f / np.spread.Y, f / np.spread.Z,
765 seed + np.seed + oct);
767 updateResults(g, persist_buf, persistence_map, bufsize);
773 if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
774 for (size_t i = 0; i != bufsize; i++)
775 result[i] = result[i] * np.scale + np.offset;
782 void Noise::updateResults(float g, float *gmap,
783 float *persistence_map, size_t bufsize)
785 // This looks very ugly, but it is 50-70% faster than having
786 // conditional statements inside the loop
787 if (np.flags & NOISE_FLAG_ABSVALUE) {
788 if (persistence_map) {
789 for (size_t i = 0; i != bufsize; i++) {
790 result[i] += gmap[i] * fabs(gradient_buf[i]);
791 gmap[i] *= persistence_map[i];
794 for (size_t i = 0; i != bufsize; i++)
795 result[i] += g * fabs(gradient_buf[i]);
798 if (persistence_map) {
799 for (size_t i = 0; i != bufsize; i++) {
800 result[i] += gmap[i] * gradient_buf[i];
801 gmap[i] *= persistence_map[i];
804 for (size_t i = 0; i != bufsize; i++)
805 result[i] += g * gradient_buf[i];