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