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 u32 bound = max - min + 1;
120 return range(bound) + min;
124 void PcgRandom::bytes(void *out, size_t len)
126 u8 *outb = (u8 *)out;
131 if (bytes_left == 0) {
132 bytes_left = sizeof(u32);
144 s32 PcgRandom::randNormalDist(s32 min, s32 max, int num_trials)
147 for (int i = 0; i != num_trials; i++)
148 accum += range(min, max);
149 return ((float)accum / num_trials) + 0.5f;
152 ///////////////////////////////////////////////////////////////////////////////
154 float noise2d(int x, int y, int seed)
156 int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y
157 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
159 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
160 return 1.f - (float)n / 0x40000000;
164 float noise3d(int x, int y, int z, int seed)
166 int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y + NOISE_MAGIC_Z * z
167 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
169 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
170 return 1.f - (float)n / 0x40000000;
174 inline float dotProduct(float vx, float vy, float wx, float wy)
176 return vx * wx + vy * wy;
180 inline float linearInterpolation(float v0, float v1, float t)
182 return v0 + (v1 - v0) * t;
186 inline float biLinearInterpolation(
187 float v00, float v10,
188 float v01, float v11,
191 float tx = easeCurve(x);
192 float ty = easeCurve(y);
195 v00 * (1 - tx) * (1 - ty) +
196 v10 * tx * (1 - ty) +
197 v01 * (1 - tx) * ty +
201 float u = linearInterpolation(v00, v10, tx);
202 float v = linearInterpolation(v01, v11, tx);
203 return linearInterpolation(u, v, ty);
207 inline float biLinearInterpolationNoEase(
208 float v00, float v10,
209 float v01, float v11,
212 float u = linearInterpolation(v00, v10, x);
213 float v = linearInterpolation(v01, v11, x);
214 return linearInterpolation(u, v, y);
218 float triLinearInterpolation(
219 float v000, float v100, float v010, float v110,
220 float v001, float v101, float v011, float v111,
221 float x, float y, float z)
223 float tx = easeCurve(x);
224 float ty = easeCurve(y);
225 float tz = easeCurve(z);
228 v000 * (1 - tx) * (1 - ty) * (1 - tz) +
229 v100 * tx * (1 - ty) * (1 - tz) +
230 v010 * (1 - tx) * ty * (1 - tz) +
231 v110 * tx * ty * (1 - tz) +
232 v001 * (1 - tx) * (1 - ty) * tz +
233 v101 * tx * (1 - ty) * tz +
234 v011 * (1 - tx) * ty * tz +
238 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, tx, ty);
239 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, tx, ty);
240 return linearInterpolation(u, v, tz);
243 float triLinearInterpolationNoEase(
244 float v000, float v100, float v010, float v110,
245 float v001, float v101, float v011, float v111,
246 float x, float y, float z)
248 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, x, y);
249 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, x, y);
250 return linearInterpolation(u, v, z);
255 float noise2d_gradient(float x, float y, int seed)
257 // Calculate the integer coordinates
258 int x0 = (x > 0.0 ? (int)x : (int)x - 1);
259 int y0 = (y > 0.0 ? (int)y : (int)y - 1);
260 // Calculate the remaining part of the coordinates
261 float xl = x - (float)x0;
262 float yl = y - (float)y0;
263 // Calculate random cosine lookup table indices for the integer corners.
264 // They are looked up as unit vector gradients from the lookup table.
265 int n00 = (int)((noise2d(x0, y0, seed)+1)*8);
266 int n10 = (int)((noise2d(x0+1, y0, seed)+1)*8);
267 int n01 = (int)((noise2d(x0, y0+1, seed)+1)*8);
268 int n11 = (int)((noise2d(x0+1, y0+1, seed)+1)*8);
269 // Make a dot product for the gradients and the positions, to get the values
270 float s = dotProduct(cos_lookup[n00], cos_lookup[(n00+12)%16], xl, yl);
271 float u = dotProduct(-cos_lookup[n10], cos_lookup[(n10+12)%16], 1.-xl, yl);
272 float v = dotProduct(cos_lookup[n01], -cos_lookup[(n01+12)%16], xl, 1.-yl);
273 float w = dotProduct(-cos_lookup[n11], -cos_lookup[(n11+12)%16], 1.-xl, 1.-yl);
274 // Interpolate between the values
275 return biLinearInterpolation(s,u,v,w,xl,yl);
280 float noise2d_gradient(float x, float y, int seed, bool eased)
282 // Calculate the integer coordinates
285 // Calculate the remaining part of the coordinates
286 float xl = x - (float)x0;
287 float yl = y - (float)y0;
288 // Get values for corners of square
289 float v00 = noise2d(x0, y0, seed);
290 float v10 = noise2d(x0+1, y0, seed);
291 float v01 = noise2d(x0, y0+1, seed);
292 float v11 = noise2d(x0+1, y0+1, seed);
295 return biLinearInterpolation(v00, v10, v01, v11, xl, yl);
297 return biLinearInterpolationNoEase(v00, v10, v01, v11, xl, yl);
301 float noise3d_gradient(float x, float y, float z, int seed, bool eased)
303 // Calculate the integer coordinates
307 // Calculate the remaining part of the coordinates
308 float xl = x - (float)x0;
309 float yl = y - (float)y0;
310 float zl = z - (float)z0;
311 // Get values for corners of cube
312 float v000 = noise3d(x0, y0, z0, seed);
313 float v100 = noise3d(x0 + 1, y0, z0, seed);
314 float v010 = noise3d(x0, y0 + 1, z0, seed);
315 float v110 = noise3d(x0 + 1, y0 + 1, z0, seed);
316 float v001 = noise3d(x0, y0, z0 + 1, seed);
317 float v101 = noise3d(x0 + 1, y0, z0 + 1, seed);
318 float v011 = noise3d(x0, y0 + 1, z0 + 1, seed);
319 float v111 = noise3d(x0 + 1, y0 + 1, z0 + 1, seed);
322 return triLinearInterpolation(
323 v000, v100, v010, v110,
324 v001, v101, v011, v111,
327 return triLinearInterpolationNoEase(
328 v000, v100, v010, v110,
329 v001, v101, v011, v111,
335 float noise2d_perlin(float x, float y, int seed,
336 int octaves, float persistence, bool eased)
341 for (int i = 0; i < octaves; i++)
343 a += g * noise2d_gradient(x * f, y * f, seed + i, eased);
351 float noise2d_perlin_abs(float x, float y, int seed,
352 int octaves, float persistence, bool eased)
357 for (int i = 0; i < octaves; i++) {
358 a += g * fabs(noise2d_gradient(x * f, y * f, seed + i, eased));
366 float noise3d_perlin(float x, float y, float z, int seed,
367 int octaves, float persistence, bool eased)
372 for (int i = 0; i < octaves; i++) {
373 a += g * noise3d_gradient(x * f, y * f, z * f, seed + i, eased);
381 float noise3d_perlin_abs(float x, float y, float z, int seed,
382 int octaves, float persistence, bool eased)
387 for (int i = 0; i < octaves; i++) {
388 a += g * fabs(noise3d_gradient(x * f, y * f, z * f, seed + i, eased));
396 float contour(float v)
405 ///////////////////////// [ New noise ] ////////////////////////////
408 float NoisePerlin2D(NoiseParams *np, float x, float y, int seed)
418 for (size_t i = 0; i < np->octaves; i++) {
419 float noiseval = noise2d_gradient(x * f, y * f, seed + i,
420 np->flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED));
422 if (np->flags & NOISE_FLAG_ABSVALUE)
423 noiseval = fabs(noiseval);
430 return np->offset + a * np->scale;
434 float NoisePerlin3D(NoiseParams *np, float x, float y, float z, int seed)
445 for (size_t i = 0; i < np->octaves; i++) {
446 float noiseval = noise3d_gradient(x * f, y * f, z * f, seed + i,
447 np->flags & NOISE_FLAG_EASED);
449 if (np->flags & NOISE_FLAG_ABSVALUE)
450 noiseval = fabs(noiseval);
457 return np->offset + a * np->scale;
461 Noise::Noise(NoiseParams *np_, int seed, int sx, int sy, int sz)
463 memcpy(&np, np_, sizeof(np));
469 this->persist_buf = NULL;
470 this->gradient_buf = NULL;
479 delete[] gradient_buf;
480 delete[] persist_buf;
486 void Noise::allocBuffers()
495 this->noise_buf = NULL;
496 resizeNoiseBuf(sz > 1);
498 delete[] gradient_buf;
499 delete[] persist_buf;
503 size_t bufsize = sx * sy * sz;
504 this->persist_buf = NULL;
505 this->gradient_buf = new float[bufsize];
506 this->result = new float[bufsize];
507 } catch (std::bad_alloc &e) {
508 throw InvalidNoiseParamsException();
513 void Noise::setSize(int sx, int sy, int sz)
523 void Noise::setSpreadFactor(v3f spread)
525 this->np.spread = spread;
527 resizeNoiseBuf(sz > 1);
531 void Noise::setOctaves(int octaves)
533 this->np.octaves = octaves;
535 resizeNoiseBuf(sz > 1);
539 void Noise::resizeNoiseBuf(bool is3d)
541 //maximum possible spread value factor
542 float ofactor = (np.lacunarity > 1.0) ?
543 pow(np.lacunarity, np.octaves - 1) :
546 // noise lattice point count
547 // (int)(sz * spread * ofactor) is # of lattice points crossed due to length
548 float num_noise_points_x = sx * ofactor / np.spread.X;
549 float num_noise_points_y = sy * ofactor / np.spread.Y;
550 float num_noise_points_z = sz * ofactor / np.spread.Z;
552 // protect against obviously invalid parameters
553 if (num_noise_points_x > 1000000000.f ||
554 num_noise_points_y > 1000000000.f ||
555 num_noise_points_z > 1000000000.f)
556 throw InvalidNoiseParamsException();
558 // + 2 for the two initial endpoints
559 // + 1 for potentially crossing a boundary due to offset
560 size_t nlx = (size_t)ceil(num_noise_points_x) + 3;
561 size_t nly = (size_t)ceil(num_noise_points_y) + 3;
562 size_t nlz = is3d ? (size_t)ceil(num_noise_points_z) + 3 : 1;
566 noise_buf = new float[nlx * nly * nlz];
567 } catch (std::bad_alloc &e) {
568 throw InvalidNoiseParamsException();
574 * NB: This algorithm is not optimal in terms of space complexity. The entire
575 * integer lattice of noise points could be done as 2 lines instead, and for 3D,
576 * 2 lines + 2 planes.
577 * However, this would require the noise calls to be interposed with the
578 * interpolation loops, which may trash the icache, leading to lower overall
580 * Another optimization that could save half as many noise calls is to carry over
581 * values from the previous noise lattice as midpoints in the new lattice for the
584 #define idx(x, y) ((y) * nlx + (x))
585 void Noise::gradientMap2D(
587 float step_x, float step_y,
590 float v00, v01, v10, v11, u, v, orig_u;
591 int index, i, j, x0, y0, noisex, noisey;
594 bool eased = np.flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED);
595 Interp2dFxn interpolate = eased ?
596 biLinearInterpolation : biLinearInterpolationNoEase;
604 //calculate noise point lattice
605 nlx = (int)(u + sx * step_x) + 2;
606 nly = (int)(v + sy * step_y) + 2;
608 for (j = 0; j != nly; j++)
609 for (i = 0; i != nlx; i++)
610 noise_buf[index++] = noise2d(x0 + i, y0 + j, seed);
612 //calculate interpolations
615 for (j = 0; j != sy; j++) {
616 v00 = noise_buf[idx(0, noisey)];
617 v10 = noise_buf[idx(1, noisey)];
618 v01 = noise_buf[idx(0, noisey + 1)];
619 v11 = noise_buf[idx(1, noisey + 1)];
623 for (i = 0; i != sx; i++) {
624 gradient_buf[index++] = interpolate(v00, v10, v01, v11, u, v);
632 v10 = noise_buf[idx(noisex + 1, noisey)];
633 v11 = noise_buf[idx(noisex + 1, noisey + 1)];
647 #define idx(x, y, z) ((z) * nly * nlx + (y) * nlx + (x))
648 void Noise::gradientMap3D(
649 float x, float y, float z,
650 float step_x, float step_y, float step_z,
653 float v000, v010, v100, v110;
654 float v001, v011, v101, v111;
655 float u, v, w, orig_u, orig_v;
656 int index, i, j, k, x0, y0, z0, noisex, noisey, noisez;
659 Interp3dFxn interpolate = (np.flags & NOISE_FLAG_EASED) ?
660 triLinearInterpolation : triLinearInterpolationNoEase;
671 //calculate noise point lattice
672 nlx = (int)(u + sx * step_x) + 2;
673 nly = (int)(v + sy * step_y) + 2;
674 nlz = (int)(w + sz * step_z) + 2;
676 for (k = 0; k != nlz; k++)
677 for (j = 0; j != nly; j++)
678 for (i = 0; i != nlx; i++)
679 noise_buf[index++] = noise3d(x0 + i, y0 + j, z0 + k, seed);
681 //calculate interpolations
685 for (k = 0; k != sz; k++) {
688 for (j = 0; j != sy; j++) {
689 v000 = noise_buf[idx(0, noisey, noisez)];
690 v100 = noise_buf[idx(1, noisey, noisez)];
691 v010 = noise_buf[idx(0, noisey + 1, noisez)];
692 v110 = noise_buf[idx(1, noisey + 1, noisez)];
693 v001 = noise_buf[idx(0, noisey, noisez + 1)];
694 v101 = noise_buf[idx(1, noisey, noisez + 1)];
695 v011 = noise_buf[idx(0, noisey + 1, noisez + 1)];
696 v111 = noise_buf[idx(1, noisey + 1, noisez + 1)];
700 for (i = 0; i != sx; i++) {
701 gradient_buf[index++] = interpolate(
702 v000, v100, v010, v110,
703 v001, v101, v011, v111,
712 v100 = noise_buf[idx(noisex + 1, noisey, noisez)];
713 v110 = noise_buf[idx(noisex + 1, noisey + 1, noisez)];
716 v101 = noise_buf[idx(noisex + 1, noisey, noisez + 1)];
717 v111 = noise_buf[idx(noisex + 1, noisey + 1, noisez + 1)];
738 float *Noise::perlinMap2D(float x, float y, float *persistence_map)
740 float f = 1.0, g = 1.0;
741 size_t bufsize = sx * sy;
746 memset(result, 0, sizeof(float) * bufsize);
748 if (persistence_map) {
750 persist_buf = new float[bufsize];
751 for (size_t i = 0; i != bufsize; i++)
752 persist_buf[i] = 1.0;
755 for (size_t oct = 0; oct < np.octaves; oct++) {
756 gradientMap2D(x * f, y * f,
757 f / np.spread.X, f / np.spread.Y,
758 seed + np.seed + oct);
760 updateResults(g, persist_buf, persistence_map, bufsize);
766 if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
767 for (size_t i = 0; i != bufsize; i++)
768 result[i] = result[i] * np.scale + np.offset;
775 float *Noise::perlinMap3D(float x, float y, float z, float *persistence_map)
777 float f = 1.0, g = 1.0;
778 size_t bufsize = sx * sy * sz;
784 memset(result, 0, sizeof(float) * bufsize);
786 if (persistence_map) {
788 persist_buf = new float[bufsize];
789 for (size_t i = 0; i != bufsize; i++)
790 persist_buf[i] = 1.0;
793 for (size_t oct = 0; oct < np.octaves; oct++) {
794 gradientMap3D(x * f, y * f, z * f,
795 f / np.spread.X, f / np.spread.Y, f / np.spread.Z,
796 seed + np.seed + oct);
798 updateResults(g, persist_buf, persistence_map, bufsize);
804 if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
805 for (size_t i = 0; i != bufsize; i++)
806 result[i] = result[i] * np.scale + np.offset;
813 void Noise::updateResults(float g, float *gmap,
814 float *persistence_map, size_t bufsize)
816 // This looks very ugly, but it is 50-70% faster than having
817 // conditional statements inside the loop
818 if (np.flags & NOISE_FLAG_ABSVALUE) {
819 if (persistence_map) {
820 for (size_t i = 0; i != bufsize; i++) {
821 result[i] += gmap[i] * fabs(gradient_buf[i]);
822 gmap[i] *= persistence_map[i];
825 for (size_t i = 0; i != bufsize; i++)
826 result[i] += g * fabs(gradient_buf[i]);
829 if (persistence_map) {
830 for (size_t i = 0; i != bufsize; i++) {
831 result[i] += gmap[i] * gradient_buf[i];
832 gmap[i] *= persistence_map[i];
835 for (size_t i = 0; i != bufsize; i++)
836 result[i] += g * gradient_buf[i];