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