noise2d.cpp

Engine/source/util/noise2d.cpp

More...

Public Functions

clamp(F32 f, F32 m)

Detailed Description

Public Functions

clamp(F32 f, F32 m)

  1
  2//-----------------------------------------------------------------------------
  3// Copyright (c) 2012 GarageGames, LLC
  4//
  5// Permission is hereby granted, free of charge, to any person obtaining a copy
  6// of this software and associated documentation files (the "Software"), to
  7// deal in the Software without restriction, including without limitation the
  8// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
  9// sell copies of the Software, and to permit persons to whom the Software is
 10// furnished to do so, subject to the following conditions:
 11//
 12// The above copyright notice and this permission notice shall be included in
 13// all copies or substantial portions of the Software.
 14//
 15// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 16// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 17// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 18// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 19// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 20// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
 21// IN THE SOFTWARE.
 22//-----------------------------------------------------------------------------
 23
 24#include "util/noise2d.h"
 25#include "core/util/tVector.h"
 26
 27//--------------------------------------
 28Noise2D::Noise2D()
 29{
 30   dMemset(mPermutation, 0, sizeof(mPermutation));
 31   dMemset(mGradient, 0, sizeof(mGradient));
 32   mSeed   = 0;
 33}
 34
 35Noise2D::~Noise2D()
 36{
 37}
 38
 39
 40//--------------------------------------
 41void Noise2D::normalize(F32 v[2])
 42{
 43   F32 s;
 44
 45   s = mSqrt(v[0] * v[0] + v[1] * v[1]);
 46   v[0] = v[0] / s;
 47   v[1] = v[1] / s;
 48}
 49
 50
 51//--------------------------------------
 52void Noise2D::setSeed(U32 seed)
 53{
 54   if (mSeed == seed)
 55      return;
 56   mSeed = seed;
 57   mRandom.setSeed(mSeed);
 58
 59   S32 i, j, k;
 60
 61   for (i = 0 ; i < SIZE ; i++) {
 62      mPermutation[i] = i;
 63
 64      for (j = 0 ; j < 2 ; j++)
 65         mGradient[i][j] = mRandom.randF( -1.0f, 1.0f );
 66      normalize(mGradient[i]);
 67   }
 68
 69   while (--i) {
 70      k = mPermutation[i];
 71      j = mRandom.randI(0, SIZE-1);
 72      mPermutation[i] = mPermutation[j];
 73      mPermutation[j] = k;
 74   }
 75
 76   // extend the size of the arrays x2 to get rid of a bunch of MODs
 77   // we'd have to do later in the code
 78   for (i = 0 ; i < SIZE + 2 ; i++) {
 79      mPermutation[SIZE + i] = mPermutation[i];
 80      for (j = 0 ; j < 2 ; j++)
 81         mGradient[SIZE + i][j] = mGradient[i][j];
 82   }
 83}
 84
 85
 86//--------------------------------------
 87U32 Noise2D::getSeed()
 88{
 89   return mSeed;
 90}
 91
 92
 93inline F32 Noise2D::lerp(F32 t, F32 a, F32 b)
 94{
 95   return a + t * (b - a);
 96}
 97
 98
 99inline F32 Noise2D::curve(F32 t)
