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)
127 u8 *outb = (u8 *)out;
129 size_t len_alignment = (uintptr_t)out % sizeof(u32);
131 len -= len_alignment;
133 while (len_alignment--) {
140 size_t len_dwords = len / sizeof(u32);
141 while (len_dwords--) {
143 *(u32 *)outb = next();
147 size_t len_remaining = len % sizeof(u32);
150 while (len_remaining--) {
159 s32 PcgRandom::randNormalDist(s32 min, s32 max, int num_trials)
162 for (int i = 0; i != num_trials; i++)
163 accum += range(min, max);
164 return ((float)accum / num_trials) + 0.5f;
167 ///////////////////////////////////////////////////////////////////////////////
169 float noise2d(int x, int y, int seed)
171 int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y
172 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
174 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
175 return 1.f - (float)n / 0x40000000;
179 float noise3d(int x, int y, int z, int seed)
181 int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y + NOISE_MAGIC_Z * z
182 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
184 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
185 return 1.f - (float)n / 0x40000000;
189 inline float dotProduct(float vx, float vy, float wx, float wy)
191 return vx * wx + vy * wy;
195 inline float linearInterpolation(float v0, float v1, float t)
197 return v0 + (v1 - v0) * t;
201 inline float biLinearInterpolation(
202 float v00, float v10,
203 float v01, float v11,
206 float tx = easeCurve(x);
207 float ty = easeCurve(y);
210 v00 * (1 - tx) * (1 - ty) +
211 v10 * tx * (1 - ty) +
212 v01 * (1 - tx) * ty +
216 float u = linearInterpolation(v00, v10, tx);
217 float v = linearInterpolation(v01, v11, tx);
218 return linearInterpolation(u, v, ty);
222 inline float biLinearInterpolationNoEase(
223 float v00, float v10,
224 float v01, float v11,
227 float u = linearInterpolation(v00, v10, x);
228 float v = linearInterpolation(v01, v11, x);
229 return linearInterpolation(u, v, y);
233 float triLinearInterpolation(
234 float v000, float v100, float v010, float v110,
235 float v001, float v101, float v011, float v111,
236 float x, float y, float z)
238 float tx = easeCurve(x);
239 float ty = easeCurve(y);
240 float tz = easeCurve(z);
243 v000 * (1 - tx) * (1 - ty) * (1 - tz) +
244 v100 * tx * (1 - ty) * (1 - tz) +
245 v010 * (1 - tx) * ty * (1 - tz) +
246 v110 * tx * ty * (1 - tz) +
247 v001 * (1 - tx) * (1 - ty) * tz +
248 v101 * tx * (1 - ty) * tz +
249 v011 * (1 - tx) * ty * tz +
253 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, tx, ty);
254 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, tx, ty);
255 return linearInterpolation(u, v, tz);
258 float triLinearInterpolationNoEase(
259 float v000, float v100, float v010, float v110,
260 float v001, float v101, float v011, float v111,
261 float x, float y, float z)
263 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, x, y);
264 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, x, y);
265 return linearInterpolation(u, v, z);
270 float noise2d_gradient(float x, float y, int seed)
272 // Calculate the integer coordinates
273 int x0 = (x > 0.0 ? (int)x : (int)x - 1);
274 int y0 = (y > 0.0 ? (int)y : (int)y - 1);
275 // Calculate the remaining part of the coordinates
276 float xl = x - (float)x0;
277 float yl = y - (float)y0;
278 // Calculate random cosine lookup table indices for the integer corners.
279 // They are looked up as unit vector gradients from the lookup table.
280 int n00 = (int)((noise2d(x0, y0, seed)+1)*8);
281 int n10 = (int)((noise2d(x0+1, y0, seed)+1)*8);
282 int n01 = (int)((noise2d(x0, y0+1, seed)+1)*8);
283 int n11 = (int)((noise2d(x0+1, y0+1, seed)+1)*8);
284 // Make a dot product for the gradients and the positions, to get the values
285 float s = dotProduct(cos_lookup[n00], cos_lookup[(n00+12)%16], xl, yl);
286 float u = dotProduct(-cos_lookup[n10], cos_lookup[(n10+12)%16], 1.-xl, yl);
287 float v = dotProduct(cos_lookup[n01], -cos_lookup[(n01+12)%16], xl, 1.-yl);
288 float w = dotProduct(-cos_lookup[n11], -cos_lookup[(n11+12)%16], 1.-xl, 1.-yl);
289 // Interpolate between the values
290 return biLinearInterpolation(s,u,v,w,xl,yl);
295 float noise2d_gradient(float x, float y, int seed, bool eased)
297 // Calculate the integer coordinates
300 // Calculate the remaining part of the coordinates
301 float xl = x - (float)x0;
302 float yl = y - (float)y0;
303 // Get values for corners of square
304 float v00 = noise2d(x0, y0, seed);
305 float v10 = noise2d(x0+1, y0, seed);
306 float v01 = noise2d(x0, y0+1, seed);
307 float v11 = noise2d(x0+1, y0+1, seed);
310 return biLinearInterpolation(v00, v10, v01, v11, xl, yl);
312 return biLinearInterpolationNoEase(v00, v10, v01, v11, xl, yl);
316 float noise3d_gradient(float x, float y, float z, int seed, bool eased)
318 // Calculate the integer coordinates
322 // Calculate the remaining part of the coordinates
323 float xl = x - (float)x0;
324 float yl = y - (float)y0;
325 float zl = z - (float)z0;
326 // Get values for corners of cube
327 float v000 = noise3d(x0, y0, z0, seed);
328 float v100 = noise3d(x0 + 1, y0, z0, seed);
329 float v010 = noise3d(x0, y0 + 1, z0, seed);
330 float v110 = noise3d(x0 + 1, y0 + 1, z0, seed);
331 float v001 = noise3d(x0, y0, z0 + 1, seed);
332 float v101 = noise3d(x0 + 1, y0, z0 + 1, seed);
333 float v011 = noise3d(x0, y0 + 1, z0 + 1, seed);
334 float v111 = noise3d(x0 + 1, y0 + 1, z0 + 1, seed);
337 return triLinearInterpolation(
338 v000, v100, v010, v110,
339 v001, v101, v011, v111,
342 return triLinearInterpolationNoEase(
343 v000, v100, v010, v110,
344 v001, v101, v011, v111,
350 float noise2d_perlin(float x, float y, int seed,
351 int octaves, float persistence, bool eased)
356 for (int i = 0; i < octaves; i++)
358 a += g * noise2d_gradient(x * f, y * f, seed + i, eased);
366 float noise2d_perlin_abs(float x, float y, int seed,
367 int octaves, float persistence, bool eased)
372 for (int i = 0; i < octaves; i++) {
373 a += g * fabs(noise2d_gradient(x * f, y * f, seed + i, eased));
381 float noise3d_perlin(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 * noise3d_gradient(x * f, y * f, z * f, seed + i, eased);
396 float noise3d_perlin_abs(float x, float y, float z, int seed,
397 int octaves, float persistence, bool eased)
402 for (int i = 0; i < octaves; i++) {
403 a += g * fabs(noise3d_gradient(x * f, y * f, z * f, seed + i, eased));
411 float contour(float v)
420 ///////////////////////// [ New noise ] ////////////////////////////
423 float NoisePerlin2D(NoiseParams *np, float x, float y, int seed)
433 for (size_t i = 0; i < np->octaves; i++) {
434 float noiseval = noise2d_gradient(x * f, y * f, seed + i,
435 np->flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED));
437 if (np->flags & NOISE_FLAG_ABSVALUE)
438 noiseval = fabs(noiseval);
445 return np->offset + a * np->scale;
449 float NoisePerlin3D(NoiseParams *np, float x, float y, float z, int seed)
460 for (size_t i = 0; i < np->octaves; i++) {
461 float noiseval = noise3d_gradient(x * f, y * f, z * f, seed + i,
462 np->flags & NOISE_FLAG_EASED);
464 if (np->flags & NOISE_FLAG_ABSVALUE)
465 noiseval = fabs(noiseval);
472 return np->offset + a * np->scale;
476 Noise::Noise(NoiseParams *np_, int seed, int sx, int sy, int sz)
478 memcpy(&np, np_, sizeof(np));
484 this->persist_buf = NULL;
485 this->gradient_buf = NULL;
494 delete[] gradient_buf;
495 delete[] persist_buf;
501 void Noise::allocBuffers()
503 this->noise_buf = NULL;
504 resizeNoiseBuf(sz > 1);
506 delete[] gradient_buf;
507 delete[] persist_buf;
511 size_t bufsize = sx * sy * sz;
512 this->persist_buf = NULL;
513 this->gradient_buf = new float[bufsize];
514 this->result = new float[bufsize];
515 } catch (std::bad_alloc &e) {
516 throw InvalidNoiseParamsException();
521 void Noise::setSize(int sx, int sy, int sz)
531 void Noise::setSpreadFactor(v3f spread)
533 this->np.spread = spread;
535 resizeNoiseBuf(sz > 1);
539 void Noise::setOctaves(int octaves)
541 this->np.octaves = octaves;
543 resizeNoiseBuf(sz > 1);
547 void Noise::resizeNoiseBuf(bool is3d)
552 //maximum possible spread value factor
553 ofactor = pow(np.lacunarity, np.octaves - 1);
555 //noise lattice point count
556 //(int)(sz * spread * ofactor) is # of lattice points crossed due to length
557 // + 2 for the two initial endpoints
558 // + 1 for potentially crossing a boundary due to offset
559 nlx = (int)ceil(sx * ofactor / np.spread.X) + 3;
560 nly = (int)ceil(sy * ofactor / np.spread.Y) + 3;
561 nlz = is3d ? (int)ceil(sz * ofactor / np.spread.Z) + 3 : 1;
565 noise_buf = new float[nlx * nly * nlz];
566 } catch (std::bad_alloc &e) {
567 throw InvalidNoiseParamsException();
573 * NB: This algorithm is not optimal in terms of space complexity. The entire
574 * integer lattice of noise points could be done as 2 lines instead, and for 3D,
575 * 2 lines + 2 planes.
576 * However, this would require the noise calls to be interposed with the
577 * interpolation loops, which may trash the icache, leading to lower overall
579 * Another optimization that could save half as many noise calls is to carry over
580 * values from the previous noise lattice as midpoints in the new lattice for the
583 #define idx(x, y) ((y) * nlx + (x))
584 void Noise::gradientMap2D(
586 float step_x, float step_y,
589 float v00, v01, v10, v11, u, v, orig_u;
590 int index, i, j, x0, y0, noisex, noisey;
593 bool eased = np.flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED);
594 Interp2dFxn interpolate = eased ?
595 biLinearInterpolation : biLinearInterpolationNoEase;
603 //calculate noise point lattice
604 nlx = (int)(u + sx * step_x) + 2;
605 nly = (int)(v + sy * step_y) + 2;
607 for (j = 0; j != nly; j++)
608 for (i = 0; i != nlx; i++)
609 noise_buf[index++] = noise2d(x0 + i, y0 + j, seed);
611 //calculate interpolations
614 for (j = 0; j != sy; j++) {
615 v00 = noise_buf[idx(0, noisey)];
616 v10 = noise_buf[idx(1, noisey)];
617 v01 = noise_buf[idx(0, noisey + 1)];
618 v11 = noise_buf[idx(1, noisey + 1)];
622 for (i = 0; i != sx; i++) {
623 gradient_buf[index++] = interpolate(v00, v10, v01, v11, u, v);
631 v10 = noise_buf[idx(noisex + 1, noisey)];
632 v11 = noise_buf[idx(noisex + 1, noisey + 1)];
646 #define idx(x, y, z) ((z) * nly * nlx + (y) * nlx + (x))
647 void Noise::gradientMap3D(
648 float x, float y, float z,
649 float step_x, float step_y, float step_z,
652 float v000, v010, v100, v110;
653 float v001, v011, v101, v111;
654 float u, v, w, orig_u, orig_v;
655 int index, i, j, k, x0, y0, z0, noisex, noisey, noisez;
658 Interp3dFxn interpolate = (np.flags & NOISE_FLAG_EASED) ?
659 triLinearInterpolation : triLinearInterpolationNoEase;
670 //calculate noise point lattice
671 nlx = (int)(u + sx * step_x) + 2;
672 nly = (int)(v + sy * step_y) + 2;
673 nlz = (int)(w + sz * step_z) + 2;
675 for (k = 0; k != nlz; k++)
676 for (j = 0; j != nly; j++)
677 for (i = 0; i != nlx; i++)
678 noise_buf[index++] = noise3d(x0 + i, y0 + j, z0 + k, seed);
680 //calculate interpolations
684 for (k = 0; k != sz; k++) {
687 for (j = 0; j != sy; j++) {
688 v000 = noise_buf[idx(0, noisey, noisez)];
689 v100 = noise_buf[idx(1, noisey, noisez)];
690 v010 = noise_buf[idx(0, noisey + 1, noisez)];
691 v110 = noise_buf[idx(1, noisey + 1, noisez)];
692 v001 = noise_buf[idx(0, noisey, noisez + 1)];
693 v101 = noise_buf[idx(1, noisey, noisez + 1)];
694 v011 = noise_buf[idx(0, noisey + 1, noisez + 1)];
695 v111 = noise_buf[idx(1, noisey + 1, noisez + 1)];
699 for (i = 0; i != sx; i++) {
700 gradient_buf[index++] = interpolate(
701 v000, v100, v010, v110,
702 v001, v101, v011, v111,
711 v100 = noise_buf[idx(noisex + 1, noisey, noisez)];
712 v110 = noise_buf[idx(noisex + 1, noisey + 1, noisez)];
715 v101 = noise_buf[idx(noisex + 1, noisey, noisez + 1)];
716 v111 = noise_buf[idx(noisex + 1, noisey + 1, noisez + 1)];
737 float *Noise::perlinMap2D(float x, float y, float *persistence_map)
739 float f = 1.0, g = 1.0;
740 size_t bufsize = sx * sy;
745 memset(result, 0, sizeof(float) * bufsize);
747 if (persistence_map) {
749 persist_buf = new float[bufsize];
750 for (size_t i = 0; i != bufsize; i++)
751 persist_buf[i] = 1.0;
754 for (size_t oct = 0; oct < np.octaves; oct++) {
755 gradientMap2D(x * f, y * f,
756 f / np.spread.X, f / np.spread.Y,
757 seed + np.seed + oct);
759 updateResults(g, persist_buf, persistence_map, bufsize);
765 if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
766 for (size_t i = 0; i != bufsize; i++)
767 result[i] = result[i] * np.scale + np.offset;
774 float *Noise::perlinMap3D(float x, float y, float z, float *persistence_map)
776 float f = 1.0, g = 1.0;
777 size_t bufsize = sx * sy * sz;
783 memset(result, 0, sizeof(float) * bufsize);
785 if (persistence_map) {
787 persist_buf = new float[bufsize];
788 for (size_t i = 0; i != bufsize; i++)
789 persist_buf[i] = 1.0;
792 for (size_t oct = 0; oct < np.octaves; oct++) {
793 gradientMap3D(x * f, y * f, z * f,
794 f / np.spread.X, f / np.spread.Y, f / np.spread.Z,
795 seed + np.seed + oct);
797 updateResults(g, persist_buf, persistence_map, bufsize);
803 if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
804 for (size_t i = 0; i != bufsize; i++)
805 result[i] = result[i] * np.scale + np.offset;
812 void Noise::updateResults(float g, float *gmap,
813 float *persistence_map, size_t bufsize)
815 // This looks very ugly, but it is 50-70% faster than having
816 // conditional statements inside the loop
817 if (np.flags & NOISE_FLAG_ABSVALUE) {
818 if (persistence_map) {
819 for (size_t i = 0; i != bufsize; i++) {
820 result[i] += gmap[i] * fabs(gradient_buf[i]);
821 gmap[i] *= persistence_map[i];
824 for (size_t i = 0; i != bufsize; i++)
825 result[i] += g * fabs(gradient_buf[i]);
828 if (persistence_map) {
829 for (size_t i = 0; i != bufsize; i++) {
830 result[i] += gmap[i] * gradient_buf[i];
831 gmap[i] *= persistence_map[i];
834 for (size_t i = 0; i != bufsize; i++)
835 result[i] += g * gradient_buf[i];