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         // + 2 for the two initial endpoints
 524         // + 1 for potentially crossing a boundary due to offset
 525         size_t nlx = (size_t)std::ceil(num_noise_points_x) + 3;
 526         size_t nly = (size_t)std::ceil(num_noise_points_y) + 3;
 527         size_t nlz = is3d ? (size_t)std::ceil(num_noise_points_z) + 3 : 1;
 528
 529         delete[] noise_buf;
 530         try {
 531                 noise_buf = new float[nlx * nly * nlz];
 532         } catch (std::bad_alloc &e) {
 533                 throw InvalidNoiseParamsException();
 534         }
 535 }
 536
 537
 538 /*
 539  * NB:  This algorithm is not optimal in terms of space complexity.  The entire
 540  * integer lattice of noise points could be done as 2 lines instead, and for 3D,
 541  * 2 lines + 2 planes.
 542  * However, this would require the noise calls to be interposed with the
 543  * interpolation loops, which may trash the icache, leading to lower overall
 544  * performance.
 545  * Another optimization that could save half as many noise calls is to carry over
 546  * values from the previous noise lattice as midpoints in the new lattice for the
 547  * next octave.
 548  */
 549 #define idx(x, y) ((y) * nlx + (x))
 550 void Noise::gradientMap2D(
 551                 float x, float y,
 552                 float step_x, float step_y,
 553                 s32 seed)
 554 {
 555         float v00, v01, v10, v11, u, v, orig_u;
 556         u32 index, i, j, noisex, noisey;
 557         u32 nlx, nly;
 558         s32 x0, y0;
 559
 560         bool eased = np.flags & (NOISE_FLAG_DEFAULTS | NOISE_FLAG_EASED);
 561         Interp2dFxn interpolate = eased ?
 562                 biLinearInterpolation : biLinearInterpolationNoEase;
 563
 564         x0 = std::floor(x);
 565         y0 = std::floor(y);
 566         u = x - (float)x0;
 567         v = y - (float)y0;
 568         orig_u = u;
 569
 570         //calculate noise point lattice
 571         nlx = (u32)(u + sx * step_x) + 2;
 572         nly = (u32)(v + sy * step_y) + 2;
 573         index = 0;
 574         for (j = 0; j != nly; j++)
 575                 for (i = 0; i != nlx; i++)
 576                         noise_buf[index++] = noise2d(x0 + i, y0 + j, seed);
 577
 578         //calculate interpolations
 579         index  = 0;
 580         noisey = 0;
 581         for (j = 0; j != sy; j++) {
 582                 v00 = noise_buf[idx(0, noisey)];
 583                 v10 = noise_buf[idx(1, noisey)];
 584                 v01 = noise_buf[idx(0, noisey + 1)];
 585                 v11 = noise_buf[idx(1, noisey + 1)];
 586
 587                 u = orig_u;
 588                 noisex = 0;
 589                 for (i = 0; i != sx; i++) {
 590                         gradient_buf[index++] = interpolate(v00, v10, v01, v11, u, v);
 591
 592                         u += step_x;
 593                         if (u >= 1.0) {
 594                                 u -= 1.0;
 595                                 noisex++;
 596                                 v00 = v10;
 597                                 v01 = v11;
 598                                 v10 = noise_buf[idx(noisex + 1, noisey)];
 599                                 v11 = noise_buf[idx(noisex + 1, noisey + 1)];
 600                         }
 601                 }
 602
 603                 v += step_y;
 604                 if (v >= 1.0) {
 605                         v -= 1.0;
 606                         noisey++;
 607                 }
 608         }
 609 }
 610 #undef idx
 611
 612
 613 #define idx(x, y, z) ((z) * nly * nlx + (y) * nlx + (x))
 614 void Noise::gradientMap3D(
 615                 float x, float y, float z,
 616                 float step_x, float step_y, float step_z,
 617                 s32 seed)
 618 {
 619         float v000, v010, v100, v110;
 620         float v001, v011, v101, v111;
 621         float u, v, w, orig_u, orig_v;
 622         u32 index, i, j, k, noisex, noisey, noisez;
 623         u32 nlx, nly, nlz;
 624         s32 x0, y0, z0;
 625
 626         Interp3dFxn interpolate = (np.flags & NOISE_FLAG_EASED) ?
 627                 triLinearInterpolation : triLinearInterpolationNoEase;
 628
 629         x0 = std::floor(x);
 630         y0 = std::floor(y);
 631         z0 = std::floor(z);
 632         u = x - (float)x0;
 633         v = y - (float)y0;
 634         w = z - (float)z0;
 635         orig_u = u;
 636         orig_v = v;
 637
 638         //calculate noise point lattice
 639         nlx = (u32)(u + sx * step_x) + 2;
 640         nly = (u32)(v + sy * step_y) + 2;
 641         nlz = (u32)(w + sz * step_z) + 2;
 642         index = 0;
 643         for (k = 0; k != nlz; k++)
 644                 for (j = 0; j != nly; j++)
 645                         for (i = 0; i != nlx; i++)
 646                                 noise_buf[index++] = noise3d(x0 + i, y0 + j, z0 + k, seed);
 647
 648         //calculate interpolations
 649         index  = 0;
 650         noisey = 0;
 651         noisez = 0;
 652         for (k = 0; k != sz; k++) {
 653                 v = orig_v;
 654                 noisey = 0;
 655                 for (j = 0; j != sy; j++) {
 656                         v000 = noise_buf[idx(0, noisey,     noisez)];
 657                         v100 = noise_buf[idx(1, noisey,     noisez)];
 658                         v010 = noise_buf[idx(0, noisey + 1, noisez)];
 659                         v110 = noise_buf[idx(1, noisey + 1, noisez)];
 660                         v001 = noise_buf[idx(0, noisey,     noisez + 1)];
 661                         v101 = noise_buf[idx(1, noisey,     noisez + 1)];
 662                         v011 = noise_buf[idx(0, noisey + 1, noisez + 1)];
 663                         v111 = noise_buf[idx(1, noisey + 1, noisez + 1)];
 664
 665                         u = orig_u;
 666                         noisex = 0;
 667                         for (i = 0; i != sx; i++) {
 668                                 gradient_buf[index++] = interpolate(
 669                                         v000, v100, v010, v110,
 670                                         v001, v101, v011, v111,
 671                                         u, v, w);
 672
 673                                 u += step_x;
 674                                 if (u >= 1.0) {
 675                                         u -= 1.0;
 676                                         noisex++;
 677                                         v000 = v100;
 678                                         v010 = v110;
 679                                         v100 = noise_buf[idx(noisex + 1, noisey,     noisez)];
 680                                         v110 = noise_buf[idx(noisex + 1, noisey + 1, noisez)];
 681                                         v001 = v101;
 682                                         v011 = v111;
 683                                         v101 = noise_buf[idx(noisex + 1, noisey,     noisez + 1)];
 684                                         v111 = noise_buf[idx(noisex + 1, noisey + 1, noisez + 1)];
 685                                 }
 686                         }
 687
 688                         v += step_y;
 689                         if (v >= 1.0) {
 690                                 v -= 1.0;
 691                                 noisey++;
 692                         }
 693                 }
 694
 695                 w += step_z;
 696                 if (w >= 1.0) {
 697                         w -= 1.0;
 698                         noisez++;
 699                 }
 700         }
 701 }
 702 #undef idx
 703
 704
 705 float *Noise::perlinMap2D(float x, float y, float *persistence_map)
 706 {
 707         float f = 1.0, g = 1.0;
 708         size_t bufsize = sx * sy;
 709
 710         x /= np.spread.X;
 711         y /= np.spread.Y;
 712
 713         memset(result, 0, sizeof(float) * bufsize);
 714
 715         if (persistence_map) {
 716                 if (!persist_buf)
 717                         persist_buf = new float[bufsize];
 718                 for (size_t i = 0; i != bufsize; i++)
 719                         persist_buf[i] = 1.0;
 720         }
 721
 722         for (size_t oct = 0; oct < np.octaves; oct++) {
 723                 gradientMap2D(x * f, y * f,
 724                         f / np.spread.X, f / np.spread.Y,
 725                         seed + np.seed + oct);
 726
 727                 updateResults(g, persist_buf, persistence_map, bufsize);
 728
 729                 f *= np.lacunarity;
 730                 g *= np.persist;
 731         }
 732
 733         if (std::fabs(np.offset - 0.f) > 0.00001 || std::fabs(np.scale - 1.f) > 0.00001) {
 734                 for (size_t i = 0; i != bufsize; i++)
 735                         result[i] = result[i] * np.scale + np.offset;
 736         }
 737
 738         return result;
 739 }
 740
 741
 742 float *Noise::perlinMap3D(float x, float y, float z, float *persistence_map)
 743 {
 744         float f = 1.0, g = 1.0;
 745         size_t bufsize = sx * sy * sz;
 746
 747         x /= np.spread.X;
 748         y /= np.spread.Y;
 749         z /= np.spread.Z;
 750
 751         memset(result, 0, sizeof(float) * bufsize);
 752
 753         if (persistence_map) {
 754                 if (!persist_buf)
 755                         persist_buf = new float[bufsize];
 756                 for (size_t i = 0; i != bufsize; i++)
 757                         persist_buf[i] = 1.0;
 758         }
 759
 760         for (size_t oct = 0; oct < np.octaves; oct++) {
 761                 gradientMap3D(x * f, y * f, z * f,
 762                         f / np.spread.X, f / np.spread.Y, f / np.spread.Z,
 763                         seed + np.seed + oct);
 764
 765                 updateResults(g, persist_buf, persistence_map, bufsize);
 766
 767                 f *= np.lacunarity;
 768                 g *= np.persist;
 769         }
 770
 771         if (std::fabs(np.offset - 0.f) > 0.00001 || std::fabs(np.scale - 1.f) > 0.00001) {
 772                 for (size_t i = 0; i != bufsize; i++)
 773                         result[i] = result[i] * np.scale + np.offset;
 774         }
 775
 776         return result;
 777 }
 778
 779
 780 void Noise::updateResults(float g, float *gmap,
 781         const float *persistence_map, size_t bufsize)
 782 {
 783         // This looks very ugly, but it is 50-70% faster than having
 784         // conditional statements inside the loop
 785         if (np.flags & NOISE_FLAG_ABSVALUE) {
 786                 if (persistence_map) {
 787                         for (size_t i = 0; i != bufsize; i++) {
 788                                 result[i] += gmap[i] * std::fabs(gradient_buf[i]);
 789                                 gmap[i] *= persistence_map[i];
 790                         }
 791                 } else {
 792                         for (size_t i = 0; i != bufsize; i++)
 793                                 result[i] += g * std::fabs(gradient_buf[i]);
 794                 }
 795         } else {
 796                 if (persistence_map) {
 797                         for (size_t i = 0; i != bufsize; i++) {
 798                                 result[i] += gmap[i] * gradient_buf[i];
 799                                 gmap[i] *= persistence_map[i];
 800                         }
 801                 } else {
 802                         for (size_t i = 0; i != bufsize; i++)
 803                                 result[i] += g * gradient_buf[i];
 804                 }
 805         }
 806 }