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 noise2d_perlin_abs(float x, float y, s32 seed,
 316         int octaves, float persistence, bool eased)
 317 {
 318         float a = 0;
 319         float f = 1.0;
 320         float g = 1.0;
 321         for (int i = 0; i < octaves; i++) {
 322                 a += g * std::fabs(noise2d_gradient(x * f, y * f, seed + i, eased));
 323                 f *= 2.0;
 324                 g *= persistence;
 325         }
 326         return a;
 327 }
 328
 329
 330 float noise3d_perlin(float x, float y, float z, s32 seed,
 331         int octaves, float persistence, bool eased)
 332 {
 333         float a = 0;
 334         float f = 1.0;
 335         float g = 1.0;
 336         for (int i = 0; i < octaves; i++) {
 337                 a += g * noise3d_gradient(x * f, y * f, z * f, seed + i, eased);
 338                 f *= 2.0;
 339                 g *= persistence;
 340         }
 341         return a;
 342 }
 343
 344
 345 float noise3d_perlin_abs(float x, float y, float z, s32 seed,
 346         int octaves, float persistence, bool eased)
 347 {
 348         float a = 0;
 349         float f = 1.0;
 350         float g = 1.0;
 351         for (int i = 0; i < octaves; i++) {
 352                 a += g * std::fabs(noise3d_gradient(x * f, y * f, z * f, seed + i, eased));
 353                 f *= 2.0;
 354                 g *= persistence;
 355         }
 356         return a;
 357 }
 358
 359
 360 float contour(float v)
 361 {
 362         v = std::fabs(v);
 363         if (v >= 1.0)
 364                 return 0.0;
 365         return (1.0 - v);
 366 }
 367
 368
 369 ///////////////////////// [ New noise ] ////////////////////////////
 370
 371
 372 float NoisePerlin2D(NoiseParams *np, float x, float y, s32 seed)
 373 {
 374         float a = 0;
 375         float f = 1.0;
 376         float g = 1.0;
 377
 378         x /= np->spread.X;
 379         y /= np->spread.Y;
 380         seed += np->seed;
 381
 382         for (size_t i = 0; i < np->octaves; i++) {
 383                 float noiseval = noise2d_gradient(x * f, y * f, seed + i,
 384                         np->flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED));
 385
 386                 if (np->flags & NOISE_FLAG_ABSVALUE)
 387                         noiseval = std::fabs(noiseval);
 388
 389                 a += g * noiseval;
 390                 f *= np->lacunarity;
 391                 g *= np->persist;
 392         }
 393
 394         return np->offset + a * np->scale;
 395 }
 396
 397
 398 float NoisePerlin3D(NoiseParams *np, float x, float y, float z, s32 seed)
 399 {
 400         float a = 0;
 401         float f = 1.0;
 402         float g = 1.0;
 403
 404         x /= np->spread.X;
 405         y /= np->spread.Y;
 406         z /= np->spread.Z;
 407         seed += np->seed;
 408
 409         for (size_t i = 0; i < np->octaves; i++) {
 410                 float noiseval = noise3d_gradient(x * f, y * f, z * f, seed + i,
 411                         np->flags & NOISE_FLAG_EASED);
 412
 413                 if (np->flags & NOISE_FLAG_ABSVALUE)
 414                         noiseval = std::fabs(noiseval);
 415
 416                 a += g * noiseval;
 417                 f *= np->lacunarity;
 418                 g *= np->persist;
 419         }
 420
 421         return np->offset + a * np->scale;
 422 }
 423
 424
 425 Noise::Noise(NoiseParams *np_, s32 seed, u32 sx, u32 sy, u32 sz)
 426 {
 427         memcpy(&np, np_, sizeof(np));
 428         this->seed = seed;
 429         this->sx   = sx;
 430         this->sy   = sy;
 431         this->sz   = sz;
 432
 433         allocBuffers();
 434 }
 435
 436
 437 Noise::~Noise()
 438 {
 439         delete[] gradient_buf;
 440         delete[] persist_buf;
 441         delete[] noise_buf;
 442         delete[] result;
 443 }
 444
 445
 446 void Noise::allocBuffers()
 447 {
 448         if (sx < 1)
 449                 sx = 1;
 450         if (sy < 1)
 451                 sy = 1;
 452         if (sz < 1)
 453                 sz = 1;
 454
 455         this->noise_buf = NULL;
 456         resizeNoiseBuf(sz > 1);
 457
 458         delete[] gradient_buf;
 459         delete[] persist_buf;
 460         delete[] result;
 461
 462         try {
 463                 size_t bufsize = sx * sy * sz;
 464                 this->persist_buf  = NULL;
 465                 this->gradient_buf = new float[bufsize];
 466                 this->result       = new float[bufsize];
 467         } catch (std::bad_alloc &e) {
 468                 throw InvalidNoiseParamsException();
 469         }
 470 }
 471
 472
 473 void Noise::setSize(u32 sx, u32 sy, u32 sz)
 474 {
 475         this->sx = sx;
 476         this->sy = sy;
 477         this->sz = sz;
 478
 479         allocBuffers();
 480 }
 481
 482
 483 void Noise::setSpreadFactor(v3f spread)
 484 {
 485         this->np.spread = spread;
 486
 487         resizeNoiseBuf(sz > 1);
 488 }
 489
 490
 491 void Noise::setOctaves(int octaves)
 492 {
 493         this->np.octaves = octaves;
 494
 495         resizeNoiseBuf(sz > 1);
 496 }
 497
 498
 499 void Noise::resizeNoiseBuf(bool is3d)
 500 {
 501         // Maximum possible spread value factor
 502         float ofactor = (np.lacunarity > 1.0) ?
 503                 pow(np.lacunarity, np.octaves - 1) :
 504                 np.lacunarity;
 505
 506         // Noise lattice point count
 507         // (int)(sz * spread * ofactor) is # of lattice points crossed due to length
 508         float num_noise_points_x = sx * ofactor / np.spread.X;
 509         float num_noise_points_y = sy * ofactor / np.spread.Y;
 510         float num_noise_points_z = sz * ofactor / np.spread.Z;
 511
 512         // Protect against obviously invalid parameters
 513         if (num_noise_points_x > 1000000000.f ||
 514                         num_noise_points_y > 1000000000.f ||
 515                         num_noise_points_z > 1000000000.f)
 516                 throw InvalidNoiseParamsException();
 517
 518         // Protect against an octave having a spread < 1, causing broken noise values
 519         if (np.spread.X / ofactor < 1.0f ||
 520                         np.spread.Y / ofactor < 1.0f ||
 521                         np.spread.Z / ofactor < 1.0f) {
 522                 errorstream << "A noise parameter has too many octaves: "
 523                         << np.octaves << " octaves" << std::endl;
 524                 throw InvalidNoiseParamsException("A noise parameter has too many octaves");
 525         }
 526
 527         // + 2 for the two initial endpoints
 528         // + 1 for potentially crossing a boundary due to offset
 529         size_t nlx = (size_t)std::ceil(num_noise_points_x) + 3;
 530         size_t nly = (size_t)std::ceil(num_noise_points_y) + 3;
 531         size_t nlz = is3d ? (size_t)std::ceil(num_noise_points_z) + 3 : 1;
 532
 533         delete[] noise_buf;
 534         try {
 535                 noise_buf = new float[nlx * nly * nlz];
 536         } catch (std::bad_alloc &e) {
 537                 throw InvalidNoiseParamsException();
 538         }
 539 }
 540
 541
 542 /*
 543  * NB:  This algorithm is not optimal in terms of space complexity.  The entire
 544  * integer lattice of noise points could be done as 2 lines instead, and for 3D,
 545  * 2 lines + 2 planes.
 546  * However, this would require the noise calls to be interposed with the
 547  * interpolation loops, which may trash the icache, leading to lower overall
 548  * performance.
 549  * Another optimization that could save half as many noise calls is to carry over
 550  * values from the previous noise lattice as midpoints in the new lattice for the
 551  * next octave.
 552  */
 553 #define idx(x, y) ((y) * nlx + (x))
 554 void Noise::gradientMap2D(
 555                 float x, float y,
 556                 float step_x, float step_y,
 557                 s32 seed)
 558 {
 559         float v00, v01, v10, v11, u, v, orig_u;
 560         u32 index, i, j, noisex, noisey;
 561         u32 nlx, nly;
 562         s32 x0, y0;
 563
 564         bool eased = np.flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED);
 565         Interp2dFxn interpolate = eased ?
 566                 biLinearInterpolation : biLinearInterpolationNoEase;
 567
 568         x0 = std::floor(x);
 569         y0 = std::floor(y);
 570         u = x - (float)x0;
 571         v = y - (float)y0;
 572         orig_u = u;
 573
 574         //calculate noise point lattice
 575         nlx = (u32)(u + sx * step_x) + 2;
 576         nly = (u32)(v + sy * step_y) + 2;
 577         index = 0;
 578         for (j = 0; j != nly; j++)
 579                 for (i = 0; i != nlx; i++)
 580                         noise_buf[index++] = noise2d(x0 + i, y0 + j, seed);
 581
 582         //calculate interpolations
 583         index  = 0;
 584         noisey = 0;
 585         for (j = 0; j != sy; j++) {
 586                 v00 = noise_buf[idx(0, noisey)];
 587                 v10 = noise_buf[idx(1, noisey)];
 588                 v01 = noise_buf[idx(0, noisey + 1)];
 589                 v11 = noise_buf[idx(1, noisey + 1)];
 590
 591                 u = orig_u;
 592                 noisex = 0;
 593                 for (i = 0; i != sx; i++) {
 594                         gradient_buf[index++] = interpolate(v00, v10, v01, v11, u, v);
 595
 596                         u += step_x;
 597                         if (u >= 1.0) {
 598                                 u -= 1.0;
 599                                 noisex++;
 600                                 v00 = v10;
 601                                 v01 = v11;
 602                                 v10 = noise_buf[idx(noisex + 1, noisey)];
 603                                 v11 = noise_buf[idx(noisex + 1, noisey + 1)];
 604                         }
 605                 }
 606
 607                 v += step_y;
 608                 if (v >= 1.0) {
 609                         v -= 1.0;
 610                         noisey++;
 611                 }
 612         }
 613 }
 614 #undef idx
 615
 616
 617 #define idx(x, y, z) ((z) * nly * nlx + (y) * nlx + (x))
 618 void Noise::gradientMap3D(
 619                 float x, float y, float z,
 620                 float step_x, float step_y, float step_z,
 621                 s32 seed)
 622 {
 623         float v000, v010, v100, v110;
 624         float v001, v011, v101, v111;
 625         float u, v, w, orig_u, orig_v;
 626         u32 index, i, j, k, noisex, noisey, noisez;
 627         u32 nlx, nly, nlz;
 628         s32 x0, y0, z0;
 629
 630         Interp3dFxn interpolate = (np.flags & NOISE_FLAG_EASED) ?
 631                 triLinearInterpolation : triLinearInterpolationNoEase;
 632
 633         x0 = std::floor(x);
 634         y0 = std::floor(y);
 635         z0 = std::floor(z);
 636         u = x - (float)x0;
 637         v = y - (float)y0;
 638         w = z - (float)z0;
 639         orig_u = u;
 640         orig_v = v;
 641
 642         //calculate noise point lattice
 643         nlx = (u32)(u + sx * step_x) + 2;
 644         nly = (u32)(v + sy * step_y) + 2;
 645         nlz = (u32)(w + sz * step_z) + 2;
 646         index = 0;
 647         for (k = 0; k != nlz; k++)
 648                 for (j = 0; j != nly; j++)
 649                         for (i = 0; i != nlx; i++)
 650                                 noise_buf[index++] = noise3d(x0 + i, y0 + j, z0 + k, seed);
 651
 652         //calculate interpolations
 653         index  = 0;
 654         noisey = 0;
 655         noisez = 0;
 656         for (k = 0; k != sz; k++) {
 657                 v = orig_v;
 658                 noisey = 0;
 659                 for (j = 0; j != sy; j++) {
 660                         v000 = noise_buf[idx(0, noisey,     noisez)];
 661                         v100 = noise_buf[idx(1, noisey,     noisez)];
 662                         v010 = noise_buf[idx(0, noisey + 1, noisez)];
 663                         v110 = noise_buf[idx(1, noisey + 1, noisez)];
 664                         v001 = noise_buf[idx(0, noisey,     noisez + 1)];
 665                         v101 = noise_buf[idx(1, noisey,     noisez + 1)];
 666                         v011 = noise_buf[idx(0, noisey + 1, noisez + 1)];
 667                         v111 = noise_buf[idx(1, noisey + 1, noisez + 1)];
 668
 669                         u = orig_u;
 670                         noisex = 0;
 671                         for (i = 0; i != sx; i++) {
 672                                 gradient_buf[index++] = interpolate(
 673                                         v000, v100, v010, v110,
 674                                         v001, v101, v011, v111,
 675                                         u, v, w);
 676
 677                                 u += step_x;
 678                                 if (u >= 1.0) {
 679                                         u -= 1.0;
 680                                         noisex++;
 681                                         v000 = v100;
 682                                         v010 = v110;
 683                                         v100 = noise_buf[idx(noisex + 1, noisey,     noisez)];
 684                                         v110 = noise_buf[idx(noisex + 1, noisey + 1, noisez)];
 685                                         v001 = v101;
 686                                         v011 = v111;
 687                                         v101 = noise_buf[idx(noisex + 1, noisey,     noisez + 1)];
 688                                         v111 = noise_buf[idx(noisex + 1, noisey + 1, noisez + 1)];
 689                                 }
 690                         }
 691
 692                         v += step_y;
 693                         if (v >= 1.0) {
 694                                 v -= 1.0;
 695                                 noisey++;
 696                         }
 697                 }
 698
 699                 w += step_z;
 700                 if (w >= 1.0) {
 701                         w -= 1.0;
 702                         noisez++;
 703                 }
 704         }
 705 }
 706 #undef idx
 707
 708
 709 float *Noise::perlinMap2D(float x, float y, float *persistence_map)
 710 {
 711         float f = 1.0, g = 1.0;
 712         size_t bufsize = sx * sy;
 713
 714         x /= np.spread.X;
 715         y /= np.spread.Y;
 716
 717         memset(result, 0, sizeof(float) * bufsize);
 718
 719         if (persistence_map) {
 720                 if (!persist_buf)
 721                         persist_buf = new float[bufsize];
 722                 for (size_t i = 0; i != bufsize; i++)
 723                         persist_buf[i] = 1.0;
 724         }
 725
 726         for (size_t oct = 0; oct < np.octaves; oct++) {
 727                 gradientMap2D(x * f, y * f,
 728                         f / np.spread.X, f / np.spread.Y,
 729                         seed + np.seed + oct);
 730
 731                 updateResults(g, persist_buf, persistence_map, bufsize);
 732
 733                 f *= np.lacunarity;
 734                 g *= np.persist;
 735         }
 736
 737         if (std::fabs(np.offset - 0.f) > 0.00001 || std::fabs(np.scale - 1.f) > 0.00001) {
 738                 for (size_t i = 0; i != bufsize; i++)
 739                         result[i] = result[i] * np.scale + np.offset;
 740         }
 741
 742         return result;
 743 }
 744
 745
 746 float *Noise::perlinMap3D(float x, float y, float z, float *persistence_map)
 747 {
 748         float f = 1.0, g = 1.0;
 749         size_t bufsize = sx * sy * sz;
 750
 751         x /= np.spread.X;
 752         y /= np.spread.Y;
 753         z /= np.spread.Z;
 754
 755         memset(result, 0, sizeof(float) * bufsize);
 756
 757         if (persistence_map) {
 758                 if (!persist_buf)
 759                         persist_buf = new float[bufsize];
 760                 for (size_t i = 0; i != bufsize; i++)
 761                         persist_buf[i] = 1.0;
 762         }
 763
 764         for (size_t oct = 0; oct < np.octaves; oct++) {
 765                 gradientMap3D(x * f, y * f, z * f,
 766                         f / np.spread.X, f / np.spread.Y, f / np.spread.Z,
 767                         seed + np.seed + oct);
 768
 769                 updateResults(g, persist_buf, persistence_map, bufsize);
 770
 771                 f *= np.lacunarity;
 772                 g *= np.persist;
 773         }
 774
 775         if (std::fabs(np.offset - 0.f) > 0.00001 || std::fabs(np.scale - 1.f) > 0.00001) {
 776                 for (size_t i = 0; i != bufsize; i++)
 777                         result[i] = result[i] * np.scale + np.offset;
 778         }
 779
 780         return result;
 781 }
 782
 783
 784 void Noise::updateResults(float g, float *gmap,
 785         const float *persistence_map, size_t bufsize)
 786 {
 787         // This looks very ugly, but it is 50-70% faster than having
 788         // conditional statements inside the loop
 789         if (np.flags & NOISE_FLAG_ABSVALUE) {
 790                 if (persistence_map) {
 791                         for (size_t i = 0; i != bufsize; i++) {
 792                                 result[i] += gmap[i] * std::fabs(gradient_buf[i]);
 793                                 gmap[i] *= persistence_map[i];
 794                         }
 795                 } else {
 796                         for (size_t i = 0; i != bufsize; i++)
 797                                 result[i] += g * std::fabs(gradient_buf[i]);
 798                 }
 799         } else {
 800                 if (persistence_map) {
 801                         for (size_t i = 0; i != bufsize; i++) {
 802                                 result[i] += gmap[i] * gradient_buf[i];
 803                                 gmap[i] *= persistence_map[i];
 804                         }
 805                 } else {
 806                         for (size_t i = 0; i != bufsize; i++)
 807                                 result[i] += g * gradient_buf[i];
 808                 }
 809         }
 810 }