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);
197 v00 * (1 - tx) * (1 - ty) +
198 v10 * tx * (1 - ty) +
199 v01 * (1 - tx) * ty +
203 float u = linearInterpolation(v00, v10, tx);
204 float v = linearInterpolation(v01, v11, tx);
205 return linearInterpolation(u, v, ty);
209 inline float biLinearInterpolationNoEase(
210 float v00, float v10,
211 float v01, float v11,
214 float u = linearInterpolation(v00, v10, x);
215 float v = linearInterpolation(v01, v11, x);
216 return linearInterpolation(u, v, y);
220 float triLinearInterpolation(
221 float v000, float v100, float v010, float v110,
222 float v001, float v101, float v011, float v111,
223 float x, float y, float z)
225 float tx = easeCurve(x);
226 float ty = easeCurve(y);
227 float tz = easeCurve(z);
230 v000 * (1 - tx) * (1 - ty) * (1 - tz) +
231 v100 * tx * (1 - ty) * (1 - tz) +
232 v010 * (1 - tx) * ty * (1 - tz) +
233 v110 * tx * ty * (1 - tz) +
234 v001 * (1 - tx) * (1 - ty) * tz +
235 v101 * tx * (1 - ty) * tz +
236 v011 * (1 - tx) * ty * tz +
240 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, tx, ty);
241 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, tx, ty);
242 return linearInterpolation(u, v, tz);
245 float triLinearInterpolationNoEase(
246 float v000, float v100, float v010, float v110,
247 float v001, float v101, float v011, float v111,
248 float x, float y, float z)
250 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, x, y);
251 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, x, y);
252 return linearInterpolation(u, v, z);
257 float noise2d_gradient(float x, float y, int seed)
259 // Calculate the integer coordinates
260 int x0 = (x > 0.0 ? (int)x : (int)x - 1);
261 int y0 = (y > 0.0 ? (int)y : (int)y - 1);
262 // Calculate the remaining part of the coordinates
263 float xl = x - (float)x0;
264 float yl = y - (float)y0;
265 // Calculate random cosine lookup table indices for the integer corners.
266 // They are looked up as unit vector gradients from the lookup table.
267 int n00 = (int)((noise2d(x0, y0, seed)+1)*8);
268 int n10 = (int)((noise2d(x0+1, y0, seed)+1)*8);
269 int n01 = (int)((noise2d(x0, y0+1, seed)+1)*8);
270 int n11 = (int)((noise2d(x0+1, y0+1, seed)+1)*8);
271 // Make a dot product for the gradients and the positions, to get the values
272 float s = dotProduct(cos_lookup[n00], cos_lookup[(n00+12)%16], xl, yl);
273 float u = dotProduct(-cos_lookup[n10], cos_lookup[(n10+12)%16], 1.-xl, yl);
274 float v = dotProduct(cos_lookup[n01], -cos_lookup[(n01+12)%16], xl, 1.-yl);
275 float w = dotProduct(-cos_lookup[n11], -cos_lookup[(n11+12)%16], 1.-xl, 1.-yl);
276 // Interpolate between the values
277 return biLinearInterpolation(s,u,v,w,xl,yl);
282 float noise2d_gradient(float x, float y, int seed, bool eased)
284 // Calculate the integer coordinates
287 // Calculate the remaining part of the coordinates
288 float xl = x - (float)x0;
289 float yl = y - (float)y0;
290 // Get values for corners of square
291 float v00 = noise2d(x0, y0, seed);
292 float v10 = noise2d(x0+1, y0, seed);
293 float v01 = noise2d(x0, y0+1, seed);
294 float v11 = noise2d(x0+1, y0+1, seed);
297 return biLinearInterpolation(v00, v10, v01, v11, xl, yl);
299 return biLinearInterpolationNoEase(v00, v10, v01, v11, xl, yl);
303 float noise3d_gradient(float x, float y, float z, int seed, bool eased)
305 // Calculate the integer coordinates
309 // Calculate the remaining part of the coordinates
310 float xl = x - (float)x0;
311 float yl = y - (float)y0;
312 float zl = z - (float)z0;
313 // Get values for corners of cube
314 float v000 = noise3d(x0, y0, z0, seed);
315 float v100 = noise3d(x0 + 1, y0, z0, seed);
316 float v010 = noise3d(x0, y0 + 1, z0, seed);
317 float v110 = noise3d(x0 + 1, y0 + 1, z0, seed);
318 float v001 = noise3d(x0, y0, z0 + 1, seed);
319 float v101 = noise3d(x0 + 1, y0, z0 + 1, seed);
320 float v011 = noise3d(x0, y0 + 1, z0 + 1, seed);
321 float v111 = noise3d(x0 + 1, y0 + 1, z0 + 1, seed);
324 return triLinearInterpolation(
325 v000, v100, v010, v110,
326 v001, v101, v011, v111,
329 return triLinearInterpolationNoEase(
330 v000, v100, v010, v110,
331 v001, v101, v011, v111,
337 float noise2d_perlin(float x, float y, int seed,
338 int octaves, float persistence, bool eased)
343 for (int i = 0; i < octaves; i++)
345 a += g * noise2d_gradient(x * f, y * f, seed + i, eased);
353 float noise2d_perlin_abs(float x, float y, int seed,
354 int octaves, float persistence, bool eased)
359 for (int i = 0; i < octaves; i++) {
360 a += g * fabs(noise2d_gradient(x * f, y * f, seed + i, eased));
368 float noise3d_perlin(float x, float y, float z, int seed,
369 int octaves, float persistence, bool eased)
374 for (int i = 0; i < octaves; i++) {
375 a += g * noise3d_gradient(x * f, y * f, z * f, seed + i, eased);
383 float noise3d_perlin_abs(float x, float y, float z, int seed,
384 int octaves, float persistence, bool eased)
389 for (int i = 0; i < octaves; i++) {
390 a += g * fabs(noise3d_gradient(x * f, y * f, z * f, seed + i, eased));
398 float contour(float v)
407 ///////////////////////// [ New noise ] ////////////////////////////
410 float NoisePerlin2D(NoiseParams *np, float x, float y, int seed)
420 for (size_t i = 0; i < np->octaves; i++) {
421 float noiseval = noise2d_gradient(x * f, y * f, seed + i,
422 np->flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED));
424 if (np->flags & NOISE_FLAG_ABSVALUE)
425 noiseval = fabs(noiseval);
432 return np->offset + a * np->scale;
436 float NoisePerlin3D(NoiseParams *np, float x, float y, float z, int seed)
447 for (size_t i = 0; i < np->octaves; i++) {
448 float noiseval = noise3d_gradient(x * f, y * f, z * f, seed + i,
449 np->flags & NOISE_FLAG_EASED);
451 if (np->flags & NOISE_FLAG_ABSVALUE)
452 noiseval = fabs(noiseval);
459 return np->offset + a * np->scale;
463 Noise::Noise(NoiseParams *np_, int seed, u32 sx, u32 sy, u32 sz)
465 memcpy(&np, np_, sizeof(np));
471 this->persist_buf = NULL;
472 this->gradient_buf = NULL;
481 delete[] gradient_buf;
482 delete[] persist_buf;
488 void Noise::allocBuffers()
497 this->noise_buf = NULL;
498 resizeNoiseBuf(sz > 1);
500 delete[] gradient_buf;
501 delete[] persist_buf;
505 size_t bufsize = sx * sy * sz;
506 this->persist_buf = NULL;
507 this->gradient_buf = new float[bufsize];
508 this->result = new float[bufsize];
509 } catch (std::bad_alloc &e) {
510 throw InvalidNoiseParamsException();
515 void Noise::setSize(u32 sx, u32 sy, u32 sz)
525 void Noise::setSpreadFactor(v3f spread)
527 this->np.spread = spread;
529 resizeNoiseBuf(sz > 1);
533 void Noise::setOctaves(int octaves)
535 this->np.octaves = octaves;
537 resizeNoiseBuf(sz > 1);
541 void Noise::resizeNoiseBuf(bool is3d)
543 //maximum possible spread value factor
544 float ofactor = (np.lacunarity > 1.0) ?
545 pow(np.lacunarity, np.octaves - 1) :
548 // noise lattice point count
549 // (int)(sz * spread * ofactor) is # of lattice points crossed due to length
550 float num_noise_points_x = sx * ofactor / np.spread.X;
551 float num_noise_points_y = sy * ofactor / np.spread.Y;
552 float num_noise_points_z = sz * ofactor / np.spread.Z;
554 // protect against obviously invalid parameters
555 if (num_noise_points_x > 1000000000.f ||
556 num_noise_points_y > 1000000000.f ||
557 num_noise_points_z > 1000000000.f)
558 throw InvalidNoiseParamsException();
560 // + 2 for the two initial endpoints
561 // + 1 for potentially crossing a boundary due to offset
562 size_t nlx = (size_t)ceil(num_noise_points_x) + 3;
563 size_t nly = (size_t)ceil(num_noise_points_y) + 3;
564 size_t nlz = is3d ? (size_t)ceil(num_noise_points_z) + 3 : 1;
568 noise_buf = new float[nlx * nly * nlz];
569 } catch (std::bad_alloc &e) {
570 throw InvalidNoiseParamsException();
576 * NB: This algorithm is not optimal in terms of space complexity. The entire
577 * integer lattice of noise points could be done as 2 lines instead, and for 3D,
578 * 2 lines + 2 planes.
579 * However, this would require the noise calls to be interposed with the
580 * interpolation loops, which may trash the icache, leading to lower overall
582 * Another optimization that could save half as many noise calls is to carry over
583 * values from the previous noise lattice as midpoints in the new lattice for the
586 #define idx(x, y) ((y) * nlx + (x))
587 void Noise::gradientMap2D(
589 float step_x, float step_y,
592 float v00, v01, v10, v11, u, v, orig_u;
593 u32 index, i, j, noisex, noisey;
597 bool eased = np.flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED);
598 Interp2dFxn interpolate = eased ?
599 biLinearInterpolation : biLinearInterpolationNoEase;
607 //calculate noise point lattice
608 nlx = (u32)(u + sx * step_x) + 2;
609 nly = (u32)(v + sy * step_y) + 2;
611 for (j = 0; j != nly; j++)
612 for (i = 0; i != nlx; i++)
613 noise_buf[index++] = noise2d(x0 + i, y0 + j, seed);
615 //calculate interpolations
618 for (j = 0; j != sy; j++) {
619 v00 = noise_buf[idx(0, noisey)];
620 v10 = noise_buf[idx(1, noisey)];
621 v01 = noise_buf[idx(0, noisey + 1)];
622 v11 = noise_buf[idx(1, noisey + 1)];
626 for (i = 0; i != sx; i++) {
627 gradient_buf[index++] = interpolate(v00, v10, v01, v11, u, v);
635 v10 = noise_buf[idx(noisex + 1, noisey)];
636 v11 = noise_buf[idx(noisex + 1, noisey + 1)];
650 #define idx(x, y, z) ((z) * nly * nlx + (y) * nlx + (x))
651 void Noise::gradientMap3D(
652 float x, float y, float z,
653 float step_x, float step_y, float step_z,
656 float v000, v010, v100, v110;
657 float v001, v011, v101, v111;
658 float u, v, w, orig_u, orig_v;
659 u32 index, i, j, k, noisex, noisey, noisez;
663 Interp3dFxn interpolate = (np.flags & NOISE_FLAG_EASED) ?
664 triLinearInterpolation : triLinearInterpolationNoEase;
675 //calculate noise point lattice
676 nlx = (u32)(u + sx * step_x) + 2;
677 nly = (u32)(v + sy * step_y) + 2;
678 nlz = (u32)(w + sz * step_z) + 2;
680 for (k = 0; k != nlz; k++)
681 for (j = 0; j != nly; j++)
682 for (i = 0; i != nlx; i++)
683 noise_buf[index++] = noise3d(x0 + i, y0 + j, z0 + k, seed);
685 //calculate interpolations
689 for (k = 0; k != sz; k++) {
692 for (j = 0; j != sy; j++) {
693 v000 = noise_buf[idx(0, noisey, noisez)];
694 v100 = noise_buf[idx(1, noisey, noisez)];
695 v010 = noise_buf[idx(0, noisey + 1, noisez)];
696 v110 = noise_buf[idx(1, noisey + 1, noisez)];
697 v001 = noise_buf[idx(0, noisey, noisez + 1)];
698 v101 = noise_buf[idx(1, noisey, noisez + 1)];
699 v011 = noise_buf[idx(0, noisey + 1, noisez + 1)];
700 v111 = noise_buf[idx(1, noisey + 1, noisez + 1)];
704 for (i = 0; i != sx; i++) {
705 gradient_buf[index++] = interpolate(
706 v000, v100, v010, v110,
707 v001, v101, v011, v111,
716 v100 = noise_buf[idx(noisex + 1, noisey, noisez)];
717 v110 = noise_buf[idx(noisex + 1, noisey + 1, noisez)];
720 v101 = noise_buf[idx(noisex + 1, noisey, noisez + 1)];
721 v111 = noise_buf[idx(noisex + 1, noisey + 1, noisez + 1)];
742 float *Noise::perlinMap2D(float x, float y, float *persistence_map)
744 float f = 1.0, g = 1.0;
745 size_t bufsize = sx * sy;
750 memset(result, 0, sizeof(float) * bufsize);
752 if (persistence_map) {
754 persist_buf = new float[bufsize];
755 for (size_t i = 0; i != bufsize; i++)
756 persist_buf[i] = 1.0;
759 for (size_t oct = 0; oct < np.octaves; oct++) {
760 gradientMap2D(x * f, y * f,
761 f / np.spread.X, f / np.spread.Y,
762 seed + np.seed + oct);
764 updateResults(g, persist_buf, persistence_map, bufsize);
770 if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
771 for (size_t i = 0; i != bufsize; i++)
772 result[i] = result[i] * np.scale + np.offset;
779 float *Noise::perlinMap3D(float x, float y, float z, float *persistence_map)
781 float f = 1.0, g = 1.0;
782 size_t bufsize = sx * sy * sz;
788 memset(result, 0, sizeof(float) * bufsize);
790 if (persistence_map) {
792 persist_buf = new float[bufsize];
793 for (size_t i = 0; i != bufsize; i++)
794 persist_buf[i] = 1.0;
797 for (size_t oct = 0; oct < np.octaves; oct++) {
798 gradientMap3D(x * f, y * f, z * f,
799 f / np.spread.X, f / np.spread.Y, f / np.spread.Z,
800 seed + np.seed + oct);
802 updateResults(g, persist_buf, persistence_map, bufsize);
808 if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
809 for (size_t i = 0; i != bufsize; i++)
810 result[i] = result[i] * np.scale + np.offset;
817 void Noise::updateResults(float g, float *gmap,
818 float *persistence_map, size_t bufsize)
820 // This looks very ugly, but it is 50-70% faster than having
821 // conditional statements inside the loop
822 if (np.flags & NOISE_FLAG_ABSVALUE) {
823 if (persistence_map) {
824 for (size_t i = 0; i != bufsize; i++) {
825 result[i] += gmap[i] * fabs(gradient_buf[i]);
826 gmap[i] *= persistence_map[i];
829 for (size_t i = 0; i != bufsize; i++)
830 result[i] += g * fabs(gradient_buf[i]);
833 if (persistence_map) {
834 for (size_t i = 0; i != bufsize; i++) {
835 result[i] += gmap[i] * gradient_buf[i];
836 gmap[i] *= persistence_map[i];
839 for (size_t i = 0; i != bufsize; i++)
840 result[i] += g * gradient_buf[i];