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