1 /* $Id: tgp.cpp 26056 2013-11-22 21:50:43Z rubidium $ */
4 * This file is part of OpenTTD.
5 * OpenTTD is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, version 2.
6 * OpenTTD is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
7 * See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with OpenTTD. If not, see <http://www.gnu.org/licenses/>.
10 /** @file tgp.cpp OTTD Perlin Noise Landscape Generator, aka TerraGenesis Perlin */
14 #include "clear_map.h"
17 #include "core/random_func.hpp"
18 #include "landscape_type.h"
20 #include "safeguards.h"
24 * Quickie guide to Perlin Noise
25 * Perlin noise is a predictable pseudo random number sequence. By generating
26 * it in 2 dimensions, it becomes a useful random map that, for a given seed
27 * and starting X & Y, is entirely predictable. On the face of it, that may not
28 * be useful. However, it means that if you want to replay a map in a different
29 * terrain, or just vary the sea level, you just re-run the generator with the
30 * same seed. The seed is an int32, and is randomised on each run of New Game.
31 * The Scenario Generator does not randomise the value, so that you can
32 * experiment with one terrain until you are happy, or click "Random" for a new
35 * Perlin Noise is a series of "octaves" of random noise added together. By
36 * reducing the amplitude of the noise with each octave, the first octave of
37 * noise defines the main terrain sweep, the next the ripples on that, and the
38 * next the ripples on that. I use 6 octaves, with the amplitude controlled by
39 * a power ratio, usually known as a persistence or p value. This I vary by the
40 * smoothness selection, as can be seen in the table below. The closer to 1,
41 * the more of that octave is added. Each octave is however raised to the power
42 * of its position in the list, so the last entry in the "smooth" row, 0.35, is
43 * raised to the power of 6, so can only add 0.001838... of the amplitude to
46 * In other words; the first p value sets the general shape of the terrain, the
47 * second sets the major variations to that, ... until finally the smallest
50 * Usefully, this routine is totally scaleable; so when 32bpp comes along, the
51 * terrain can be as bumpy as you like! It is also infinitely expandable; a
52 * single random seed terrain continues in X & Y as far as you care to
53 * calculate. In theory, we could use just one seed value, but randomly select
54 * where in the Perlin XY space we use for the terrain. Personally I prefer
55 * using a simple (0, 0) to (X, Y), with a varying seed.
58 * Other things i have had to do: mountainous wasn't mountainous enough, and
59 * since we only have 0..15 heights available, I add a second generated map
60 * (with a modified seed), onto the original. This generally raises the
61 * terrain, which then needs scaling back down. Overall effect is a general
64 * However, the values on the top of mountains are then almost guaranteed to go
65 * too high, so large flat plateaus appeared at height 15. To counter this, I
66 * scale all heights above 12 to proportion up to 15. It still makes the
67 * mountains have flattish tops, rather than craggy peaks, but at least they
68 * aren't smooth as glass.
71 * For a full discussion of Perlin Noise, please visit:
72 * http://freespace.virgin.net/hugo.elias/models/m_perlin.htm
77 * The algorithm as described in the above link suggests to compute each tile height
78 * as composition of several noise waves. Some of them are computed directly by
79 * noise(x, y) function, some are calculated using linear approximation. Our
80 * first implementation of perlin_noise_2D() used 4 noise(x, y) calls plus
81 * 3 linear interpolations. It was called 6 times for each tile. This was a bit
84 * The following implementation uses optimized algorithm that should produce
85 * the same quality result with much less computations, but more memory accesses.
86 * The overall speedup should be 300% to 800% depending on CPU and memory speed.
88 * I will try to explain it on the example below:
90 * Have a map of 4 x 4 tiles, our simplified noise generator produces only two
91 * values -1 and +1, use 3 octaves with wave length 1, 2 and 4, with amplitudes
92 * 3, 2, 1. Original algorithm produces:
94 * h00 = lerp(lerp(-3, 3, 0/4), lerp(3, -3, 0/4), 0/4) + lerp(lerp(-2, 2, 0/2), lerp( 2, -2, 0/2), 0/2) + -1 = lerp(-3.0, 3.0, 0/4) + lerp(-2, 2, 0/2) + -1 = -3.0 + -2 + -1 = -6.0
95 * h01 = lerp(lerp(-3, 3, 1/4), lerp(3, -3, 1/4), 0/4) + lerp(lerp(-2, 2, 1/2), lerp( 2, -2, 1/2), 0/2) + 1 = lerp(-1.5, 1.5, 0/4) + lerp( 0, 0, 0/2) + 1 = -1.5 + 0 + 1 = -0.5
96 * h02 = lerp(lerp(-3, 3, 2/4), lerp(3, -3, 2/4), 0/4) + lerp(lerp( 2, -2, 0/2), lerp(-2, 2, 0/2), 0/2) + -1 = lerp( 0, 0, 0/4) + lerp( 2, -2, 0/2) + -1 = 0 + 2 + -1 = 1.0
97 * h03 = lerp(lerp(-3, 3, 3/4), lerp(3, -3, 3/4), 0/4) + lerp(lerp( 2, -2, 1/2), lerp(-2, 2, 1/2), 0/2) + 1 = lerp( 1.5, -1.5, 0/4) + lerp( 0, 0, 0/2) + 1 = 1.5 + 0 + 1 = 2.5
99 * h10 = lerp(lerp(-3, 3, 0/4), lerp(3, -3, 0/4), 1/4) + lerp(lerp(-2, 2, 0/2), lerp( 2, -2, 0/2), 1/2) + 1 = lerp(-3.0, 3.0, 1/4) + lerp(-2, 2, 1/2) + 1 = -1.5 + 0 + 1 = -0.5
100 * h11 = lerp(lerp(-3, 3, 1/4), lerp(3, -3, 1/4), 1/4) + lerp(lerp(-2, 2, 1/2), lerp( 2, -2, 1/2), 1/2) + -1 = lerp(-1.5, 1.5, 1/4) + lerp( 0, 0, 1/2) + -1 = -0.75 + 0 + -1 = -1.75
101 * h12 = lerp(lerp(-3, 3, 2/4), lerp(3, -3, 2/4), 1/4) + lerp(lerp( 2, -2, 0/2), lerp(-2, 2, 0/2), 1/2) + 1 = lerp( 0, 0, 1/4) + lerp( 2, -2, 1/2) + 1 = 0 + 0 + 1 = 1.0
102 * h13 = lerp(lerp(-3, 3, 3/4), lerp(3, -3, 3/4), 1/4) + lerp(lerp( 2, -2, 1/2), lerp(-2, 2, 1/2), 1/2) + -1 = lerp( 1.5, -1.5, 1/4) + lerp( 0, 0, 1/2) + -1 = 0.75 + 0 + -1 = -0.25
107 * 1) we need to allocate a bit more tiles: (size_x + 1) * (size_y + 1) = (5 * 5):
109 * 2) setup corner values using amplitude 3
116 * 3a) interpolate values in the middle
117 * { -3.0 X 0.0 X 3.0 }
119 * { 0.0 X 0.0 X 0.0 }
121 * { 3.0 X 0.0 X -3.0 }
123 * 3b) add patches with amplitude 2 to them
124 * { -5.0 X 2.0 X 1.0 }
126 * { 2.0 X -2.0 X 2.0 }
128 * { 1.0 X 2.0 X -5.0 }
130 * 4a) interpolate values in the middle
131 * { -5.0 -1.5 2.0 1.5 1.0 }
132 * { -1.5 -0.75 0.0 0.75 1.5 }
133 * { 2.0 0.0 -2.0 0.0 2.0 }
134 * { 1.5 0.75 0.0 -0.75 -1.5 }
135 * { 1.0 1.5 2.0 -1.5 -5.0 }
137 * 4b) add patches with amplitude 1 to them
138 * { -6.0 -0.5 1.0 2.5 0.0 }
139 * { -0.5 -1.75 1.0 -0.25 2.5 }
140 * { 1.0 1.0 -3.0 1.0 1.0 }
141 * { 2.5 -0.25 1.0 -1.75 -0.5 }
142 * { 0.0 2.5 1.0 -0.5 -6.0 }
148 * As you can see above, each noise function was called just once. Therefore
149 * we don't need to use noise function that calculates the noise from x, y and
150 * some prime. The same quality result we can obtain using standard Random()
155 /** Fixed point type for heights */
156 typedef int16 height_t
;
157 static const int height_decimal_bits
= 4;
159 /** Fixed point array for amplitudes (and percent values) */
160 typedef int amplitude_t
;
161 static const int amplitude_decimal_bits
= 10;
163 /** Height map - allocated array of heights (MapSizeX() + 1) x (MapSizeY() + 1) */
166 height_t
*h
; //< array of heights
167 /* Even though the sizes are always positive, there are many cases where
168 * X and Y need to be signed integers due to subtractions. */
169 int dim_x
; //< height map size_x MapSizeX() + 1
170 int total_size
; //< height map total size
171 int size_x
; //< MapSizeX()
172 int size_y
; //< MapSizeY()
175 * Height map accessor
176 * @param x X position
177 * @param y Y position
178 * @return height as fixed point number
180 inline height_t
&height(uint x
, uint y
)
182 return h
[x
+ y
* dim_x
];
186 /** Global height map instance */
187 static HeightMap _height_map
= {nullptr, 0, 0, 0, 0};
189 /** Conversion: int to height_t */
190 #define I2H(i) ((i) << height_decimal_bits)
191 /** Conversion: height_t to int */
192 #define H2I(i) ((i) >> height_decimal_bits)
194 /** Conversion: int to amplitude_t */
195 #define I2A(i) ((i) << amplitude_decimal_bits)
196 /** Conversion: amplitude_t to int */
197 #define A2I(i) ((i) >> amplitude_decimal_bits)
199 /** Conversion: amplitude_t to height_t */
200 #define A2H(a) ((a) >> (amplitude_decimal_bits - height_decimal_bits))
203 /** Walk through all items of _height_map.h */
204 #define FOR_ALL_TILES_IN_HEIGHT(h) for (h = _height_map.h; h < &_height_map.h[_height_map.total_size]; h++)
206 /** Maximum number of TGP noise frequencies. */
207 static const int MAX_TGP_FREQUENCIES
= 10;
209 /** Desired water percentage (100% == 1024) - indexed by _settings_game.difficulty.quantity_sea_lakes */
210 static const amplitude_t _water_percent
[4] = {70, 170, 270, 420};
213 * Gets the maximum allowed height while generating a map based on
214 * mapsize, terraintype, and the maximum height level.
215 * @return The maximum height for the map generation.
216 * @note Values should never be lower than 3 since the minimum snowline height is 2.
218 static height_t
TGPGetMaxHeight()
221 * Desired maximum height - indexed by:
222 * - _settings_game.difficulty.terrain_type
223 * - min(MapLogX(), MapLogY()) - MIN_MAP_SIZE_BITS
225 * It is indexed by map size as well as terrain type since the map size limits the height of
226 * a usable mountain. For example, on a 64x64 map a 24 high single peak mountain (as if you
227 * raised land 24 times in the center of the map) will leave only a ring of about 10 tiles
228 * around the mountain to build on. On a 4096x4096 map, it won't cover any major part of the map.
230 static const int max_height
[5][MAX_MAP_SIZE_BITS
- MIN_MAP_SIZE_BITS
+ 1] = {
231 /* 64 128 256 512 1024 2048 4096 */
232 { 3, 3, 3, 3, 4, 5, 7 }, ///< Very flat
233 { 5, 7, 8, 9, 14, 19, 31 }, ///< Flat
234 { 8, 9, 10, 15, 23, 37, 61 }, ///< Hilly
235 { 10, 11, 17, 19, 49, 63, 73 }, ///< Mountainous
236 { 12, 19, 25, 31, 67, 75, 87 }, ///< Alpinist
239 int max_height_from_table
= max_height
[_settings_game
.difficulty
.terrain_type
][min(MapLogX(), MapLogY()) - MIN_MAP_SIZE_BITS
];
240 return I2H(min(max_height_from_table
, _settings_game
.construction
.max_heightlevel
));
244 * Get the amplitude associated with the currently selected
245 * smoothness and maximum height level.
246 * @param frequency The frequency to get the amplitudes for
247 * @return The amplitudes to apply to the map.
249 static amplitude_t
GetAmplitude(int frequency
)
251 /* Base noise amplitudes (multiplied by 1024) and indexed by "smoothness setting" and log2(frequency). */
252 static const amplitude_t amplitudes
[][7] = {
253 /* lowest frequency ...... highest (every corner) */
254 {16000, 5600, 1968, 688, 240, 16, 16}, ///< Very smooth
255 {24000, 12800, 6400, 2700, 1024, 128, 16}, ///< Smooth
256 {32000, 19200, 12800, 8000, 3200, 256, 64}, ///< Rough
257 {48000, 24000, 19200, 16000, 8000, 512, 320}, ///< Very rough
260 * Extrapolation factors for ranges before the table.
261 * The extrapolation is needed to account for the higher map heights. They need larger
262 * areas with a particular gradient so that we are able to create maps without too
263 * many steep slopes up to the wanted height level. It's definitely not perfect since
264 * it will bring larger rectangles with similar slopes which makes the rectangular
265 * behaviour of TGP more noticable. However, these height differentiations cannot
266 * happen over much smaller areas; we basically double the "range" to give a similar
267 * slope for every doubling of map height.
269 static const double extrapolation_factors
[] = { 3.3, 2.8, 2.3, 1.8 };
271 int smoothness
= _settings_game
.game_creation
.tgen_smoothness
;
273 /* Get the table index, and return that value if possible. */
274 int index
= frequency
- MAX_TGP_FREQUENCIES
+ lengthof(amplitudes
[smoothness
]);
275 amplitude_t amplitude
= amplitudes
[smoothness
][max(0, index
)];
276 if (index
>= 0) return amplitude
;
278 /* We need to extrapolate the amplitude. */
279 double extrapolation_factor
= extrapolation_factors
[smoothness
];
280 int height_range
= I2H(16);
282 amplitude
= (amplitude_t
)(extrapolation_factor
* (double)amplitude
);
287 return Clamp((TGPGetMaxHeight() - height_range
) / height_range
, 0, 1) * amplitude
;
291 * Check if a X/Y set are within the map.
292 * @param x coordinate x
293 * @param y coordinate y
294 * @return true if within the map
296 static inline bool IsValidXY(int x
, int y
)
298 return x
>= 0 && x
< _height_map
.size_x
&& y
>= 0 && y
< _height_map
.size_y
;
303 * Allocate array of (MapSizeX()+1)*(MapSizeY()+1) heights and init the _height_map structure members
304 * @return true on success
306 static inline bool AllocHeightMap()
310 _height_map
.size_x
= MapSizeX();
311 _height_map
.size_y
= MapSizeY();
313 /* Allocate memory block for height map row pointers */
314 _height_map
.total_size
= (_height_map
.size_x
+ 1) * (_height_map
.size_y
+ 1);
315 _height_map
.dim_x
= _height_map
.size_x
+ 1;
316 _height_map
.h
= CallocT
<height_t
>(_height_map
.total_size
);
318 /* Iterate through height map and initialise values. */
319 FOR_ALL_TILES_IN_HEIGHT(h
) *h
= 0;
324 /** Free height map */
325 static inline void FreeHeightMap()
328 _height_map
.h
= nullptr;
332 * Generates new random height in given amplitude (generated numbers will range from - amplitude to + amplitude)
333 * @param rMax Limit of result
334 * @return generated height
336 static inline height_t
RandomHeight(amplitude_t rMax
)
338 /* Spread height into range -rMax..+rMax */
339 return A2H(RandomRange(2 * rMax
+ 1) - rMax
);
343 * Base Perlin noise generator - fills height map with raw Perlin noise.
345 * This runs several iterations with increasing precision; the last iteration looks at areas
346 * of 1 by 1 tiles, the second to last at 2 by 2 tiles and the initial 2**MAX_TGP_FREQUENCIES
347 * by 2**MAX_TGP_FREQUENCIES tiles.
349 static void HeightMapGenerate()
351 /* Trying to apply noise to uninitialized height map */
352 assert(_height_map
.h
!= nullptr);
354 int start
= max(MAX_TGP_FREQUENCIES
- (int)min(MapLogX(), MapLogY()), 0);
357 for (int frequency
= start
; frequency
< MAX_TGP_FREQUENCIES
; frequency
++) {
358 const amplitude_t amplitude
= GetAmplitude(frequency
);
360 /* Ignore zero amplitudes; it means our map isn't height enough for this
361 * amplitude, so ignore it and continue with the next set of amplitude. */
362 if (amplitude
== 0) continue;
364 const int step
= 1 << (MAX_TGP_FREQUENCIES
- frequency
- 1);
367 /* This is first round, we need to establish base heights with step = size_min */
368 for (int y
= 0; y
<= _height_map
.size_y
; y
+= step
) {
369 for (int x
= 0; x
<= _height_map
.size_x
; x
+= step
) {
370 height_t height
= (amplitude
> 0) ? RandomHeight(amplitude
) : 0;
371 _height_map
.height(x
, y
) = height
;
378 /* It is regular iteration round.
379 * Interpolate height values at odd x, even y tiles */
380 for (int y
= 0; y
<= _height_map
.size_y
; y
+= 2 * step
) {
381 for (int x
= 0; x
<= _height_map
.size_x
- 2 * step
; x
+= 2 * step
) {
382 height_t h00
= _height_map
.height(x
+ 0 * step
, y
);
383 height_t h02
= _height_map
.height(x
+ 2 * step
, y
);
384 height_t h01
= (h00
+ h02
) / 2;
385 _height_map
.height(x
+ 1 * step
, y
) = h01
;
389 /* Interpolate height values at odd y tiles */
390 for (int y
= 0; y
<= _height_map
.size_y
- 2 * step
; y
+= 2 * step
) {
391 for (int x
= 0; x
<= _height_map
.size_x
; x
+= step
) {
392 height_t h00
= _height_map
.height(x
, y
+ 0 * step
);
393 height_t h20
= _height_map
.height(x
, y
+ 2 * step
);
394 height_t h10
= (h00
+ h20
) / 2;
395 _height_map
.height(x
, y
+ 1 * step
) = h10
;
399 /* Add noise for next higher frequency (smaller steps) */
400 for (int y
= 0; y
<= _height_map
.size_y
; y
+= step
) {
401 for (int x
= 0; x
<= _height_map
.size_x
; x
+= step
) {
402 _height_map
.height(x
, y
) += RandomHeight(amplitude
);
408 /** Returns min, max and average height from height map */
409 static void HeightMapGetMinMaxAvg(height_t
*min_ptr
, height_t
*max_ptr
, height_t
*avg_ptr
)
411 height_t h_min
, h_max
, h_avg
, *h
;
413 h_min
= h_max
= _height_map
.height(0, 0);
415 /* Get h_min, h_max and accumulate heights into h_accu */
416 FOR_ALL_TILES_IN_HEIGHT(h
) {
417 if (*h
< h_min
) h_min
= *h
;
418 if (*h
> h_max
) h_max
= *h
;
422 /* Get average height */
423 h_avg
= (height_t
)(h_accu
/ (_height_map
.size_x
* _height_map
.size_y
));
425 /* Return required results */
426 if (min_ptr
!= nullptr) *min_ptr
= h_min
;
427 if (max_ptr
!= nullptr) *max_ptr
= h_max
;
428 if (avg_ptr
!= nullptr) *avg_ptr
= h_avg
;
431 /** Dill histogram and return pointer to its base point - to the count of zero heights */
432 static int *HeightMapMakeHistogram(height_t h_min
, height_t h_max
, int *hist_buf
)
434 int *hist
= hist_buf
- h_min
;
437 /* Count the heights and fill the histogram */
438 FOR_ALL_TILES_IN_HEIGHT(h
) {
446 /** Applies sine wave redistribution onto height map */
447 static void HeightMapSineTransform(height_t h_min
, height_t h_max
)
451 FOR_ALL_TILES_IN_HEIGHT(h
) {
454 if (*h
< h_min
) continue;
456 /* Transform height into 0..1 space */
457 fheight
= (double)(*h
- h_min
) / (double)(h_max
- h_min
);
458 /* Apply sine transform depending on landscape type */
459 switch (_settings_game
.game_creation
.landscape
) {
462 /* Move and scale 0..1 into -1..+1 */
463 fheight
= 2 * fheight
- 1;
465 fheight
= sin(fheight
* M_PI_2
);
466 /* Transform it back from -1..1 into 0..1 space */
467 fheight
= 0.5 * (fheight
+ 1);
472 /* Arctic terrain needs special height distribution.
473 * Redistribute heights to have more tiles at highest (75%..100%) range */
474 double sine_upper_limit
= 0.75;
475 double linear_compression
= 2;
476 if (fheight
>= sine_upper_limit
) {
477 /* Over the limit we do linear compression up */
478 fheight
= 1.0 - (1.0 - fheight
) / linear_compression
;
480 double m
= 1.0 - (1.0 - sine_upper_limit
) / linear_compression
;
481 /* Get 0..sine_upper_limit into -1..1 */
482 fheight
= 2.0 * fheight
/ sine_upper_limit
- 1.0;
483 /* Sine wave transform */
484 fheight
= sin(fheight
* M_PI_2
);
485 /* Get -1..1 back to 0..(1 - (1 - sine_upper_limit) / linear_compression) == 0.0..m */
486 fheight
= 0.5 * (fheight
+ 1.0) * m
;
493 /* Desert terrain needs special height distribution.
494 * Half of tiles should be at lowest (0..25%) heights */
495 double sine_lower_limit
= 0.5;
496 double linear_compression
= 2;
497 if (fheight
<= sine_lower_limit
) {
498 /* Under the limit we do linear compression down */
499 fheight
= fheight
/ linear_compression
;
501 double m
= sine_lower_limit
/ linear_compression
;
502 /* Get sine_lower_limit..1 into -1..1 */
503 fheight
= 2.0 * ((fheight
- sine_lower_limit
) / (1.0 - sine_lower_limit
)) - 1.0;
504 /* Sine wave transform */
505 fheight
= sin(fheight
* M_PI_2
);
506 /* Get -1..1 back to (sine_lower_limit / linear_compression)..1.0 */
507 fheight
= 0.5 * ((1.0 - m
) * fheight
+ (1.0 + m
));
516 /* Transform it back into h_min..h_max space */
517 *h
= (height_t
)(fheight
* (h_max
- h_min
) + h_min
);
518 if (*h
< 0) *h
= I2H(0);
519 if (*h
>= h_max
) *h
= h_max
- 1;
524 * Additional map variety is provided by applying different curve maps
525 * to different parts of the map. A randomized low resolution grid contains
526 * which curve map to use on each part of the make. This filtered non-linearly
527 * to smooth out transitions between curves, so each tile could have between
528 * 100% of one map applied or 25% of four maps.
530 * The curve maps define different land styles, i.e. lakes, low-lands, hills
531 * and mountain ranges, although these are dependent on the landscape style
534 * The level parameter dictates the resolution of the grid. A low resolution
535 * grid will result in larger continuous areas of a land style, a higher
536 * resolution grid splits the style into smaller areas.
537 * @param level Rough indication of the size of the grid sections to style. Small level means large grid sections.
539 static void HeightMapCurves(uint level
)
541 height_t mh
= TGPGetMaxHeight() - I2H(1); // height levels above sea level only
543 /** Basically scale height X to height Y. Everything in between is interpolated. */
544 struct control_point_t
{
545 height_t x
; ///< The height to scale from.
546 height_t y
; ///< The height to scale to.
548 /* Scaled curve maps; value is in height_ts. */
549 #define F(fraction) ((height_t)(fraction * mh))
550 const control_point_t curve_map_1
[] = { { F(0.0), F(0.0) }, { F(0.8), F(0.13) }, { F(1.0), F(0.4) } };
551 const control_point_t curve_map_2
[] = { { F(0.0), F(0.0) }, { F(0.53), F(0.13) }, { F(0.8), F(0.27) }, { F(1.0), F(0.6) } };
552 const control_point_t curve_map_3
[] = { { F(0.0), F(0.0) }, { F(0.53), F(0.27) }, { F(0.8), F(0.57) }, { F(1.0), F(0.8) } };
553 const control_point_t curve_map_4
[] = { { F(0.0), F(0.0) }, { F(0.4), F(0.3) }, { F(0.7), F(0.8) }, { F(0.92), F(0.99) }, { F(1.0), F(0.99) } };
556 /** Helper structure to index the different curve maps. */
557 struct control_point_list_t
{
558 size_t length
; ///< The length of the curve map.
559 const control_point_t
*list
; ///< The actual curve map.
561 const control_point_list_t curve_maps
[] = {
562 { lengthof(curve_map_1
), curve_map_1
},
563 { lengthof(curve_map_2
), curve_map_2
},
564 { lengthof(curve_map_3
), curve_map_3
},
565 { lengthof(curve_map_4
), curve_map_4
},
568 height_t ht
[lengthof(curve_maps
)];
569 MemSetT(ht
, 0, lengthof(ht
));
571 /* Set up a grid to choose curve maps based on location; attempt to get a somewhat square grid */
572 float factor
= sqrt((float)_height_map
.size_x
/ (float)_height_map
.size_y
);
573 uint sx
= Clamp((int)(((1 << level
) * factor
) + 0.5), 1, 128);
574 uint sy
= Clamp((int)(((1 << level
) / factor
) + 0.5), 1, 128);
575 byte
*c
= AllocaM(byte
, sx
* sy
);
577 for (uint i
= 0; i
< sx
* sy
; i
++) {
578 c
[i
] = Random() % lengthof(curve_maps
);
582 for (int x
= 0; x
< _height_map
.size_x
; x
++) {
584 /* Get our X grid positions and bi-linear ratio */
585 float fx
= (float)(sx
* x
) / _height_map
.size_x
+ 1.0f
;
588 float xr
= 2.0f
* (fx
- x1
) - 1.0f
;
589 xr
= sin(xr
* M_PI_2
);
590 xr
= sin(xr
* M_PI_2
);
591 xr
= 0.5f
* (xr
+ 1.0f
);
592 float xri
= 1.0f
- xr
;
599 for (int y
= 0; y
< _height_map
.size_y
; y
++) {
601 /* Get our Y grid position and bi-linear ratio */
602 float fy
= (float)(sy
* y
) / _height_map
.size_y
+ 1.0f
;
605 float yr
= 2.0f
* (fy
- y1
) - 1.0f
;
606 yr
= sin(yr
* M_PI_2
);
607 yr
= sin(yr
* M_PI_2
);
608 yr
= 0.5f
* (yr
+ 1.0f
);
609 float yri
= 1.0f
- yr
;
616 uint corner_a
= c
[x1
+ sx
* y1
];
617 uint corner_b
= c
[x1
+ sx
* y2
];
618 uint corner_c
= c
[x2
+ sx
* y1
];
619 uint corner_d
= c
[x2
+ sx
* y2
];
621 /* Bitmask of which curve maps are chosen, so that we do not bother
622 * calculating a curve which won't be used. */
623 uint corner_bits
= 0;
624 corner_bits
|= 1 << corner_a
;
625 corner_bits
|= 1 << corner_b
;
626 corner_bits
|= 1 << corner_c
;
627 corner_bits
|= 1 << corner_d
;
629 height_t
*h
= &_height_map
.height(x
, y
);
631 /* Do not touch sea level */
632 if (*h
< I2H(1)) continue;
634 /* Only scale above sea level */
637 /* Apply all curve maps that are used on this tile. */
638 for (uint t
= 0; t
< lengthof(curve_maps
); t
++) {
639 if (!HasBit(corner_bits
, t
)) continue;
642 const control_point_t
*cm
= curve_maps
[t
].list
;
643 for (uint i
= 0; i
< curve_maps
[t
].length
- 1; i
++) {
644 const control_point_t
&p1
= cm
[i
];
645 const control_point_t
&p2
= cm
[i
+ 1];
647 if (*h
>= p1
.x
&& *h
< p2
.x
) {
648 ht
[t
] = p1
.y
+ (*h
- p1
.x
) * (p2
.y
- p1
.y
) / (p2
.x
- p1
.x
);
656 /* Apply interpolation of curve map results. */
657 *h
= (height_t
)((ht
[corner_a
] * yri
+ ht
[corner_b
] * yr
) * xri
+ (ht
[corner_c
] * yri
+ ht
[corner_d
] * yr
) * xr
);
659 /* Readd sea level */
665 /** Adjusts heights in height map to contain required amount of water tiles */
666 static void HeightMapAdjustWaterLevel(amplitude_t water_percent
, height_t h_max_new
)
668 height_t h_min
, h_max
, h_avg
, h_water_level
;
669 int64 water_tiles
, desired_water_tiles
;
673 HeightMapGetMinMaxAvg(&h_min
, &h_max
, &h_avg
);
675 /* Allocate histogram buffer and clear its cells */
676 int *hist_buf
= CallocT
<int>(h_max
- h_min
+ 1);
678 hist
= HeightMapMakeHistogram(h_min
, h_max
, hist_buf
);
680 /* How many water tiles do we want? */
681 desired_water_tiles
= A2I(((int64
)water_percent
) * (int64
)(_height_map
.size_x
* _height_map
.size_y
));
683 /* Raise water_level and accumulate values from histogram until we reach required number of water tiles */
684 for (h_water_level
= h_min
, water_tiles
= 0; h_water_level
< h_max
; h_water_level
++) {
685 water_tiles
+= hist
[h_water_level
];
686 if (water_tiles
>= desired_water_tiles
) break;
689 /* We now have the proper water level value.
690 * Transform the height map into new (normalized) height map:
691 * values from range: h_min..h_water_level will become negative so it will be clamped to 0
692 * values from range: h_water_level..h_max are transformed into 0..h_max_new
693 * where h_max_new is depending on terrain type and map size.
695 FOR_ALL_TILES_IN_HEIGHT(h
) {
696 /* Transform height from range h_water_level..h_max into 0..h_max_new range */
697 *h
= (height_t
)(((int)h_max_new
) * (*h
- h_water_level
) / (h_max
- h_water_level
)) + I2H(1);
698 /* Make sure all values are in the proper range (0..h_max_new) */
699 if (*h
< 0) *h
= I2H(0);
700 if (*h
>= h_max_new
) *h
= h_max_new
- 1;
706 static double perlin_coast_noise_2D(const double x
, const double y
, const double p
, const int prime
);
709 * This routine sculpts in from the edge a random amount, again a Perlin
710 * sequence, to avoid the rigid flat-edge slopes that were present before. The
711 * Perlin noise map doesn't know where we are going to slice across, and so we
712 * often cut straight through high terrain. The smoothing routine makes it
713 * legal, gradually increasing up from the edge to the original terrain height.
714 * By cutting parts of this away, it gives a far more irregular edge to the
715 * map-edge. Sometimes it works beautifully with the existing sea & lakes, and
716 * creates a very realistic coastline. Other times the variation is less, and
717 * the map-edge shows its cliff-like roots.
719 * This routine may be extended to randomly sculpt the height of the terrain
720 * near the edge. This will have the coast edge at low level (1-3), rising in
721 * smoothed steps inland to about 15 tiles in. This should make it look as
722 * though the map has been built for the map size, rather than a slice through
725 * Please note that all the small numbers; 53, 101, 167, etc. are small primes
726 * to help give the perlin noise a bit more of a random feel.
728 static void HeightMapCoastLines(uint8 water_borders
)
730 int smallest_size
= min(_settings_game
.game_creation
.map_x
, _settings_game
.game_creation
.map_y
);
731 const int margin
= 4;
736 /* Lower to sea level */
737 for (y
= 0; y
<= _height_map
.size_y
; y
++) {
738 if (HasBit(water_borders
, BORDER_NE
)) {
740 max_x
= abs((perlin_coast_noise_2D(_height_map
.size_y
- y
, y
, 0.9, 53) + 0.25) * 5 + (perlin_coast_noise_2D(y
, y
, 0.35, 179) + 1) * 12);
741 max_x
= max((smallest_size
* smallest_size
/ 64) + max_x
, (smallest_size
* smallest_size
/ 64) + margin
- max_x
);
742 if (smallest_size
< 8 && max_x
> 5) max_x
/= 1.5;
743 for (x
= 0; x
< max_x
; x
++) {
744 _height_map
.height(x
, y
) = 0;
748 if (HasBit(water_borders
, BORDER_SW
)) {
750 max_x
= abs((perlin_coast_noise_2D(_height_map
.size_y
- y
, y
, 0.85, 101) + 0.3) * 6 + (perlin_coast_noise_2D(y
, y
, 0.45, 67) + 0.75) * 8);
751 max_x
= max((smallest_size
* smallest_size
/ 64) + max_x
, (smallest_size
* smallest_size
/ 64) + margin
- max_x
);
752 if (smallest_size
< 8 && max_x
> 5) max_x
/= 1.5;
753 for (x
= _height_map
.size_x
; x
> (_height_map
.size_x
- 1 - max_x
); x
--) {
754 _height_map
.height(x
, y
) = 0;
759 /* Lower to sea level */
760 for (x
= 0; x
<= _height_map
.size_x
; x
++) {
761 if (HasBit(water_borders
, BORDER_NW
)) {
763 max_y
= abs((perlin_coast_noise_2D(x
, _height_map
.size_y
/ 2, 0.9, 167) + 0.4) * 5 + (perlin_coast_noise_2D(x
, _height_map
.size_y
/ 3, 0.4, 211) + 0.7) * 9);
764 max_y
= max((smallest_size
* smallest_size
/ 64) + max_y
, (smallest_size
* smallest_size
/ 64) + margin
- max_y
);
765 if (smallest_size
< 8 && max_y
> 5) max_y
/= 1.5;
766 for (y
= 0; y
< max_y
; y
++) {
767 _height_map
.height(x
, y
) = 0;
771 if (HasBit(water_borders
, BORDER_SE
)) {
773 max_y
= abs((perlin_coast_noise_2D(x
, _height_map
.size_y
/ 3, 0.85, 71) + 0.25) * 6 + (perlin_coast_noise_2D(x
, _height_map
.size_y
/ 3, 0.35, 193) + 0.75) * 12);
774 max_y
= max((smallest_size
* smallest_size
/ 64) + max_y
, (smallest_size
* smallest_size
/ 64) + margin
- max_y
);
775 if (smallest_size
< 8 && max_y
> 5) max_y
/= 1.5;
776 for (y
= _height_map
.size_y
; y
> (_height_map
.size_y
- 1 - max_y
); y
--) {
777 _height_map
.height(x
, y
) = 0;
783 /** Start at given point, move in given direction, find and Smooth coast in that direction */
784 static void HeightMapSmoothCoastInDirection(int org_x
, int org_y
, int dir_x
, int dir_y
)
786 const int max_coast_dist_from_edge
= 35;
787 const int max_coast_Smooth_depth
= 35;
790 int ed
; // coast distance from edge
793 height_t h_prev
= I2H(1);
796 assert(IsValidXY(org_x
, org_y
));
798 /* Search for the coast (first non-water tile) */
799 for (x
= org_x
, y
= org_y
, ed
= 0; IsValidXY(x
, y
) && ed
< max_coast_dist_from_edge
; x
+= dir_x
, y
+= dir_y
, ed
++) {
801 if (_height_map
.height(x
, y
) >= I2H(1)) break;
803 /* Coast found in the neighborhood? */
804 if (IsValidXY(x
+ dir_y
, y
+ dir_x
) && _height_map
.height(x
+ dir_y
, y
+ dir_x
) > 0) break;
806 /* Coast found in the neighborhood on the other side */
807 if (IsValidXY(x
- dir_y
, y
- dir_x
) && _height_map
.height(x
- dir_y
, y
- dir_x
) > 0) break;
810 /* Coast found or max_coast_dist_from_edge has been reached.
811 * Soften the coast slope */
812 for (depth
= 0; IsValidXY(x
, y
) && depth
<= max_coast_Smooth_depth
; depth
++, x
+= dir_x
, y
+= dir_y
) {
813 h
= _height_map
.height(x
, y
);
814 h
= min(h
, h_prev
+ (4 + depth
)); // coast softening formula
815 _height_map
.height(x
, y
) = h
;
820 /** Smooth coasts by modulating height of tiles close to map edges with cosine of distance from edge */
821 static void HeightMapSmoothCoasts(uint8 water_borders
)
824 /* First Smooth NW and SE coasts (y close to 0 and y close to size_y) */
825 for (x
= 0; x
< _height_map
.size_x
; x
++) {
826 if (HasBit(water_borders
, BORDER_NW
)) HeightMapSmoothCoastInDirection(x
, 0, 0, 1);
827 if (HasBit(water_borders
, BORDER_SE
)) HeightMapSmoothCoastInDirection(x
, _height_map
.size_y
- 1, 0, -1);
829 /* First Smooth NE and SW coasts (x close to 0 and x close to size_x) */
830 for (y
= 0; y
< _height_map
.size_y
; y
++) {
831 if (HasBit(water_borders
, BORDER_NE
)) HeightMapSmoothCoastInDirection(0, y
, 1, 0);
832 if (HasBit(water_borders
, BORDER_SW
)) HeightMapSmoothCoastInDirection(_height_map
.size_x
- 1, y
, -1, 0);
837 * This routine provides the essential cleanup necessary before OTTD can
838 * display the terrain. When generated, the terrain heights can jump more than
839 * one level between tiles. This routine smooths out those differences so that
840 * the most it can change is one level. When OTTD can support cliffs, this
841 * routine may not be necessary.
843 static void HeightMapSmoothSlopes(height_t dh_max
)
845 for (int y
= 0; y
<= (int)_height_map
.size_y
; y
++) {
846 for (int x
= 0; x
<= (int)_height_map
.size_x
; x
++) {
847 height_t h_max
= min(_height_map
.height(x
> 0 ? x
- 1 : x
, y
), _height_map
.height(x
, y
> 0 ? y
- 1 : y
)) + dh_max
;
848 if (_height_map
.height(x
, y
) > h_max
) _height_map
.height(x
, y
) = h_max
;
851 for (int y
= _height_map
.size_y
; y
>= 0; y
--) {
852 for (int x
= _height_map
.size_x
; x
>= 0; x
--) {
853 height_t h_max
= min(_height_map
.height(x
< _height_map
.size_x
? x
+ 1 : x
, y
), _height_map
.height(x
, y
< _height_map
.size_y
? y
+ 1 : y
)) + dh_max
;
854 if (_height_map
.height(x
, y
) > h_max
) _height_map
.height(x
, y
) = h_max
;
860 * Height map terraform post processing:
861 * - water level adjusting
864 * - height histogram redistribution by sine wave transform
866 static void HeightMapNormalize()
868 int sea_level_setting
= _settings_game
.difficulty
.quantity_sea_lakes
;
869 const amplitude_t water_percent
= sea_level_setting
!= (int)CUSTOM_SEA_LEVEL_NUMBER_DIFFICULTY
? _water_percent
[sea_level_setting
] : _settings_game
.game_creation
.custom_sea_level
* 1024 / 100;
870 const height_t h_max_new
= TGPGetMaxHeight();
871 const height_t roughness
= 7 + 3 * _settings_game
.game_creation
.tgen_smoothness
;
873 HeightMapAdjustWaterLevel(water_percent
, h_max_new
);
875 byte water_borders
= _settings_game
.construction
.freeform_edges
? _settings_game
.game_creation
.water_borders
: 0xF;
876 if (water_borders
== BORDERS_RANDOM
) water_borders
= GB(Random(), 0, 4);
878 HeightMapCoastLines(water_borders
);
879 HeightMapSmoothSlopes(roughness
);
881 HeightMapSmoothCoasts(water_borders
);
882 HeightMapSmoothSlopes(roughness
);
884 HeightMapSineTransform(I2H(1), h_max_new
);
886 if (_settings_game
.game_creation
.variety
> 0) {
887 HeightMapCurves(_settings_game
.game_creation
.variety
);
890 HeightMapSmoothSlopes(I2H(1));
894 * The Perlin Noise calculation using large primes
895 * The initial number is adjusted by two values; the generation_seed, and the
896 * passed parameter; prime.
897 * prime is used to allow the perlin noise generator to create useful random
898 * numbers from slightly different series.
900 static double int_noise(const long x
, const long y
, const int prime
)
902 long n
= x
+ y
* prime
+ _settings_game
.game_creation
.generation_seed
;
906 /* Pseudo-random number generator, using several large primes */
907 return 1.0 - (double)((n
* (n
* n
* 15731 + 789221) + 1376312589) & 0x7fffffff) / 1073741824.0;
912 * This routine determines the interpolated value between a and b
914 static inline double linear_interpolate(const double a
, const double b
, const double x
)
916 return a
+ x
* (b
- a
);
921 * This routine returns the smoothed interpolated noise for an x and y, using
922 * the values from the surrounding positions.
924 static double interpolated_noise(const double x
, const double y
, const int prime
)
926 const int integer_X
= (int)x
;
927 const int integer_Y
= (int)y
;
929 const double fractional_X
= x
- (double)integer_X
;
930 const double fractional_Y
= y
- (double)integer_Y
;
932 const double v1
= int_noise(integer_X
, integer_Y
, prime
);
933 const double v2
= int_noise(integer_X
+ 1, integer_Y
, prime
);
934 const double v3
= int_noise(integer_X
, integer_Y
+ 1, prime
);
935 const double v4
= int_noise(integer_X
+ 1, integer_Y
+ 1, prime
);
937 const double i1
= linear_interpolate(v1
, v2
, fractional_X
);
938 const double i2
= linear_interpolate(v3
, v4
, fractional_X
);
940 return linear_interpolate(i1
, i2
, fractional_Y
);
945 * This is a similar function to the main perlin noise calculation, but uses
946 * the value p passed as a parameter rather than selected from the predefined
947 * sequences. as you can guess by its title, i use this to create the indented
948 * coastline, which is just another perlin sequence.
950 static double perlin_coast_noise_2D(const double x
, const double y
, const double p
, const int prime
)
954 for (int i
= 0; i
< 6; i
++) {
955 const double frequency
= (double)(1 << i
);
956 const double amplitude
= pow(p
, (double)i
);
958 total
+= interpolated_noise((x
* frequency
) / 64.0, (y
* frequency
) / 64.0, prime
) * amplitude
;
965 /** A small helper function to initialize the terrain */
966 static void TgenSetTileHeight(TileIndex tile
, int height
)
968 SetTileHeight(tile
, height
);
970 /* Only clear the tiles within the map area. */
971 if (IsInnerTile(tile
)) {
972 MakeClear(tile
, CLEAR_GRASS
, 3);
977 * The main new land generator using Perlin noise. Desert landscape is handled
978 * different to all others to give a desert valley between two high mountains.
979 * Clearly if a low height terrain (flat/very flat) is chosen, then the tropic
980 * areas wont be high enough, and there will be very little tropic on the map.
981 * Thus Tropic works best on Hilly or Mountainous.
983 void GenerateTerrainPerlin()
985 if (!AllocHeightMap()) return;
986 GenerateWorldSetAbortCallback(FreeHeightMap
);
990 IncreaseGeneratingWorldProgress(GWP_LANDSCAPE
);
992 HeightMapNormalize();
994 IncreaseGeneratingWorldProgress(GWP_LANDSCAPE
);
996 /* First make sure the tiles at the north border are void tiles if needed. */
997 if (_settings_game
.construction
.freeform_edges
) {
998 for (int y
= 0; y
< _height_map
.size_y
- 1; y
++) MakeVoid(_height_map
.size_x
* y
);
999 for (int x
= 0; x
< _height_map
.size_x
; x
++) MakeVoid(x
);
1002 int max_height
= H2I(TGPGetMaxHeight());
1004 /* Transfer height map into OTTD map */
1005 for (int y
= 0; y
< _height_map
.size_y
; y
++) {
1006 for (int x
= 0; x
< _height_map
.size_x
; x
++) {
1007 TgenSetTileHeight(TileXY(x
, y
), Clamp(H2I(_height_map
.height(x
, y
)), 0, max_height
));
1011 IncreaseGeneratingWorldProgress(GWP_LANDSCAPE
);
1014 GenerateWorldSetAbortCallback(nullptr);