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