1 /*
2 * Modification History
3 *
4 * 2006-September-26 Jason Rohrer
5 * Switched to cosine interpolation.
6 * Added optimizations discovered with profiler. Reduced running time by 18%.
7 */
8
9
10
11 #include "landscape.h"
12
13 #include <math.h>
14 #include <stdlib.h>
15
16 #include "minorGems/util/SimpleVector.h"
17
18
19
landscape(double inX,double inY,double inT,double inBaseFrequency,double inRoughness,int inDetail)20 double landscape( double inX, double inY, double inT,
21 double inBaseFrequency, double inRoughness,
22 int inDetail ) {
23
24 if( inDetail < 0 ) {
25 return 0.0;
26 }
27 else {
28 // frequency of octave
29 double frequency = inBaseFrequency * pow( 2, inDetail );
30 double amplitude = pow( inRoughness, inDetail );
31 return amplitude * noise4d( inX * frequency,
32 inY * frequency,
33 // index different planes of noise
34 inDetail,
35 inT )
36 + landscape( inX, inY, inT, inBaseFrequency, inRoughness,
37 inDetail - 1 );
38 }
39 }
40
41
42
variableRoughnessLandscape(double inX,double inY,double inT,double inBaseFrequency,double inRoughnessChangeFrequency,double inMinRoughness,double inMaxRoughness,int inDetail)43 double variableRoughnessLandscape( double inX, double inY, double inT,
44 double inBaseFrequency,
45 double inRoughnessChangeFrequency,
46 double inMinRoughness,
47 double inMaxRoughness,
48 int inDetail ) {
49
50 double roughnessFreqX = inX * inRoughnessChangeFrequency;
51 double roughnessFreqY = inY * inRoughnessChangeFrequency;
52
53 // use low-frequency noise 4d to select landscape roughness
54 // between 0 and 1
55 double roughness =
56 ( noise4d( 6, roughnessFreqX, roughnessFreqY, inT ) + 1 ) / 2;
57
58 // move roughness into specified range
59 roughness =
60 roughness * ( inMaxRoughness - inMinRoughness ) +
61 inMinRoughness;
62
63 return landscape( inX, inY, inT, inBaseFrequency, roughness, inDetail );
64 }
65
66
67
getRandomPoints(double inStartX,double inEndX,double inStartY,double inEndY,double inT,double inSampleStepSize,double inDensity,double ** outXCoordinates,double ** outYCoordinates)68 int getRandomPoints( double inStartX, double inEndX,
69 double inStartY, double inEndY,
70 double inT,
71 double inSampleStepSize,
72 double inDensity,
73 double **outXCoordinates,
74 double **outYCoordinates ) {
75
76 SimpleVector<double> *xCoordinates = new SimpleVector<double>();
77 SimpleVector<double> *yCoordinates = new SimpleVector<double>();
78
79 // discretize startX and start Y so that sample grid for differently-placed
80 // windows always meshes
81 // use ceil to ensure that starting points are always inside the
82 // inStart/inEnd bounds
83 double discretizedStartX =
84 inSampleStepSize * ceil( inStartX / inSampleStepSize );
85 double discretizedStartY =
86 inSampleStepSize * ceil( inStartY / inSampleStepSize );
87
88 // put a point wherever we have a zero-crossing
89 double lastSample = 1;
90
91 for( double x=discretizedStartX; x<=inEndX; x+=inSampleStepSize ) {
92 for( double y=discretizedStartY; y<=inEndY; y+=inSampleStepSize ) {
93 double landscapeSample =
94 variableRoughnessLandscape(
95 30 * x + 1000, 30 * y + 1000, inT + 1000,
96 0.01, 0.001, 0.25, 0.65, 0 );
97
98 // shift landscape up to reduce chance of zero-crossing
99 landscapeSample = (1-inDensity) * 0.5 + 0.5 + landscapeSample ;
100
101 if( landscapeSample < 0 &&
102 lastSample >= 0 ||
103 landscapeSample >= 0 &&
104 landscapeSample < 0 ) {
105
106 // sign change
107
108 // hit
109 xCoordinates->push_back( x );
110 yCoordinates->push_back( y );
111 }
112
113 lastSample = landscapeSample;
114 }
115 }
116
117 *outXCoordinates = xCoordinates->getElementArray();
118 *outYCoordinates = yCoordinates->getElementArray();
119
120 int numPoints = xCoordinates->size();
121
122 delete xCoordinates;
123 delete yCoordinates;
124
125 return numPoints;
126 }
127
128
129
130 /**
131 * Computes a 32-bit random number.
132 * Use linear congruential method.
133 *
134 * @param inSeed the seed to use.
135 */
136 // this is the readable version of the funcion
137 // it has been turned into a set of macros below
random32_readable(unsigned long inSeed)138 inline unsigned long random32_readable( unsigned long inSeed ) {
139 // this is the true hot-spot of the entire landscape function
140 // thus, optimization is warranted.
141
142 // multiplier = 3141592621
143 // use hex to avoid warnings
144 //unsigned long multiplier = 0xBB40E62D;
145 //unsigned long increment = 1;
146
147 // better:
148 // unsigned long multiplier = 196314165
149 // unsigned long increment = 907633515
150
151 // this will automatically be mod-ed by 2^32 because of the limit
152 // of the unsigned long type
153 // return multiplier * inSeed + increment;
154 //return 0xBB40E62D * inSeed + 1;
155 //return 196314165 * inSeed + 907633515;
156
157 //int n = ( inSeed << 13 ) ^ inSeed;
158 //return n * (n * n * 15731 + 789221) + 1376312589;
159
160 //const unsigned long Num1 = (inSeed * 0xFEA09B9DLU) + 1;
161 //const unsigned long Num2 = ((inSeed * 0xB89C8895LU) + 1) >> 16;
162 //return Num1 ^ Num2;
163
164 /*
165 unsigned int rseed=(inSeed*15064013)^(inSeed*99991+604322121)^(inSeed*45120321)^(inSeed*5034121+13);
166
167 const unsigned long Num1 = (inSeed * 0xFEA09B9DLU) + 1;
168
169 const unsigned long Num2 = ((inSeed * 0xB89C8895LU) + 1) >> 16;
170
171 rseed *= Num1 ^ Num2;
172
173 return rseed;
174 */
175
176 const unsigned long Num1 = (inSeed * 0xFEA09B9DLU) + 1;
177 const unsigned long Num2 = ((inSeed^Num1) * 0x9C129511LU) + 1;
178 const unsigned long Num3 = (inSeed * 0x2512CFB8LU) + 1;
179 const unsigned long Num4 = ((inSeed^Num3) * 0xB89C8895LU) + 1;
180 const unsigned long Num5 = (inSeed * 0x6BF962C1LU) + 1;
181 const unsigned long Num6 = ((inSeed^Num5) * 0x4BF962C1LU) + 1;
182
183 return Num2 ^ (Num4 >> 11) ^ (Num6 >> 22);
184 }
185
186
187 // faster as a set of macros
188 #define Num1( inSeed ) \
189 ( ( inSeed * 0xFEA09B9DLU ) + 1 )
190
191 #define Num2( inSeed ) \
192 ( ( ( inSeed ^ Num1( inSeed ) ) * 0x9C129511LU ) + 1 )
193
194 #define Num3( inSeed ) \
195 ( ( inSeed * 0x2512CFB8LU ) + 1 )
196
197 #define Num4( inSeed ) \
198 ( ( ( inSeed ^ Num3( inSeed ) ) * 0xB89C8895LU ) + 1 )
199
200 #define Num5( inSeed ) \
201 ( ( inSeed * 0x6BF962C1LU ) + 1 )
202
203 #define Num6( inSeed ) \
204 ( ( ( inSeed ^ Num5( inSeed ) ) * 0x4BF962C1LU ) + 1 )
205
206 #define random32( inSeed ) \
207 ( Num2( inSeed ) ^ (Num4( inSeed ) >> 11) ^ (Num6( inSeed ) >> 22) )
208
209
210
211
212
213
214
215 #define invMaxIntAsDouble 2.32830643708e-10
216 // 1/(x/2) = 2*(1/x)
217 //double invHalfMaxIntAsDouble = 2 * invMaxIntAsDouble;
218 // 2.32830643708e-10
219 //+ 2.32830643708e-10
220 //-------------------
221 // 4.65661287416e-10
222 #define invHalfMaxIntAsDouble 4.65661287416e-10
223
224
225
226 #define mixFour( x, y, z, t ) ( x ^ (y * 57) ^ (z * 131) ^ (t * 2383) )
227
228
229
230 /**
231 * Maps 4d integer coordinates into a [-1..1] noise space.
232 *
233 * @param x, y, z, t the 4d coordinates.
234 *
235 * @return a random value in the range [-1..1]
236 */
237 // keep readable version around for reference
238 // it has been replaced by a macro below
noise4dInt32_readable(unsigned long x,unsigned long y,unsigned long z,unsigned long t)239 inline double noise4dInt32_readable( unsigned long x,
240 unsigned long y,
241 unsigned long z,
242 unsigned long t ) {
243
244
245 //double maxIntAsDouble = 4294967295.0;
246
247 // modular addition automatic
248 // multiply x, y, z, and t by distinct primes to
249 // avoid correllations.
250 // using xor ( ^ ) here seems to avoid correllations that show
251 // up when using addition.
252
253 // mix x, y, z, and t
254 unsigned long randomSeed =
255 x ^
256 y * 57 ^
257 z * 131 ^
258 t * 2383;
259
260 // a random value between 0 and max unsigned long
261 unsigned long randomValue = random32( randomSeed );
262
263 // a random value between 0 and 2
264 double zeroTwoValue = randomValue * invHalfMaxIntAsDouble;
265
266 // a random value between -1 and 1
267 return zeroTwoValue - 1;
268 }
269
270
271 // noise4dInt32 function call itself was the slowest spot in code
272 // (found with profiler)
273 // turn into a set of macros
274
275 // matches original parameter format
276 #define noise4dInt32( x, y, z, t ) \
277 random32( mixFour( x, y, z, t ) ) \
278 * invHalfMaxIntAsDouble - 1
279
280 // problem: now that random32 is a macro, we are passing the unevaluated
281 // expression, ( x ^ (y * 57) ^ (z * 131) ^ (t * 2383) ), down into it.
282 // it is being evaluated 6 times within the depths of the random32 macro
283
284 // thus, we need to provide a new format where the caller can precompute
285 // the mix for us. This is even faster.
286 #define noise1dInt32( precomputedMix ) \
287 random32( precomputedMix ) \
288 * invHalfMaxIntAsDouble - 1
289
290
291
292
293 /*
294 * The following functions (blendNoiseNd) do 4d linear interpolation
295 * one dimension at a time.
296 *
297 * The end result is 8 calls to blendNoise1d (and 16 calls to noise4dInt32).
298 *
299 * This method was inspired by the functional implementations---I am
300 * decomposing a complicated problem into sub-problems that are easier
301 * to solve.
302 */
303
304
305 // faster than f * b + (1-f) * a
306 // one less multiply
307 #define linearInterpolation( t, a, b ) ( a + t * ( b - a ) )
308
309
310 /**
311 * Blends 4d discrete (integer-parameter) noise function along one dimension
312 * with 3 fixed integer parameters.
313 */
blendNoise1d(double x,unsigned long y,unsigned long z,unsigned long t)314 inline double blendNoise1d( double x,
315 unsigned long y,
316 unsigned long z,
317 unsigned long t ) {
318
319 double floorX = floor( x );
320 unsigned long floorIntX = (unsigned long)floorX;
321
322 if( floorX == x ) {
323 unsigned long precomputedMix = mixFour( floorIntX, y, z, t );
324
325 return noise1dInt32( precomputedMix );
326 }
327 else {
328 unsigned long ceilIntX = floorIntX + 1;
329
330 // cosine interpolation
331 // from http://freespace.virgin.net/hugo.elias/models/m_perlin.htm
332 double ft = ( x - floorX ) * M_PI;
333 double f = ( 1 - cos( ft ) ) * .5;
334
335
336 // need to pre-store intermediate values because noise4dInt32 is a
337 // macro
338 // thus, we end up calling the noise1dInt32 function instead
339
340 unsigned long precomputedMix = mixFour( floorIntX, y, z, t );
341 double valueAtFloor = noise1dInt32( precomputedMix );
342
343 precomputedMix = mixFour( ceilIntX, y, z, t );
344 double valueAtCeiling = noise1dInt32( precomputedMix );
345
346 return linearInterpolation( f, valueAtFloor, valueAtCeiling );
347 }
348 }
349
350
351
352 /**
353 * Blends 4d discrete (integer-parameter) noise function along 2 dimensions
354 * with 2 fixed integer parameters.
355 */
blendNoise2d(double x,double y,unsigned long z,unsigned long t)356 double blendNoise2d( double x,
357 double y,
358 unsigned long z,
359 unsigned long t ) {
360
361 double floorY = floor( y );
362 unsigned long floorIntY = (unsigned long)floorY;
363
364 if( floorY == y ) {
365 return blendNoise1d( x, floorIntY, z, t );
366 }
367 else {
368 unsigned long ceilIntY = floorIntY + 1;
369
370 // cosine interpolation
371 // from http://freespace.virgin.net/hugo.elias/models/m_perlin.htm
372 double ft = ( y - floorY ) * M_PI;
373 double f = ( 1 - cos( ft ) ) * .5;
374
375 return ( f ) * blendNoise1d( x, ceilIntY, z, t ) +
376 ( 1 - f ) * blendNoise1d( x, floorIntY, z, t );
377 }
378 }
379
380
381
382 /**
383 * Blends 4d discrete (integer-parameter) noise function along 3 dimensions
384 * with 1 fixed integer parameters.
385 */
blendNoise3d(double x,double y,double z,unsigned long t)386 double blendNoise3d( double x,
387 double y,
388 double z,
389 unsigned long t ) {
390
391 double floorZ = floor( z );
392 unsigned long floorIntZ = (unsigned long)floorZ;
393
394 if( floorZ == z ) {
395 return blendNoise2d( x, y, floorIntZ, t );
396 }
397 else {
398 unsigned long ceilIntZ = floorIntZ + 1;
399
400 // cosine interpolation
401 // from http://freespace.virgin.net/hugo.elias/models/m_perlin.htm
402 double ft = ( z - floorZ ) * M_PI;
403 double f = ( 1 - cos( ft ) ) * .5;
404
405 return ( f ) * blendNoise2d( x, y, ceilIntZ, t ) +
406 ( 1 - f ) * blendNoise2d( x, y, floorIntZ, t );
407 }
408 }
409
410
411
412 /**
413 * Blends 4d discrete (integer-parameter) noise function along 4 dimensions.
414 */
noise4d(double x,double y,double z,double t)415 double noise4d( double x,
416 double y,
417 double z,
418 double t ) {
419
420 double floorT = floor( t );
421 unsigned long floorIntT = (unsigned long)floorT;
422
423 if( floorT == t ) {
424 return blendNoise3d( x, y, z, floorIntT );
425 }
426 else {
427 unsigned long ceilIntT = floorIntT + 1;
428
429 // cosine interpolation
430 // from http://freespace.virgin.net/hugo.elias/models/m_perlin.htm
431 double ft = ( t - floorT ) * M_PI;
432 double f = ( 1 - cos( ft ) ) * .5;
433
434 return ( f ) * blendNoise3d( x, y, z, ceilIntT ) +
435 ( 1 - f ) * blendNoise3d( x, y, z, floorIntT );
436 }
437 }
438
439
440
441
442
443