src/noise.cpp

   1 /*
   2 Minetest
   3 Copyright (C) 2010-2013 celeron55, Perttu Ahola <celeron55@gmail.com>
   4 Copyright (C) 2010-2013 kwolekr, Ryan Kwolek <kwolekr@minetest.net>
   5
   6 This program is free software; you can redistribute it and/or modify
   7 it under the terms of the GNU Lesser General Public License as published by
   8 the Free Software Foundation; either version 2.1 of the License, or
   9 (at your option) any later version.
  10
  11 This program is distributed in the hope that it will be useful,
  12 but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 GNU Lesser General Public License for more details.
  15
  16 You should have received a copy of the GNU Lesser General Public License along
  17 with this program; if not, write to the Free Software Foundation, Inc.,
  18 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
  19 */
  20
  21 #include <math.h>
  22 #include "noise.h"
  23 #include <iostream>
  24 #include <string.h> // memset
  25 #include "debug.h"
  26 #include "util/numeric.h"
  27
  28 #define NOISE_MAGIC_X    1619
  29 #define NOISE_MAGIC_Y    31337
  30 #define NOISE_MAGIC_Z    52591
  31 #define NOISE_MAGIC_SEED 1013
  32
  33 typedef float (*Interp3dFxn)(
  34                 float v000, float v100, float v010, float v110,
  35                 float v001, float v101, float v011, float v111,
  36                 float x, float y, float z);
  37
  38 float cos_lookup[16] = {
  39         1.0,  0.9238,  0.7071,  0.3826, 0, -0.3826, -0.7071, -0.9238,
  40         1.0, -0.9238, -0.7071, -0.3826, 0,  0.3826,  0.7071,  0.9238
  41 };
  42
  43
  44 ///////////////////////////////////////////////////////////////////////////////
  45
  46
  47 //noise poly:  p(n) = 60493n^3 + 19990303n + 137612589
  48 float noise2d(int x, int y, int seed)
  49 {
  50         int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y
  51                         + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
  52         n = (n >> 13) ^ n;
  53         n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
  54         return 1.f - (float)n / 0x40000000;
  55 }
  56
  57
  58 float noise3d(int x, int y, int z, int seed)
  59 {
  60         int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y + NOISE_MAGIC_Z * z
  61                         + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
  62         n = (n >> 13) ^ n;
  63         n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
  64         return 1.f - (float)n / 0x40000000;
  65 }
  66
  67
  68 float dotProduct(float vx, float vy, float wx, float wy)
  69 {
  70         return vx * wx + vy * wy;
  71 }
  72
  73
  74 inline float linearInterpolation(float v0, float v1, float t)
  75 {
  76     return v0 + (v1 - v0) * t;
  77 }
  78
  79
  80 float biLinearInterpolation(
  81                 float v00, float v10,
  82                 float v01, float v11,
  83                 float x, float y)
  84 {
  85     float tx = easeCurve(x);
  86     float ty = easeCurve(y);
  87     float u = linearInterpolation(v00, v10, tx);
  88     float v = linearInterpolation(v01, v11, tx);
  89     return linearInterpolation(u, v, ty);
  90 }
  91
  92
  93 float biLinearInterpolationNoEase(
  94                 float x0y0, float x1y0,
  95                 float x0y1, float x1y1,
  96                 float x, float y)
  97 {
  98     float u = linearInterpolation(x0y0, x1y0, x);
  99     float v = linearInterpolation(x0y1, x1y1, x);
 100     return linearInterpolation(u, v, y);
 101 }
 102
 103 /*
 104 float triLinearInterpolation(
 105                 float v000, float v100, float v010, float v110,
 106                 float v001, float v101, float v011, float v111,
 107                 float x, float y, float z)
 108 {
 109         float u = biLinearInterpolation(v000, v100, v010, v110, x, y);
 110         float v = biLinearInterpolation(v001, v101, v011, v111, x, y);
 111         return linearInterpolation(u, v, z);
 112 }
 113
 114
 115 float triLinearInterpolationNoEase(
 116                 float v000, float v100, float v010, float v110,
 117                 float v001, float v101, float v011, float v111,
 118                 float x, float y, float z)
 119 {
 120         float u = biLinearInterpolationNoEase(v000, v100, v010, v110, x, y);
 121         float v = biLinearInterpolationNoEase(v001, v101, v011, v111, x, y);
 122         return linearInterpolation(u, v, z);
 123 }
 124 */
 125
 126
 127 float triLinearInterpolation(
 128                 float v000, float v100, float v010, float v110,
 129                 float v001, float v101, float v011, float v111,
 130                 float x, float y, float z)
 131 {
 132         float tx = easeCurve(x);
 133         float ty = easeCurve(y);
 134         float tz = easeCurve(z);
 135
 136         return (
 137                 v000 * (1 - tx) * (1 - ty) * (1 - tz) +
 138                 v100 * tx * (1 - ty) * (1 - tz) +
 139                 v010 * (1 - tx) * ty * (1 - tz) +
 140                 v110 * tx * ty * (1 - tz) +
 141                 v001 * (1 - tx) * (1 - ty) * tz +
 142                 v101 * tx * (1 - ty) * tz +
 143                 v011 * (1 - tx) * ty * tz +
 144                 v111 * tx * ty * tz
 145         );
 146 }
 147
 148 float triLinearInterpolationNoEase(
 149                 float v000, float v100, float v010, float v110,
 150                 float v001, float v101, float v011, float v111,
 151                 float x, float y, float z)
 152 {
 153         float tx = x;
 154         float ty = y;
 155         float tz = z;
 156         return (
 157                 v000 * (1 - tx) * (1 - ty) * (1 - tz) +
 158                 v100 * tx * (1 - ty) * (1 - tz) +
 159                 v010 * (1 - tx) * ty * (1 - tz) +
 160                 v110 * tx * ty * (1 - tz) +
 161                 v001 * (1 - tx) * (1 - ty) * tz +
 162                 v101 * tx * (1 - ty) * tz +
 163                 v011 * (1 - tx) * ty * tz +
 164                 v111 * tx * ty * tz
 165         );
 166 }
 167
 168
 169 #if 0
 170 float noise2d_gradient(float x, float y, int seed)
 171 {
 172         // Calculate the integer coordinates
 173         int x0 = (x > 0.0 ? (int)x : (int)x - 1);
 174         int y0 = (y > 0.0 ? (int)y : (int)y - 1);
 175         // Calculate the remaining part of the coordinates
 176         float xl = x - (float)x0;
 177         float yl = y - (float)y0;
 178         // Calculate random cosine lookup table indices for the integer corners.
 179         // They are looked up as unit vector gradients from the lookup table.
 180         int n00 = (int)((noise2d(x0, y0, seed)+1)*8);
 181         int n10 = (int)((noise2d(x0+1, y0, seed)+1)*8);
 182         int n01 = (int)((noise2d(x0, y0+1, seed)+1)*8);
 183         int n11 = (int)((noise2d(x0+1, y0+1, seed)+1)*8);
 184         // Make a dot product for the gradients and the positions, to get the values
 185         float s = dotProduct(cos_lookup[n00], cos_lookup[(n00+12)%16], xl, yl);
 186         float u = dotProduct(-cos_lookup[n10], cos_lookup[(n10+12)%16], 1.-xl, yl);
 187         float v = dotProduct(cos_lookup[n01], -cos_lookup[(n01+12)%16], xl, 1.-yl);
 188         float w = dotProduct(-cos_lookup[n11], -cos_lookup[(n11+12)%16], 1.-xl, 1.-yl);
 189         // Interpolate between the values
 190         return biLinearInterpolation(s,u,v,w,xl,yl);
 191 }
 192 #endif
 193
 194
 195 float noise2d_gradient(float x, float y, int seed)
 196 {
 197         // Calculate the integer coordinates
 198         int x0 = myfloor(x);
 199         int y0 = myfloor(y);
 200         // Calculate the remaining part of the coordinates
 201         float xl = x - (float)x0;
 202         float yl = y - (float)y0;
 203         // Get values for corners of square
 204         float v00 = noise2d(x0, y0, seed);
 205         float v10 = noise2d(x0+1, y0, seed);
 206         float v01 = noise2d(x0, y0+1, seed);
 207         float v11 = noise2d(x0+1, y0+1, seed);
 208         // Interpolate
 209         return biLinearInterpolation(v00, v10, v01, v11, xl, yl);
 210 }
 211
 212
 213 float noise3d_gradient(float x, float y, float z, int seed)
 214 {
 215         // Calculate the integer coordinates
 216         int x0 = myfloor(x);
 217         int y0 = myfloor(y);
 218         int z0 = myfloor(z);
 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);
 232         // Interpolate
 233         return triLinearInterpolationNoEase(
 234                 v000, v100, v010, v110,
 235                 v001, v101, v011, v111,
 236                 xl, yl, zl);
 237 }
 238
 239
 240 float noise2d_perlin(float x, float y, int seed,
 241                 int octaves, float persistence)
 242 {
 243         float a = 0;
 244         float f = 1.0;
 245         float g = 1.0;
 246         for (int i = 0; i < octaves; i++)
 247         {
 248                 a += g * noise2d_gradient(x * f, y * f, seed + i);
 249                 f *= 2.0;
 250                 g *= persistence;
 251         }
 252         return a;
 253 }
 254
 255
 256 float noise2d_perlin_abs(float x, float y, int seed,
 257                 int octaves, float persistence)
 258 {
 259         float a = 0;
 260         float f = 1.0;
 261         float g = 1.0;
 262         for (int i = 0; i < octaves; i++)
 263         {
 264                 a += g * fabs(noise2d_gradient(x * f, y * f, seed + i));
 265                 f *= 2.0;
 266                 g *= persistence;
 267         }
 268         return a;
 269 }
 270
 271
 272 float noise3d_perlin(float x, float y, float z, int seed,
 273                 int octaves, float persistence)
 274 {
 275         float a = 0;
 276         float f = 1.0;
 277         float g = 1.0;
 278         for (int i = 0; i < octaves; i++)
 279         {
 280                 a += g * noise3d_gradient(x * f, y * f, z * f, seed + i);
 281                 f *= 2.0;
 282                 g *= persistence;
 283         }
 284         return a;
 285 }
 286
 287
 288 float noise3d_perlin_abs(float x, float y, float z, int seed,
 289                 int octaves, float persistence)
 290 {
 291         float a = 0;
 292         float f = 1.0;
 293         float g = 1.0;
 294         for (int i = 0; i < octaves; i++)
 295         {
 296                 a += g * fabs(noise3d_gradient(x * f, y * f, z * f, seed + i));
 297                 f *= 2.0;
 298                 g *= persistence;
 299         }
 300         return a;
 301 }
 302
 303
 304 // -1->0, 0->1, 1->0
 305 float contour(float v)
 306 {
 307         v = fabs(v);
 308         if(v >= 1.0)
 309                 return 0.0;
 310         return (1.0 - v);
 311 }
 312
 313
 314 ///////////////////////// [ New perlin stuff ] ////////////////////////////
 315
 316
 317 Noise::Noise(NoiseParams *np, int seed, int sx, int sy, int sz)
 318 {
 319         this->np   = np;
 320         this->seed = seed;
 321         this->sx   = sx;
 322         this->sy   = sy;
 323         this->sz   = sz;
 324
 325         this->noisebuf = NULL;
 326         resizeNoiseBuf(sz > 1);
 327
 328         this->buf    = new float[sx * sy * sz];
 329         this->result = new float[sx * sy * sz];
 330 }
 331
 332
 333 Noise::~Noise()
 334 {
 335         delete[] buf;
 336         delete[] result;
 337         delete[] noisebuf;
 338 }
 339
 340
 341 void Noise::setSize(int sx, int sy, int sz)
 342 {
 343         this->sx = sx;
 344         this->sy = sy;
 345         this->sz = sz;
 346
 347         this->noisebuf = NULL;
 348         resizeNoiseBuf(sz > 1);
 349
 350         delete[] buf;
 351         delete[] result;
 352         this->buf    = new float[sx * sy * sz];
 353         this->result = new float[sx * sy * sz];
 354 }
 355
 356
 357 void Noise::setSpreadFactor(v3f spread)
 358 {
 359         this->np->spread = spread;
 360
 361         resizeNoiseBuf(sz > 1);
 362 }
 363
 364
 365 void Noise::setOctaves(int octaves)
 366 {
 367         this->np->octaves = octaves;
 368
 369         resizeNoiseBuf(sz > 1);
 370 }
 371
 372
 373 void Noise::resizeNoiseBuf(bool is3d)
 374 {
 375         int nlx, nly, nlz;
 376         float ofactor;
 377
 378         //maximum possible spread value factor
 379         ofactor = (float)(1 << (np->octaves - 1));
 380
 381         //noise lattice point count
 382         //(int)(sz * spread * ofactor) is # of lattice points crossed due to length
 383         // + 2 for the two initial endpoints
 384         // + 1 for potentially crossing a boundary due to offset
 385         nlx = (int)(sx * ofactor / np->spread.X) + 3;
 386         nly = (int)(sy * ofactor / np->spread.Y) + 3;
 387         nlz = is3d ? (int)(sz * ofactor / np->spread.Z) + 3 : 1;
 388
 389         if (noisebuf)
 390                 delete[] noisebuf;
 391         noisebuf = new float[nlx * nly * nlz];
 392 }
 393
 394
 395 /*
 396  * NB:  This algorithm is not optimal in terms of space complexity.  The entire
 397  * integer lattice of noise points could be done as 2 lines instead, and for 3D,
 398  * 2 lines + 2 planes.
 399  * However, this would require the noise calls to be interposed with the
 400  * interpolation loops, which may trash the icache, leading to lower overall
 401  * performance.
 402  * Another optimization that could save half as many noise calls is to carry over
 403  * values from the previous noise lattice as midpoints in the new lattice for the
 404  * next octave.
 405  */
 406 #define idx(x, y) ((y) * nlx + (x))
 407 void Noise::gradientMap2D(
 408                 float x, float y,
 409                 float step_x, float step_y,
 410                 int seed)
 411 {
 412         float v00, v01, v10, v11, u, v, orig_u;
 413         int index, i, j, x0, y0, noisex, noisey;
 414         int nlx, nly;
 415
 416         x0 = floor(x);
 417         y0 = floor(y);
 418         u = x - (float)x0;
 419         v = y - (float)y0;
 420         orig_u = u;
 421
 422         //calculate noise point lattice
 423         nlx = (int)(u + sx * step_x) + 2;
 424         nly = (int)(v + sy * step_y) + 2;
 425         index = 0;
 426         for (j = 0; j != nly; j++)
 427                 for (i = 0; i != nlx; i++)
 428                         noisebuf[index++] = noise2d(x0 + i, y0 + j, seed);
 429
 430         //calculate interpolations
 431         index  = 0;
 432         noisey = 0;
 433         for (j = 0; j != sy; j++) {
 434                 v00 = noisebuf[idx(0, noisey)];
 435                 v10 = noisebuf[idx(1, noisey)];
 436                 v01 = noisebuf[idx(0, noisey + 1)];
 437                 v11 = noisebuf[idx(1, noisey + 1)];
 438
 439                 u = orig_u;
 440                 noisex = 0;
 441                 for (i = 0; i != sx; i++) {
 442                         buf[index++] = biLinearInterpolation(v00, v10, v01, v11, u, v);
 443                         u += step_x;
 444                         if (u >= 1.0) {
 445                                 u -= 1.0;
 446                                 noisex++;
 447                                 v00 = v10;
 448                                 v01 = v11;
 449                                 v10 = noisebuf[idx(noisex + 1, noisey)];
 450                                 v11 = noisebuf[idx(noisex + 1, noisey + 1)];
 451                         }
 452                 }
 453
 454                 v += step_y;
 455                 if (v >= 1.0) {
 456                         v -= 1.0;
 457                         noisey++;
 458                 }
 459         }
 460 }
 461 #undef idx
 462
 463
 464 #define idx(x, y, z) ((z) * nly * nlx + (y) * nlx + (x))
 465 void Noise::gradientMap3D(
 466                 float x, float y, float z,
 467                 float step_x, float step_y, float step_z,
 468                 int seed, bool eased)
 469 {
 470         float v000, v010, v100, v110;
 471         float v001, v011, v101, v111;
 472         float u, v, w, orig_u, orig_v;
 473         int index, i, j, k, x0, y0, z0, noisex, noisey, noisez;
 474         int nlx, nly, nlz;
 475
 476         Interp3dFxn interpolate = eased ?
 477                 triLinearInterpolation : triLinearInterpolationNoEase;
 478
 479         x0 = floor(x);
 480         y0 = floor(y);
 481         z0 = floor(z);
 482         u = x - (float)x0;
 483         v = y - (float)y0;
 484         w = z - (float)z0;
 485         orig_u = u;
 486         orig_v = v;
 487
 488         //calculate noise point lattice
 489         nlx = (int)(u + sx * step_x) + 2;
 490         nly = (int)(v + sy * step_y) + 2;
 491         nlz = (int)(w + sz * step_z) + 2;
 492         index = 0;
 493         for (k = 0; k != nlz; k++)
 494                 for (j = 0; j != nly; j++)
 495                         for (i = 0; i != nlx; i++)
 496                                 noisebuf[index++] = noise3d(x0 + i, y0 + j, z0 + k, seed);
 497
 498         //calculate interpolations
 499         index  = 0;
 500         noisey = 0;
 501         noisez = 0;
 502         for (k = 0; k != sz; k++) {
 503                 v = orig_v;
 504                 noisey = 0;
 505                 for (j = 0; j != sy; j++) {
 506                         v000 = noisebuf[idx(0, noisey,     noisez)];
 507                         v100 = noisebuf[idx(1, noisey,     noisez)];
 508                         v010 = noisebuf[idx(0, noisey + 1, noisez)];
 509                         v110 = noisebuf[idx(1, noisey + 1, noisez)];
 510                         v001 = noisebuf[idx(0, noisey,     noisez + 1)];
 511                         v101 = noisebuf[idx(1, noisey,     noisez + 1)];
 512                         v011 = noisebuf[idx(0, noisey + 1, noisez + 1)];
 513                         v111 = noisebuf[idx(1, noisey + 1, noisez + 1)];
 514
 515                         u = orig_u;
 516                         noisex = 0;
 517                         for (i = 0; i != sx; i++) {
 518                                 buf[index++] = interpolate(
 519                                                                         v000, v100, v010, v110,
 520                                                                         v001, v101, v011, v111,
 521                                                                         u, v, w);
 522
 523                                 u += step_x;
 524                                 if (u >= 1.0) {
 525                                         u -= 1.0;
 526                                         noisex++;
 527                                         v000 = v100;
 528                                         v010 = v110;
 529                                         v100 = noisebuf[idx(noisex + 1, noisey,     noisez)];
 530                                         v110 = noisebuf[idx(noisex + 1, noisey + 1, noisez)];
 531                                         v001 = v101;
 532                                         v011 = v111;
 533                                         v101 = noisebuf[idx(noisex + 1, noisey,     noisez + 1)];
 534                                         v111 = noisebuf[idx(noisex + 1, noisey + 1, noisez + 1)];
 535                                 }
 536                         }
 537
 538                         v += step_y;
 539                         if (v >= 1.0) {
 540                                 v -= 1.0;
 541                                 noisey++;
 542                         }
 543                 }
 544
 545                 w += step_z;
 546                 if (w >= 1.0) {
 547                         w -= 1.0;
 548                         noisez++;
 549                 }
 550         }
 551 }
 552 #undef idx
 553
 554
 555 float *Noise::perlinMap2D(float x, float y)
 556 {
 557         float f = 1.0, g = 1.0;
 558         size_t bufsize = sx * sy;
 559
 560         x /= np->spread.X;
 561         y /= np->spread.Y;
 562
 563         memset(result, 0, sizeof(float) * bufsize);
 564
 565         for (int oct = 0; oct < np->octaves; oct++) {
 566                 gradientMap2D(x * f, y * f,
 567                         f / np->spread.X, f / np->spread.Y,
 568                         seed + np->seed + oct);
 569
 570                 for (size_t i = 0; i != bufsize; i++)
 571                         result[i] += g * buf[i];
 572
 573                 f *= 2.0;
 574                 g *= np->persist;
 575         }
 576
 577         return result;
 578 }
 579
 580
 581 float *Noise::perlinMap2DModulated(float x, float y, float *persist_map)
 582 {
 583         float f = 1.0;
 584         size_t bufsize = sx * sy;
 585
 586         x /= np->spread.X;
 587         y /= np->spread.Y;
 588
 589         memset(result, 0, sizeof(float) * bufsize);
 590
 591         float *g = new float[bufsize];
 592         for (size_t i = 0; i != bufsize; i++)
 593                 g[i] = 1.0;
 594
 595         for (int oct = 0; oct < np->octaves; oct++) {
 596                 gradientMap2D(x * f, y * f,
 597                         f / np->spread.X, f / np->spread.Y,
 598                         seed + np->seed + oct);
 599
 600                 for (size_t i = 0; i != bufsize; i++) {
 601                         result[i] += g[i] * buf[i];
 602                         g[i] *= persist_map[i];
 603                 }
 604
 605                 f *= 2.0;
 606         }
 607
 608         delete[] g;
 609         return result;
 610 }
 611
 612
 613 float *Noise::perlinMap3D(float x, float y, float z, bool eased)
 614 {
 615         float f = 1.0, g = 1.0;
 616         size_t bufsize = sx * sy * sz;
 617
 618         x /= np->spread.X;
 619         y /= np->spread.Y;
 620         z /= np->spread.Z;
 621
 622         memset(result, 0, sizeof(float) * bufsize);
 623
 624         for (int oct = 0; oct < np->octaves; oct++) {
 625                 gradientMap3D(x * f, y * f, z * f,
 626                         f / np->spread.X, f / np->spread.Y, f / np->spread.Z,
 627                         seed + np->seed + oct, eased);
 628
 629                 for (size_t i = 0; i != bufsize; i++)
 630                         result[i] += g * buf[i];
 631
 632                 f *= 2.0;
 633                 g *= np->persist;
 634         }
 635
 636         return result;
 637 }
 638
 639
 640 void Noise::transformNoiseMap()
 641 {
 642         size_t i = 0;
 643
 644         for (int z = 0; z != sz; z++)
 645         for (int y = 0; y != sy; y++)
 646         for (int x = 0; x != sx; x++) {
 647                 result[i] = result[i] * np->scale + np->offset;
 648                 i++;
 649         }
 650 }