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