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