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 "nstring.h"
51 
52 #include "SDL.h"
53 #include "SDL_thread.h"
54 #include "threadpool.h"
55 
56 #include "log.h"
57 #include "rng.h"
58 #include "nfile.h"
59 
60 
61 #define SIMPLEX_SCALE 0.5f
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 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[NOISE_MAX_OCTAVES]; /**< Not sure. */
81 };
82 
83 
84 /**
85  * @brief Threading stuff.
86  */
87 typedef struct thread_args_ {
88    int z; /**< Z level working on. */
89    float zoom; /**< Zoom level of detail. */
90    int n; /**< Number of layers to generate. */
91    int h; /**< Height. */
92    int w; /**< Width. */
93    perlin_data_t *noise; /**< Parent noise. */
94    int octaves; /**< Octave parameters. */
95    float *max; /**< Maximum value. */
96    float *nebula; /**< Nebula loading into. */
97 } thread_args;
98 
99 
100 /*
101  * prototypes
102  */
103 /* normalizing. */
104 static void normalize3( float f[3] );
105 static void normalize2( float f[2] );
106 /* noise processing. */
107 static float lattice3( perlin_data_t *pdata, int ix, float fx,
108       int iy, float fy, int iz, float fz );
109 static float lattice2( perlin_data_t *pdata, int ix, float fx, int iy, float fy );
110 static float lattice1( perlin_data_t *pdata, int ix, float fx );
111 /*Threading */
112 static int noise_genNebulaMap_thread( void *data );
113 
114 
115 /**
116  * @brief Not sure what it does.
117  */
lattice3(perlin_data_t * pdata,int ix,float fx,int iy,float fy,int iz,float fz)118 static float lattice3( perlin_data_t *pdata, int ix, float fx, int iy, float fy,
119 int iz, float fz )
120 {
121    int nIndex;
122    float value;
123 
124    nIndex = 0;
125    nIndex = pdata->map[(nIndex + ix) & 0xFF];
126    nIndex = pdata->map[(nIndex + iy) & 0xFF];
127    nIndex = pdata->map[(nIndex + iz) & 0xFF];
128 
129    value  = pdata->buffer[nIndex][0] * fx;
130    value += pdata->buffer[nIndex][1] * fy;
131    value += pdata->buffer[nIndex][2] * fz;
132 
133    return value;
134 }
135 
136 
137 /**
138  * @brief Not sure what it does.
139  */
lattice2(perlin_data_t * pdata,int ix,float fx,int iy,float fy)140 static float lattice2( perlin_data_t *pdata, int ix, float fx, int iy, float fy )
141 {
142    int nIndex;
143    float value;
144 
145    nIndex = 0;
146    nIndex = pdata->map[(nIndex + ix) & 0xFF];
147    nIndex = pdata->map[(nIndex + iy) & 0xFF];
148 
149    value  = pdata->buffer[nIndex][0] * fx;
150    value += pdata->buffer[nIndex][1] * fy;
151 
152    return value;
153 }
154 
155 
156 /**
157  * @brief Not sure what it does.
158  */
lattice1(perlin_data_t * pdata,int ix,float fx)159 static float lattice1( perlin_data_t *pdata, int ix, float fx )
160 {
161    int nIndex;
162    float value;
163 
164    nIndex = 0;
165    nIndex = pdata->map[(nIndex + ix) & 0xFF];
166 
167    value  = pdata->buffer[nIndex][0] * fx;
168 
169    return value;
170 }
171 
172 
173 #define SWAP(a, b, t)      (t) = (a); (a) = (b); (b) = (t) /**< Swaps two values. */
174 #define FLOOR(a)           ((int)(a) - ((a) < 0 && (a) != (int)(a))) /**< Limits to 0. */
175 #define CUBIC(a)           ( (a) * (a) * (3 - 2*(a)) ) /**< Does cubic filtering. */
176 
177 
178 /**
179  * @brief Normalizes a 3d vector.
180  *
181  *    @param f Vector to normalize.
182  */
normalize3(float f[3])183 static void normalize3( float f[3] )
184 {
185    float magnitude;
186 
187    magnitude = 1. / sqrtf(f[0]*f[0] + f[1]*f[1] + f[2]*f[2]);
188    f[0] *= magnitude;
189    f[1] *= magnitude;
190    f[2] *= magnitude;
191 }
192 
193 
194 /**
195  * @brief Normalizes a 2d vector.
196  *
197  *    @param f Vector to normalize.
198  */
normalize2(float f[2])199 static void normalize2( float f[2] )
200 {
201    float magnitude;
202 
203    magnitude = 1. / sqrtf(f[0]*f[0] + f[1]*f[1]);
204    f[0] *= magnitude;
205    f[1] *= magnitude;
206 }
207 
208 
209 /**
210  * @brief Creates a new perlin noise generator.
211  *
212  *    @param dim Dimension of the noise.
213  *    @param hurst
214  *    @param lacunarity
215  */
noise_new(int dim,float hurst,float lacunarity)216 perlin_data_t* noise_new( int dim, float hurst, float lacunarity )
217 {
218    perlin_data_t *pdata;
219    int i, j;
220    unsigned char tmp;
221    float f;
222 
223    /* Create the data. */
224    pdata = calloc(sizeof(perlin_data_t),1);
225    pdata->ndim = dim;
226 
227    /* Create the buffer and map. */
228    if (dim == 3) {
229       for(i=0; i<256; i++) {
230          pdata->map[i] = (unsigned char)i;
231          pdata->buffer[i][0] = RNGF()-0.5;
232          pdata->buffer[i][1] = RNGF()-0.5;
233          pdata->buffer[i][2] = RNGF()-0.5;
234          normalize3(pdata->buffer[i]);
235       }
236    }
237    else if (dim == 2) {
238       for(i=0; i<256; i++) {
239          pdata->map[i] = (unsigned char)i;
240          pdata->buffer[i][0] = RNGF()-0.5;
241          pdata->buffer[i][1] = RNGF()-0.5;
242          normalize2(pdata->buffer[i]);
243       }
244    }
245    else {
246       for (i=0; i<256; i++) {
247          pdata->map[i] = (unsigned char)i;
248          pdata->buffer[i][0] = 1.;
249       }
250    }
251 
252    while (--i) {
253       j = RNG(0, 255);
254       SWAP(pdata->map[i], pdata->map[j], tmp);
255    }
256 
257    f = 1.;
258    pdata->H = hurst;
259    pdata->lacunarity = lacunarity;
260    for(i=0; i<NOISE_MAX_OCTAVES; i++) {
261       /*exponent[i] = powf(f, -H); */
262       pdata->exponent[i] = 1. / f;
263       f *= lacunarity;
264    }
265 
266    return pdata;
267 }
268 
269 
270 /**
271  * @brief Gets some 3D Perlin noise from the data.
272  *
273  * Somewhat optimized for speed, probably can't get optimized much more.
274  *
275  *    @param pdata Perlin data to use.
276  *    @param f Position of the noise to get.
277  */
noise_get3(perlin_data_t * pdata,float f[3])278 float noise_get3( perlin_data_t* pdata, float f[3] )
279 {
280    int n[3] __attribute__ ((aligned (32))); /* Indexes to pass to lattice function */
281    float r[3] __attribute__ ((aligned (32))); /* Remainders to pass to lattice function */
282    float w[3] __attribute__ ((aligned (32))); /* Cubic values to pass to interpolation function */
283    float value;
284    float v[8] __attribute__ ((aligned (32)));
285 
286    n[0] = (int)f[0];
287    n[1] = (int)f[1];
288    n[2] = (int)f[2];
289 
290    r[0] = f[0] - n[0];
291    r[1] = f[1] - n[1];
292    r[2] = f[2] - n[2];
293 
294    w[0] = CUBIC(r[0]);
295    w[1] = CUBIC(r[1]);
296    w[2] = CUBIC(r[2]);
297 
298    /*
299     * This is the big ugly bit in dire need of optimization
300     */
301    v[0] = lattice3(pdata, n[0],   r[0],   n[1],   r[1],   n[2],   r[2]);
302    v[1] = lattice3(pdata, n[0]+1, r[0]-1, n[1],   r[1],   n[2],   r[2]);
303    v[2] = lattice3(pdata, n[0],   r[0],   n[1]+1, r[1]-1, n[2],   r[2]);
304    v[3] = lattice3(pdata, n[0]+1, r[0]-1, n[1]+1, r[1]-1, n[2],   r[2]);
305    v[4] = lattice3(pdata, n[0],   r[0],   n[1],   r[1],   n[2]+1, r[2]-1);
306    v[5] = lattice3(pdata, n[0]+1, r[0]-1, n[1],   r[1],   n[2]+1, r[2]-1);
307    v[6] = lattice3(pdata, n[0],   r[0],   n[1]+1, r[1]-1, n[2]+1, r[2]-1);
308    v[7] = lattice3(pdata, n[0]+1, r[0]-1, n[1]+1, r[1]-1, n[2]+1, r[2]-1);
309    value = LERP(
310          LERP(
311             LERP(v[0], v[1], w[0]),
312             LERP(v[2], v[3], w[0]),
313             w[1]
314             ),
315          LERP(
316             LERP(v[4], v[5], w[0]),
317             LERP(v[6], v[7], w[0]),
318             w[1]
319             ),
320          w[2]
321          );
322 
323    return CLAMP(-0.99999f, 0.99999f, value);
324 }
325 
326 
327 /**
328  * @brief Gets some 2D Perlin noise from the data.
329  *
330  * Somewhat optimized for speed, probably can't get optimized much more.
331  *
332  *    @param pdata Perlin data to use.
333  *    @param f Position of the noise to get.
334  */
noise_get2(perlin_data_t * pdata,float f[2])335 float noise_get2( perlin_data_t* pdata, float f[2] )
336 {
337    int n[2] __attribute__ ((aligned (32))); /* Indexes to pass to lattice function */
338    float r[2] __attribute__ ((aligned (32))); /* Remainders to pass to lattice function */
339    float w[2] __attribute__ ((aligned (32))); /* Cubic values to pass to interpolation function */
340    float value __attribute__ ((aligned (32)));
341    float v[4] __attribute__ ((aligned (32)));
342 
343    n[0] = FLOOR(f[0]);
344    n[1] = FLOOR(f[1]);
345 
346    r[0] = f[0] - n[0];
347    r[1] = f[1] - n[1];
348 
349    w[0] = CUBIC(r[0]);
350    w[1] = CUBIC(r[1]);
351 
352    /*
353     * Much faster in 2d.
354     */
355    v[0] = lattice2(pdata,n[0],   r[0],   n[1],   r[1]);
356    v[1] = lattice2(pdata,n[0]+1, r[0]-1, n[1],   r[1]);
357    v[2] = lattice2(pdata,n[0],   r[0],   n[1]+1, r[1]-1);
358    v[3] = lattice2(pdata,n[0]+1, r[0]-1, n[1]+1, r[1]-1);
359    value = LERP(LERP(v[0], v[1], w[0]),
360          LERP(v[2], v[3], w[0]),
361          w[1]);
362 
363    return CLAMP(-0.99999f, 0.99999f, value);
364 }
365 
366 
367 /**
368  * @brief Gets some 1D Perlin noise from the data.
369  *
370  * Somewhat optimized for speed, probably can't get optimized much more.
371  *
372  *    @param pdata Perlin data to use.
373  *    @param f Position of the noise to get.
374  */
noise_get1(perlin_data_t * pdata,float f[1])375 float noise_get1( perlin_data_t* pdata, float f[1] )
376 {
377    int n[1] __attribute__ ((aligned (32))); /* Indexes to pass to lattice function */
378    float r[1] __attribute__ ((aligned (32))); /* Remainders to pass to lattice function */
379    float w[1] __attribute__ ((aligned (32))); /* Cubic values to pass to interpolation function */
380    float value __attribute__ ((aligned (32)));
381    float v[2] __attribute__ ((aligned (32)));
382 
383    n[0] = FLOOR(f[0]);
384 
385    r[0] = f[0] - n[0];
386 
387    w[0] = CUBIC(r[0]);
388 
389    v[0] = lattice1(pdata,n[0],   r[0]   );
390    v[1] = lattice1(pdata,n[0]+1, r[0]-1 );
391    value = LERP( v[0], v[1], w[0] );
392 
393    return CLAMP(-0.99999f, 0.99999f, value);
394 }
395 
396 
397 /**
398  * @brief Gets 3d Turbulence noise for a position.
399  *
400  *    @param pdata Perlin data to generate noise from.
401  *    @param f Position of the noise.
402  *    @param octaves Octaves to use.
403  *    @return The noise level at the position.
404  */
noise_turbulence3(perlin_data_t * pdata,float f[3],int octaves)405 float noise_turbulence3( perlin_data_t* pdata, float f[3], int octaves )
406 {
407    float tf[3];
408    /* Initialize locals */
409    float value = 0;
410    int i;
411 
412    tf[0] = f[0];
413    tf[1] = f[1];
414    tf[2] = f[2];
415 
416    /* Inner loop of spectral construction, where the fractal is built */
417    for(i=0; i<octaves; i++)
418    {
419       value += ABS(noise_get3(pdata,tf)) * pdata->exponent[i];
420       tf[0] *= pdata->lacunarity;
421       tf[1] *= pdata->lacunarity;
422       tf[2] *= pdata->lacunarity;
423    }
424 
425    return CLAMP(-0.99999f, 0.99999f, value);
426 }
427 
428 
429 /**
430  * @brief Gets 2d Turbulence noise for a position.
431  *
432  *    @param pdata Perlin data to generate noise from.
433  *    @param f Position of the noise.
434  *    @param octaves Octaves to use.
435  *    @return The noise level at the position.
436  */
noise_turbulence2(perlin_data_t * pdata,float f[2],int octaves)437 float noise_turbulence2( perlin_data_t* pdata, float f[2], int octaves )
438 {
439    float tf[2];
440    /* Initialize locals */
441    float value = 0;
442    int i;
443 
444    tf[0] = f[0];
445    tf[1] = f[1];
446 
447    /* Inner loop of spectral construction, where the fractal is built */
448    for(i=0; i<octaves; i++)
449    {
450       value += ABS(noise_get2(pdata,tf)) * pdata->exponent[i];
451       tf[0] *= pdata->lacunarity;
452       tf[1] *= pdata->lacunarity;
453    }
454 
455    return CLAMP(-0.99999f, 0.99999f, value);
456 }
457 
458 
459 /**
460  * @brief Gets 1d Turbulence noise for a position.
461  *
462  *    @param pdata Perlin data to generate noise from.
463  *    @param f Position of the noise.
464  *    @param octaves Octaves to use.
465  *    @return The noise level at the position.
466  */
noise_turbulence1(perlin_data_t * pdata,float f[1],int octaves)467 float noise_turbulence1( perlin_data_t* pdata, float f[1], int octaves )
468 {
469    float tf[1];
470    /* Initialize locals */
471    float value = 0;
472    int i;
473 
474    tf[0] = f[0];
475 
476    /* Inner loop of spectral construction, where the fractal is built */
477    for(i=0; i<octaves; i++)
478    {
479       value += ABS(noise_get1(pdata,tf)) * pdata->exponent[i];
480       tf[0] *= pdata->lacunarity;
481    }
482 
483    return CLAMP(-0.99999f, 0.99999f, value);
484 }
485 
486 
487 
488 #define NOISE_SIMPLEX_GRADIENT_1D(n,h,x) { float grad; h &= 0xF; grad=1.0f+(h & 7); if ( h & 8 ) grad = -grad; n = grad * x; }
489 
490 
491 /**
492  * @brief Gets 1D simplex noise for a position.
493  *
494  *    @param pdata Perlin data to generate noise from.
495  *    @param f Position of the noise.
496  */
noise_simplex1(perlin_data_t * pdata,float f[1])497 float noise_simplex1( perlin_data_t* pdata, float f[1] )
498 {
499    int i0   = (int)FLOOR( f[0]*SIMPLEX_SCALE );
500    int i1   = i0+1;
501    float x0 = f[0]*SIMPLEX_SCALE - i0;
502    float x1 = x0 - 1.0f;
503    float t0 = 1.0f - x0*x0;
504    float t1 = 1.0f - x1*x1;
505    float n0,n1;
506 
507    t0    = t0*t0;
508    t1    = t1*t1;
509    i0    = pdata->map[i0&0xFF];
510    NOISE_SIMPLEX_GRADIENT_1D( n0, i0, x0 );
511    n0   *= t0*t0;
512    i1    = pdata->map[i1&0xFF];
513    NOISE_SIMPLEX_GRADIENT_1D( n1, i1, x1 );
514    n1   *= t1*t1;
515 
516    return 0.25f * (n0+n1);
517 
518 }
519 
520 
521 /**
522  * @brief Frees some noise data.
523  *
524  *    @param noise Noise data to free.
525  */
noise_delete(perlin_data_t * pdata)526 void noise_delete( perlin_data_t* pdata )
527 {
528    free(pdata);
529 }
530 
531 
532 /**
533  * @brief Generates radar interference.
534  *
535  *    @param w Width to generate.
536  *    @param h Height to generate.
537  *    @param rug Rugosity of the interference.
538  *    @return The map generated.
539  */
noise_genRadarInt(const int w,const int h,float rug)540 float* noise_genRadarInt( const int w, const int h, float rug )
541 {
542    int x, y;
543    float f[2];
544    float hurst;
545    float lacunarity;
546    perlin_data_t* noise;
547    float *map;
548    float value;
549 
550    /* pretty default values */
551    hurst       = NOISE_DEFAULT_HURST;
552    lacunarity  = NOISE_DEFAULT_LACUNARITY;
553 
554    /* create noise and data */
555    noise       = noise_new( 2, hurst, lacunarity );
556    map         = malloc(sizeof(float)*w*h);
557    if (map == NULL) {
558       noise_delete( noise );
559       WARN("Out of memory!");
560       return NULL;
561    }
562 
563    /* Start to create the nebula */
564    for (y=0; y<h; y++) {
565 
566       f[1] = rug * (float)y / (float)h;
567       for (x=0; x<w; x++) {
568 
569          f[0] = rug * (float)x / (float)w;
570 
571          /* Get the 2d noise. */
572          value = noise_get2( noise, f );
573 
574          /* Set the value to [0,1]. */
575          map[y*w + x] = (value + 1.) / 2.;
576       }
577    }
578 
579    /* Clean up */
580    noise_delete( noise );
581 
582    /* Results */
583    return map;
584 }
585 
586 
587 /**
588  * @brief Thread worker for generating nebula stuff.
589  *
590  *    @param data Data to pass.
591  */
noise_genNebulaMap_thread(void * data)592 static int noise_genNebulaMap_thread( void *data )
593 {
594    thread_args *args = (thread_args*) data;
595    float f[3];
596    float value;
597    int y, x;
598    float max;
599 
600    /* Generate the layer. */
601    max = 0;
602    f[2] = args->zoom * (float)args->z / (float)args->n;
603 
604    for (y=0; y<args->h; y++) {
605       f[1] = args->zoom * (float)y / (float)args->h;
606 
607       for (x=0; x<args->w; x++) {
608          f[0] = args->zoom * (float)x / (float)args->w;
609 
610          value = noise_turbulence3( args->noise, f, args->octaves );
611          if (max < value)
612             max = value;
613 
614          args->nebula[args->z * args->w * args->h + y * args->w + x] = value;
615        }
616    }
617 
618    /* Set up output. */
619    *args->max = max;
620 
621    /* Clean up. */
622    free( args );
623    return 0;
624 }
625 
626 
627 /**
628  * @brief Generates a 3d nebula map.
629  *
630  *    @param w Width of the map.
631  *    @param h Height of the map.
632  *    @param n Number of slices of the map (2d planes).
633  *    @param rug Rugosity of the map.
634  *    @return The map generated.
635  */
noise_genNebulaMap(const int w,const int h,const int n,float rug)636 float* noise_genNebulaMap( const int w, const int h, const int n, float rug )
637 {
638    int x, y, z, i;
639    int octaves;
640    float hurst;
641    float lacunarity;
642    perlin_data_t* noise;
643    float *nebula;
644    float value;
645    float zoom;
646    float *_max;
647    float max;
648    unsigned int s;
649    thread_args *args;
650    ThreadQueue *vpool;
651 
652    /* pretty default values */
653    octaves     = 3;
654    hurst       = NOISE_DEFAULT_HURST;
655    lacunarity  = NOISE_DEFAULT_LACUNARITY;
656    zoom        = rug * ((float)h/768.)*((float)w/1024.);
657 
658    /* create noise and data */
659    noise      = noise_new( 3, hurst, lacunarity );
660    nebula     = malloc(sizeof(float)*w*h*n);
661    if (nebula == NULL) {
662       noise_delete( noise );
663       WARN("Out of memory!");
664       return NULL;
665    }
666 
667    /* Some debug information and time setting */
668    s = SDL_GetTicks();
669    DEBUG("Generating Nebula of size %dx%dx%d", w, h, n);
670 
671    /* Prepare for generation. */
672    _max        = malloc( sizeof(float) * n );
673 
674    /* Initialize vpool */
675    vpool = vpool_create();
676 
677    /* Start to create the nebula */
678    for (z=0; z<n; z++) {
679       /* Make ze arguments! */
680       args     = malloc( sizeof(thread_args) );
681       args->z  = z;
682       args->zoom = zoom;
683       args->n  = n;
684       args->h  = h;
685       args->w  = w;
686       args->noise = noise;
687       args->octaves = octaves;
688       args->max = &_max[z];
689       args->nebula = nebula;
690 
691       /* Launch ze thread. */
692       vpool_enqueue( vpool, noise_genNebulaMap_thread, args );
693    }
694 
695    /* Wait for threads to signal completion. */
696    vpool_wait( vpool );
697    max = 0.;
698    for (i=0; i<n; i++) {
699       if (_max[i]>max)
700          max = _max[i];
701    }
702 
703    /* Post filtering */
704    value = 1. - max;
705    for (z=0; z<n; z++)
706       for (y=0; y<h; y++)
707          for (x=0; x<w; x++)
708             nebula[z*w*h + y*w + x] += value;
709 
710    /* Clean up */
711    noise_delete( noise );
712    free(_max);
713 
714    /* Results */
715    DEBUG("Nebula Generated in %d ms", SDL_GetTicks() - s );
716    return nebula;
717 }
718 
719 
720 /**
721  * @brief Generates tiny nebula puffs
722  *
723  *    @param w Width of the puff to generate.
724  *    @param h Height of the puff to generate.
725  *    @param rug Rugosity of the puff.
726  *    @return The puff generated.
727  */
noise_genNebulaPuffMap(const int w,const int h,float rug)728 float* noise_genNebulaPuffMap( const int w, const int h, float rug )
729 {
730    int x,y, hw,hh;
731    float d;
732    float f[2];
733    int octaves;
734    float hurst;
735    float lacunarity;
736    perlin_data_t* noise;
737    float *nebula;
738    float value;
739    float zoom;
740    float max;
741 
742    /* pretty default values */
743    octaves     = 3;
744    hurst       = NOISE_DEFAULT_HURST;
745    lacunarity  = NOISE_DEFAULT_LACUNARITY;
746    zoom        = rug;
747 
748    /* create noise and data */
749    noise       = noise_new( 2, hurst, lacunarity );
750    nebula      = malloc(sizeof(float)*w*h);
751    if (nebula == NULL) {
752       noise_delete( noise );
753       WARN("Out of memory!");
754       return NULL;
755    }
756 
757    /* Start to create the nebula */
758    max   = 0.;
759    hw    = w/2;
760    hh    = h/2;
761    d     = (float)MIN(hw,hh);
762    for (y=0; y<h; y++) {
763 
764       f[1] = zoom * (float)y / (float)h;
765       for (x=0; x<w; x++) {
766 
767          f[0] = zoom * (float)x / (float)w;
768 
769          /* Get the 2d noise. */
770          value = noise_turbulence2( noise, f, octaves );
771 
772          /* Make value also depend on distance from center */
773          value *= (d - 1. - sqrtf( (float)((x-hw)*(x-hw) + (y-hh)*(y-hh)) )) / d;
774          if (value < 0.)
775             value = 0.;
776 
777          /* Cap at maximum. */
778          if (max < value)
779             max = value;
780 
781          /* Set the value. */
782          nebula[y*w + x] = value;
783       }
784    }
785 
786    /* Clean up */
787    noise_delete( noise );
788 
789    /* Results */
790    return nebula;
791 }
792 
793 
794