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"
33 #define NOISE_MAGIC_X 1619
34 #define NOISE_MAGIC_Y 31337
35 #define NOISE_MAGIC_Z 52591
36 #define NOISE_MAGIC_SEED 1013
38 typedef float (*Interp3dFxn)(
39 float v000, float v100, float v010, float v110,
40 float v001, float v101, float v011, float v111,
41 float x, float y, float z);
43 float cos_lookup[16] = {
44 1.0, 0.9238, 0.7071, 0.3826, 0, -0.3826, -0.7071, -0.9238,
45 1.0, -0.9238, -0.7071, -0.3826, 0, 0.3826, 0.7071, 0.9238
49 ///////////////////////////////////////////////////////////////////////////////
52 //noise poly: p(n) = 60493n^3 + 19990303n + 137612589
53 float noise2d(int x, int y, int seed)
55 int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y
56 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
58 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
59 return 1.f - (float)n / 0x40000000;
63 float noise3d(int x, int y, int z, int seed)
65 int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y + NOISE_MAGIC_Z * z
66 + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
68 n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
69 return 1.f - (float)n / 0x40000000;
73 float dotProduct(float vx, float vy, float wx, float wy)
75 return vx * wx + vy * wy;
79 inline float linearInterpolation(float v0, float v1, float t)
81 return v0 + (v1 - v0) * t;
85 float biLinearInterpolation(
90 float tx = easeCurve(x);
91 float ty = easeCurve(y);
92 float u = linearInterpolation(v00, v10, tx);
93 float v = linearInterpolation(v01, v11, tx);
94 return linearInterpolation(u, v, ty);
98 float biLinearInterpolationNoEase(
99 float x0y0, float x1y0,
100 float x0y1, float x1y1,
103 float u = linearInterpolation(x0y0, x1y0, x);
104 float v = linearInterpolation(x0y1, x1y1, x);
105 return linearInterpolation(u, v, y);
109 float triLinearInterpolation(
110 float v000, float v100, float v010, float v110,
111 float v001, float v101, float v011, float v111,
112 float x, float y, float z)
114 float u = biLinearInterpolation(v000, v100, v010, v110, x, y);
115 float v = biLinearInterpolation(v001, v101, v011, v111, x, y);
116 return linearInterpolation(u, v, z);
120 float triLinearInterpolationNoEase(
121 float v000, float v100, float v010, float v110,
122 float v001, float v101, float v011, float v111,
123 float x, float y, float z)
125 float u = biLinearInterpolationNoEase(v000, v100, v010, v110, x, y);
126 float v = biLinearInterpolationNoEase(v001, v101, v011, v111, x, y);
127 return linearInterpolation(u, v, z);
132 float triLinearInterpolation(
133 float v000, float v100, float v010, float v110,
134 float v001, float v101, float v011, float v111,
135 float x, float y, float z)
137 float tx = easeCurve(x);
138 float ty = easeCurve(y);
139 float tz = easeCurve(z);
142 v000 * (1 - tx) * (1 - ty) * (1 - tz) +
143 v100 * tx * (1 - ty) * (1 - tz) +
144 v010 * (1 - tx) * ty * (1 - tz) +
145 v110 * tx * ty * (1 - tz) +
146 v001 * (1 - tx) * (1 - ty) * tz +
147 v101 * tx * (1 - ty) * tz +
148 v011 * (1 - tx) * ty * tz +
153 float triLinearInterpolationNoEase(
154 float v000, float v100, float v010, float v110,
155 float v001, float v101, float v011, float v111,
156 float x, float y, float z)
162 v000 * (1 - tx) * (1 - ty) * (1 - tz) +
163 v100 * tx * (1 - ty) * (1 - tz) +
164 v010 * (1 - tx) * ty * (1 - tz) +
165 v110 * tx * ty * (1 - tz) +
166 v001 * (1 - tx) * (1 - ty) * tz +
167 v101 * tx * (1 - ty) * tz +
168 v011 * (1 - tx) * ty * tz +
175 float noise2d_gradient(float x, float y, int seed)
177 // Calculate the integer coordinates
178 int x0 = (x > 0.0 ? (int)x : (int)x - 1);
179 int y0 = (y > 0.0 ? (int)y : (int)y - 1);
180 // Calculate the remaining part of the coordinates
181 float xl = x - (float)x0;
182 float yl = y - (float)y0;
183 // Calculate random cosine lookup table indices for the integer corners.
184 // They are looked up as unit vector gradients from the lookup table.
185 int n00 = (int)((noise2d(x0, y0, seed)+1)*8);
186 int n10 = (int)((noise2d(x0+1, y0, seed)+1)*8);
187 int n01 = (int)((noise2d(x0, y0+1, seed)+1)*8);
188 int n11 = (int)((noise2d(x0+1, y0+1, seed)+1)*8);
189 // Make a dot product for the gradients and the positions, to get the values
190 float s = dotProduct(cos_lookup[n00], cos_lookup[(n00+12)%16], xl, yl);
191 float u = dotProduct(-cos_lookup[n10], cos_lookup[(n10+12)%16], 1.-xl, yl);
192 float v = dotProduct(cos_lookup[n01], -cos_lookup[(n01+12)%16], xl, 1.-yl);
193 float w = dotProduct(-cos_lookup[n11], -cos_lookup[(n11+12)%16], 1.-xl, 1.-yl);
194 // Interpolate between the values
195 return biLinearInterpolation(s,u,v,w,xl,yl);
200 float noise2d_gradient(float x, float y, int seed)
202 // Calculate the integer coordinates
205 // Calculate the remaining part of the coordinates
206 float xl = x - (float)x0;
207 float yl = y - (float)y0;
208 // Get values for corners of square
209 float v00 = noise2d(x0, y0, seed);
210 float v10 = noise2d(x0+1, y0, seed);
211 float v01 = noise2d(x0, y0+1, seed);
212 float v11 = noise2d(x0+1, y0+1, seed);
214 return biLinearInterpolation(v00, v10, v01, v11, xl, yl);
218 float noise3d_gradient(float x, float y, float z, int seed)
220 // Calculate the integer coordinates
224 // Calculate the remaining part of the coordinates
225 float xl = x - (float)x0;
226 float yl = y - (float)y0;
227 float zl = z - (float)z0;
228 // Get values for corners of cube
229 float v000 = noise3d(x0, y0, z0, seed);
230 float v100 = noise3d(x0 + 1, y0, z0, seed);
231 float v010 = noise3d(x0, y0 + 1, z0, seed);
232 float v110 = noise3d(x0 + 1, y0 + 1, z0, seed);
233 float v001 = noise3d(x0, y0, z0 + 1, seed);
234 float v101 = noise3d(x0 + 1, y0, z0 + 1, seed);
235 float v011 = noise3d(x0, y0 + 1, z0 + 1, seed);
236 float v111 = noise3d(x0 + 1, y0 + 1, z0 + 1, seed);
238 return triLinearInterpolationNoEase(
239 v000, v100, v010, v110,
240 v001, v101, v011, v111,
245 float noise2d_perlin(float x, float y, int seed,
246 int octaves, float persistence)
251 for (int i = 0; i < octaves; i++)
253 a += g * noise2d_gradient(x * f, y * f, seed + i);
261 float noise2d_perlin_abs(float x, float y, int seed,
262 int octaves, float persistence)
267 for (int i = 0; i < octaves; i++)
269 a += g * fabs(noise2d_gradient(x * f, y * f, seed + i));
277 float noise3d_perlin(float x, float y, float z, int seed,
278 int octaves, float persistence)
283 for (int i = 0; i < octaves; i++)
285 a += g * noise3d_gradient(x * f, y * f, z * f, seed + i);
293 float noise3d_perlin_abs(float x, float y, float z, int seed,
294 int octaves, float persistence)
299 for (int i = 0; i < octaves; i++)
301 a += g * fabs(noise3d_gradient(x * f, y * f, z * f, seed + i));
309 float contour(float v)
318 ///////////////////////// [ New perlin stuff ] ////////////////////////////
321 Noise::Noise(NoiseParams *np, int seed, int sx, int sy, int sz)
329 this->noisebuf = NULL;
330 resizeNoiseBuf(sz > 1);
332 this->buf = new float[sx * sy * sz];
333 this->result = new float[sx * sy * sz];
345 void Noise::setSize(int sx, int sy, int sz)
351 this->noisebuf = NULL;
352 resizeNoiseBuf(sz > 1);
356 this->buf = new float[sx * sy * sz];
357 this->result = new float[sx * sy * sz];
361 void Noise::setSpreadFactor(v3f spread)
363 this->np->spread = spread;
365 resizeNoiseBuf(sz > 1);
369 void Noise::setOctaves(int octaves)
371 this->np->octaves = octaves;
373 resizeNoiseBuf(sz > 1);
377 void Noise::resizeNoiseBuf(bool is3d)
382 //maximum possible spread value factor
383 ofactor = (float)(1 << (np->octaves - 1));
385 //noise lattice point count
386 //(int)(sz * spread * ofactor) is # of lattice points crossed due to length
387 // + 2 for the two initial endpoints
388 // + 1 for potentially crossing a boundary due to offset
389 nlx = (int)(sx * ofactor / np->spread.X) + 3;
390 nly = (int)(sy * ofactor / np->spread.Y) + 3;
391 nlz = is3d ? (int)(sz * ofactor / np->spread.Z) + 3 : 1;
395 noisebuf = new float[nlx * nly * nlz];
400 * NB: This algorithm is not optimal in terms of space complexity. The entire
401 * integer lattice of noise points could be done as 2 lines instead, and for 3D,
402 * 2 lines + 2 planes.
403 * However, this would require the noise calls to be interposed with the
404 * interpolation loops, which may trash the icache, leading to lower overall
406 * Another optimization that could save half as many noise calls is to carry over
407 * values from the previous noise lattice as midpoints in the new lattice for the
410 #define idx(x, y) ((y) * nlx + (x))
411 void Noise::gradientMap2D(
413 float step_x, float step_y,
416 float v00, v01, v10, v11, u, v, orig_u;
417 int index, i, j, x0, y0, noisex, noisey;
426 //calculate noise point lattice
427 nlx = (int)(u + sx * step_x) + 2;
428 nly = (int)(v + sy * step_y) + 2;
430 for (j = 0; j != nly; j++)
431 for (i = 0; i != nlx; i++)
432 noisebuf[index++] = noise2d(x0 + i, y0 + j, seed);
434 //calculate interpolations
437 for (j = 0; j != sy; j++) {
438 v00 = noisebuf[idx(0, noisey)];
439 v10 = noisebuf[idx(1, noisey)];
440 v01 = noisebuf[idx(0, noisey + 1)];
441 v11 = noisebuf[idx(1, noisey + 1)];
445 for (i = 0; i != sx; i++) {
446 buf[index++] = biLinearInterpolation(v00, v10, v01, v11, u, v);
453 v10 = noisebuf[idx(noisex + 1, noisey)];
454 v11 = noisebuf[idx(noisex + 1, noisey + 1)];
468 #define idx(x, y, z) ((z) * nly * nlx + (y) * nlx + (x))
469 void Noise::gradientMap3D(
470 float x, float y, float z,
471 float step_x, float step_y, float step_z,
472 int seed, bool eased)
474 float v000, v010, v100, v110;
475 float v001, v011, v101, v111;
476 float u, v, w, orig_u, orig_v;
477 int index, i, j, k, x0, y0, z0, noisex, noisey, noisez;
480 Interp3dFxn interpolate = eased ?
481 triLinearInterpolation : triLinearInterpolationNoEase;
492 //calculate noise point lattice
493 nlx = (int)(u + sx * step_x) + 2;
494 nly = (int)(v + sy * step_y) + 2;
495 nlz = (int)(w + sz * step_z) + 2;
497 for (k = 0; k != nlz; k++)
498 for (j = 0; j != nly; j++)
499 for (i = 0; i != nlx; i++)
500 noisebuf[index++] = noise3d(x0 + i, y0 + j, z0 + k, seed);
502 //calculate interpolations
506 for (k = 0; k != sz; k++) {
509 for (j = 0; j != sy; j++) {
510 v000 = noisebuf[idx(0, noisey, noisez)];
511 v100 = noisebuf[idx(1, noisey, noisez)];
512 v010 = noisebuf[idx(0, noisey + 1, noisez)];
513 v110 = noisebuf[idx(1, noisey + 1, noisez)];
514 v001 = noisebuf[idx(0, noisey, noisez + 1)];
515 v101 = noisebuf[idx(1, noisey, noisez + 1)];
516 v011 = noisebuf[idx(0, noisey + 1, noisez + 1)];
517 v111 = noisebuf[idx(1, noisey + 1, noisez + 1)];
521 for (i = 0; i != sx; i++) {
522 buf[index++] = interpolate(
523 v000, v100, v010, v110,
524 v001, v101, v011, v111,
533 v100 = noisebuf[idx(noisex + 1, noisey, noisez)];
534 v110 = noisebuf[idx(noisex + 1, noisey + 1, noisez)];
537 v101 = noisebuf[idx(noisex + 1, noisey, noisez + 1)];
538 v111 = noisebuf[idx(noisex + 1, noisey + 1, noisez + 1)];
559 float *Noise::perlinMap2D(float x, float y)
561 float f = 1.0, g = 1.0;
562 size_t bufsize = sx * sy;
567 memset(result, 0, sizeof(float) * bufsize);
569 for (int oct = 0; oct < np->octaves; oct++) {
570 gradientMap2D(x * f, y * f,
571 f / np->spread.X, f / np->spread.Y,
572 seed + np->seed + oct);
574 for (size_t i = 0; i != bufsize; i++)
575 result[i] += g * buf[i];
585 float *Noise::perlinMap2DModulated(float x, float y, float *persist_map)
588 size_t bufsize = sx * sy;
593 memset(result, 0, sizeof(float) * bufsize);
595 float *g = new float[bufsize];
596 for (size_t i = 0; i != bufsize; i++)
599 for (int oct = 0; oct < np->octaves; oct++) {
600 gradientMap2D(x * f, y * f,
601 f / np->spread.X, f / np->spread.Y,
602 seed + np->seed + oct);
604 for (size_t i = 0; i != bufsize; i++) {
605 result[i] += g[i] * buf[i];
606 g[i] *= persist_map[i];
617 float *Noise::perlinMap3D(float x, float y, float z, bool eased)
619 float f = 1.0, g = 1.0;
620 size_t bufsize = sx * sy * sz;
626 memset(result, 0, sizeof(float) * bufsize);
628 for (int oct = 0; oct < np->octaves; oct++) {
629 gradientMap3D(x * f, y * f, z * f,
630 f / np->spread.X, f / np->spread.Y, f / np->spread.Z,
631 seed + np->seed + oct, eased);
633 for (size_t i = 0; i != bufsize; i++)
634 result[i] += g * buf[i];
644 void Noise::transformNoiseMap()
648 for (int z = 0; z != sz; z++)
649 for (int y = 0; y != sy; y++)
650 for (int x = 0; x != sx; x++) {
651 result[i] = result[i] * np->scale + np->offset;