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