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