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 ///////////////////////////////////////////////////////////////////////////////
66 float noise2d(int x, int y, int seed)
68 int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y
69 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
71 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
72 return 1.f - (float)n / 0x40000000;
76 float noise3d(int x, int y, int z, int seed)
78 int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y + NOISE_MAGIC_Z * z
79 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
81 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
82 return 1.f - (float)n / 0x40000000;
86 inline float dotProduct(float vx, float vy, float wx, float wy)
88 return vx * wx + vy * wy;
92 inline float linearInterpolation(float v0, float v1, float t)
94 return v0 + (v1 - v0) * t;
98 inline float biLinearInterpolation(
100 float v01, float v11,
103 float tx = easeCurve(x);
104 float ty = easeCurve(y);
107 v00 * (1 - tx) * (1 - ty) +
108 v10 * tx * (1 - ty) +
109 v01 * (1 - tx) * ty +
113 float u = linearInterpolation(v00, v10, tx);
114 float v = linearInterpolation(v01, v11, tx);
115 return linearInterpolation(u, v, ty);
119 inline float biLinearInterpolationNoEase(
120 float v00, float v10,
121 float v01, float v11,
124 float u = linearInterpolation(v00, v10, x);
125 float v = linearInterpolation(v01, v11, x);
126 return linearInterpolation(u, v, y);
130 float triLinearInterpolation(
131 float v000, float v100, float v010, float v110,
132 float v001, float v101, float v011, float v111,
133 float x, float y, float z)
135 float tx = easeCurve(x);
136 float ty = easeCurve(y);
137 float tz = easeCurve(z);
140 v000 * (1 - tx) * (1 - ty) * (1 - tz) +
141 v100 * tx * (1 - ty) * (1 - tz) +
142 v010 * (1 - tx) * ty * (1 - tz) +
143 v110 * tx * ty * (1 - tz) +
144 v001 * (1 - tx) * (1 - ty) * tz +
145 v101 * tx * (1 - ty) * tz +
146 v011 * (1 - tx) * ty * tz +
150 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, tx, ty);
151 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, tx, ty);
152 return linearInterpolation(u, v, tz);
155 float triLinearInterpolationNoEase(
156 float v000, float v100, float v010, float v110,
157 float v001, float v101, float v011, float v111,
158 float x, float y, float z)
160 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, x, y);
161 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, x, y);
162 return linearInterpolation(u, v, z);
167 float noise2d_gradient(float x, float y, int seed)
169 // Calculate the integer coordinates
170 int x0 = (x > 0.0 ? (int)x : (int)x - 1);
171 int y0 = (y > 0.0 ? (int)y : (int)y - 1);
172 // Calculate the remaining part of the coordinates
173 float xl = x - (float)x0;
174 float yl = y - (float)y0;
175 // Calculate random cosine lookup table indices for the integer corners.
176 // They are looked up as unit vector gradients from the lookup table.
177 int n00 = (int)((noise2d(x0, y0, seed)+1)*8);
178 int n10 = (int)((noise2d(x0+1, y0, seed)+1)*8);
179 int n01 = (int)((noise2d(x0, y0+1, seed)+1)*8);
180 int n11 = (int)((noise2d(x0+1, y0+1, seed)+1)*8);
181 // Make a dot product for the gradients and the positions, to get the values
182 float s = dotProduct(cos_lookup[n00], cos_lookup[(n00+12)%16], xl, yl);
183 float u = dotProduct(-cos_lookup[n10], cos_lookup[(n10+12)%16], 1.-xl, yl);
184 float v = dotProduct(cos_lookup[n01], -cos_lookup[(n01+12)%16], xl, 1.-yl);
185 float w = dotProduct(-cos_lookup[n11], -cos_lookup[(n11+12)%16], 1.-xl, 1.-yl);
186 // Interpolate between the values
187 return biLinearInterpolation(s,u,v,w,xl,yl);
192 float noise2d_gradient(float x, float y, int seed, bool eased)
194 // Calculate the integer coordinates
197 // Calculate the remaining part of the coordinates
198 float xl = x - (float)x0;
199 float yl = y - (float)y0;
200 // Get values for corners of square
201 float v00 = noise2d(x0, y0, seed);
202 float v10 = noise2d(x0+1, y0, seed);
203 float v01 = noise2d(x0, y0+1, seed);
204 float v11 = noise2d(x0+1, y0+1, seed);
207 return biLinearInterpolation(v00, v10, v01, v11, xl, yl);
209 return biLinearInterpolationNoEase(v00, v10, v01, v11, xl, yl);
213 float noise3d_gradient(float x, float y, float z, int seed, bool eased)
215 // Calculate the integer coordinates
219 // Calculate the remaining part of the coordinates
220 float xl = x - (float)x0;
221 float yl = y - (float)y0;
222 float zl = z - (float)z0;
223 // Get values for corners of cube
224 float v000 = noise3d(x0, y0, z0, seed);
225 float v100 = noise3d(x0 + 1, y0, z0, seed);
226 float v010 = noise3d(x0, y0 + 1, z0, seed);
227 float v110 = noise3d(x0 + 1, y0 + 1, z0, seed);
228 float v001 = noise3d(x0, y0, z0 + 1, seed);
229 float v101 = noise3d(x0 + 1, y0, z0 + 1, seed);
230 float v011 = noise3d(x0, y0 + 1, z0 + 1, seed);
231 float v111 = noise3d(x0 + 1, y0 + 1, z0 + 1, seed);
234 return triLinearInterpolation(
235 v000, v100, v010, v110,
236 v001, v101, v011, v111,
239 return triLinearInterpolationNoEase(
240 v000, v100, v010, v110,
241 v001, v101, v011, v111,
247 float noise2d_perlin(float x, float y, int seed,
248 int octaves, float persistence, bool eased)
253 for (int i = 0; i < octaves; i++)
255 a += g * noise2d_gradient(x * f, y * f, seed + i, eased);
263 float noise2d_perlin_abs(float x, float y, int seed,
264 int octaves, float persistence, bool eased)
269 for (int i = 0; i < octaves; i++) {
270 a += g * fabs(noise2d_gradient(x * f, y * f, seed + i, eased));
278 float noise3d_perlin(float x, float y, float z, int seed,
279 int octaves, float persistence, bool eased)
284 for (int i = 0; i < octaves; i++) {
285 a += g * noise3d_gradient(x * f, y * f, z * f, seed + i, eased);
293 float noise3d_perlin_abs(float x, float y, float z, int seed,
294 int octaves, float persistence, bool eased)
299 for (int i = 0; i < octaves; i++) {
300 a += g * fabs(noise3d_gradient(x * f, y * f, z * f, seed + i, eased));
308 float contour(float v)
317 ///////////////////////// [ New noise ] ////////////////////////////
320 float NoisePerlin2D(NoiseParams *np, float x, float y, int seed)
330 for (size_t i = 0; i < np->octaves; i++) {
331 float noiseval = noise2d_gradient(x * f, y * f, seed + i,
332 np->flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED));
334 if (np->flags & NOISE_FLAG_ABSVALUE)
335 noiseval = fabs(noiseval);
342 return np->offset + a * np->scale;
346 float NoisePerlin3D(NoiseParams *np, float x, float y, float z, int seed)
357 for (size_t i = 0; i < np->octaves; i++) {
358 float noiseval = noise3d_gradient(x * f, y * f, z * f, seed + i,
359 np->flags & NOISE_FLAG_EASED);
361 if (np->flags & NOISE_FLAG_ABSVALUE)
362 noiseval = fabs(noiseval);
369 return np->offset + a * np->scale;
373 Noise::Noise(NoiseParams *np_, int seed, int sx, int sy, int sz)
375 memcpy(&np, np_, sizeof(np));
381 this->persist_buf = NULL;
382 this->gradient_buf = NULL;
391 delete[] gradient_buf;
392 delete[] persist_buf;
398 void Noise::allocBuffers()
400 this->noise_buf = NULL;
401 resizeNoiseBuf(sz > 1);
403 delete[] gradient_buf;
404 delete[] persist_buf;
408 size_t bufsize = sx * sy * sz;
409 this->persist_buf = NULL;
410 this->gradient_buf = new float[bufsize];
411 this->result = new float[bufsize];
412 } catch (std::bad_alloc &e) {
413 throw InvalidNoiseParamsException();
418 void Noise::setSize(int sx, int sy, int sz)
428 void Noise::setSpreadFactor(v3f spread)
430 this->np.spread = spread;
432 resizeNoiseBuf(sz > 1);
436 void Noise::setOctaves(int octaves)
438 this->np.octaves = octaves;
440 resizeNoiseBuf(sz > 1);
444 void Noise::resizeNoiseBuf(bool is3d)
449 //maximum possible spread value factor
450 ofactor = pow(np.lacunarity, np.octaves - 1);
452 //noise lattice point count
453 //(int)(sz * spread * ofactor) is # of lattice points crossed due to length
454 // + 2 for the two initial endpoints
455 // + 1 for potentially crossing a boundary due to offset
456 nlx = (int)ceil(sx * ofactor / np.spread.X) + 3;
457 nly = (int)ceil(sy * ofactor / np.spread.Y) + 3;
458 nlz = is3d ? (int)ceil(sz * ofactor / np.spread.Z) + 3 : 1;
462 noise_buf = new float[nlx * nly * nlz];
463 } catch (std::bad_alloc &e) {
464 throw InvalidNoiseParamsException();
470 * NB: This algorithm is not optimal in terms of space complexity. The entire
471 * integer lattice of noise points could be done as 2 lines instead, and for 3D,
472 * 2 lines + 2 planes.
473 * However, this would require the noise calls to be interposed with the
474 * interpolation loops, which may trash the icache, leading to lower overall
476 * Another optimization that could save half as many noise calls is to carry over
477 * values from the previous noise lattice as midpoints in the new lattice for the
480 #define idx(x, y) ((y) * nlx + (x))
481 void Noise::gradientMap2D(
483 float step_x, float step_y,
486 float v00, v01, v10, v11, u, v, orig_u;
487 int index, i, j, x0, y0, noisex, noisey;
490 bool eased = np.flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED);
491 Interp2dFxn interpolate = eased ?
492 biLinearInterpolation : biLinearInterpolationNoEase;
500 //calculate noise point lattice
501 nlx = (int)(u + sx * step_x) + 2;
502 nly = (int)(v + sy * step_y) + 2;
504 for (j = 0; j != nly; j++)
505 for (i = 0; i != nlx; i++)
506 noise_buf[index++] = noise2d(x0 + i, y0 + j, seed);
508 //calculate interpolations
511 for (j = 0; j != sy; j++) {
512 v00 = noise_buf[idx(0, noisey)];
513 v10 = noise_buf[idx(1, noisey)];
514 v01 = noise_buf[idx(0, noisey + 1)];
515 v11 = noise_buf[idx(1, noisey + 1)];
519 for (i = 0; i != sx; i++) {
520 gradient_buf[index++] = interpolate(v00, v10, v01, v11, u, v);
528 v10 = noise_buf[idx(noisex + 1, noisey)];
529 v11 = noise_buf[idx(noisex + 1, noisey + 1)];
543 #define idx(x, y, z) ((z) * nly * nlx + (y) * nlx + (x))
544 void Noise::gradientMap3D(
545 float x, float y, float z,
546 float step_x, float step_y, float step_z,
549 float v000, v010, v100, v110;
550 float v001, v011, v101, v111;
551 float u, v, w, orig_u, orig_v;
552 int index, i, j, k, x0, y0, z0, noisex, noisey, noisez;
555 Interp3dFxn interpolate = (np.flags & NOISE_FLAG_EASED) ?
556 triLinearInterpolation : triLinearInterpolationNoEase;
567 //calculate noise point lattice
568 nlx = (int)(u + sx * step_x) + 2;
569 nly = (int)(v + sy * step_y) + 2;
570 nlz = (int)(w + sz * step_z) + 2;
572 for (k = 0; k != nlz; k++)
573 for (j = 0; j != nly; j++)
574 for (i = 0; i != nlx; i++)
575 noise_buf[index++] = noise3d(x0 + i, y0 + j, z0 + k, seed);
577 //calculate interpolations
581 for (k = 0; k != sz; k++) {
584 for (j = 0; j != sy; j++) {
585 v000 = noise_buf[idx(0, noisey, noisez)];
586 v100 = noise_buf[idx(1, noisey, noisez)];
587 v010 = noise_buf[idx(0, noisey + 1, noisez)];
588 v110 = noise_buf[idx(1, noisey + 1, noisez)];
589 v001 = noise_buf[idx(0, noisey, noisez + 1)];
590 v101 = noise_buf[idx(1, noisey, noisez + 1)];
591 v011 = noise_buf[idx(0, noisey + 1, noisez + 1)];
592 v111 = noise_buf[idx(1, noisey + 1, noisez + 1)];
596 for (i = 0; i != sx; i++) {
597 gradient_buf[index++] = interpolate(
598 v000, v100, v010, v110,
599 v001, v101, v011, v111,
608 v100 = noise_buf[idx(noisex + 1, noisey, noisez)];
609 v110 = noise_buf[idx(noisex + 1, noisey + 1, noisez)];
612 v101 = noise_buf[idx(noisex + 1, noisey, noisez + 1)];
613 v111 = noise_buf[idx(noisex + 1, noisey + 1, noisez + 1)];
634 float *Noise::perlinMap2D(float x, float y, float *persistence_map)
636 float f = 1.0, g = 1.0;
637 size_t bufsize = sx * sy;
642 memset(result, 0, sizeof(float) * bufsize);
644 if (persistence_map) {
646 persist_buf = new float[bufsize];
647 for (size_t i = 0; i != bufsize; i++)
648 persist_buf[i] = 1.0;
651 for (size_t oct = 0; oct < np.octaves; oct++) {
652 gradientMap2D(x * f, y * f,
653 f / np.spread.X, f / np.spread.Y,
654 seed + np.seed + oct);
656 updateResults(g, persist_buf, persistence_map, bufsize);
662 if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
663 for (size_t i = 0; i != bufsize; i++)
664 result[i] = result[i] * np.scale + np.offset;
671 float *Noise::perlinMap3D(float x, float y, float z, float *persistence_map)
673 float f = 1.0, g = 1.0;
674 size_t bufsize = sx * sy * sz;
680 memset(result, 0, sizeof(float) * bufsize);
682 if (persistence_map) {
684 persist_buf = new float[bufsize];
685 for (size_t i = 0; i != bufsize; i++)
686 persist_buf[i] = 1.0;
689 for (size_t oct = 0; oct < np.octaves; oct++) {
690 gradientMap3D(x * f, y * f, z * f,
691 f / np.spread.X, f / np.spread.Y, f / np.spread.Z,
692 seed + np.seed + oct);
694 updateResults(g, persist_buf, persistence_map, bufsize);
700 if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
701 for (size_t i = 0; i != bufsize; i++)
702 result[i] = result[i] * np.scale + np.offset;
709 void Noise::updateResults(float g, float *gmap,
710 float *persistence_map, size_t bufsize)
712 // This looks very ugly, but it is 50-70% faster than having
713 // conditional statements inside the loop
714 if (np.flags & NOISE_FLAG_ABSVALUE) {
715 if (persistence_map) {
716 for (size_t i = 0; i != bufsize; i++) {
717 result[i] += gmap[i] * fabs(gradient_buf[i]);
718 gmap[i] *= persistence_map[i];
721 for (size_t i = 0; i != bufsize; i++)
722 result[i] += g * fabs(gradient_buf[i]);
725 if (persistence_map) {
726 for (size_t i = 0; i != bufsize; i++) {
727 result[i] += gmap[i] * gradient_buf[i];
728 gmap[i] *= persistence_map[i];
731 for (size_t i = 0; i != bufsize; i++)
732 result[i] += g * gradient_buf[i];