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 "exceptions.h"
34 #define NOISE_MAGIC_X 1619
35 #define NOISE_MAGIC_Y 31337
36 #define NOISE_MAGIC_Z 52591
37 #define NOISE_MAGIC_SEED 1013
39 typedef float (*Interp3dFxn)(
40 float v000, float v100, float v010, float v110,
41 float v001, float v101, float v011, float v111,
42 float x, float y, float z);
44 float cos_lookup[16] = {
45 1.0, 0.9238, 0.7071, 0.3826, 0, -0.3826, -0.7071, -0.9238,
46 1.0, -0.9238, -0.7071, -0.3826, 0, 0.3826, 0.7071, 0.9238
50 ///////////////////////////////////////////////////////////////////////////////
53 //noise poly: p(n) = 60493n^3 + 19990303n + 137612589
54 float noise2d(int x, int y, int seed)
56 int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y
57 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
59 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
60 return 1.f - (float)n / 0x40000000;
64 float noise3d(int x, int y, int z, int seed)
66 int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y + NOISE_MAGIC_Z * z
67 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
69 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
70 return 1.f - (float)n / 0x40000000;
74 float dotProduct(float vx, float vy, float wx, float wy)
76 return vx * wx + vy * wy;
80 inline float linearInterpolation(float v0, float v1, float t)
82 return v0 + (v1 - v0) * t;
86 float biLinearInterpolation(
91 float tx = easeCurve(x);
92 float ty = easeCurve(y);
93 float u = linearInterpolation(v00, v10, tx);
94 float v = linearInterpolation(v01, v11, tx);
95 return linearInterpolation(u, v, ty);
99 float biLinearInterpolationNoEase(
100 float x0y0, float x1y0,
101 float x0y1, float x1y1,
104 float u = linearInterpolation(x0y0, x1y0, x);
105 float v = linearInterpolation(x0y1, x1y1, x);
106 return linearInterpolation(u, v, y);
110 float triLinearInterpolation(
111 float v000, float v100, float v010, float v110,
112 float v001, float v101, float v011, float v111,
113 float x, float y, float z)
115 float u = biLinearInterpolation(v000, v100, v010, v110, x, y);
116 float v = biLinearInterpolation(v001, v101, v011, v111, x, y);
117 return linearInterpolation(u, v, z);
121 float triLinearInterpolationNoEase(
122 float v000, float v100, float v010, float v110,
123 float v001, float v101, float v011, float v111,
124 float x, float y, float z)
126 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, x, y);
127 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, x, y);
128 return linearInterpolation(u, v, z);
133 float triLinearInterpolation(
134 float v000, float v100, float v010, float v110,
135 float v001, float v101, float v011, float v111,
136 float x, float y, float z)
138 float tx = easeCurve(x);
139 float ty = easeCurve(y);
140 float tz = easeCurve(z);
143 v000 * (1 - tx) * (1 - ty) * (1 - tz) +
144 v100 * tx * (1 - ty) * (1 - tz) +
145 v010 * (1 - tx) * ty * (1 - tz) +
146 v110 * tx * ty * (1 - tz) +
147 v001 * (1 - tx) * (1 - ty) * tz +
148 v101 * tx * (1 - ty) * tz +
149 v011 * (1 - tx) * ty * tz +
154 float triLinearInterpolationNoEase(
155 float v000, float v100, float v010, float v110,
156 float v001, float v101, float v011, float v111,
157 float x, float y, float z)
163 v000 * (1 - tx) * (1 - ty) * (1 - tz) +
164 v100 * tx * (1 - ty) * (1 - tz) +
165 v010 * (1 - tx) * ty * (1 - tz) +
166 v110 * tx * ty * (1 - tz) +
167 v001 * (1 - tx) * (1 - ty) * tz +
168 v101 * tx * (1 - ty) * tz +
169 v011 * (1 - tx) * ty * tz +
176 float noise2d_gradient(float x, float y, int seed)
178 // Calculate the integer coordinates
179 int x0 = (x > 0.0 ? (int)x : (int)x - 1);
180 int y0 = (y > 0.0 ? (int)y : (int)y - 1);
181 // Calculate the remaining part of the coordinates
182 float xl = x - (float)x0;
183 float yl = y - (float)y0;
184 // Calculate random cosine lookup table indices for the integer corners.
185 // They are looked up as unit vector gradients from the lookup table.
186 int n00 = (int)((noise2d(x0, y0, seed)+1)*8);
187 int n10 = (int)((noise2d(x0+1, y0, seed)+1)*8);
188 int n01 = (int)((noise2d(x0, y0+1, seed)+1)*8);
189 int n11 = (int)((noise2d(x0+1, y0+1, seed)+1)*8);
190 // Make a dot product for the gradients and the positions, to get the values
191 float s = dotProduct(cos_lookup[n00], cos_lookup[(n00+12)%16], xl, yl);
192 float u = dotProduct(-cos_lookup[n10], cos_lookup[(n10+12)%16], 1.-xl, yl);
193 float v = dotProduct(cos_lookup[n01], -cos_lookup[(n01+12)%16], xl, 1.-yl);
194 float w = dotProduct(-cos_lookup[n11], -cos_lookup[(n11+12)%16], 1.-xl, 1.-yl);
195 // Interpolate between the values
196 return biLinearInterpolation(s,u,v,w,xl,yl);
201 float noise2d_gradient(float x, float y, int seed)
203 // Calculate the integer coordinates
206 // Calculate the remaining part of the coordinates
207 float xl = x - (float)x0;
208 float yl = y - (float)y0;
209 // Get values for corners of square
210 float v00 = noise2d(x0, y0, seed);
211 float v10 = noise2d(x0+1, y0, seed);
212 float v01 = noise2d(x0, y0+1, seed);
213 float v11 = noise2d(x0+1, y0+1, seed);
215 return biLinearInterpolation(v00, v10, v01, v11, xl, yl);
219 float noise3d_gradient(float x, float y, float z, int seed, bool eased)
221 // Calculate the integer coordinates
225 // Calculate the remaining part of the coordinates
226 float xl = x - (float)x0;
227 float yl = y - (float)y0;
228 float zl = z - (float)z0;
229 // Get values for corners of cube
230 float v000 = noise3d(x0, y0, z0, seed);
231 float v100 = noise3d(x0 + 1, y0, z0, seed);
232 float v010 = noise3d(x0, y0 + 1, z0, seed);
233 float v110 = noise3d(x0 + 1, y0 + 1, z0, seed);
234 float v001 = noise3d(x0, y0, z0 + 1, seed);
235 float v101 = noise3d(x0 + 1, y0, z0 + 1, seed);
236 float v011 = noise3d(x0, y0 + 1, z0 + 1, seed);
237 float v111 = noise3d(x0 + 1, y0 + 1, z0 + 1, seed);
240 return triLinearInterpolation(
241 v000, v100, v010, v110,
242 v001, v101, v011, v111,
245 return triLinearInterpolationNoEase(
246 v000, v100, v010, v110,
247 v001, v101, v011, v111,
253 float noise2d_perlin(float x, float y, int seed,
254 int octaves, float persistence)
259 for (int i = 0; i < octaves; i++)
261 a += g * noise2d_gradient(x * f, y * f, seed + i);
269 float noise2d_perlin_abs(float x, float y, int seed,
270 int octaves, float persistence)
275 for (int i = 0; i < octaves; i++)
277 a += g * fabs(noise2d_gradient(x * f, y * f, seed + i));
285 float noise3d_perlin(float x, float y, float z, int seed,
286 int octaves, float persistence, bool eased)
291 for (int i = 0; i < octaves; i++)
293 a += g * noise3d_gradient(x * f, y * f, z * f, seed + i, eased);
301 float noise3d_perlin_abs(float x, float y, float z, int seed,
302 int octaves, float persistence, bool eased)
307 for (int i = 0; i < octaves; i++)
309 a += g * fabs(noise3d_gradient(x * f, y * f, z * f, seed + i, eased));
317 float contour(float v)
326 ///////////////////////// [ New perlin stuff ] ////////////////////////////
329 Noise::Noise(NoiseParams *np, int seed, int sx, int sy, int sz)
337 this->noisebuf = NULL;
338 resizeNoiseBuf(sz > 1);
341 this->buf = new float[sx * sy * sz];
342 this->result = new float[sx * sy * sz];
343 } catch (std::bad_alloc &e) {
344 throw InvalidNoiseParamsException();
357 void Noise::setSize(int sx, int sy, int sz)
363 this->noisebuf = NULL;
364 resizeNoiseBuf(sz > 1);
369 this->buf = new float[sx * sy * sz];
370 this->result = new float[sx * sy * sz];
371 } catch (std::bad_alloc &e) {
372 throw InvalidNoiseParamsException();
377 void Noise::setSpreadFactor(v3f spread)
379 this->np->spread = spread;
381 resizeNoiseBuf(sz > 1);
385 void Noise::setOctaves(int octaves)
387 this->np->octaves = octaves;
389 resizeNoiseBuf(sz > 1);
393 void Noise::resizeNoiseBuf(bool is3d)
398 //maximum possible spread value factor
399 ofactor = (float)(1 << (np->octaves - 1));
401 //noise lattice point count
402 //(int)(sz * spread * ofactor) is # of lattice points crossed due to length
403 // + 2 for the two initial endpoints
404 // + 1 for potentially crossing a boundary due to offset
405 nlx = (int)(sx * ofactor / np->spread.X) + 3;
406 nly = (int)(sy * ofactor / np->spread.Y) + 3;
407 nlz = is3d ? (int)(sz * ofactor / np->spread.Z) + 3 : 1;
412 noisebuf = new float[nlx * nly * nlz];
413 } catch (std::bad_alloc &e) {
414 throw InvalidNoiseParamsException();
420 * NB: This algorithm is not optimal in terms of space complexity. The entire
421 * integer lattice of noise points could be done as 2 lines instead, and for 3D,
422 * 2 lines + 2 planes.
423 * However, this would require the noise calls to be interposed with the
424 * interpolation loops, which may trash the icache, leading to lower overall
426 * Another optimization that could save half as many noise calls is to carry over
427 * values from the previous noise lattice as midpoints in the new lattice for the
430 #define idx(x, y) ((y) * nlx + (x))
431 void Noise::gradientMap2D(
433 float step_x, float step_y,
436 float v00, v01, v10, v11, u, v, orig_u;
437 int index, i, j, x0, y0, noisex, noisey;
446 //calculate noise point lattice
447 nlx = (int)(u + sx * step_x) + 2;
448 nly = (int)(v + sy * step_y) + 2;
450 for (j = 0; j != nly; j++)
451 for (i = 0; i != nlx; i++)
452 noisebuf[index++] = noise2d(x0 + i, y0 + j, seed);
454 //calculate interpolations
457 for (j = 0; j != sy; j++) {
458 v00 = noisebuf[idx(0, noisey)];
459 v10 = noisebuf[idx(1, noisey)];
460 v01 = noisebuf[idx(0, noisey + 1)];
461 v11 = noisebuf[idx(1, noisey + 1)];
465 for (i = 0; i != sx; i++) {
466 buf[index++] = biLinearInterpolation(v00, v10, v01, v11, u, v);
473 v10 = noisebuf[idx(noisex + 1, noisey)];
474 v11 = noisebuf[idx(noisex + 1, noisey + 1)];
488 #define idx(x, y, z) ((z) * nly * nlx + (y) * nlx + (x))
489 void Noise::gradientMap3D(
490 float x, float y, float z,
491 float step_x, float step_y, float step_z,
492 int seed, bool eased)
494 float v000, v010, v100, v110;
495 float v001, v011, v101, v111;
496 float u, v, w, orig_u, orig_v;
497 int index, i, j, k, x0, y0, z0, noisex, noisey, noisez;
500 Interp3dFxn interpolate = eased ?
501 triLinearInterpolation : triLinearInterpolationNoEase;
512 //calculate noise point lattice
513 nlx = (int)(u + sx * step_x) + 2;
514 nly = (int)(v + sy * step_y) + 2;
515 nlz = (int)(w + sz * step_z) + 2;
517 for (k = 0; k != nlz; k++)
518 for (j = 0; j != nly; j++)
519 for (i = 0; i != nlx; i++)
520 noisebuf[index++] = noise3d(x0 + i, y0 + j, z0 + k, seed);
522 //calculate interpolations
526 for (k = 0; k != sz; k++) {
529 for (j = 0; j != sy; j++) {
530 v000 = noisebuf[idx(0, noisey, noisez)];
531 v100 = noisebuf[idx(1, noisey, noisez)];
532 v010 = noisebuf[idx(0, noisey + 1, noisez)];
533 v110 = noisebuf[idx(1, noisey + 1, noisez)];
534 v001 = noisebuf[idx(0, noisey, noisez + 1)];
535 v101 = noisebuf[idx(1, noisey, noisez + 1)];
536 v011 = noisebuf[idx(0, noisey + 1, noisez + 1)];
537 v111 = noisebuf[idx(1, noisey + 1, noisez + 1)];
541 for (i = 0; i != sx; i++) {
542 buf[index++] = interpolate(
543 v000, v100, v010, v110,
544 v001, v101, v011, v111,
553 v100 = noisebuf[idx(noisex + 1, noisey, noisez)];
554 v110 = noisebuf[idx(noisex + 1, noisey + 1, noisez)];
557 v101 = noisebuf[idx(noisex + 1, noisey, noisez + 1)];
558 v111 = noisebuf[idx(noisex + 1, noisey + 1, noisez + 1)];
579 float *Noise::perlinMap2D(float x, float y)
581 float f = 1.0, g = 1.0;
582 size_t bufsize = sx * sy;
587 memset(result, 0, sizeof(float) * bufsize);
589 for (int oct = 0; oct < np->octaves; oct++) {
590 gradientMap2D(x * f, y * f,
591 f / np->spread.X, f / np->spread.Y,
592 seed + np->seed + oct);
594 for (size_t i = 0; i != bufsize; i++)
595 result[i] += g * buf[i];
605 float *Noise::perlinMap2DModulated(float x, float y, float *persist_map)
608 size_t bufsize = sx * sy;
613 memset(result, 0, sizeof(float) * bufsize);
615 float *g = new float[bufsize];
616 for (size_t i = 0; i != bufsize; i++)
619 for (int oct = 0; oct < np->octaves; oct++) {
620 gradientMap2D(x * f, y * f,
621 f / np->spread.X, f / np->spread.Y,
622 seed + np->seed + oct);
624 for (size_t i = 0; i != bufsize; i++) {
625 result[i] += g[i] * buf[i];
626 g[i] *= persist_map[i];
637 float *Noise::perlinMap3D(float x, float y, float z, bool eased)
639 float f = 1.0, g = 1.0;
640 size_t bufsize = sx * sy * sz;
646 memset(result, 0, sizeof(float) * bufsize);
648 for (int oct = 0; oct < np->octaves; oct++) {
649 gradientMap3D(x * f, y * f, z * f,
650 f / np->spread.X, f / np->spread.Y, f / np->spread.Z,
651 seed + np->seed + oct, eased);
653 for (size_t i = 0; i != bufsize; i++)
654 result[i] += g * buf[i];
664 void Noise::transformNoiseMap()
668 for (int z = 0; z != sz; z++)
669 for (int y = 0; y != sy; y++)
670 for (int x = 0; x != sx; x++) {
671 result[i] = result[i] * np->scale + np->offset;