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 perlin stuff ] ////////////////////////////
320 Noise::Noise(NoiseParams *np_, int seed, int sx, int sy, int sz)
322 memcpy(&np, np_, sizeof(np));
328 this->persist_buf = NULL;
329 this->gradient_buf = NULL;
332 if (np.flags & NOISE_FLAG_DEFAULTS) {
333 // By default, only 2d noise is eased.
335 np.flags |= NOISE_FLAG_EASED;
344 delete[] gradient_buf;
345 delete[] persist_buf;
351 void Noise::allocBuffers()
353 this->noise_buf = NULL;
354 resizeNoiseBuf(sz > 1);
356 delete[] gradient_buf;
357 delete[] persist_buf;
361 size_t bufsize = sx * sy * sz;
362 this->persist_buf = NULL;
363 this->gradient_buf = new float[bufsize];
364 this->result = new float[bufsize];
365 } catch (std::bad_alloc &e) {
366 throw InvalidNoiseParamsException();
371 void Noise::setSize(int sx, int sy, int sz)
381 void Noise::setSpreadFactor(v3f spread)
383 this->np.spread = spread;
385 resizeNoiseBuf(sz > 1);
389 void Noise::setOctaves(int octaves)
391 this->np.octaves = octaves;
393 resizeNoiseBuf(sz > 1);
397 void Noise::resizeNoiseBuf(bool is3d)
402 //maximum possible spread value factor
403 ofactor = pow(np.lacunarity, np.octaves - 1);
405 //noise lattice point count
406 //(int)(sz * spread * ofactor) is # of lattice points crossed due to length
407 // + 2 for the two initial endpoints
408 // + 1 for potentially crossing a boundary due to offset
409 nlx = (int)ceil(sx * ofactor / np.spread.X) + 3;
410 nly = (int)ceil(sy * ofactor / np.spread.Y) + 3;
411 nlz = is3d ? (int)ceil(sz * ofactor / np.spread.Z) + 3 : 1;
415 noise_buf = new float[nlx * nly * nlz];
416 } catch (std::bad_alloc &e) {
417 throw InvalidNoiseParamsException();
423 * NB: This algorithm is not optimal in terms of space complexity. The entire
424 * integer lattice of noise points could be done as 2 lines instead, and for 3D,
425 * 2 lines + 2 planes.
426 * However, this would require the noise calls to be interposed with the
427 * interpolation loops, which may trash the icache, leading to lower overall
429 * Another optimization that could save half as many noise calls is to carry over
430 * values from the previous noise lattice as midpoints in the new lattice for the
433 #define idx(x, y) ((y) * nlx + (x))
434 void Noise::gradientMap2D(
436 float step_x, float step_y,
439 float v00, v01, v10, v11, u, v, orig_u;
440 int index, i, j, x0, y0, noisex, noisey;
443 Interp2dFxn interpolate = (np.flags & NOISE_FLAG_EASED) ?
444 biLinearInterpolation : biLinearInterpolationNoEase;
452 //calculate noise point lattice
453 nlx = (int)(u + sx * step_x) + 2;
454 nly = (int)(v + sy * step_y) + 2;
456 for (j = 0; j != nly; j++)
457 for (i = 0; i != nlx; i++)
458 noise_buf[index++] = noise2d(x0 + i, y0 + j, seed);
460 //calculate interpolations
463 for (j = 0; j != sy; j++) {
464 v00 = noise_buf[idx(0, noisey)];
465 v10 = noise_buf[idx(1, noisey)];
466 v01 = noise_buf[idx(0, noisey + 1)];
467 v11 = noise_buf[idx(1, noisey + 1)];
471 for (i = 0; i != sx; i++) {
472 gradient_buf[index++] = interpolate(v00, v10, v01, v11, u, v);
480 v10 = noise_buf[idx(noisex + 1, noisey)];
481 v11 = noise_buf[idx(noisex + 1, noisey + 1)];
495 #define idx(x, y, z) ((z) * nly * nlx + (y) * nlx + (x))
496 void Noise::gradientMap3D(
497 float x, float y, float z,
498 float step_x, float step_y, float step_z,
501 float v000, v010, v100, v110;
502 float v001, v011, v101, v111;
503 float u, v, w, orig_u, orig_v;
504 int index, i, j, k, x0, y0, z0, noisex, noisey, noisez;
507 Interp3dFxn interpolate = (np.flags & NOISE_FLAG_EASED) ?
508 triLinearInterpolation : triLinearInterpolationNoEase;
519 //calculate noise point lattice
520 nlx = (int)(u + sx * step_x) + 2;
521 nly = (int)(v + sy * step_y) + 2;
522 nlz = (int)(w + sz * step_z) + 2;
524 for (k = 0; k != nlz; k++)
525 for (j = 0; j != nly; j++)
526 for (i = 0; i != nlx; i++)
527 noise_buf[index++] = noise3d(x0 + i, y0 + j, z0 + k, seed);
529 //calculate interpolations
533 for (k = 0; k != sz; k++) {
536 for (j = 0; j != sy; j++) {
537 v000 = noise_buf[idx(0, noisey, noisez)];
538 v100 = noise_buf[idx(1, noisey, noisez)];
539 v010 = noise_buf[idx(0, noisey + 1, noisez)];
540 v110 = noise_buf[idx(1, noisey + 1, noisez)];
541 v001 = noise_buf[idx(0, noisey, noisez + 1)];
542 v101 = noise_buf[idx(1, noisey, noisez + 1)];
543 v011 = noise_buf[idx(0, noisey + 1, noisez + 1)];
544 v111 = noise_buf[idx(1, noisey + 1, noisez + 1)];
548 for (i = 0; i != sx; i++) {
549 gradient_buf[index++] = interpolate(
550 v000, v100, v010, v110,
551 v001, v101, v011, v111,
560 v100 = noise_buf[idx(noisex + 1, noisey, noisez)];
561 v110 = noise_buf[idx(noisex + 1, noisey + 1, noisez)];
564 v101 = noise_buf[idx(noisex + 1, noisey, noisez + 1)];
565 v111 = noise_buf[idx(noisex + 1, noisey + 1, noisez + 1)];
586 float *Noise::perlinMap2D(float x, float y, float *persistence_map)
588 float f = 1.0, g = 1.0;
589 size_t bufsize = sx * sy;
594 memset(result, 0, sizeof(float) * bufsize);
596 if (persistence_map) {
598 persist_buf = new float[bufsize];
599 for (size_t i = 0; i != bufsize; i++)
600 persist_buf[i] = 1.0;
603 for (size_t oct = 0; oct < np.octaves; oct++) {
604 gradientMap2D(x * f, y * f,
605 f / np.spread.X, f / np.spread.Y,
606 seed + np.seed + oct);
608 updateResults(g, persist_buf, persistence_map, bufsize);
614 if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
615 for (size_t i = 0; i != bufsize; i++)
616 result[i] = result[i] * np.scale + np.offset;
623 float *Noise::perlinMap3D(float x, float y, float z, float *persistence_map)
625 float f = 1.0, g = 1.0;
626 size_t bufsize = sx * sy * sz;
632 memset(result, 0, sizeof(float) * bufsize);
634 if (persistence_map) {
636 persist_buf = new float[bufsize];
637 for (size_t i = 0; i != bufsize; i++)
638 persist_buf[i] = 1.0;
641 for (size_t oct = 0; oct < np.octaves; oct++) {
642 gradientMap3D(x * f, y * f, z * f,
643 f / np.spread.X, f / np.spread.Y, f / np.spread.Z,
644 seed + np.seed + oct);
646 updateResults(g, persist_buf, persistence_map, bufsize);
652 if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
653 for (size_t i = 0; i != bufsize; i++)
654 result[i] = result[i] * np.scale + np.offset;
661 void Noise::updateResults(float g, float *gmap,
662 float *persistence_map, size_t bufsize)
664 // This looks very ugly, but it is 50-70% faster than having
665 // conditional statements inside the loop
666 if (np.flags & NOISE_FLAG_ABSVALUE) {
667 if (persistence_map) {
668 for (size_t i = 0; i != bufsize; i++) {
669 result[i] += gmap[i] * fabs(gradient_buf[i]);
670 gmap[i] *= persistence_map[i];
673 for (size_t i = 0; i != bufsize; i++)
674 result[i] += g * fabs(gradient_buf[i]);
677 if (persistence_map) {
678 for (size_t i = 0; i != bufsize; i++) {
679 result[i] += gmap[i] * gradient_buf[i];
680 gmap[i] *= persistence_map[i];
683 for (size_t i = 0; i != bufsize; i++)
684 result[i] += g * gradient_buf[i];