100{
101   return t * t * (3.0f - 2.0f * t);
102}
103
104
105inline F32 clamp(F32 f, F32 m)
106{
107   while (f > m)
108      f -= m;
109   while (f < 0.0f)
110      f += m;
111   return f;
112}
113
114
115//--------------------------------------
116void Noise2D::fBm( Vector<F32> *dst, U32 size, U32 interval, F32 h, F32 octaves )
117{
118   interval = getMin(U32(128), getMax(U32(1), interval));
119   F32 H = getMin(1.0f, getMax(0.0f, h));
120   octaves = getMin(5.0f, getMax(1.0f, octaves));
121   F32 lacunarity = 2.0f;
122
123   F32 exponent_array[32];
124
125   U32 shift = getBinLog2( size );
126
127   // precompute and store spectral weights
128   // seize required memory for exponent_array
129   F32 frequency = 1.0;
130   for (U32 i=0; i<=octaves; i++)
131   {
132      // compute weight for each frequency
133      exponent_array[i] = mPow( frequency, -H );
134      frequency *= lacunarity;
135   }
136
137   // initialize dst
138   for (S32 k=0; k < (size*size); k++)
139      (*dst)[k] = 0.0f;
140
141   F32 scale = 1.0f / (F32)size * interval;
142   for (S32 o=0; o<octaves; o++)
143   {
144      F32 exp = exponent_array[o];
145      for (S32 y=0; y<size; y++)
146      {
147         F32 fy = (F32)y * scale;
148         for (S32 x=0; x<size; x++)
149         {
150            F32 fx = (F32)x * scale;
151            F32 noise = getValue(fx, fy, interval);
152            (*dst)[x + (y << shift)] += noise * exp;
153         }
154      }
155      scale    *= lacunarity;
156      interval  = (U32)(interval * lacunarity);
157   }
158}
159
160
161//--------------------------------------
162void Noise2D::rigidMultiFractal(Vector<F32> *dst, Vector<F32> *sig, U32 size, U32 interval, F32 h, F32 octaves)
163{
164   interval = getMin(U32(128), getMax(U32(1), interval));
165   F32 H = getMin(1.0f, getMax(0.0f, h));
166   octaves = getMin(5.0f, getMax(1.0f, octaves));
167   F32 lacunarity = 2.0f;
168   F32 offset     = 1.0f;
169   F32 gain       = 2.0f;
170
171   U32 shift = getBinLog2( size );
172
173   F32 exponent_array[32];
174
175   // precompute and store spectral weights
176   // seize required memory for exponent_array
177   F32 frequency = 1.0;
178   for (U32 i=0; i<=octaves; i++)
179   {
180      // compute weight for each frequency
181      exponent_array[i] = mPow( frequency, -H );
182      frequency *= lacunarity;
183   }
184
185   F32 scale = 1.0f / (F32)size * interval;
186
187   //--------------------------------------
188   // compute first octave
189   for (S32 y=0; y<size; y++)
190   {
191      F32 fy = (F32)y * scale;
192      for (S32 x=0; x<size; x++)
193      {
194         F32 fx = (F32)x * scale;
195
196         F32 signal = mFabs(getValue(fx,fy,interval));   // get absolute value of signal (this creates the ridges)
197         //signal = mSqrt(signal);
198         signal = offset - signal;  // invert and translate (note that "offset" should be ~= 1.0)
199         signal *= signal + 0.1;          // square the signal, to increase "sharpness" of ridges
200
201         // assign initial values
202         (*dst)[x + (y << shift)] = signal;
203         (*sig)[x + (y << shift)] = signal;
204      }
205   }
206
207   //--------------------------------------
208   // compute remaining octaves
209   for (S32 o=1; o<octaves; o++)
210   {
211      // increase the frequency
212      scale    *= lacunarity;
213      interval  = (U32)(interval * lacunarity);
214      F32 exp   = exponent_array[o];
215
216      for (S32 y=0; y<size; y++)
217      {
218         F32 fy = (F32)y * scale;
219         for (S32 x=0; x<size; x++)
220         {
221            F32 fx = (F32)x * scale;
222            U32 index  = x + (y << shift);
223            F32 result = (*dst)[index];
224            F32 signal = (*sig)[index];
225
226            // weight successive contributions by previous signal
227            F32 weight = mClampF(signal * gain, 0.0f, 1.0f);
228
229            signal = mFabs(getValue( fx, fy, interval ));
230
231            signal = offset - signal;
232            signal *= signal + 0.2;
233            // weight the contribution
234            signal *= weight;
235            result += signal * exp;
236
237            (*dst)[index] = result;
238            (*sig)[index] = signal;
239         }
240      }
241   }
242   
243   for (S32 k=0; k < (size*size); k++)
244      (*dst)[k] = ((*dst)[k]-1.0f)/2.0f;
245}
246
247bool Noise2D::erodeHydraulic( Vector<F32> *src, Vector<F32> *dst, U32 iterations, U32 size )
248{
249   // early out if there is nothing to do
250   if (iterations == 0 )
251   {
252      *dst = *src;
253      return true;
254   }
255
256   F32 fmin, fmax;
257   getMinMax( src, &fmin, &fmax, size);
258
259   U32 shift = getBinLog2( size );
260   U32 mask = size - 1;
261
262
263// currently using SCRATCH_3 for debugging -- Rick
264   Vector<F32> scratch = *src;
265   U32 *o = (U32*)scratch.address();
266   Vector<F32> a = *src;
267   Vector<F32> b = *src;
268   Vector<F32> c = *src;
269
270   for (S32 k=0; k < (size*size); k++)
271      c[k] = 0.0f;
272
273   for (S32 i=0; i<iterations; i++)
274   {
275      b = a;
276
277      for (S32 y=0; y<size; y++)
278      {
279         for (S32 x=0; x<size; x++)
280         {
281            U32 srcOffset = (x + (y << shift));
282            F32 height    = a[srcOffset];
283            o[srcOffset]  = srcOffset;
284            for (S32 y1=y-1; y1 <= y+1; y1++)
285            {
286               F32 maxDelta = 0.0f;
287               S32 ywrap = (y1 & mask);
288               for (S32 x1=x-1; x1 <= x+1; x1++)
289               {
290                  if (x1 != x && y1 != y)
291                  {
292                     U32 adjOffset  = ((x1 & mask) + (ywrap << shift));
293                     F32 &adjHeight = a[adjOffset];
294                     F32 delta   = height - adjHeight;
295                     if (x1 != x || y1 != y)
296                        delta *= 1.414213562f;    // compensate for diagonals
297                     if (delta > maxDelta)
298                     {
299                        maxDelta = delta;
300                        o[srcOffset] = adjOffset;
301                     }
302                  }
303               }
304            }
305         }
306      }
307      for (S32 j=0; j < (size*size); j++)
308      {
309         F32 &s = a[j];
310         F32 &d = b[ o[j] ];
311         F32 delta = s - d;
312         if (delta > 0.0f)
313         {
314            F32 alt = (s-fmin) / (fmax-fmin);
315            F32 amt = delta * (0.1f * (1.0f-alt));
316            s -= amt;
317            d += amt;
318         }
319      }
320// debug only
321      for (S32 k=0; k < (size*size); k++)
322         c[k] += b[k] - a[k];
323
324      Vector<F32> tmp = a;
325      a = b;
326      b = tmp;
327   }
328   *dst = b;
329   //*dst = *c;
330
331   return true;
332}
333
334
335
336bool Noise2D::erodeThermal(Vector<F32> *src, Vector<F32> *dst, F32 slope, F32 materialLoss, U32 iterations, U32 size, U32 squareSize, F32 maxHeight )
337{
338   // early out if there is nothing to do
339   if (iterations == 0 )
340   {
341      *dst = *src;
342      return true;
343   }
344
345
346   F32 fmin, fmax;
347   getMinMax(src, &fmin, &fmax, size);
348
349   Vector<F32> a = *src;
350   // Heightfield *b = getScratch(1);
351   Vector<F32> r;
352   r.setSize( size * size );
353   //dMemset( r.address(), 0, r.memSize() );
354
355   F32 conservation = 1.0f - mClampF(materialLoss, 0.0f,  100.0f)/100.0f;
356   slope            = mClampF(conservation, 0.0f, 89.0f);                  // clamp to 0-89 degrees
357
358   F32 talusConst = mTan(mDegToRad(slope)) * squareSize; // in world units
359   talusConst = talusConst * (fmax-fmin) / maxHeight;     // scale to current height units
360   F32 p = 0.1f;
361
362   U32 mask = size - 1;
363   U32 shift = getBinLog2( size );
364
365   for (U32 i=0; i<iterations; i++)
366   {
367      // clear out the rubble accumulation field
368      dMemset( r.address(), 0, r.memSize() );
369
370      for (S32 y=0; y<size; y++)
371      {
372         for (S32 x=0; x<size; x++)
373         {
374            F32 *height    = &a[ x + ( y << shift )];
375            F32 *dstHeight = &r[ x + ( y << shift )];
376
377            // for each height look at the immediate surrounding heights
378            // if any are higher than talusConst erode on me
379            for (S32 y1=y-1; y1 <= y+1; y1++)
380            {
381               S32 ywrap = (y1 & mask);
382               for (S32 x1=x-1; x1 <= x+1; x1++)
383               {
384                  if (x1 != x && y1 != y)
385                  {
386                     S32 adjOffset = ((x1 & mask) + (ywrap << shift));
387                     F32 adjHeight = a[adjOffset];
388                     F32 delta     = adjHeight - *height;
389                     if (delta > talusConst)
390                     {
391                        F32 rubble    = p * (delta - talusConst);
392                        r[adjOffset] -= rubble;
393                        *dstHeight   += rubble * conservation;
394                     }
395                  }
396               }
397            }
398         }
399      }
400      for (S32 k=0; k < (size*size); k++)
401         a[k] += r[k];
402   }
403   *dst = a;
404   return true;
405}
406
407void Noise2D::getMinMax( Vector<F32> *src, F32 *fmin, F32 *fmax, U32 size )
408{
409   if (!src)
410      return;
411
412   F32 *p = (*src).address();
413   *fmin = *p;
414   *fmax = *p;
415   for (S32 i=0; i < (size*size); i++, p++)
416   {
417      if (*fmin > *p) *fmin = *p;
418      if (*fmax < *p) *fmax = *p;
419   }
420}
421
422//--------------------------------------
423F32 Noise2D::turbulence(F32 x, F32 y, F32 freq)
424{
425   F32 t, x2, y2;
426
427   for ( t = 0.0f ; freq >= 3.0f ; freq /= 2.0f)
428   {
429      x2 = freq * x;
430      y2 = freq * y;
431      t += mFabs(getValue(x2, y2, (S32)freq)) / freq;
432   }
433   return t;
434}
435
436
437//--------------------------------------
438inline void Noise2D::setup(F32 t, S32 &b0, S32 &b1, F32 &r0, F32 &r1)
439{
440   // find the bounding integers of u
441   b0 = S32(t) & SIZE_MASK;
442   b1 = (b0+1) & SIZE_MASK;
443
444   // seperate the fractional components
445   r0 = t - (S32)t;
446   r1 = r0 - 1.0f;
447}
448
449inline F32 Noise2D::dot(const F32 *q, F32 rx, F32 ry)
450{
451   return (rx * q[0] + ry * q[1] );
452}
453
454
455
456//--------------------------------------
457F32 Noise2D::getValue(F32 x, F32 y, S32 interval)
458{
459   S32 bx0, bx1, by0, by1;
460   F32 rx0, rx1, ry0, ry1;
461
462   // Imagine having a square of the type
463   //  p0---p1    Where p0 = (bx0, by0)   +----> U
464   //  |(u,v)|          p1 = (bx1, by0)   |
465   //  |     |          p2 = (bx0, by1)   |    Coordinate System
466   //  p2---p3          p3 = (bx1, by1)   V
467   // The u, v point in 2D texture space is bounded by this rectangle.
468
469   // Goal, determine the scalar at the points p0, p1, p2, p3.
470   // Then the scalar of the point (u, v) will be found by linear interpolation.
471
472   // First step:  Get the 2D coordinates of the points p0, p1, p2, p3.
473   // We also need vectors pointing from each point in the square above and
474   // ending at the (u,v) coordinate located inside the square.
475   // The vector (rx0, ry0) goes from P0 to the (u,v) coordinate.
476   // The vector (rx1, ry0) goes from P1 to the (u,v) coordinate.
477   // The vector (rx0, ry1) goes from P2 to the (u,v) coordinate.
478   // The vector (rx1, ry1) goes from P3 to the (u,v) coordinate.
479
480   setup(x, bx0, bx1, rx0, rx1);
481   setup(y, by0, by1, ry0, ry1);
482
483   // Make sure the box corners fall within the interval
484   // so that the final output will wrap on itself
485   bx0 = bx0 % interval;
486   bx1 = bx1 % interval;
487   by0 = by0 % interval;
488   by1 = by1 % interval;
489
490   S32 i = mPermutation[ bx0 ];
491   S32 j = mPermutation[ bx1 ];
492
493   S32 b00 = mPermutation[ i + by0 ];
494   S32 b10 = mPermutation[ j + by0 ];
495   S32 b01 = mPermutation[ i + by1 ];
496   S32 b11 = mPermutation[ j + by1 ];
497
498   // Next, calculate the dropoff component about the point p0.
499   F32 sx = curve(rx0);
500   F32 sy = curve(ry0);
501
502   // Now, for each point in the square shown above, calculate the dot
503   // product of the gradiant vector and the vector going from each square
504   // corner point to the (u,v) point inside the square.
505   F32 u = dot(mGradient[ b00 ], rx0,ry0);
506   F32 v = dot(mGradient[ b10 ], rx1,ry0);
507
508   // Interpolation along the X axis.
509   F32 a = lerp(sx, u, v);
510
511   u = dot(mGradient[ b01 ], rx0,ry1);
512   v = dot(mGradient[ b11 ], rx1,ry1);
513
514   // Interpolation along the Y axis.
515   F32 b = lerp(sx, u, v);
516
517   // Final Interpolation
518   return lerp(sy, a, b);
519}
520
521
522