src/perlin.c

   1 /*
   2 * libtcod 1.3.1
   3 * Copyright (c) 2007,2008 J.C.Wilk
   4 * All rights reserved.
   5 *
   6 * Redistribution and use in source and binary forms, with or without
   7 * modification, are permitted provided that the following conditions are met:
   8 *     * Redistributions of source code must retain the above copyright
   9 *       notice, this list of conditions and the following disclaimer.
  10 *     * Redistributions in binary form must reproduce the above copyright
  11 *       notice, this list of conditions and the following disclaimer in the
  12 *       documentation and/or other materials provided with the distribution.
  13 *     * The name of J.C.Wilk may not be used to endorse or promote products
  14 *       derived from this software without specific prior written permission.
  15 *
  16 * THIS SOFTWARE IS PROVIDED BY J.C.WILK ``AS IS'' AND ANY
  17 * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  18 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  19 * DISCLAIMED. IN NO EVENT SHALL J.C.WILK BE LIABLE FOR ANY
  20 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
  21 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  22 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
  23 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  24 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  25 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  26 */
  27
  28 /**
  29  * @file perlin.c
  30  *
  31  * @brief Handles creating noise based on perlin noise.
  32  *
  33  * Code tries to handle basically 2d/3d cases, without much genericness
  34  *  because it needs to be pretty fast.  Originally sped up the code from
  35  *  about 20 seconds to 8 seconds per Nebula image with the manual loop
  36  *  unrolling.
  37  *
  38  * @note Tried to optimize a while back with SSE and the works, but because
  39  *       of the nature of how it's implemented in non-linear fashion it just
  40  *       wound up complicating the code without actually making it faster.
  41  */
  42
  43
  44 #include "perlin.h"
  45
  46 #include "naev.h"
  47
  48 #include <math.h>
  49 #include <stdlib.h>
  50 #include <string.h>
  51
  52 #include "SDL.h"
  53
  54 #include "log.h"
  55 #include "rng.h"
  56 #include "nfile.h"
  57
  58
  59 #define TCOD_NOISE_MAX_OCTAVES            4 /**< Default octaves for noise. */
  60 #define TCOD_NOISE_DEFAULT_HURST          0.5 /**< Default hurst for noise. */
  61 #define TCOD_NOISE_DEFAULT_LACUNARITY     2. /**< Default lacunarity for noise. */
  62
  63
  64 /**
  65  * @brief Linearly Interpolates x between a and b.
  66  */
  67 #define LERP(a, b, x)      ( a + x * (b - a) )
  68
  69
  70 /**
  71  * @brief Structure used for generating noise.
  72  */
  73 typedef struct perlin_data_s {
  74    int ndim; /**< Dimension of the noise. */
  75    unsigned char map[256]; /**< Randomized map of indexes into buffer */
  76    float buffer[256][3];   /**< Random 256 x 3 buffer */
  77    /* fractal stuff */
  78    float H; /**< Not sure. */
  79    float lacunarity; /**< Not sure. */
  80    float exponent[TCOD_NOISE_MAX_OCTAVES]; /**< Not sure. */
  81 } perlin_data_t; /**< Internal perlin noise data. */
  82
  83
  84 /*
  85  * prototypes
  86  */
  87 /* perlin data handling. */
  88 static perlin_data_t* TCOD_noise_new( int dim, float hurst, float lacunarity );
  89 static void TCOD_noise_delete( perlin_data_t* noise );
  90 /* normalizing. */
  91 static void normalize3( float f[3] );
  92 static void normalize2( float f[2] );
  93 /* noise processing. */
  94 static float lattice3( perlin_data_t *pdata, int ix, float fx,
  95       int iy, float fy, int iz, float fz );
  96 static float lattice2( perlin_data_t *pdata, int ix, float fx, int iy, float fy );
  97 /* basic perlin noise */
  98 static float TCOD_noise_get3( perlin_data_t* pdata, float f[3] );
  99 static float TCOD_noise_get2( perlin_data_t* pdata, float f[2] );
 100 /* turbulence */
 101 static float TCOD_noise_turbulence3( perlin_data_t* noise, float f[3], int octaves );
 102 static float TCOD_noise_turbulence2( perlin_data_t* noise, float f[2], int octaves );
 103
 104
 105 /**
 106  * @brief Not sure what it does.
 107  */
 108 static __inline float lattice3( perlin_data_t *pdata, int ix, float fx, int iy, float fy,
 109 int iz, float fz )
 110 {
 111    int nIndex;
 112    float value;
 113
 114    nIndex = 0;
 115    nIndex = pdata->map[(nIndex + ix) & 0xFF];
 116    nIndex = pdata->map[(nIndex + iy) & 0xFF];
 117    nIndex = pdata->map[(nIndex + iz) & 0xFF];
 118
 119    value  = pdata->buffer[nIndex][0] * fx;
 120    value += pdata->buffer[nIndex][1] * fy;
 121    value += pdata->buffer[nIndex][2] * fz;
 122
 123    return value;
 124 }
 125
 126
 127 /**
 128  * @brief Not sure what it does.
 129  */
 130 static float lattice2( perlin_data_t *pdata, int ix, float fx, int iy, float fy )
 131 {
 132    int nIndex;
 133    float value;
 134
 135    nIndex = 0;
 136    nIndex = pdata->map[(nIndex + ix) & 0xFF];
 137    nIndex = pdata->map[(nIndex + iy) & 0xFF];
 138
 139    value  = pdata->buffer[nIndex][0] * fx;
 140    value += pdata->buffer[nIndex][1] * fy;
 141
 142    return value;
 143 }
 144
 145
 146 #define SWAP(a, b, t)      t = a; a = b; b = t /**< Swaps two values. */
 147 #define FLOOR(a) ((int)a - (a < 0 && a != (int)a)) /**< Limits to 0. */
 148 #define CUBIC(a)  ( a * a * (3 - 2*a) ) /**< Does cubic filtering. */
 149
 150
 151 /**
 152  * @brief Normalizes a 3d vector.
 153  *
 154  *    @param f Vector to normalize.
 155  */
 156 static void normalize3( float f[3] )
 157 {
 158    float magnitude;
 159
 160    magnitude = 1. / sqrtf(f[0]*f[0] + f[1]*f[1] + f[2]*f[2]);
 161    f[0] *= magnitude;
 162    f[1] *= magnitude;
 163    f[2] *= magnitude;
 164 }
 165
 166
 167 /**
 168  * @brief Normalizes a 2d vector.
 169  *
 170  *    @param f Vector to normalize.
 171  */
 172 static void normalize2( float f[2] )
 173 {
 174    float magnitude;
 175
 176    magnitude = 1. / sqrtf(f[0]*f[0] + f[1]*f[1]);
 177    f[0] *= magnitude;
 178    f[1] *= magnitude;
 179 }
 180
 181
 182 /**
 183  * @brief Creates a new perlin noise generator.
 184  *
 185  *    @param dim Dimension of the noise.
 186  *    @param hurst
 187  *    @param lacunarity
 188  */
 189 static perlin_data_t* TCOD_noise_new( int dim, float hurst, float lacunarity )
 190 {
 191    perlin_data_t *pdata;
 192    int i, j;
 193    unsigned char tmp;
 194    float f;
 195
 196    /* Create the data. */
 197    pdata = calloc(sizeof(perlin_data_t),1);
 198    pdata->ndim = dim;
 199
 200    /* Create the buffer and map. */
 201    if (dim == 3) {
 202       for(i=0; i<256; i++) {
 203          pdata->map[i] = (unsigned char)i;
 204          pdata->buffer[i][0] = RNGF()-0.5;
 205          pdata->buffer[i][1] = RNGF()-0.5;
 206          pdata->buffer[i][2] = RNGF()-0.5;
 207          normalize3(pdata->buffer[i]);
 208       }
 209    }
 210    else if (dim == 2) {
 211       for(i=0; i<256; i++) {
 212          pdata->map[i] = (unsigned char)i;
 213          pdata->buffer[i][0] = RNGF()-0.5;
 214          pdata->buffer[i][1] = RNGF()-0.5;
 215          normalize2(pdata->buffer[i]);
 216       }
 217    }
 218    else {
 219       i = 0;
 220    }
 221
 222    while(--i) {
 223       j = RNG(0, 255);
 224       SWAP(pdata->map[i], pdata->map[j], tmp);
 225    }
 226
 227    f = 1.;
 228    pdata->H = hurst;
 229    pdata->lacunarity = lacunarity;
 230    for(i=0; i<TCOD_NOISE_MAX_OCTAVES; i++) {
 231       /*exponent[i] = powf(f, -H); */
 232       pdata->exponent[i] = 1. / f;
 233       f *= lacunarity;
 234    }
 235
 236    return pdata;
 237 }
 238
 239
 240 /**
 241  * @brief Gets some 3D Perlin noise from the data.
 242  *
 243  * Somewhat optimized for speed, probably can't get optimized much more.
 244  *
 245  *    @param pdata Perlin data to use.
 246  *    @param f Position of the noise to get.
 247  */
 248 static float TCOD_noise_get3( perlin_data_t* pdata, float f[3] )
 249 {
 250    int n[3] __attribute__ ((aligned (32))); /* Indexes to pass to lattice function */
 251    float r[3] __attribute__ ((aligned (32))); /* Remainders to pass to lattice function */
 252    float w[3] __attribute__ ((aligned (32))); /* Cubic values to pass to interpolation function */
 253    float value;
 254    float v[8] __attribute__ ((aligned (32)));
 255
 256    n[0] = (int)f[0];
 257    n[1] = (int)f[1];
 258    n[2] = (int)f[2];
 259
 260    r[0] = f[0] - n[0];
 261    r[1] = f[1] - n[1];
 262    r[2] = f[2] - n[2];
 263
 264    w[0] = CUBIC(r[0]);
 265    w[1] = CUBIC(r[1]);
 266    w[2] = CUBIC(r[2]);
 267
 268    /*
 269     * This is the big ugly bit in dire need of optimization
 270     */
 271    v[0] = lattice3(pdata, n[0],   r[0],   n[1],   r[1],   n[2],   r[2]);
 272    v[1] = lattice3(pdata, n[0]+1, r[0]-1, n[1],   r[1],   n[2],   r[2]);
 273    v[2] = lattice3(pdata, n[0],   r[0],   n[1]+1, r[1]-1, n[2],   r[2]);
 274    v[3] = lattice3(pdata, n[0]+1, r[0]-1, n[1]+1, r[1]-1, n[2],   r[2]);
 275    v[4] = lattice3(pdata, n[0],   r[0],   n[1],   r[1],   n[2]+1, r[2]-1);
 276    v[5] = lattice3(pdata, n[0]+1, r[0]-1, n[1],   r[1],   n[2]+1, r[2]-1);
 277    v[6] = lattice3(pdata, n[0],   r[0],   n[1]+1, r[1]-1, n[2]+1, r[2]-1);
 278    v[7] = lattice3(pdata, n[0]+1, r[0]-1, n[1]+1, r[1]-1, n[2]+1, r[2]-1);
 279    value = LERP(
 280          LERP(
 281             LERP(v[0], v[1], w[0]),
 282             LERP(v[2], v[3], w[0]),
 283             w[1]
 284             ),
 285          LERP(
 286             LERP(v[4], v[5], w[0]),
 287             LERP(v[6], v[7], w[0]),
 288             w[1]
 289             ),
 290          w[2]
 291          );
 292
 293    return CLAMP(-0.99999f, 0.99999f, value);
 294 }
 295
 296
 297 /**
 298  * @brief Gets some 2D Perlin noise from the data.
 299  *
 300  * Somewhat optimized for speed, probably can't get optimized much more.
 301  *
 302  *    @param pdata Perlin data to use.
 303  *    @param f Position of the noise to get.
 304  */
 305 static float TCOD_noise_get2( perlin_data_t* pdata, float f[2] )
 306 {
 307    int n[2] __attribute__ ((aligned (32))); /* Indexes to pass to lattice function */
 308    float r[2] __attribute__ ((aligned (32))); /* Remainders to pass to lattice function */
 309    float w[2] __attribute__ ((aligned (32))); /* Cubic values to pass to interpolation function */
 310    float value __attribute__ ((aligned (32)));
 311    float v[4] __attribute__ ((aligned (32)));
 312
 313    n[0] = FLOOR(f[0]);
 314    n[1] = FLOOR(f[1]);
 315
 316    r[0] = f[0] - n[0];
 317    r[1] = f[1] - n[1];
 318
 319    w[0] = CUBIC(r[0]);
 320    w[1] = CUBIC(r[1]);
 321
 322    /*
 323     * Much faster in 2d.
 324     */
 325    v[0] = lattice2(pdata,n[0],   r[0],   n[1],   r[1]);
 326    v[1] = lattice2(pdata,n[0]+1, r[0]-1, n[1],   r[1]);
 327    v[2] = lattice2(pdata,n[0],   r[0],   n[1]+1, r[1]-1);
 328    v[3] = lattice2(pdata,n[0]+1, r[0]-1, n[1]+1, r[1]-1);
 329    value = LERP(LERP(v[0], v[1], w[0]),
 330          LERP(v[2], v[3], w[0]),
 331          w[1]);
 332
 333    return CLAMP(-0.99999f, 0.99999f, value);
 334 }
 335
 336
 337 /**
 338  * @brief Gets 3d Turbulence noise for a position.
 339  *
 340  *    @param noise Perlin data to generate noise from.
 341  *    @param f Position of the noise.
 342  *    @param octaves Octaves to use.
 343  *    @return The noise level at the position.
 344  */
 345 static float TCOD_noise_turbulence3( perlin_data_t* noise, float f[3], int octaves )
 346 {
 347    float tf[3];
 348    perlin_data_t *pdata=(perlin_data_t *)noise;
 349    /* Initialize locals */
 350    float value = 0;
 351    int i;
 352
 353    tf[0] = f[0];
 354    tf[1] = f[1];
 355    tf[2] = f[2];
 356
 357    /* Inner loop of spectral construction, where the fractal is built */
 358    for(i=0; i<octaves; i++)
 359    {
 360       value += ABS(TCOD_noise_get3(noise,tf)) * pdata->exponent[i];
 361       tf[0] *= pdata->lacunarity;
 362       tf[1] *= pdata->lacunarity;
 363       tf[2] *= pdata->lacunarity;
 364    }
 365
 366    return CLAMP(-0.99999f, 0.99999f, value);
 367 }
 368
 369
 370 /**
 371  * @brief Gets 2d Turbulence noise for a position.
 372  *
 373  *    @param noise Perlin data to generate noise from.
 374  *    @param f Position of the noise.
 375  *    @param octaves Octaves to use.
 376  *    @return The noise level at the position.
 377  */
 378 static float TCOD_noise_turbulence2( perlin_data_t* noise, float f[2], int octaves )
 379 {
 380    float tf[2];
 381    perlin_data_t *pdata=(perlin_data_t *)noise;
 382    /* Initialize locals */
 383    float value = 0;
 384    int i;
 385
 386    tf[0] = f[0];
 387    tf[1] = f[1];
 388
 389    /* Inner loop of spectral construction, where the fractal is built */
 390    for(i=0; i<octaves; i++)
 391    {
 392       value += ABS(TCOD_noise_get2(noise,tf)) * pdata->exponent[i];
 393       tf[0] *= pdata->lacunarity;
 394       tf[1] *= pdata->lacunarity;
 395    }
 396
 397    return CLAMP(-0.99999f, 0.99999f, value);
 398 }
 399
 400
 401 /**
 402  * @brief Frees some noise data.
 403  *
 404  *    @param noise Noise data to free.
 405  */
 406 void TCOD_noise_delete( perlin_data_t* noise )
 407 {
 408    free(noise);
 409 }
 410
 411
 412 /**
 413  * @brief Generates radar interference.
 414  *
 415  *    @param w Width to generate.
 416  *    @param h Height to generate.
 417  *    @param rug Rugosity of the interference.
 418  *    @return The map generated.
 419  */
 420 float* noise_genRadarInt( const int w, const int h, float rug )
 421 {
 422    int x, y;
 423    float f[2];
 424    int octaves;
 425    float hurst;
 426    float lacunarity;
 427    perlin_data_t* noise;
 428    float *map;
 429    float value;
 430
 431    /* pretty default values */
 432    octaves     = 3;
 433    hurst       = TCOD_NOISE_DEFAULT_HURST;
 434    lacunarity  = TCOD_NOISE_DEFAULT_LACUNARITY;
 435
 436    /* create noise and data */
 437    noise       = TCOD_noise_new( 2, hurst, lacunarity );
 438    map         = malloc(sizeof(float)*w*h);
 439    if (map == NULL) {
 440       WARN("Out of memory!");
 441       return NULL;
 442    }
 443
 444    /* Start to create the nebula */
 445    for (y=0; y<h; y++) {
 446
 447       f[1] = rug * (float)y / (float)h;
 448       for (x=0; x<w; x++) {
 449
 450          f[0] = rug * (float)x / (float)w;
 451
 452          /* Get the 2d noise. */
 453          value = TCOD_noise_get2( noise, f );
 454
 455          /* Set the value to [0,1]. */
 456          map[y*w + x] = (value + 1.) / 2.;
 457       }
 458    }
 459
 460    /* Clean up */
 461    TCOD_noise_delete( noise );
 462
 463    /* Results */
 464    return map;
 465 }
 466
 467
 468 /**
 469  * @brief Generates a 3d nebula map.
 470  *
 471  *    @param w Width of the map.
 472  *    @param h Height of the map.
 473  *    @param n Number of slices of the map (2d planes).
 474  *    @param rug Rugosity of the map.
 475  *    @return The map generated.
 476  */
 477 float* noise_genNebulaMap( const int w, const int h, const int n, float rug )
 478 {
 479    int x, y, z;
 480    float f[3];
 481    int octaves;
 482    float hurst;
 483    float lacunarity;
 484    perlin_data_t* noise;
 485    float *nebula;
 486    float value;
 487    float zoom;
 488    float max;
 489    unsigned int *t, s;
 490
 491    /* pretty default values */
 492    octaves     = 3;
 493    hurst       = TCOD_NOISE_DEFAULT_HURST;
 494    lacunarity  = TCOD_NOISE_DEFAULT_LACUNARITY;
 495    zoom        = rug * ((float)h/768.)*((float)w/1024.);
 496
 497    /* create noise and data */
 498    noise       = TCOD_noise_new( 3, hurst, lacunarity );
 499    nebula     = malloc(sizeof(float)*w*h*n);
 500    if (nebula == NULL) {
 501       WARN("Out of memory!");
 502       return NULL;
 503    }
 504
 505    /* Some debug information and time setting */
 506    s = SDL_GetTicks();
 507    t = malloc(sizeof(unsigned int)*n);
 508    DEBUG("Generating Nebula of size %dx%dx%d", w, h, n);
 509
 510    /* Start to create the nebula */
 511    max = 0.;
 512    for (z=0; z<n; z++) {
 513
 514       f[2] = zoom * (float)z / (float)n;
 515
 516       for (y=0; y<h; y++) {
 517
 518          f[1] = zoom * (float)y / (float)h;
 519
 520          for (x=0; x<w; x++) {
 521
 522             f[0] = zoom * (float)x / (float)w;
 523
 524             value = TCOD_noise_turbulence3( noise, f, octaves );
 525             if (max < value) max = value;
 526
 527             nebula[z*w*h + y*w + x] = value;
 528          }
 529       }
 530
 531       /* More time magic debug */
 532       t[z] = SDL_GetTicks();
 533       DEBUG("   Layer %d/%d generated in %d ms", z+1, n,
 534             (z>0) ? t[z] - t[z-1] : t[z] - s );
 535    }
 536
 537    /* Post filtering */
 538    value = 1. - max;
 539    for (z=0; z<n; z++)
 540       for (y=0; y<h; y++)
 541          for (x=0; x<w; x++)
 542             nebula[z*w*h + y*w + x] += value;
 543
 544    /* Clean up */
 545    TCOD_noise_delete( noise );
 546
 547    /* Results */
 548    DEBUG("Nebula Generated in %d ms", SDL_GetTicks() - s );
 549    return nebula;
 550 }
 551
 552
 553 /**
 554  * @brief Generates tiny nebula puffs
 555  *
 556  *    @param w Width of the puff to generate.
 557  *    @param h Height of the puff to generate.
 558  *    @param rug Rugosity of the puff.
 559  *    @return The puff generated.
 560  */
 561 float* noise_genNebulaPuffMap( const int w, const int h, float rug )
 562 {
 563    int x,y, hw,hh;
 564    float d;
 565    float f[2];
 566    int octaves;
 567    float hurst;
 568    float lacunarity;
 569    perlin_data_t* noise;
 570    float *nebula;
 571    float value;
 572    float zoom;
 573    float max;
 574
 575    /* pretty default values */
 576    octaves     = 3;
 577    hurst       = TCOD_NOISE_DEFAULT_HURST;
 578    lacunarity  = TCOD_NOISE_DEFAULT_LACUNARITY;
 579    zoom        = rug;
 580
 581    /* create noise and data */
 582    noise       = TCOD_noise_new( 2, hurst, lacunarity );
 583    nebula     = malloc(sizeof(float)*w*h);
 584    if (nebula == NULL) {
 585       WARN("Out of memory!");
 586       return NULL;
 587    }
 588
 589    /* Start to create the nebula */
 590    max   = 0.;
 591    hw    = w/2;
 592    hh    = h/2;
 593    d     = (float)MIN(hw,hh);
 594    for (y=0; y<h; y++) {
 595
 596       f[1] = zoom * (float)y / (float)h;
 597       for (x=0; x<w; x++) {
 598
 599          f[0] = zoom * (float)x / (float)w;
 600
 601          /* Get the 2d noise. */
 602          value = TCOD_noise_turbulence2( noise, f, octaves );
 603
 604          /* Make value also depend on distance from center */
 605          value *= (d - 1. - sqrtf( (float)((x-hw)*(x-hw) + (y-hh)*(y-hh)) )) / d;
 606          if (value < 0.)
 607             value = 0.;
 608
 609          /* Cap at maximum. */
 610          if (max < value)
 611             max = value;
 612
 613          /* Set the value. */
 614          nebula[y*w + x] = value;
 615       }
 616    }
 617
 618    /* Clean up */
 619    TCOD_noise_delete( noise );
 620
 621    /* Results */
 622    return nebula;
 623 }
 624
 625