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