src/noise.cpp

   1 /*
   2  * Minetest
   3  * Copyright (C) 2010-2014 celeron55, Perttu Ahola <celeron55@gmail.com>
   4  * Copyright (C) 2010-2014 kwolekr, Ryan Kwolek <kwolekr@minetest.net>
   5  * All rights reserved.
   6  *
   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.
  14  *
  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.
  24  */
  25
  26 #include <math.h>
  27 #include "noise.h"
  28 #include <iostream>
  29 #include <string.h> // memset
  30 #include "debug.h"
  31 #include "util/numeric.h"
  32 #include "util/string.h"
  33 #include "exceptions.h"
  34
  35 #define NOISE_MAGIC_X    1619
  36 #define NOISE_MAGIC_Y    31337
  37 #define NOISE_MAGIC_Z    52591
  38 #define NOISE_MAGIC_SEED 1013
  39
  40 typedef float (*Interp2dFxn)(
  41                 float v00, float v10, float v01, float v11,
  42                 float x, float y);
  43
  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);
  48
  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
  52 };
  53
  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},
  60         {NULL,          0}
  61 };
  62
  63 ///////////////////////////////////////////////////////////////////////////////
  64
  65 PcgRandom::PcgRandom(u64 state, u64 seq)
  66 {
  67         seed(state, seq);
  68 }
  69
  70 void PcgRandom::seed(u64 state, u64 seq)
  71 {
  72         m_state = 0U;
  73         m_inc = (seq << 1u) | 1u;
  74         next();
  75         m_state += state;
  76         next();
  77 }
  78
  79
  80 u32 PcgRandom::next()
  81 {
  82         u64 oldstate = m_state;
  83         m_state = oldstate * 6364136223846793005ULL + m_inc;
  84
  85         u32 xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;
  86         u32 rot = oldstate >> 59u;
  87         return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));
  88 }
  89
  90
  91 u32 PcgRandom::range(u32 bound)
  92 {
  93         // If the bound is 0, we cover the whole RNG's range
  94         if (bound == 0)
  95                 return next();
  96         /*
  97         If the bound is not a multiple of the RNG's range, it may cause bias,
  98         e.g. a RNG has a range from 0 to 3 and we take want a number 0 to 2.
  99         Using rand() % 3, the number 0 would be twice as likely to appear.
 100         With a very large RNG range, the effect becomes less prevalent but
 101         still present.  This can be solved by modifying the range of the RNG
 102         to become a multiple of bound by dropping values above the a threshhold.
 103         In our example, threshhold == 4 - 3 = 1 % 3 == 1, so reject 0, thus
 104         making the range 3 with no bias.
 105
 106         This loop looks dangerous, but will always terminate due to the
 107         RNG's property of uniformity.
 108         */
 109         u32 threshhold = -bound % bound;
 110         u32 r;
 111
 112         while ((r = next()) < threshhold)
 113                 ;
 114
 115         return r % bound;
 116 }
 117
 118
 119 s32 PcgRandom::range(s32 min, s32 max)
 120 {
 121         if (max < min)
 122                 throw PrngException("Invalid range (max < min)");
 123
 124         u32 bound = max - min + 1;
 125         return range(bound) + min;
 126 }
 127
 128
 129 void PcgRandom::bytes(void *out, size_t len)
 130 {
 131         u8 *outb = (u8 *)out;
 132         int bytes_left = 0;
 133         u32 r;
 134
 135         while (len--) {
 136                 if (bytes_left == 0) {
 137                         bytes_left = sizeof(u32);
 138                         r = next();
 139                 }
 140
 141                 *outb = r & 0xFF;
 142                 outb++;
 143                 bytes_left--;
 144                 r >>= CHAR_BIT;
 145         }
 146 }
 147
 148
 149 s32 PcgRandom::randNormalDist(s32 min, s32 max, int num_trials)
 150 {
 151         s32 accum = 0;
 152         for (int i = 0; i != num_trials; i++)
 153                 accum += range(min, max);
 154         return myround((float)accum / num_trials);
 155 }
 156
 157 ///////////////////////////////////////////////////////////////////////////////
 158
 159 float noise2d(int x, int y, int seed)
 160 {
 161         int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y
 162                         + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
 163         n = (n >> 13) ^ n;
 164         n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
 165         return 1.f - (float)n / 0x40000000;
 166 }
 167
 168
 169 float noise3d(int x, int y, int z, int seed)
 170 {
 171         int n = (NOISE_MAGIC_X * x + NOISE_MAGIC_Y * y + NOISE_MAGIC_Z * z
 172                         + NOISE_MAGIC_SEED * seed) & 0x7fffffff;
 173         n = (n >> 13) ^ n;
 174         n = (n * (n * n * 60493 + 19990303) + 1376312589) & 0x7fffffff;
 175         return 1.f - (float)n / 0x40000000;
 176 }
 177
 178
 179 inline float dotProduct(float vx, float vy, float wx, float wy)
 180 {
 181         return vx * wx + vy * wy;
 182 }
 183
 184
 185 inline float linearInterpolation(float v0, float v1, float t)
 186 {
 187         return v0 + (v1 - v0) * t;
 188 }
 189
 190
 191 inline float biLinearInterpolation(
 192         float v00, float v10,
 193         float v01, float v11,
 194         float x, float y)
 195 {
 196         float tx = easeCurve(x);
 197         float ty = easeCurve(y);
 198         float u = linearInterpolation(v00, v10, tx);
 199         float v = linearInterpolation(v01, v11, tx);
 200         return linearInterpolation(u, v, ty);
 201 }
 202
 203
 204 inline float biLinearInterpolationNoEase(
 205         float v00, float v10,
 206         float v01, float v11,
 207         float x, float y)
 208 {
 209         float u = linearInterpolation(v00, v10, x);
 210         float v = linearInterpolation(v01, v11, x);
 211         return linearInterpolation(u, v, y);
 212 }
 213
 214
 215 float triLinearInterpolation(
 216         float v000, float v100, float v010, float v110,
 217         float v001, float v101, float v011, float v111,
 218         float x, float y, float z)
 219 {
 220         float tx = easeCurve(x);
 221         float ty = easeCurve(y);
 222         float tz = easeCurve(z);
 223         float u = biLinearInterpolationNoEase(v000, v100, v010, v110, tx, ty);
 224         float v = biLinearInterpolationNoEase(v001, v101, v011, v111, tx, ty);
 225         return linearInterpolation(u, v, tz);
 226 }
 227
 228 float triLinearInterpolationNoEase(
 229         float v000, float v100, float v010, float v110,
 230         float v001, float v101, float v011, float v111,
 231         float x, float y, float z)
 232 {
 233         float u = biLinearInterpolationNoEase(v000, v100, v010, v110, x, y);
 234         float v = biLinearInterpolationNoEase(v001, v101, v011, v111, x, y);
 235         return linearInterpolation(u, v, z);
 236 }
 237
 238 float noise2d_gradient(float x, float y, int seed, bool eased)
 239 {
 240         // Calculate the integer coordinates
 241         int x0 = myfloor(x);
 242         int y0 = myfloor(y);
 243         // Calculate the remaining part of the coordinates
 244         float xl = x - (float)x0;
 245         float yl = y - (float)y0;
 246         // Get values for corners of square
 247         float v00 = noise2d(x0, y0, seed);
 248         float v10 = noise2d(x0+1, y0, seed);
 249         float v01 = noise2d(x0, y0+1, seed);
 250         float v11 = noise2d(x0+1, y0+1, seed);
 251         // Interpolate
 252         if (eased)
 253                 return biLinearInterpolation(v00, v10, v01, v11, xl, yl);
 254         else
 255                 return biLinearInterpolationNoEase(v00, v10, v01, v11, xl, yl);
 256 }
 257
 258
 259 float noise3d_gradient(float x, float y, float z, int seed, bool eased)
 260 {
 261         // Calculate the integer coordinates
 262         int x0 = myfloor(x);
 263         int y0 = myfloor(y);
 264         int z0 = myfloor(z);
 265         // Calculate the remaining part of the coordinates
 266         float xl = x - (float)x0;
 267         float yl = y - (float)y0;
 268         float zl = z - (float)z0;
 269         // Get values for corners of cube
 270         float v000 = noise3d(x0,     y0,     z0,     seed);
 271         float v100 = noise3d(x0 + 1, y0,     z0,     seed);
 272         float v010 = noise3d(x0,     y0 + 1, z0,     seed);
 273         float v110 = noise3d(x0 + 1, y0 + 1, z0,     seed);
 274         float v001 = noise3d(x0,     y0,     z0 + 1, seed);
 275         float v101 = noise3d(x0 + 1, y0,     z0 + 1, seed);
 276         float v011 = noise3d(x0,     y0 + 1, z0 + 1, seed);
 277         float v111 = noise3d(x0 + 1, y0 + 1, z0 + 1, seed);
 278         // Interpolate
 279         if (eased) {
 280                 return triLinearInterpolation(
 281                         v000, v100, v010, v110,
 282                         v001, v101, v011, v111,
 283                         xl, yl, zl);
 284         } else {
 285                 return triLinearInterpolationNoEase(
 286                         v000, v100, v010, v110,
 287                         v001, v101, v011, v111,
 288                         xl, yl, zl);
 289         }
 290 }
 291
 292
 293 float noise2d_perlin(float x, float y, int seed,
 294         int octaves, float persistence, bool eased)
 295 {
 296         float a = 0;
 297         float f = 1.0;
 298         float g = 1.0;
 299         for (int i = 0; i < octaves; i++)
 300         {
 301                 a += g * noise2d_gradient(x * f, y * f, seed + i, eased);
 302                 f *= 2.0;
 303                 g *= persistence;
 304         }
 305         return a;
 306 }
 307
 308
 309 float noise2d_perlin_abs(float x, float y, int seed,
 310         int octaves, float persistence, bool eased)
 311 {
 312         float a = 0;
 313         float f = 1.0;
 314         float g = 1.0;
 315         for (int i = 0; i < octaves; i++) {
 316                 a += g * fabs(noise2d_gradient(x * f, y * f, seed + i, eased));
 317                 f *= 2.0;
 318                 g *= persistence;
 319         }
 320         return a;
 321 }
 322
 323
 324 float noise3d_perlin(float x, float y, float z, int seed,
 325         int octaves, float persistence, bool eased)
 326 {
 327         float a = 0;
 328         float f = 1.0;
 329         float g = 1.0;
 330         for (int i = 0; i < octaves; i++) {
 331                 a += g * noise3d_gradient(x * f, y * f, z * f, seed + i, eased);
 332                 f *= 2.0;
 333                 g *= persistence;
 334         }
 335         return a;
 336 }
 337
 338
 339 float noise3d_perlin_abs(float x, float y, float z, int seed,
 340         int octaves, float persistence, bool eased)
 341 {
 342         float a = 0;
 343         float f = 1.0;
 344         float g = 1.0;
 345         for (int i = 0; i < octaves; i++) {
 346                 a += g * fabs(noise3d_gradient(x * f, y * f, z * f, seed + i, eased));
 347                 f *= 2.0;
 348                 g *= persistence;
 349         }
 350         return a;
 351 }
 352
 353
 354 float contour(float v)
 355 {
 356         v = fabs(v);
 357         if (v >= 1.0)
 358                 return 0.0;
 359         return (1.0 - v);
 360 }
 361
 362
 363 ///////////////////////// [ New noise ] ////////////////////////////
 364
 365
 366 float NoisePerlin2D(NoiseParams *np, float x, float y, int seed)
 367 {
 368         float a = 0;
 369         float f = 1.0;
 370         float g = 1.0;
 371
 372         x /= np->spread.X;
 373         y /= np->spread.Y;
 374         seed += np->seed;
 375
 376         for (size_t i = 0; i < np->octaves; i++) {
 377                 float noiseval = noise2d_gradient(x * f, y * f, seed + i,
 378                         np->flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED));
 379
 380                 if (np->flags & NOISE_FLAG_ABSVALUE)
 381                         noiseval = fabs(noiseval);
 382
 383                 a += g * noiseval;
 384                 f *= np->lacunarity;
 385                 g *= np->persist;
 386         }
 387
 388         return np->offset + a * np->scale;
 389 }
 390
 391
 392 float NoisePerlin3D(NoiseParams *np, float x, float y, float z, int seed)
 393 {
 394         float a = 0;
 395         float f = 1.0;
 396         float g = 1.0;
 397
 398         x /= np->spread.X;
 399         y /= np->spread.Y;
 400         z /= np->spread.Z;
 401         seed += np->seed;
 402
 403         for (size_t i = 0; i < np->octaves; i++) {
 404                 float noiseval = noise3d_gradient(x * f, y * f, z * f, seed + i,
 405                         np->flags & NOISE_FLAG_EASED);
 406
 407                 if (np->flags & NOISE_FLAG_ABSVALUE)
 408                         noiseval = fabs(noiseval);
 409
 410                 a += g * noiseval;
 411                 f *= np->lacunarity;
 412                 g *= np->persist;
 413         }
 414
 415         return np->offset + a * np->scale;
 416 }
 417
 418
 419 Noise::Noise(NoiseParams *np_, int seed, u32 sx, u32 sy, u32 sz)
 420 {
 421         memcpy(&np, np_, sizeof(np));
 422         this->seed = seed;
 423         this->sx   = sx;
 424         this->sy   = sy;
 425         this->sz   = sz;
 426
 427         this->persist_buf  = NULL;
 428         this->gradient_buf = NULL;
 429         this->result       = NULL;
 430
 431         allocBuffers();
 432 }
 433
 434
 435 Noise::~Noise()
 436 {
 437         delete[] gradient_buf;
 438         delete[] persist_buf;
 439         delete[] noise_buf;
 440         delete[] result;
 441 }
 442
 443
 444 void Noise::allocBuffers()
 445 {
 446         if (sx < 1)
 447                 sx = 1;
 448         if (sy < 1)
 449                 sy = 1;
 450         if (sz < 1)
 451                 sz = 1;
 452
 453         this->noise_buf = NULL;
 454         resizeNoiseBuf(sz > 1);
 455
 456         delete[] gradient_buf;
 457         delete[] persist_buf;
 458         delete[] result;
 459
 460         try {
 461                 size_t bufsize = sx * sy * sz;
 462                 this->persist_buf  = NULL;
 463                 this->gradient_buf = new float[bufsize];
 464                 this->result       = new float[bufsize];
 465         } catch (std::bad_alloc &e) {
 466                 throw InvalidNoiseParamsException();
 467         }
 468 }
 469
 470
 471 void Noise::setSize(u32 sx, u32 sy, u32 sz)
 472 {
 473         this->sx = sx;
 474         this->sy = sy;
 475         this->sz = sz;
 476
 477         allocBuffers();
 478 }
 479
 480
 481 void Noise::setSpreadFactor(v3f spread)
 482 {
 483         this->np.spread = spread;
 484
 485         resizeNoiseBuf(sz > 1);
 486 }
 487
 488
 489 void Noise::setOctaves(int octaves)
 490 {
 491         this->np.octaves = octaves;
 492
 493         resizeNoiseBuf(sz > 1);
 494 }
 495
 496
 497 void Noise::resizeNoiseBuf(bool is3d)
 498 {
 499         //maximum possible spread value factor
 500         float ofactor = (np.lacunarity > 1.0) ?
 501                 pow(np.lacunarity, np.octaves - 1) :
 502                 np.lacunarity;
 503
 504         // noise lattice point count
 505         // (int)(sz * spread * ofactor) is # of lattice points crossed due to length
 506         float num_noise_points_x = sx * ofactor / np.spread.X;
 507         float num_noise_points_y = sy * ofactor / np.spread.Y;
 508         float num_noise_points_z = sz * ofactor / np.spread.Z;
 509
 510         // protect against obviously invalid parameters
 511         if (num_noise_points_x > 1000000000.f ||
 512                 num_noise_points_y > 1000000000.f ||
 513                 num_noise_points_z > 1000000000.f)
 514                 throw InvalidNoiseParamsException();
 515
 516         // + 2 for the two initial endpoints
 517         // + 1 for potentially crossing a boundary due to offset
 518         size_t nlx = (size_t)ceil(num_noise_points_x) + 3;
 519         size_t nly = (size_t)ceil(num_noise_points_y) + 3;
 520         size_t nlz = is3d ? (size_t)ceil(num_noise_points_z) + 3 : 1;
 521
 522         delete[] noise_buf;
 523         try {
 524                 noise_buf = new float[nlx * nly * nlz];
 525         } catch (std::bad_alloc &e) {
 526                 throw InvalidNoiseParamsException();
 527         }
 528 }
 529
 530
 531 /*
 532  * NB:  This algorithm is not optimal in terms of space complexity.  The entire
 533  * integer lattice of noise points could be done as 2 lines instead, and for 3D,
 534  * 2 lines + 2 planes.
 535  * However, this would require the noise calls to be interposed with the
 536  * interpolation loops, which may trash the icache, leading to lower overall
 537  * performance.
 538  * Another optimization that could save half as many noise calls is to carry over
 539  * values from the previous noise lattice as midpoints in the new lattice for the
 540  * next octave.
 541  */
 542 #define idx(x, y) ((y) * nlx + (x))
 543 void Noise::gradientMap2D(
 544                 float x, float y,
 545                 float step_x, float step_y,
 546                 int seed)
 547 {
 548         float v00, v01, v10, v11, u, v, orig_u;
 549         u32 index, i, j, noisex, noisey;
 550         u32 nlx, nly;
 551         s32 x0, y0;
 552
 553         bool eased = np.flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED);
 554         Interp2dFxn interpolate = eased ?
 555                 biLinearInterpolation : biLinearInterpolationNoEase;
 556
 557         x0 = floor(x);
 558         y0 = floor(y);
 559         u = x - (float)x0;
 560         v = y - (float)y0;
 561         orig_u = u;
 562
 563         //calculate noise point lattice
 564         nlx = (u32)(u + sx * step_x) + 2;
 565         nly = (u32)(v + sy * step_y) + 2;
 566         index = 0;
 567         for (j = 0; j != nly; j++)
 568                 for (i = 0; i != nlx; i++)
 569                         noise_buf[index++] = noise2d(x0 + i, y0 + j, seed);
 570
 571         //calculate interpolations
 572         index  = 0;
 573         noisey = 0;
 574         for (j = 0; j != sy; j++) {
 575                 v00 = noise_buf[idx(0, noisey)];
 576                 v10 = noise_buf[idx(1, noisey)];
 577                 v01 = noise_buf[idx(0, noisey + 1)];
 578                 v11 = noise_buf[idx(1, noisey + 1)];
 579
 580                 u = orig_u;
 581                 noisex = 0;
 582                 for (i = 0; i != sx; i++) {
 583                         gradient_buf[index++] = interpolate(v00, v10, v01, v11, u, v);
 584
 585                         u += step_x;
 586                         if (u >= 1.0) {
 587                                 u -= 1.0;
 588                                 noisex++;
 589                                 v00 = v10;
 590                                 v01 = v11;
 591                                 v10 = noise_buf[idx(noisex + 1, noisey)];
 592                                 v11 = noise_buf[idx(noisex + 1, noisey + 1)];
 593                         }
 594                 }
 595
 596                 v += step_y;
 597                 if (v >= 1.0) {
 598                         v -= 1.0;
 599                         noisey++;
 600                 }
 601         }
 602 }
 603 #undef idx
 604
 605
 606 #define idx(x, y, z) ((z) * nly * nlx + (y) * nlx + (x))
 607 void Noise::gradientMap3D(
 608                 float x, float y, float z,
 609                 float step_x, float step_y, float step_z,
 610                 int seed)
 611 {
 612         float v000, v010, v100, v110;
 613         float v001, v011, v101, v111;
 614         float u, v, w, orig_u, orig_v;
 615         u32 index, i, j, k, noisex, noisey, noisez;
 616         u32 nlx, nly, nlz;
 617         s32 x0, y0, z0;
 618
 619         Interp3dFxn interpolate = (np.flags & NOISE_FLAG_EASED) ?
 620                 triLinearInterpolation : triLinearInterpolationNoEase;
 621
 622         x0 = floor(x);
 623         y0 = floor(y);
 624         z0 = floor(z);
 625         u = x - (float)x0;
 626         v = y - (float)y0;
 627         w = z - (float)z0;
 628         orig_u = u;
 629         orig_v = v;
 630
 631         //calculate noise point lattice
 632         nlx = (u32)(u + sx * step_x) + 2;
 633         nly = (u32)(v + sy * step_y) + 2;
 634         nlz = (u32)(w + sz * step_z) + 2;
 635         index = 0;
 636         for (k = 0; k != nlz; k++)
 637                 for (j = 0; j != nly; j++)
 638                         for (i = 0; i != nlx; i++)
 639                                 noise_buf[index++] = noise3d(x0 + i, y0 + j, z0 + k, seed);
 640
 641         //calculate interpolations
 642         index  = 0;
 643         noisey = 0;
 644         noisez = 0;
 645         for (k = 0; k != sz; k++) {
 646                 v = orig_v;
 647                 noisey = 0;
 648                 for (j = 0; j != sy; j++) {
 649                         v000 = noise_buf[idx(0, noisey,     noisez)];
 650                         v100 = noise_buf[idx(1, noisey,     noisez)];
 651                         v010 = noise_buf[idx(0, noisey + 1, noisez)];
 652                         v110 = noise_buf[idx(1, noisey + 1, noisez)];
 653                         v001 = noise_buf[idx(0, noisey,     noisez + 1)];
 654                         v101 = noise_buf[idx(1, noisey,     noisez + 1)];
 655                         v011 = noise_buf[idx(0, noisey + 1, noisez + 1)];
 656                         v111 = noise_buf[idx(1, noisey + 1, noisez + 1)];
 657
 658                         u = orig_u;
 659                         noisex = 0;
 660                         for (i = 0; i != sx; i++) {
 661                                 gradient_buf[index++] = interpolate(
 662                                         v000, v100, v010, v110,
 663                                         v001, v101, v011, v111,
 664                                         u, v, w);
 665
 666                                 u += step_x;
 667                                 if (u >= 1.0) {
 668                                         u -= 1.0;
 669                                         noisex++;
 670                                         v000 = v100;
 671                                         v010 = v110;
 672                                         v100 = noise_buf[idx(noisex + 1, noisey,     noisez)];
 673                                         v110 = noise_buf[idx(noisex + 1, noisey + 1, noisez)];
 674                                         v001 = v101;
 675                                         v011 = v111;
 676                                         v101 = noise_buf[idx(noisex + 1, noisey,     noisez + 1)];
 677                                         v111 = noise_buf[idx(noisex + 1, noisey + 1, noisez + 1)];
 678                                 }
 679                         }
 680
 681                         v += step_y;
 682                         if (v >= 1.0) {
 683                                 v -= 1.0;
 684                                 noisey++;
 685                         }
 686                 }
 687
 688                 w += step_z;
 689                 if (w >= 1.0) {
 690                         w -= 1.0;
 691                         noisez++;
 692                 }
 693         }
 694 }
 695 #undef idx
 696
 697
 698 float *Noise::perlinMap2D(float x, float y, float *persistence_map)
 699 {
 700         float f = 1.0, g = 1.0;
 701         size_t bufsize = sx * sy;
 702
 703         x /= np.spread.X;
 704         y /= np.spread.Y;
 705
 706         memset(result, 0, sizeof(float) * bufsize);
 707
 708         if (persistence_map) {
 709                 if (!persist_buf)
 710                         persist_buf = new float[bufsize];
 711                 for (size_t i = 0; i != bufsize; i++)
 712                         persist_buf[i] = 1.0;
 713         }
 714
 715         for (size_t oct = 0; oct < np.octaves; oct++) {
 716                 gradientMap2D(x * f, y * f,
 717                         f / np.spread.X, f / np.spread.Y,
 718                         seed + np.seed + oct);
 719
 720                 updateResults(g, persist_buf, persistence_map, bufsize);
 721
 722                 f *= np.lacunarity;
 723                 g *= np.persist;
 724         }
 725
 726         if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
 727                 for (size_t i = 0; i != bufsize; i++)
 728                         result[i] = result[i] * np.scale + np.offset;
 729         }
 730
 731         return result;
 732 }
 733
 734
 735 float *Noise::perlinMap3D(float x, float y, float z, float *persistence_map)
 736 {
 737         float f = 1.0, g = 1.0;
 738         size_t bufsize = sx * sy * sz;
 739
 740         x /= np.spread.X;
 741         y /= np.spread.Y;
 742         z /= np.spread.Z;
 743
 744         memset(result, 0, sizeof(float) * bufsize);
 745
 746         if (persistence_map) {
 747                 if (!persist_buf)
 748                         persist_buf = new float[bufsize];
 749                 for (size_t i = 0; i != bufsize; i++)
 750                         persist_buf[i] = 1.0;
 751         }
 752
 753         for (size_t oct = 0; oct < np.octaves; oct++) {
 754                 gradientMap3D(x * f, y * f, z * f,
 755                         f / np.spread.X, f / np.spread.Y, f / np.spread.Z,
 756                         seed + np.seed + oct);
 757
 758                 updateResults(g, persist_buf, persistence_map, bufsize);
 759
 760                 f *= np.lacunarity;
 761                 g *= np.persist;
 762         }
 763
 764         if (fabs(np.offset - 0.f) > 0.00001 || fabs(np.scale - 1.f) > 0.00001) {
 765                 for (size_t i = 0; i != bufsize; i++)
 766                         result[i] = result[i] * np.scale + np.offset;
 767         }
 768
 769         return result;
 770 }
 771
 772
 773 void Noise::updateResults(float g, float *gmap,
 774         float *persistence_map, size_t bufsize)
 775 {
 776         // This looks very ugly, but it is 50-70% faster than having
 777         // conditional statements inside the loop
 778         if (np.flags & NOISE_FLAG_ABSVALUE) {
 779                 if (persistence_map) {
 780                         for (size_t i = 0; i != bufsize; i++) {
 781                                 result[i] += gmap[i] * fabs(gradient_buf[i]);
 782                                 gmap[i] *= persistence_map[i];
 783                         }
 784                 } else {
 785                         for (size_t i = 0; i != bufsize; i++)
 786                                 result[i] += g * fabs(gradient_buf[i]);
 787                 }
 788         } else {
 789                 if (persistence_map) {
 790                         for (size_t i = 0; i != bufsize; i++) {
 791                                 result[i] += gmap[i] * gradient_buf[i];
 792                                 gmap[i] *= persistence_map[i];
 793                         }
 794                 } else {
 795                         for (size_t i = 0; i != bufsize; i++)
 796                                 result[i] += g * gradient_buf[i];
 797                 }
 798         }
 799 }