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