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)
94 If the bound is not a multiple of the RNG's range, it may cause bias,
95 e.g. a RNG has a range from 0 to 3 and we take want a number 0 to 2.
96 Using rand() % 3, the number 0 would be twice as likely to appear.
97 With a very large RNG range, the effect becomes less prevalent but
98 still present. This can be solved by modifying the range of the RNG
99 to become a multiple of bound by dropping values above the a threshhold.
100 In our example, threshhold == 4 - 3 = 1 % 3 == 1, so reject 0, thus
101 making the range 3 with no bias.
103 This loop looks dangerous, but will always terminate due to the
104 RNG's property of uniformity.
106 u32 threshhold = -bound % bound;
109 while ((r = next()) < threshhold)
116 s32 PcgRandom::range(s32 min, s32 max)
119 throw PrngException("Invalid range (max < min)");
121 u32 bound = max - min + 1;
122 return range(bound) + min;
126 void PcgRandom::bytes(void *out, size_t len)
128 u8 *outb = (u8 *)out;
133 if (bytes_left == 0) {
134 bytes_left = sizeof(u32);
146 s32 PcgRandom::randNormalDist(s32 min, s32 max, int num_trials)
149 for (int i = 0; i != num_trials; i++)
150 accum += range(min, max);
151 return myround((float)accum / num_trials);
154 ///////////////////////////////////////////////////////////////////////////////
156 float noise2d(int x, int y, int seed)
158 int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y
159 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
161 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
162 return 1.f - (float)n / 0x40000000;
166 float noise3d(int x, int y, int z, int seed)
168 int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y + NOISE_MAGIC_Z * z
169 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
171 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
172 return 1.f - (float)n / 0x40000000;
176 inline float dotProduct(float vx, float vy, float wx, float wy)
178 return vx * wx + vy * wy;
182 inline float linearInterpolation(float v0, float v1, float t)
184 return v0 + (v1 - v0) * t;
188 inline float biLinearInterpolation(
189 float v00, float v10,
190 float v01, float v11,
193 float tx = easeCurve(x);
194 float ty = easeCurve(y);
195 float u = linearInterpolation(v00, v10, tx);
196 float v = linearInterpolation(v01, v11, tx);
197 return linearInterpolation(u, v, ty);
201 inline float biLinearInterpolationNoEase(
202 float v00, float v10,
203 float v01, float v11,
206 float u = linearInterpolation(v00, v10, x);
207 float v = linearInterpolation(v01, v11, x);
208 return linearInterpolation(u, v, y);
212 float triLinearInterpolation(
213 float v000, float v100, float v010, float v110,
214 float v001, float v101, float v011, float v111,
215 float x, float y, float z)
217 float tx = easeCurve(x);
218 float ty = easeCurve(y);
219 float tz = easeCurve(z);
220 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, tx, ty);
221 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, tx, ty);
222 return linearInterpolation(u, v, tz);
225 float triLinearInterpolationNoEase(
226 float v000, float v100, float v010, float v110,
227 float v001, float v101, float v011, float v111,
228 float x, float y, float z)
230 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, x, y);
231 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, x, y);
232 return linearInterpolation(u, v, z);
235 float noise2d_gradient(float x, float y, int seed, bool eased)
237 // Calculate the integer coordinates
240 // Calculate the remaining part of the coordinates
241 float xl = x - (float)x0;
242 float yl = y - (float)y0;
243 // Get values for corners of square
244 float v00 = noise2d(x0, y0, seed);
245 float v10 = noise2d(x0+1, y0, seed);
246 float v01 = noise2d(x0, y0+1, seed);
247 float v11 = noise2d(x0+1, y0+1, seed);
250 return biLinearInterpolation(v00, v10, v01, v11, xl, yl);
252 return biLinearInterpolationNoEase(v00, v10, v01, v11, xl, yl);
256 float noise3d_gradient(float x, float y, float z, int seed, bool eased)
258 // Calculate the integer coordinates
262 // Calculate the remaining part of the coordinates
263 float xl = x - (float)x0;
264 float yl = y - (float)y0;
265 float zl = z - (float)z0;
266 // Get values for corners of cube
267 float v000 = noise3d(x0, y0, z0, seed);
268 float v100 = noise3d(x0 + 1, y0, z0, seed);
269 float v010 = noise3d(x0, y0 + 1, z0, seed);
270 float v110 = noise3d(x0 + 1, y0 + 1, z0, seed);
271 float v001 = noise3d(x0, y0, z0 + 1, seed);
272 float v101 = noise3d(x0 + 1, y0, z0 + 1, seed);
273 float v011 = noise3d(x0, y0 + 1, z0 + 1, seed);
274 float v111 = noise3d(x0 + 1, y0 + 1, z0 + 1, seed);
277 return triLinearInterpolation(
278 v000, v100, v010, v110,
279 v001, v101, v011, v111,
282 return triLinearInterpolationNoEase(
283 v000, v100, v010, v110,
284 v001, v101, v011, v111,
290 float noise2d_perlin(float x, float y, int seed,
291 int octaves, float persistence, bool eased)
296 for (int i = 0; i < octaves; i++)
298 a += g * noise2d_gradient(x * f, y * f, seed + i, eased);
306 float noise2d_perlin_abs(float x, float y, int seed,
307 int octaves, float persistence, bool eased)
312 for (int i = 0; i < octaves; i++) {
313 a += g * fabs(noise2d_gradient(x * f, y * f, seed + i, eased));
321 float noise3d_perlin(float x, float y, float z, int seed,
322 int octaves, float persistence, bool eased)
327 for (int i = 0; i < octaves; i++) {
328 a += g * noise3d_gradient(x * f, y * f, z * f, seed + i, eased);
336 float noise3d_perlin_abs(float x, float y, float z, int seed,
337 int octaves, float persistence, bool eased)
342 for (int i = 0; i < octaves; i++) {
343 a += g * fabs(noise3d_gradient(x * f, y * f, z * f, seed + i, eased));
351 float contour(float v)
360 ///////////////////////// [ New noise ] ////////////////////////////
363 float NoisePerlin2D(NoiseParams *np, float x, float y, int seed)
373 for (size_t i = 0; i < np->octaves; i++) {
374 float noiseval = noise2d_gradient(x * f, y * f, seed + i,
375 np->flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED));
377 if (np->flags & NOISE_FLAG_ABSVALUE)
378 noiseval = fabs(noiseval);
385 return np->offset + a * np->scale;
389 float NoisePerlin3D(NoiseParams *np, float x, float y, float z, int seed)
400 for (size_t i = 0; i < np->octaves; i++) {
401 float noiseval = noise3d_gradient(x * f, y * f, z * f, seed + i,
402 np->flags & NOISE_FLAG_EASED);
404 if (np->flags & NOISE_FLAG_ABSVALUE)
405 noiseval = fabs(noiseval);
412 return np->offset + a * np->scale;
416 Noise::Noise(NoiseParams *np_, int seed, u32 sx, u32 sy, u32 sz)
418 memcpy(&np, np_, sizeof(np));
424 this->persist_buf = NULL;
425 this->gradient_buf = NULL;
434 delete[] gradient_buf;
435 delete[] persist_buf;
441 void Noise::allocBuffers()
450 this->noise_buf = NULL;
451 resizeNoiseBuf(sz > 1);
453 delete[] gradient_buf;
454 delete[] persist_buf;
458 size_t bufsize = sx * sy * sz;
459 this->persist_buf = NULL;
460 this->gradient_buf = new float[bufsize];
461 this->result = new float[bufsize];
462 } catch (std::bad_alloc &e) {
463 throw InvalidNoiseParamsException();
468 void Noise::setSize(u32 sx, u32 sy, u32 sz)
478 void Noise::setSpreadFactor(v3f spread)
480 this->np.spread = spread;
482 resizeNoiseBuf(sz > 1);
486 void Noise::setOctaves(int octaves)
488 this->np.octaves = octaves;
490 resizeNoiseBuf(sz > 1);
494 void Noise::resizeNoiseBuf(bool is3d)
496 //maximum possible spread value factor
497 float ofactor = (np.lacunarity > 1.0) ?
498 pow(np.lacunarity, np.octaves - 1) :
501 // noise lattice point count
502 // (int)(sz * spread * ofactor) is # of lattice points crossed due to length
503 float num_noise_points_x = sx * ofactor / np.spread.X;
504 float num_noise_points_y = sy * ofactor / np.spread.Y;
505 float num_noise_points_z = sz * ofactor / np.spread.Z;
507 // protect against obviously invalid parameters
508 if (num_noise_points_x > 1000000000.f ||
509 num_noise_points_y > 1000000000.f ||
510 num_noise_points_z > 1000000000.f)
511 throw InvalidNoiseParamsException();
513 // + 2 for the two initial endpoints
514 // + 1 for potentially crossing a boundary due to offset
515 size_t nlx = (size_t)ceil(num_noise_points_x) + 3;
516 size_t nly = (size_t)ceil(num_noise_points_y) + 3;
517 size_t nlz = is3d ? (size_t)ceil(num_noise_points_z) + 3 : 1;
521 noise_buf = new float[nlx * nly * nlz];
522 } catch (std::bad_alloc &e) {
523 throw InvalidNoiseParamsException();
529 * NB: This algorithm is not optimal in terms of space complexity. The entire
530 * integer lattice of noise points could be done as 2 lines instead, and for 3D,
531 * 2 lines + 2 planes.
532 * However, this would require the noise calls to be interposed with the
533 * interpolation loops, which may trash the icache, leading to lower overall
535 * Another optimization that could save half as many noise calls is to carry over
536 * values from the previous noise lattice as midpoints in the new lattice for the
539 #define idx(x, y) ((y) * nlx + (x))
540 void Noise::gradientMap2D(
542 float step_x, float step_y,
545 float v00, v01, v10, v11, u, v, orig_u;
546 u32 index, i, j, noisex, noisey;
550 bool eased = np.flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED);
551 Interp2dFxn interpolate = eased ?
552 biLinearInterpolation : biLinearInterpolationNoEase;
560 //calculate noise point lattice
561 nlx = (u32)(u + sx * step_x) + 2;
562 nly = (u32)(v + sy * step_y) + 2;
564 for (j = 0; j != nly; j++)
565 for (i = 0; i != nlx; i++)
566 noise_buf[index++] = noise2d(x0 + i, y0 + j, seed);
568 //calculate interpolations
571 for (j = 0; j != sy; j++) {
572 v00 = noise_buf[idx(0, noisey)];
573 v10 = noise_buf[idx(1, noisey)];
574 v01 = noise_buf[idx(0, noisey + 1)];
575 v11 = noise_buf[idx(1, noisey + 1)];
579 for (i = 0; i != sx; i++) {
580 gradient_buf[index++] = interpolate(v00, v10, v01, v11, u, v);
588 v10 = noise_buf[idx(noisex + 1, noisey)];
589 v11 = noise_buf[idx(noisex + 1, noisey + 1)];
603 #define idx(x, y, z) ((z) * nly * nlx + (y) * nlx + (x))
604 void Noise::gradientMap3D(
605 float x, float y, float z,
606 float step_x, float step_y, float step_z,
609 float v000, v010, v100, v110;
610 float v001, v011, v101, v111;
611 float u, v, w, orig_u, orig_v;
612 u32 index, i, j, k, noisex, noisey, noisez;
616 Interp3dFxn interpolate = (np.flags & NOISE_FLAG_EASED) ?
617 triLinearInterpolation : triLinearInterpolationNoEase;
628 //calculate noise point lattice
629 nlx = (u32)(u + sx * step_x) + 2;
630 nly = (u32)(v + sy * step_y) + 2;
631 nlz = (u32)(w + sz * step_z) + 2;
633 for (k = 0; k != nlz; k++)
634 for (j = 0; j != nly; j++)
635 for (i = 0; i != nlx; i++)
636 noise_buf[index++] = noise3d(x0 + i, y0 + j, z0 + k, seed);
638 //calculate interpolations
642 for (k = 0; k != sz; k++) {
645 for (j = 0; j != sy; j++) {
646 v000 = noise_buf[idx(0, noisey, noisez)];
647 v100 = noise_buf[idx(1, noisey, noisez)];
648 v010 = noise_buf[idx(0, noisey + 1, noisez)];
649 v110 = noise_buf[idx(1, noisey + 1, noisez)];
650 v001 = noise_buf[idx(0, noisey, noisez + 1)];
651 v101 = noise_buf[idx(1, noisey, noisez + 1)];
652 v011 = noise_buf[idx(0, noisey + 1, noisez + 1)];
653 v111 = noise_buf[idx(1, noisey + 1, noisez + 1)];
657 for (i = 0; i != sx; i++) {
658 gradient_buf[index++] = interpolate(
659 v000, v100, v010, v110,
660 v001, v101, v011, v111,
669 v100 = noise_buf[idx(noisex + 1, noisey, noisez)];
670 v110 = noise_buf[idx(noisex + 1, noisey + 1, noisez)];
673 v101 = noise_buf[idx(noisex + 1, noisey, noisez + 1)];
674 v111 = noise_buf[idx(noisex + 1, noisey + 1, noisez + 1)];
695 float *Noise::perlinMap2D(float x, float y, float *persistence_map)
697 float f = 1.0, g = 1.0;
698 size_t bufsize = sx * sy;
703 memset(result, 0, sizeof(float) * bufsize);
705 if (persistence_map) {
707 persist_buf = new float[bufsize];
708 for (size_t i = 0; i != bufsize; i++)
709 persist_buf[i] = 1.0;
712 for (size_t oct = 0; oct < np.octaves; oct++) {
713 gradientMap2D(x * f, y * f,
714 f / np.spread.X, f / np.spread.Y,
715 seed + np.seed + oct);
717 updateResults(g, persist_buf, persistence_map, bufsize);
723 if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
724 for (size_t i = 0; i != bufsize; i++)
725 result[i] = result[i] * np.scale + np.offset;
732 float *Noise::perlinMap3D(float x, float y, float z, float *persistence_map)
734 float f = 1.0, g = 1.0;
735 size_t bufsize = sx * sy * sz;
741 memset(result, 0, sizeof(float) * bufsize);
743 if (persistence_map) {
745 persist_buf = new float[bufsize];
746 for (size_t i = 0; i != bufsize; i++)
747 persist_buf[i] = 1.0;
750 for (size_t oct = 0; oct < np.octaves; oct++) {
751 gradientMap3D(x * f, y * f, z * f,
752 f / np.spread.X, f / np.spread.Y, f / np.spread.Z,
753 seed + np.seed + oct);
755 updateResults(g, persist_buf, persistence_map, bufsize);
761 if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
762 for (size_t i = 0; i != bufsize; i++)
763 result[i] = result[i] * np.scale + np.offset;
770 void Noise::updateResults(float g, float *gmap,
771 float *persistence_map, size_t bufsize)
773 // This looks very ugly, but it is 50-70% faster than having
774 // conditional statements inside the loop
775 if (np.flags & NOISE_FLAG_ABSVALUE) {
776 if (persistence_map) {
777 for (size_t i = 0; i != bufsize; i++) {
778 result[i] += gmap[i] * fabs(gradient_buf[i]);
779 gmap[i] *= persistence_map[i];
782 for (size_t i = 0; i != bufsize; i++)
783 result[i] += g * fabs(gradient_buf[i]);
786 if (persistence_map) {
787 for (size_t i = 0; i != bufsize; i++) {
788 result[i] += gmap[i] * gradient_buf[i];
789 gmap[i] *= persistence_map[i];
792 for (size_t i = 0; i != bufsize; i++)
793 result[i] += g * gradient_buf[i];