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 perlin stuff ] ////////////////////////////
 318
 319
 320 Noise::Noise(NoiseParams *np_, int seed, int sx, int sy, int sz)
 321 {
 322         memcpy(&np, np_, sizeof(np));
 323         this->seed = seed;
 324         this->sx   = sx;
 325         this->sy   = sy;
 326         this->sz   = sz;
 327
 328         this->persist_buf  = NULL;
 329         this->gradient_buf = NULL;
 330         this->result       = NULL;
 331
 332         if (np.flags & NOISE_FLAG_DEFAULTS) {
 333                 // By default, only 2d noise is eased.
 334                 if (sz <= 1)
 335                         np.flags |= NOISE_FLAG_EASED;
 336         }
 337
 338         allocBuffers();
 339 }
 340
 341
 342 Noise::~Noise()
 343 {
 344         delete[] gradient_buf;
 345         delete[] persist_buf;
 346         delete[] noise_buf;
 347         delete[] result;
 348 }
 349
 350
 351 void Noise::allocBuffers()
 352 {
 353         this->noise_buf = NULL;
 354         resizeNoiseBuf(sz > 1);
 355
 356         delete[] gradient_buf;
 357         delete[] persist_buf;
 358         delete[] result;
 359
 360         try {
 361                 size_t bufsize = sx * sy * sz;
 362                 this->persist_buf  = NULL;
 363                 this->gradient_buf = new float[bufsize];
 364                 this->result       = new float[bufsize];
 365         } catch (std::bad_alloc &e) {
 366                 throw InvalidNoiseParamsException();
 367         }
 368 }
 369
 370
 371 void Noise::setSize(int sx, int sy, int sz)
 372 {
 373         this->sx = sx;
 374         this->sy = sy;
 375         this->sz = sz;
 376
 377         allocBuffers();
 378 }
 379
 380
 381 void Noise::setSpreadFactor(v3f spread)
 382 {
 383         this->np.spread = spread;
 384
 385         resizeNoiseBuf(sz > 1);
 386 }
 387
 388
 389 void Noise::setOctaves(int octaves)
 390 {
 391         this->np.octaves = octaves;
 392
 393         resizeNoiseBuf(sz > 1);
 394 }
 395
 396
 397 void Noise::resizeNoiseBuf(bool is3d)
 398 {
 399         int nlx, nly, nlz;
 400         float ofactor;
 401
 402         //maximum possible spread value factor
 403         ofactor = pow(np.lacunarity, np.octaves - 1);
 404
 405         //noise lattice point count
 406         //(int)(sz * spread * ofactor) is # of lattice points crossed due to length
 407         // + 2 for the two initial endpoints
 408         // + 1 for potentially crossing a boundary due to offset
 409         nlx = (int)ceil(sx * ofactor / np.spread.X) + 3;
 410         nly = (int)ceil(sy * ofactor / np.spread.Y) + 3;
 411         nlz = is3d ? (int)ceil(sz * ofactor / np.spread.Z) + 3 : 1;
 412
 413         delete[] noise_buf;
 414         try {
 415                 noise_buf = new float[nlx * nly * nlz];
 416         } catch (std::bad_alloc &e) {
 417                 throw InvalidNoiseParamsException();
 418         }
 419 }
 420
 421
 422 /*
 423  * NB:  This algorithm is not optimal in terms of space complexity.  The entire
 424  * integer lattice of noise points could be done as 2 lines instead, and for 3D,
 425  * 2 lines + 2 planes.
 426  * However, this would require the noise calls to be interposed with the
 427  * interpolation loops, which may trash the icache, leading to lower overall
 428  * performance.
 429  * Another optimization that could save half as many noise calls is to carry over
 430  * values from the previous noise lattice as midpoints in the new lattice for the
 431  * next octave.
 432  */
 433 #define idx(x, y) ((y) * nlx + (x))
 434 void Noise::gradientMap2D(
 435                 float x, float y,
 436                 float step_x, float step_y,
 437                 int seed)
 438 {
 439         float v00, v01, v10, v11, u, v, orig_u;
 440         int index, i, j, x0, y0, noisex, noisey;
 441         int nlx, nly;
 442
 443         Interp2dFxn interpolate = (np.flags & NOISE_FLAG_EASED) ?
 444                 biLinearInterpolation : biLinearInterpolationNoEase;
 445
 446         x0 = floor(x);
 447         y0 = floor(y);
 448         u = x - (float)x0;
 449         v = y - (float)y0;
 450         orig_u = u;
 451
 452         //calculate noise point lattice
 453         nlx = (int)(u + sx * step_x) + 2;
 454         nly = (int)(v + sy * step_y) + 2;
 455         index = 0;
 456         for (j = 0; j != nly; j++)
 457                 for (i = 0; i != nlx; i++)
 458                         noise_buf[index++] = noise2d(x0 + i, y0 + j, seed);
 459
 460         //calculate interpolations
 461         index  = 0;
 462         noisey = 0;
 463         for (j = 0; j != sy; j++) {
 464                 v00 = noise_buf[idx(0, noisey)];
 465                 v10 = noise_buf[idx(1, noisey)];
 466                 v01 = noise_buf[idx(0, noisey + 1)];
 467                 v11 = noise_buf[idx(1, noisey + 1)];
 468
 469                 u = orig_u;
 470                 noisex = 0;
 471                 for (i = 0; i != sx; i++) {
 472                         gradient_buf[index++] = interpolate(v00, v10, v01, v11, u, v);
 473
 474                         u += step_x;
 475                         if (u >= 1.0) {
 476                                 u -= 1.0;
 477                                 noisex++;
 478                                 v00 = v10;
 479                                 v01 = v11;
 480                                 v10 = noise_buf[idx(noisex + 1, noisey)];
 481                                 v11 = noise_buf[idx(noisex + 1, noisey + 1)];
 482                         }
 483                 }
 484
 485                 v += step_y;
 486                 if (v >= 1.0) {
 487                         v -= 1.0;
 488                         noisey++;
 489                 }
 490         }
 491 }
 492 #undef idx
 493
 494
 495 #define idx(x, y, z) ((z) * nly * nlx + (y) * nlx + (x))
 496 void Noise::gradientMap3D(
 497                 float x, float y, float z,
 498                 float step_x, float step_y, float step_z,
 499                 int seed)
 500 {
 501         float v000, v010, v100, v110;
 502         float v001, v011, v101, v111;
 503         float u, v, w, orig_u, orig_v;
 504         int index, i, j, k, x0, y0, z0, noisex, noisey, noisez;
 505         int nlx, nly, nlz;
 506
 507         Interp3dFxn interpolate = (np.flags & NOISE_FLAG_EASED) ?
 508                 triLinearInterpolation : triLinearInterpolationNoEase;
 509
 510         x0 = floor(x);
 511         y0 = floor(y);
 512         z0 = floor(z);
 513         u = x - (float)x0;
 514         v = y - (float)y0;
 515         w = z - (float)z0;
 516         orig_u = u;
 517         orig_v = v;
 518
 519         //calculate noise point lattice
 520         nlx = (int)(u + sx * step_x) + 2;
 521         nly = (int)(v + sy * step_y) + 2;
 522         nlz = (int)(w + sz * step_z) + 2;
 523         index = 0;
 524         for (k = 0; k != nlz; k++)
 525                 for (j = 0; j != nly; j++)
 526                         for (i = 0; i != nlx; i++)
 527                                 noise_buf[index++] = noise3d(x0 + i, y0 + j, z0 + k, seed);
 528
 529         //calculate interpolations
 530         index  = 0;
 531         noisey = 0;
 532         noisez = 0;
 533         for (k = 0; k != sz; k++) {
 534                 v = orig_v;
 535                 noisey = 0;
 536                 for (j = 0; j != sy; j++) {
 537                         v000 = noise_buf[idx(0, noisey,     noisez)];
 538                         v100 = noise_buf[idx(1, noisey,     noisez)];
 539                         v010 = noise_buf[idx(0, noisey + 1, noisez)];
 540                         v110 = noise_buf[idx(1, noisey + 1, noisez)];
 541                         v001 = noise_buf[idx(0, noisey,     noisez + 1)];
 542                         v101 = noise_buf[idx(1, noisey,     noisez + 1)];
 543                         v011 = noise_buf[idx(0, noisey + 1, noisez + 1)];
 544                         v111 = noise_buf[idx(1, noisey + 1, noisez + 1)];
 545
 546                         u = orig_u;
 547                         noisex = 0;
 548                         for (i = 0; i != sx; i++) {
 549                                 gradient_buf[index++] = interpolate(
 550                                         v000, v100, v010, v110,
 551                                         v001, v101, v011, v111,
 552                                         u, v, w);
 553
 554                                 u += step_x;
 555                                 if (u >= 1.0) {
 556                                         u -= 1.0;
 557                                         noisex++;
 558                                         v000 = v100;
 559                                         v010 = v110;
 560                                         v100 = noise_buf[idx(noisex + 1, noisey,     noisez)];
 561                                         v110 = noise_buf[idx(noisex + 1, noisey + 1, noisez)];
 562                                         v001 = v101;
 563                                         v011 = v111;
 564                                         v101 = noise_buf[idx(noisex + 1, noisey,     noisez + 1)];
 565                                         v111 = noise_buf[idx(noisex + 1, noisey + 1, noisez + 1)];
 566                                 }
 567                         }
 568
 569                         v += step_y;
 570                         if (v >= 1.0) {
 571                                 v -= 1.0;
 572                                 noisey++;
 573                         }
 574                 }
 575
 576                 w += step_z;
 577                 if (w >= 1.0) {
 578                         w -= 1.0;
 579                         noisez++;
 580                 }
 581         }
 582 }
 583 #undef idx
 584
 585
 586 float *Noise::perlinMap2D(float x, float y, float *persistence_map)
 587 {
 588         float f = 1.0, g = 1.0;
 589         size_t bufsize = sx * sy;
 590
 591         x /= np.spread.X;
 592         y /= np.spread.Y;
 593
 594         memset(result, 0, sizeof(float) * bufsize);
 595
 596         if (persistence_map) {
 597                 if (!persist_buf)
 598                         persist_buf = new float[bufsize];
 599                 for (size_t i = 0; i != bufsize; i++)
 600                         persist_buf[i] = 1.0;
 601         }
 602
 603         for (size_t oct = 0; oct < np.octaves; oct++) {
 604                 gradientMap2D(x * f, y * f,
 605                         f / np.spread.X, f / np.spread.Y,
 606                         seed + np.seed + oct);
 607
 608                 updateResults(g, persist_buf, persistence_map, bufsize);
 609
 610                 f *= np.lacunarity;
 611                 g *= np.persist;
 612         }
 613
 614         return result;
 615 }
 616
 617
 618 float *Noise::perlinMap3D(float x, float y, float z, float *persistence_map)
 619 {
 620         float f = 1.0, g = 1.0;
 621         size_t bufsize = sx * sy * sz;
 622
 623         x /= np.spread.X;
 624         y /= np.spread.Y;
 625         z /= np.spread.Z;
 626
 627         memset(result, 0, sizeof(float) * bufsize);
 628
 629         if (persistence_map) {
 630                 if (!persist_buf)
 631                         persist_buf = new float[bufsize];
 632                 for (size_t i = 0; i != bufsize; i++)
 633                         persist_buf[i] = 1.0;
 634         }
 635
 636         for (size_t oct = 0; oct < np.octaves; oct++) {
 637                 gradientMap3D(x * f, y * f, z * f,
 638                         f / np.spread.X, f / np.spread.Y, f / np.spread.Z,
 639                         seed + np.seed + oct);
 640
 641                 updateResults(g, persist_buf, persistence_map, bufsize);
 642
 643                 f *= np.lacunarity;
 644                 g *= np.persist;
 645         }
 646
 647         return result;
 648 }
 649
 650
 651 void Noise::updateResults(float g, float *gmap,
 652         float *persistence_map, size_t bufsize)
 653 {
 654         // This looks very ugly, but it is 50-70% faster than having
 655         // conditional statements inside the loop
 656         if (np.flags & NOISE_FLAG_ABSVALUE) {
 657                 if (persistence_map) {
 658                         for (size_t i = 0; i != bufsize; i++) {
 659                                 result[i] += gmap[i] * fabs(gradient_buf[i]);
 660                                 gmap[i] *= persistence_map[i];
 661                         }
 662                 } else {
 663                         for (size_t i = 0; i != bufsize; i++)
 664                                 result[i] += g * fabs(gradient_buf[i]);
 665                 }
 666         } else {
 667                 if (persistence_map) {
 668                         for (size_t i = 0; i != bufsize; i++) {
 669                                 result[i] += gmap[i] * gradient_buf[i];
 670                                 gmap[i] *= persistence_map[i];
 671                         }
 672                 } else {
 673                         for (size_t i = 0; i != bufsize; i++)
 674                                 result[i] += g * gradient_buf[i];
 675                 }
 676         }
 677 }
 678
 679
 680 void Noise::transformNoiseMap()
 681 {
 682         // Because sx, sy, and sz are object members whose values may conceivably be
 683         // modified in other threads. gcc (at least) will consider the buffer size
 684         // computation as invalidated between loop comparisons, resulting in a ~2x
 685         // slowdown even with -O2.  To prevent this, store the value in a local.
 686         size_t bufsize = sx * sy * sz;
 687         for (size_t i = 0; i != bufsize; i++)
 688                 result[i] = result[i] * np.scale + np.offset;
 689 }
 690