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