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 round((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, int sx, int sy, int 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(int sx, int sy, int 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 int index, i, j, x0, y0, noisex, noisey;
596 bool eased = np.flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED);
597 Interp2dFxn interpolate = eased ?
598 biLinearInterpolation : biLinearInterpolationNoEase;
606 //calculate noise point lattice
607 nlx = (int)(u + sx * step_x) + 2;
608 nly = (int)(v + sy * step_y) + 2;
610 for (j = 0; j != nly; j++)
611 for (i = 0; i != nlx; i++)
612 noise_buf[index++] = noise2d(x0 + i, y0 + j, seed);
614 //calculate interpolations
617 for (j = 0; j != sy; j++) {
618 v00 = noise_buf[idx(0, noisey)];
619 v10 = noise_buf[idx(1, noisey)];
620 v01 = noise_buf[idx(0, noisey + 1)];
621 v11 = noise_buf[idx(1, noisey + 1)];
625 for (i = 0; i != sx; i++) {
626 gradient_buf[index++] = interpolate(v00, v10, v01, v11, u, v);
634 v10 = noise_buf[idx(noisex + 1, noisey)];
635 v11 = noise_buf[idx(noisex + 1, noisey + 1)];
649 #define idx(x, y, z) ((z) * nly * nlx + (y) * nlx + (x))
650 void Noise::gradientMap3D(
651 float x, float y, float z,
652 float step_x, float step_y, float step_z,
655 float v000, v010, v100, v110;
656 float v001, v011, v101, v111;
657 float u, v, w, orig_u, orig_v;
658 int index, i, j, k, x0, y0, z0, noisex, noisey, noisez;
661 Interp3dFxn interpolate = (np.flags & NOISE_FLAG_EASED) ?
662 triLinearInterpolation : triLinearInterpolationNoEase;
673 //calculate noise point lattice
674 nlx = (int)(u + sx * step_x) + 2;
675 nly = (int)(v + sy * step_y) + 2;
676 nlz = (int)(w + sz * step_z) + 2;
678 for (k = 0; k != nlz; k++)
679 for (j = 0; j != nly; j++)
680 for (i = 0; i != nlx; i++)
681 noise_buf[index++] = noise3d(x0 + i, y0 + j, z0 + k, seed);
683 //calculate interpolations
687 for (k = 0; k != sz; k++) {
690 for (j = 0; j != sy; j++) {
691 v000 = noise_buf[idx(0, noisey, noisez)];
692 v100 = noise_buf[idx(1, noisey, noisez)];
693 v010 = noise_buf[idx(0, noisey + 1, noisez)];
694 v110 = noise_buf[idx(1, noisey + 1, noisez)];
695 v001 = noise_buf[idx(0, noisey, noisez + 1)];
696 v101 = noise_buf[idx(1, noisey, noisez + 1)];
697 v011 = noise_buf[idx(0, noisey + 1, noisez + 1)];
698 v111 = noise_buf[idx(1, noisey + 1, noisez + 1)];
702 for (i = 0; i != sx; i++) {
703 gradient_buf[index++] = interpolate(
704 v000, v100, v010, v110,
705 v001, v101, v011, v111,
714 v100 = noise_buf[idx(noisex + 1, noisey, noisez)];
715 v110 = noise_buf[idx(noisex + 1, noisey + 1, noisez)];
718 v101 = noise_buf[idx(noisex + 1, noisey, noisez + 1)];
719 v111 = noise_buf[idx(noisex + 1, noisey + 1, noisez + 1)];
740 float *Noise::perlinMap2D(float x, float y, float *persistence_map)
742 float f = 1.0, g = 1.0;
743 size_t bufsize = sx * sy;
748 memset(result, 0, sizeof(float) * bufsize);
750 if (persistence_map) {
752 persist_buf = new float[bufsize];
753 for (size_t i = 0; i != bufsize; i++)
754 persist_buf[i] = 1.0;
757 for (size_t oct = 0; oct < np.octaves; oct++) {
758 gradientMap2D(x * f, y * f,
759 f / np.spread.X, f / np.spread.Y,
760 seed + np.seed + oct);
762 updateResults(g, persist_buf, persistence_map, bufsize);
768 if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
769 for (size_t i = 0; i != bufsize; i++)
770 result[i] = result[i] * np.scale + np.offset;
777 float *Noise::perlinMap3D(float x, float y, float z, float *persistence_map)
779 float f = 1.0, g = 1.0;
780 size_t bufsize = sx * sy * sz;
786 memset(result, 0, sizeof(float) * bufsize);
788 if (persistence_map) {
790 persist_buf = new float[bufsize];
791 for (size_t i = 0; i != bufsize; i++)
792 persist_buf[i] = 1.0;
795 for (size_t oct = 0; oct < np.octaves; oct++) {
796 gradientMap3D(x * f, y * f, z * f,
797 f / np.spread.X, f / np.spread.Y, f / np.spread.Z,
798 seed + np.seed + oct);
800 updateResults(g, persist_buf, persistence_map, bufsize);
806 if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
807 for (size_t i = 0; i != bufsize; i++)
808 result[i] = result[i] * np.scale + np.offset;
815 void Noise::updateResults(float g, float *gmap,
816 float *persistence_map, size_t bufsize)
818 // This looks very ugly, but it is 50-70% faster than having
819 // conditional statements inside the loop
820 if (np.flags & NOISE_FLAG_ABSVALUE) {
821 if (persistence_map) {
822 for (size_t i = 0; i != bufsize; i++) {
823 result[i] += gmap[i] * fabs(gradient_buf[i]);
824 gmap[i] *= persistence_map[i];
827 for (size_t i = 0; i != bufsize; i++)
828 result[i] += g * fabs(gradient_buf[i]);
831 if (persistence_map) {
832 for (size_t i = 0; i != bufsize; i++) {
833 result[i] += gmap[i] * gradient_buf[i];
834 gmap[i] *= persistence_map[i];
837 for (size_t i = 0; i != bufsize; i++)
838 result[i] += g * gradient_buf[i];