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