2 Commits

Author SHA1 Message Date
b1f4437a91 Added settings for margin, and documented in settingtypes.txt 2022-01-03 16:39:02 +01:00
a00c1cbd39 Added margin with a settable width near grid border
Elevation gets closer to -50 when approaching the border
2022-01-03 16:39:02 +01:00
16 changed files with 303 additions and 413 deletions

View File

@ -1,18 +0,0 @@
CHANGELOG
=========
## `v1.0.2` (2022-01-10)
- Use builtin logging system and appropriate loglevels
- Skip empty chunks, when generating high above ground (~20% speedup)
- Minor optimizations (turning global variables to local...)
## `v1.0.1` (2021-09-14)
- Automatically switch to `singlenode` mapgen at init time
## `v1.0` (2021-08-01)
- Rewritten pregen code (terrainlib) in pure Lua
- Optimized grid loading
- Load grid nodes on request by default
- Changed river width settings
- Added map size in settings
- Added logs

View File

@ -1,5 +1,5 @@
# Map Generator with Rivers
`mapgen_rivers v1.0.2` by Gaël de Sailly.
`mapgen_rivers v1.0.1` by Gaël de Sailly.
Semi-procedural map generator for Minetest 5.x. It aims to create realistic and nice-looking landscapes for the game, focused on river networks. It is based on algorithms modelling water flow and river erosion at a broad scale, similar to some used by researchers in Earth Sciences. It is taking some inspiration from [Fastscape](https://github.com/fastscape-lem/fastscape).
@ -13,11 +13,10 @@ It used to be composed of a Python script doing pre-generation, and a Lua mod re
License: GNU LGPLv3.0
Code: Gaël de Sailly
Flow routing algorithm concept (in `terrainlib/rivermapper.lua`): Cordonnier, G., Bovy, B., & Braun, J. (2019). A versatile, linear complexity algorithm for flow routing in topographies with depressions. Earth Surface Dynamics, 7(2), 549-562.
# Requirements
No required dependency, but [`biomegen`](https://gitlab.com/gaelysam/biomegen) recommended (provides biome system).
Mod dependencies: `default` required, and [`biomegen`](https://github.com/Gael-de-Sailly/biomegen) optional (provides biome system).
# Installation
This mod should be placed in the `mods/` directory of Minetest like any other mod.

View File

@ -1,21 +1,3 @@
local function version_is_lower(v1, v2)
local d1, c1, d2, c2
while #v1 > 0 and #v2 > 0 do
d1, c1, v1 = v1:match("^(%d*)(%D*)(.*)$")
d2, c2, v2 = v2:match("^(%d*)(%D*)(.*)$")
d1 = tonumber(d1) or -1
d2 = tonumber(d2) or -1
if d1 ~= d2 then
return d1 < d2
end
if c1 ~= c2 then
return c1 < c2
end
end
return false
end
local function fix_min_catchment(settings, is_global)
local prefix = is_global and "mapgen_rivers_" or ""
@ -39,18 +21,6 @@ local function fix_compatibility_minetest(settings)
if previous_version == "0.0" then
fix_min_catchment(settings, true)
end
if version_is_lower(previous_version, "1.0.2-dev1") then
local blocksize = tonumber(settings:get("mapgen_rivers_blocksize") or 15)
local grid_x_size = tonumber(settings:get("mapgen_rivers_grid_x_size"))
if grid_x_size then
settings:set("mapgen_rivers_map_x_size", tostring(grid_x_size * blocksize))
end
local grid_z_size = tonumber(settings:get("mapgen_rivers_grid_z_size"))
if grid_z_size then
settings:set("mapgen_rivers_map_z_size", tostring(grid_z_size * blocksize))
end
end
end
local function fix_compatibility_mapgen_rivers(settings)
@ -59,18 +29,6 @@ local function fix_compatibility_mapgen_rivers(settings)
if previous_version == "0.0" then
fix_min_catchment(settings, false)
end
if version_is_lower(previous_version, "1.0.2-dev1") then
local blocksize = tonumber(settings:get("blocksize") or 15)
local grid_x_size = tonumber(settings:get("grid_x_size"))
if grid_x_size then
settings:set("map_x_size", tostring(grid_x_size * blocksize))
end
local grid_z_size = tonumber(settings:get("grid_z_size"))
if grid_z_size then
settings:set("map_z_size", tostring(grid_z_size * blocksize))
end
end
end
return fix_compatibility_minetest, fix_compatibility_mapgen_rivers

View File

@ -6,7 +6,7 @@ local transform_quadri = dofile(modpath .. 'geometry.lua')
local sea_level = mapgen_rivers.settings.sea_level
local riverbed_slope = mapgen_rivers.settings.riverbed_slope * mapgen_rivers.settings.blocksize
local MAP_BOTTOM = -31000
local out_elev = mapgen_rivers.settings.margin_elev
-- Localize for performance
local floor, min, max = math.floor, math.min, math.max
@ -128,8 +128,8 @@ local function heightmaps(minp, maxp)
terrain_height_map[i] = terrain_height
lake_height_map[i] = lake_height
else
terrain_height_map[i] = MAP_BOTTOM
lake_height_map[i] = MAP_BOTTOM
terrain_height_map[i] = out_elev
lake_height_map[i] = out_elev
end
i = i + 1
end

View File

@ -16,11 +16,12 @@ local elevation_chill = mapgen_rivers.settings.elevation_chill
local use_distort = mapgen_rivers.settings.distort
local use_biomes = mapgen_rivers.settings.biomes
local use_biomegen_mod = use_biomes and minetest.global_exists('biomegen')
use_biomes = use_biomes and minetest.global_exists('default') and not use_biomegen_mod
use_biomes = use_biomes and not use_biomegen_mod
if use_biomegen_mod then
biomegen.set_elevation_chill(elevation_chill)
end
dofile(modpath .. 'noises.lua')
local heightmaps = dofile(modpath .. 'heightmap.lua')
@ -146,19 +147,15 @@ local function generate(minp, maxp, seed)
end
end
local c_stone = minetest.get_content_id("mapgen_stone")
local c_water = minetest.get_content_id("mapgen_water_source")
local c_rwater = minetest.get_content_id("mapgen_river_water_source")
local c_dirt, c_lawn, c_dirtsnow, c_snow, c_sand, c_ice
if use_biomes then
c_dirt = minetest.get_content_id("default:dirt")
c_lawn = minetest.get_content_id("default:dirt_with_grass")
c_dirtsnow = minetest.get_content_id("default:dirt_with_snow")
c_snow = minetest.get_content_id("default:snowblock")
c_sand = minetest.get_content_id("default:sand")
c_ice = minetest.get_content_id("default:ice")
end
local c_stone = minetest.get_content_id("default:stone")
local c_dirt = minetest.get_content_id("default:dirt")
local c_lawn = minetest.get_content_id("default:dirt_with_grass")
local c_dirtsnow = minetest.get_content_id("default:dirt_with_snow")
local c_snow = minetest.get_content_id("default:snowblock")
local c_sand = minetest.get_content_id("default:sand")
local c_water = minetest.get_content_id("default:water_source")
local c_rwater = minetest.get_content_id("default:river_water_source")
local c_ice = minetest.get_content_id("default:ice")
local vm, emin, emax = minetest.get_mapgen_object("voxelmanip")
vm:get_data(data)

View File

@ -1,3 +1,4 @@
name = mapgen_rivers
title = Map generator with realistic rivers
optional_depends = biomegen, default
depends = default
optional_depends = biomegen

80
noises.lua Normal file
View File

@ -0,0 +1,80 @@
local def_setting = mapgen_rivers.define_setting
mapgen_rivers.noise_params = {
base = def_setting('np_base', 'noise', {
offset = 0,
scale = 300,
seed = 2469,
octaves = 8,
spread = {x=2048, y=2048, z=2048},
persist = 0.6,
lacunarity = 2,
flags = "eased",
}),
distort_x = def_setting('np_distort_x', 'noise', {
offset = 0,
scale = 1,
seed = -4574,
spread = {x=64, y=32, z=64},
octaves = 3,
persistence = 0.75,
lacunarity = 2,
}),
distort_z = def_setting('np_distort_z', 'noise', {
offset = 0,
scale = 1,
seed = -7940,
spread = {x=64, y=32, z=64},
octaves = 3,
persistence = 0.75,
lacunarity = 2,
}),
distort_amplitude = def_setting('np_distort_amplitude', 'noise', {
offset = 0,
scale = 10,
seed = 676,
spread = {x=1024, y=1024, z=1024},
octaves = 5,
persistence = 0.5,
lacunarity = 2,
flags = "absvalue",
}),
heat = minetest.get_mapgen_setting_noiseparams('mg_biome_np_heat'),
heat_blend = minetest.get_mapgen_setting_noiseparams('mg_biome_np_heat_blend'),
}
-- Convert to number because Minetest API is not able to do it cleanly...
for name, np in pairs(mapgen_rivers.noise_params) do
for field, value in pairs(np) do
if field ~= 'flags' and type(value) == 'string' then
np[field] = tonumber(value) or value
elseif field == 'spread' then
for dir, v in pairs(value) do
value[dir] = tonumber(v) or v
end
end
end
end
local heat = mapgen_rivers.noise_params.heat
local base = mapgen_rivers.noise_params.base
local settings = mapgen_rivers.settings
heat.offset = heat.offset + settings.sea_level * settings.elevation_chill
base.spread.x = base.spread.x / settings.blocksize
base.spread.y = base.spread.y / settings.blocksize
base.spread.z = base.spread.z / settings.blocksize
for name, np in pairs(mapgen_rivers.noise_params) do
local lac = np.lacunarity or 2
if lac > 1 then
local omax = math.floor(math.log(math.min(np.spread.x, np.spread.y, np.spread.z)) / math.log(lac))+1
if np.octaves > omax then
minetest.log("warning", "[mapgen_rivers] Noise " .. name .. ": 'octaves' reduced to " .. omax)
np.octaves = omax
end
end
end

View File

@ -11,9 +11,8 @@ dofile(modpath .. 'load.lua')
mapgen_rivers.grid = {}
local blocksize = mapgen_rivers.settings.blocksize
local X = math.floor(mapgen_rivers.settings.map_x_size / blocksize)
local Z = math.floor(mapgen_rivers.settings.map_z_size / blocksize)
local X = mapgen_rivers.settings.grid_x_size
local Z = mapgen_rivers.settings.grid_z_size
local function offset_converter(o)
return (o + 0.5) * (1/256)
@ -81,6 +80,7 @@ local function index(x, z)
return z*X+x+1
end
local blocksize = mapgen_rivers.settings.blocksize
local min_catchment = mapgen_rivers.settings.min_catchment
local max_catchment = mapgen_rivers.settings.max_catchment

View File

@ -5,7 +5,6 @@ local blocksize = mapgen_rivers.settings.blocksize
local tectonic_speed = mapgen_rivers.settings.tectonic_speed
local np_base = table.copy(mapgen_rivers.noise_params.base)
np_base.spread = vector.divide(np_base.spread, blocksize)
local evol_params = mapgen_rivers.settings.evol_params
@ -14,14 +13,42 @@ local time_step = mapgen_rivers.settings.evol_time_step
local niter = math.ceil(time/time_step)
time_step = time / niter
local use_margin = mapgen_rivers.settings.margin
local margin_width = mapgen_rivers.settings.margin_width / blocksize
local margin_elev = mapgen_rivers.settings.margin_elev
local function margin(dem, width, elev)
local X, Y = dem.X, dem.Y
for i=1, width do
local c1 = ((i-1)/width) ^ 0.5
local c2 = (1-c1) * elev
local index = (i-1)*X + 1
for x=1, X do
dem[index] = dem[index] * c1 + c2
index = index + 1
end
index = i
for y=1, Y do
dem[index] = dem[index] * c1 + c2
index = index + X
end
index = X*(Y-i) + 1
for x=1, X do
dem[index] = dem[index] * c1 + c2
index = index + 1
end
index = X-i + 1
for y=1, Y do
dem[index] = dem[index] * c1 + c2
index = index + X
end
end
end
local function pregenerate(keep_loaded)
local grid = mapgen_rivers.grid
local size = grid.size
if size.x * size.y > 4e6 then
minetest.log("warning", "[mapgen_rivers] You are going to generate a very large grid (>4M nodes). If you experience problems, you should increase blocksize or reduce map size.")
end
local seed = tonumber(minetest.get_mapgen_setting("seed"))
np_base.seed = (np_base.seed or 0) + seed
@ -31,6 +58,10 @@ local function pregenerate(keep_loaded)
dem.X = size.x
dem.Y = size.y
if use_margin then
margin(dem, margin_width, margin_elev)
end
local model = EvolutionModel(evol_params)
model.dem = dem
local ref_dem = model:define_isostasy(dem)
@ -46,6 +77,9 @@ local function pregenerate(keep_loaded)
if i < niter then
if tectonic_step ~= 0 then
nobj_base:get_3d_map_flat({x=0, y=tectonic_step*i, z=0}, ref_dem)
if use_margin then
margin(ref_dem, margin_width, margin_elev)
end
end
model:isostasy()
end

View File

@ -1,7 +1,7 @@
local mtsettings = minetest.settings
local mgrsettings = Settings(minetest.get_worldpath() .. '/mapgen_rivers.conf')
mapgen_rivers.version = "1.0.2-dev1"
mapgen_rivers.version = "1.0.1"
local previous_version_mt = mtsettings:get("mapgen_rivers_version") or "0.0"
local previous_version_mgr = mgrsettings:get("version") or "0.0"
@ -19,33 +19,13 @@ end
mtsettings:set("mapgen_rivers_version", mapgen_rivers.version)
mgrsettings:set("version", mapgen_rivers.version)
local defaults
do
local f = io.open(mapgen_rivers.modpath .. "/settings_default.json")
defaults = minetest.parse_json(f:read("*all"))
f:close()
end
-- Convert strings to numbers in noise params because Minetest API is not able to do it cleanly...
local function clean_np(np)
for field, value in pairs(np) do
if field ~= 'flags' and type(value) == 'string' then
np[field] = tonumber(value) or value
elseif field == 'spread' then
for dir, v in pairs(value) do
value[dir] = tonumber(v) or v
end
end
end
end
function mapgen_rivers.define_setting(name, dtype, default)
if dtype == "number" or dtype == "string" then
local v = mgrsettings:get(name)
if v == nil then
v = mtsettings:get('mapgen_rivers_' .. name)
if v == nil then
v = defaults[name]
v = default
end
mgrsettings:set(name, v)
end
@ -59,7 +39,7 @@ function mapgen_rivers.define_setting(name, dtype, default)
if v == nil then
v = mtsettings:get_bool('mapgen_rivers_' .. name)
if v == nil then
v = defaults[name]
v = default
end
mgrsettings:set_bool(name, v)
end
@ -69,11 +49,10 @@ function mapgen_rivers.define_setting(name, dtype, default)
if v == nil then
v = mtsettings:get_np_group('mapgen_rivers_' .. name)
if v == nil then
v = defaults[name]
v = default
end
mgrsettings:set_np_group(name, v)
end
clean_np(v)
return v
end
end
@ -81,46 +60,36 @@ end
local def_setting = mapgen_rivers.define_setting
mapgen_rivers.settings = {
center = def_setting('center', 'bool'),
blocksize = def_setting('blocksize', 'number'),
center = def_setting('center', 'bool', true),
blocksize = def_setting('blocksize', 'number', 15),
sea_level = tonumber(minetest.get_mapgen_setting('water_level')),
min_catchment = def_setting('min_catchment', 'number'),
river_widening_power = def_setting('river_widening_power', 'number'),
riverbed_slope = def_setting('riverbed_slope', 'number'),
distort = def_setting('distort', 'bool'),
biomes = def_setting('biomes', 'bool'),
glaciers = def_setting('glaciers', 'bool'),
glacier_factor = def_setting('glacier_factor', 'number'),
elevation_chill = def_setting('elevation_chill', 'number'),
min_catchment = def_setting('min_catchment', 'number', 3600),
river_widening_power = def_setting('river_widening_power', 'number', 0.5),
riverbed_slope = def_setting('riverbed_slope', 'number', 0.4),
distort = def_setting('distort', 'bool', true),
biomes = def_setting('biomes', 'bool', true),
glaciers = def_setting('glaciers', 'bool', false),
glacier_factor = def_setting('glacier_factor', 'number', 8),
elevation_chill = def_setting('elevation_chill', 'number', 0.25),
map_x_size = def_setting('map_x_size', 'number'),
map_z_size = def_setting('map_z_size', 'number'),
grid_x_size = def_setting('grid_x_size', 'number', 1000),
grid_z_size = def_setting('grid_z_size', 'number', 1000),
margin = def_setting('margin', 'bool', true),
margin_width = def_setting('margin_width', 'number', 2000),
margin_elev = def_setting('margin_elev', 'number', -200),
evol_params = {
K = def_setting('river_erosion_coef', 'number'),
m = def_setting('river_erosion_power', 'number'),
d = def_setting('diffusive_erosion', 'number'),
compensation_radius = def_setting('compensation_radius', 'number'),
K = def_setting('river_erosion_coef', 'number', 0.5),
m = def_setting('river_erosion_power', 'number', 0.4),
d = def_setting('diffusive_erosion', 'number', 0.5),
compensation_radius = def_setting('compensation_radius', 'number', 50),
},
tectonic_speed = def_setting('tectonic_speed', 'number'),
evol_time = def_setting('evol_time', 'number'),
evol_time_step = def_setting('evol_time_step', 'number'),
tectonic_speed = def_setting('tectonic_speed', 'number', 70),
evol_time = def_setting('evol_time', 'number', 10),
evol_time_step = def_setting('evol_time_step', 'number', 1),
load_all = mtsettings:get_bool('mapgen_rivers_load_all')
}
mapgen_rivers.noise_params = {
base = def_setting("np_base", "noise"),
distort_x = def_setting("np_distort_x", "noise"),
distort_z = def_setting("np_distort_z", "noise"),
distort_amplitude = def_setting("np_distort_amplitude", "noise"),
heat = minetest.get_mapgen_setting_noiseparams('mg_biome_np_heat'),
heat_blend = minetest.get_mapgen_setting_noiseparams('mg_biome_np_heat_blend'),
}
mapgen_rivers.noise_params.heat.offset = mapgen_rivers.noise_params.heat.offset +
mapgen_rivers.settings.sea_level * mapgen_rivers.settings.elevation_chill
local function write_settings()
mgrsettings:write()
end

View File

@ -1,68 +0,0 @@
{
"version": "1.0.2",
"center": true,
"water_level": 1,
"blocksize": 15,
"min_catchment": 3600,
"river_widening_power": 0.5,
"riverbed_slope": 0.4,
"distort": true,
"biomes": true,
"glaciers": false,
"glacier_factor": 8,
"elevation_chill": 0.25,
"map_x_size": 15000,
"map_z_size": 15000,
"river_erosion_coef": 0.5,
"river_erosion_power": 0.4,
"diffusive_erosion": 0.5,
"compensation_radius": 50,
"tectonic_speed": 70,
"evol_time": 10,
"evol_time_step": 1,
"load_all": false,
"np_base": {
"offset": 0,
"scale": 300,
"seed": 2469,
"octaves": 8,
"spread": {"x": 2048, "y": 2048, "z": 2048},
"persist": 0.6,
"lacunarity": 2.0,
"flags": "eased"
},
"np_distort_x": {
"offset": 0,
"scale": 1,
"seed": -4574,
"octaves": 3,
"spread": {"x": 64, "y": 32, "z": 64},
"persist": 0.75,
"lacunarity": 2.0
},
"np_distort_z": {
"offset": 0,
"scale": 1,
"seed": -7940,
"octaves": 3,
"spread": {"x": 64, "y": 32, "z": 64},
"persist": 0.75,
"lacunarity": 2.0
},
"np_distort_amplitude": {
"offset": 0,
"scale": 10,
"seed": 676,
"octaves": 5,
"spread": {"x": 1024, "y": 1024, "z": 1024},
"persist": 0.5,
"lacunarity": 2.0,
"flags": "absvalue"
}
}

View File

@ -3,25 +3,29 @@
# Whether the map should be centered at x=0, z=0.
mapgen_rivers_center (Center map) bool true
# Every cell of the river grid will represent a square of this size.
# A lower value will result in more detailed terrain and finer computation
# of rivers, but will be slower to generate and use more resources.
#
# WARNING: Excessively low values may cause crashes at pre-generation, due to
# memory issues
# Represents horizontal map scale. Every cell of the grid will be upscaled to
# a square of this size.
# For example if the grid size is 1000x1000 and block size is 12,
# the actual size of the map will be 12000.
mapgen_rivers_blocksize (Block size) float 15.0 2.0 100.0
# X size of the map being generated
#
# X size of the river grid will be this/blocksize
# When increasing, it is recommended to increase blocksize too
mapgen_rivers_map_x_size (Map X size) int 15000 500 66000
# X size of the grid being generated
# Actual size of the map is grid_x_size * blocksize
mapgen_rivers_grid_x_size (Grid X size) int 1000 50 5000
# Z size of the map being generated
#
# Z size of the river grid will be this/blocksize
# When increasing, it is recommended to increase blocksize too
mapgen_rivers_map_z_size (Map Z size) int 15000 500 66000
# Z size of the grid being generated
# Actual size of the map is grid_z_size * blocksize
mapgen_rivers_grid_z_size (Grid Z size) int 1000 50 5000
# If margin is enabled, elevation becomes closer to a fixed value when approaching
# the edges of the map.
mapgen_rivers_margin (Margin) bool true
# Width of the transition at map borders, in nodes
mapgen_rivers_margin_width (Margin width) float 2000.0 0.0 15000.0
# Elevation toward which to converge at map borders
mapgen_rivers_margin_elev (Margin elevation) float -200.0 -31000.0 31000.0
# Minimal catchment area for a river to be drawn, in square nodes
# Lower value means bigger river density

View File

@ -1,13 +1,7 @@
-- erosion.lua
-- This is the main file of terrainlib_lua. It registers the EvolutionModel object and some of the
local function erode(model, time)
-- Apply river erosion on the model
-- Erosion model is based on the simplified version of the stream-power law Ey = K×A^m×S
-- where Ey is the vertical erosion speed, A catchment area of the river, S slope along the river, m and K local constants.
-- It is equivalent to considering a horizontal erosion wave travelling at Ex = K×A^m, and this latter approach allows much greather time steps so it is used here.
-- For each point, instead of moving upstream and see what point the erosion wave would reach, we move downstream and see from which point the erosion wave would reach the given point, then we can set the elevation.
--local tinsert = table.insert
local mmin, mmax = math.min, math.max
local dem = model.dem
local dirs = model.dirs
@ -28,9 +22,9 @@ local function erode(model, time)
if scalars then
for i=1, X*Y do
local etime = 1 / (K*rivers[i]^m) -- Inverse of erosion speed (Ex); time needed for the erosion wave to move through the river section.
local etime = 1 / (K*rivers[i]^m)
erosion_time[i] = etime
lakes[i] = mmax(lakes[i], dem[i], sea_level) -- Use lake/sea surface if higher than ground level, because rivers can not erode below.
lakes[i] = mmax(lakes[i], dem[i], sea_level)
end
else
for i=1, X*Y do
@ -45,12 +39,10 @@ local function erode(model, time)
local remaining = time
local new_elev
while true do
-- Explore downstream until we find the point 'iw' from which the erosion wave will reach 'i'
local inext = iw
local d = dirs[iw]
-- Follow the river downstream (move 'iw')
if d == 0 then -- If no flow direction, we reach the border of the map: set elevation to the latest node's elev and abort.
if d == 0 then
new_elev = lakes[iw]
break
elseif d == 1 then
@ -64,13 +56,13 @@ local function erode(model, time)
end
local etime = erosion_time[iw]
if remaining <= etime then -- We have found the node from which the erosion wave will take 'time' to arrive to 'i'.
if remaining <= etime then
local c = remaining / etime
new_elev = (1-c) * lakes[iw] + c * lakes[inext] -- Interpolate linearly between the two nodes
new_elev = (1-c) * lakes[iw] + c * lakes[inext]
break
end
remaining = remaining - etime -- If we still don't reach the target time, decrement time and move to next point.
remaining = remaining - etime
iw = inext
end
@ -79,13 +71,10 @@ local function erode(model, time)
end
local function diffuse(model, time)
-- Apply diffusion using finite differences methods
-- Adapted for small radiuses
local mmax = math.max
local dem = model.dem
local X, Y = dem.X, dem.Y
local d = model.params.d
-- 'd' is equal to 4 times the diffusion coefficient
local dmax = d
if type(d) == "table" then
dmax = -math.huge
@ -95,8 +84,6 @@ local function diffuse(model, time)
end
local diff = dmax * time
-- diff should never exceed 1 per iteration.
-- If needed, we will divide the process in enough iterations so that 'ddiff' is below 1.
local niter = math.floor(diff) + 1
local ddiff = diff / niter
@ -109,7 +96,6 @@ local function diffuse(model, time)
for x=1, X do
local iW = (x==1) and 0 or -1
local iE = (x==X) and 0 or 1
-- Laplacian Δdem × 1/4
temp[i] = (dem[i+iN]+dem[i+iE]+dem[i+iS]+dem[i+iW])*0.25 - dem[i]
i = i + 1
end
@ -119,6 +105,7 @@ local function diffuse(model, time)
dem[i] = dem[i] + temp[i]*ddiff
end
end
-- TODO Test this
end
local modpath = ""
@ -139,7 +126,6 @@ local function flow(model)
end
local function uplift(model, time)
-- Raises the terrain according to uplift rate (model.params.uplift)
local dem = model.dem
local X, Y = dem.X, dem.Y
local uplift_rate = model.params.uplift
@ -156,7 +142,6 @@ local function uplift(model, time)
end
local function noise(model, time)
-- Adds noise to the terrain according to noise depth (model.params.noise)
local random = math.random
local dem = model.dem
local noise_depth = model.params.noise * 2 * time
@ -166,12 +151,6 @@ local function noise(model, time)
end
end
-- Isostasy
-- This is the geological phenomenon that makes the lithosphere "float" over the underlying layers.
-- One of the key implications is that when a very large mass is removed from the ground, the lithosphere reacts by moving upward. This compensation only occurs at large scale (as the lithosphere is not flexible enough for small scale adjustments) so the implementation is using a very large-window Gaussian blur of the elevation array.
-- This implementation is quite simplistic, it does not do a mass balance of the lithosphere as this would introduce too many parameters. Instead, it defines a reference equilibrium elevation, and the ground will react toward this elevation (at the scale of the gaussian window).
-- A change in reference isostasy during the run can also be used to simulate tectonic forcing, like making a new mountain range appear.
local function define_isostasy(model, ref, link)
ref = ref or model.dem
if link then
@ -189,20 +168,18 @@ local function define_isostasy(model, ref, link)
return ref2
end
-- Apply isostasy
local function isostasy(model)
local dem = model.dem
local X, Y = dem.X, dem.Y
local temp = {X=X, Y=Y}
local ref = model.isostasy_ref
for i=1, X*Y do
temp[i] = ref[i] - dem[i] -- Compute the difference between the ground level and the target level
temp[i] = ref[i] - dem[i]
end
-- Blur the difference map using Gaussian blur
gaussian.gaussian_blur_approx(temp, model.params.compensation_radius, 4)
for i=1, X*Y do
dem[i] = dem[i] + temp[i] -- Apply the difference
dem[i] = dem[i] + temp[i]
end
end

View File

@ -71,6 +71,20 @@ local function box_blur_2d(map1, size, map2)
return map1
end
--[[local function gaussian_blur(map, std, tail)
local exp = math.exp
local kernel_mid = math.ceil(std*tail) + 1
local kernel_size = kernel_mid * 2 - 1
local kernel = {}
local cst1 = 1/(std*(2*math.pi)^0.5)
local cst2 = -1/(2*std^2)
for i=1, kernel_size do
kernel[i] = cst1 * exp((i-kernel_mid)^2 * cst2)
end
]]
local function gaussian_blur_approx(map, sigma, n, map2)
map2 = map2 or {}
local sizes = get_box_size(sigma, n)

View File

@ -9,26 +9,24 @@
--
-- Big thanks to them for releasing this paper under a free license ! :)
-- The algorithm here makes use of most of the paper's concepts, including the Planar Borůvka algorithm.
-- The algorithm here makes use of most of the paper's concepts, including the Planar Boruvka algorithm.
-- Only flow_local and accumulate_flow are custom algorithms.
local function flow_local_semirandom(plist)
-- Determines how water should flow at 1 node scale.
-- The straightforward approach would be "Water will flow to the lowest of the 4 neighbours", but here water flows to one of the lower neighbours, chosen randomly, but probability depends on height difference.
-- This makes rivers better follow the curvature of the topography at large scale, and be less biased by pure N/E/S/W directions.
-- 'plist': array of downward height differences (0 if upward)
local sum = 0
for i=1, #plist do
sum = sum + plist[i] -- Sum of probabilities
sum = sum + plist[i]
end
--for _, p in ipairs(plist) do
--sum = sum + p
--end
if sum == 0 then
return 0
end
local r = math.random() * sum
for i=1, #plist do
local p = plist[i]
--for i, p in ipairs(plist) do
if r < p then
return i
end
@ -37,14 +35,11 @@ local function flow_local_semirandom(plist)
return 0
end
-- Maybe implement more flow methods in the future?
local flow_methods = {
semirandom = flow_local_semirandom,
}
-- Applies all steps of the flow routing, to calculate flow direction for every node, and lake surface elevation.
-- It's quite a hard piece of code, but we will go step by step and explain what's going on, so stay with me and... let's goooooooo!
local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are optional tables to reuse for memory optimization, they may contain any data.
local function flow_routing(dem, dirs, lakes, method)
method = method or 'semirandom'
local flow_local = flow_methods[method] or flow_local_semirandom
@ -52,6 +47,7 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
lakes = lakes or {}
-- Localize for performance
--local tinsert = table.insert
local tremove = table.remove
local mmax = math.max
@ -66,15 +62,11 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
dirs2[i] = 0
end
----------------------------------------
-- STEP 1: Find local flow directions --
----------------------------------------
-- Use the local flow function and fill the flow direction tables
local singular = {}
for y=1, Y do
for x=1, X do
local zi = dem[i]
local plist = { -- Get the height difference of the 4 neighbours (and 0 if uphill)
local plist = {
y<Y and mmax(zi-dem[i+X], 0) or 0, -- Southward
x<X and mmax(zi-dem[i+1], 0) or 0, -- Eastward
y>1 and mmax(zi-dem[i-X], 0) or 0, -- Northward
@ -82,10 +74,8 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
}
local d = flow_local(plist)
-- 'dirs': Direction toward which water flow
-- 'dirs2': Directions from which water comes
dirs[i] = d
if d == 0 then -- If water can't flow from this node, add it to the list of singular nodes that will be resolved later
if d == 0 then
singular[#singular+1] = i
elseif d == 1 then
dirs2[i+X] = dirs2[i+X] + 1
@ -100,39 +90,30 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
end
end
--------------------------------------
-- STEP 2: Compute basins and links --
--------------------------------------
-- Now water can flow until it reaches a singular node (which is in most cases the bottom of a depression)
-- We will calculate the drainage basin of every singular node (all the nodes from which the water will flow in this singular node, directly or indirectly), make an adjacency list of basins, and find the lowest pass between each pair of adjacent basins (they are potential lake outlets)
-- Compute basins and links
local nbasins = #singular
local basin_id = {}
local links = {}
local basin_links
-- Function to analyse a link between two nodes
local function add_link(i1, i2, b1, isY)
-- i1, i2: coordinates of two nodes
-- b1: basin that contains i1
-- isY: whether the link is in Y direction
local b2
-- Note that basin number #0 represents the outside of the map; or if the coordinate is inside the map, means that the basin number is uninitialized.
if i2 == 0 then -- If outside the map
if i2 == 0 then
b2 = 0
else
b2 = basin_id[i2]
if b2 == 0 then -- If basin of i2 is not already computed, skip
if b2 == 0 then
return
end
end
if b2 ~= b1 then -- If these two nodes don't belong to the same basin, we have found a link between two adjacent basins
local elev = i2 == 0 and dem[i1] or mmax(dem[i1], dem[i2]) -- Elevation of the highest of the two sides of the link (or only i1 if b2 is map outside)
if b2 ~= b1 then
local elev = i2 == 0 and dem[i1] or mmax(dem[i1], dem[i2])
local l2 = basin_links[b2]
if not l2 then
l2 = {}
basin_links[b2] = l2
end
if not l2.elev or l2.elev > elev then -- If this link is lower than the lowest registered link between these two basins, register it as the new lowest pass
if not l2.elev or l2.elev > elev then
l2.elev = elev
l2.i = mmax(i1,i2)
l2.is_y = isY
@ -145,48 +126,51 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
for i=1, X*Y do
basin_id[i] = 0
end
--for ib, s in ipairs(singular) do
for ib=1, nbasins do
-- Here we will recursively search upstream from the singular node to determine its drainage basin
local queue = {singular[ib]} -- Start with the singular node, then this queue will be filled with water donors neighbours
--local s = singular[ib]
local queue = {singular[ib]}
basin_links = {}
links[#links+1] = basin_links
--tinsert(links, basin_links)
while #queue > 0 do
local i = tremove(queue)
basin_id[i] = ib
local d = dirs2[i] -- Get the directions water is coming from
local d = dirs2[i]
-- Iterate through the 4 directions
if d >= 8 then -- River coming from the East
if d >= 8 then -- River coming from East
d = d - 8
queue[#queue+1] = i+1
-- If no river is coming from the East, we might be at the limit of two basins, thus we need to test adjacency.
--tinsert(queue, i+X)
elseif i%X > 0 then
add_link(i, i+1, ib, false)
else -- If the eastern neighbour is outside the map
else
add_link(i, 0, ib, false)
end
if d >= 4 then -- River coming from the South
if d >= 4 then -- River coming from South
d = d - 4
queue[#queue+1] = i+X
--tinsert(queue, i+1)
elseif i <= X*(Y-1) then
add_link(i, i+X, ib, true)
else
add_link(i, 0, ib, true)
end
if d >= 2 then -- River coming from the West
if d >= 2 then -- River coming from West
d = d - 2
queue[#queue+1] = i-1
--tinsert(queue, i-X)
elseif i%X ~= 1 then
add_link(i, i-1, ib, false)
else
add_link(i, 0, ib, false)
end
if d >= 1 then -- River coming from the North
if d >= 1 then -- River coming from North
queue[#queue+1] = i-X
--tinsert(queue, i-1)
elseif i > X then
add_link(i, i-X, ib, true)
else
@ -202,7 +186,7 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
nlinks[i] = 0
end
-- Iterate through pairs of adjacent basins, and make the links reciprocal
--for ib1, blinks in ipairs(links) do
for ib1=1, #links do
for ib2, link in pairs(links[ib1]) do
if ib2 < ib1 then
@ -213,15 +197,6 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
end
end
-----------------------------------------------------
-- STEP 3: Compute minimal spanning tree of basins --
-----------------------------------------------------
-- We've got an adjacency list of basins with the elevation of their links.
-- We will build a minimal spanning tree of the basins (where costs are the elevation of the links). As demonstrated by Cordonnier et al., this finds the outlets of the basins, where water would naturally flow. This does not tell in which direction water is flowing, however.
-- We will use a version of Borůvka's algorithm, with Mareš' optimizations to approach linear complexity (see paper).
-- The concept of Borůvka's algorithm is to take elements and merge them with their lowest neighbour, until all elements are merged.
-- Mareš' optimizations mainly consist in skipping elements that have over 8 links, until extra links are removed when other elements are merged.
-- Note that for this step we are only working on basins, not grid nodes.
local lowlevel = {}
for i, n in pairs(nlinks) do
if n <= 8 then
@ -231,8 +206,6 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
local basin_graph = {}
for n=1, nbasins do
-- Iterate in lowlevel but its contents may change during the loop
-- 'next' called with only one argument always returns an element if table is not empty
local b1, lnk1 = next(lowlevel)
lowlevel[b1] = nil
@ -240,7 +213,6 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
local lowest = math.huge
local lnk1 = links[b1]
local i = 0
-- Look for lowest link
for bn, bdata in pairs(lnk1) do
i = i + 1
if bdata.elev < lowest then
@ -249,7 +221,7 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
end
end
-- Add link to the graph, in both directions
-- Add link to the graph
local bound = lnk1[b2]
local bb1, bb2 = bound[1], bound[2]
if not basin_graph[bb1] then
@ -267,7 +239,6 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
lnk1[b2] = nil
lnk2[b1] = nil
nlinks[b2] = nlinks[b2] - 1
-- When the number of links is changing, we need to check whether the basin can be added to / removed from 'lowlevel'
if nlinks[b2] == 8 then
lowlevel[b2] = lnk2
end
@ -276,32 +247,25 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
local lnkn = links[bn]
lnkn[b1] = nil
if lnkn[b2] then -- If bassin bn is also linked to b2
nlinks[bn] = nlinks[bn] - 1 -- Then bassin bn is losing a link because it keeps only one link toward b1/b2 after the merge
if lnkn[b2] then
nlinks[bn] = nlinks[bn] - 1
if nlinks[bn] == 8 then
lowlevel[bn] = lnkn
end
else -- If bn was linked to b1 but not to b2
nlinks[b2] = nlinks[b2] + 1 -- Then b2 is gaining a link to bn because of the merge
else
nlinks[b2] = nlinks[b2] + 1
if nlinks[b2] == 9 then
lowlevel[b2] = nil
end
end
if not lnkn[b2] or lnkn[b2].elev > bdata.elev then -- If the link b1-bn will become the new lowest link between b2 and bn, redirect the link to b2
if not lnkn[b2] or lnkn[b2].elev > bdata.elev then
lnkn[b2] = bdata
lnk2[bn] = bdata
end
end
end
--------------------------------------------------------------
-- STEP 4: Orient basin graph, and grid nodes inside basins --
--------------------------------------------------------------
-- We will finally solve those freaking singular nodes.
-- To orient the basin graph, we will consider that the ultimate basin water should flow into is the map outside (basin #0). We will start from it and recursively walk upstream to the neighbouring basins, using only links that are in the minimal spanning tree. This gives the flow direction of the links, and thus, the outlet of every basin.
-- This will also give lake elevation, which is the highest link encountered between map outside and the given basin on the spanning tree.
-- And within each basin, we need to modify flow directions to connect the singular node to the outlet.
local queue = {[0] = -math.huge}
local basin_lake = {}
for n=1, nbasins do
@ -309,17 +273,15 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
end
local reverse = {3, 4, 1, 2, [0]=0}
for n=1, nbasins do
local b1, elev1 = next(queue) -- Pop from queue
local b1, elev1 = next(queue)
queue[b1] = nil
basin_lake[b1] = elev1
-- Iterate through b1's neighbours (according to the spanning tree)
for b2, bound in pairs(basin_graph[b1]) do
-- Make b2 flow into b1
local i = bound.i -- Get the coordinate of the link (which is the basin's outlet)
local dir = bound.is_y and 3 or 4 -- And get the direction (S/E/N/W)
local i = bound.i
local dir = bound.is_y and 3 or 4
if basin_id[i] ~= b2 then
dir = dir - 2
-- Coordinate 'i' refers to the side of the link with the highest X/Y position. In case it is in the wrong basin, take the other side by decrementing by one row/column.
if bound.is_y then
i = i - X
else
@ -329,12 +291,8 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
dir = 0
end
-- Use the flow directions computed in STEP 2 to find the route from the outlet position to the singular node, and reverse this route to make the singular node flow into the outlet
-- This can make the river flow uphill, which may seem unnatural, but it can only happen below a lake (because outlet elevation defines lake surface elevation)
repeat
-- Assign i's direction to 'dir', and get i's former direction
dir, dirs[i] = dirs[i], dir
-- Move i by following its former flow direction (downstream)
if dir == 1 then
i = i + X
elseif dir == 2 then
@ -344,36 +302,26 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
elseif dir == 4 then
i = i - 1
end
-- Reverse the flow direction for the next node, which will flow into i
dir = reverse[dir]
until dir == 0 -- Stop when reaching the singular node
-- Add basin b2 into the queue, and keep the highest link elevation, that will define the elevation of the lake in b2
until dir == 0
-- Add b2 into the queue
queue[b2] = mmax(elev1, bound.elev)
-- Remove b1 from b2's neighbours to avoid coming back to b1
basin_graph[b2][b1] = nil
end
basin_graph[b1] = nil
end
-- Every node will be assigned the lake elevation of the basin it belongs to.
-- If lake elevation is lower than ground elevation, it simply means that there is no lake here.
for i=1, X*Y do
lakes[i] = basin_lake[basin_id[i]]
end
-- That's it!
return dirs, lakes
end
local function accumulate(dirs, waterq)
-- Calculates the river flow by determining the surface of the catchment area for every node
-- This means: how many nodes will give their water to that given node, directly or indirectly?
-- This is obtained by following rivers downstream and summing up the flow of every tributary, starting with a value of 1 at the sources.
-- This will give non-zero values for every node but only large values will be considered to be rivers.
waterq = waterq or {}
local X, Y = dirs.X, dirs.Y
--local tinsert = table.insert
local ndonors = {}
local waterq = {X=X, Y=Y}
@ -382,7 +330,7 @@ local function accumulate(dirs, waterq)
waterq[i] = 1
end
-- Calculate the number of direct donors
--for i1, dir in ipairs(dirs) do
for i1=1, X*Y do
local i2
local dir = dirs[i1]
@ -401,12 +349,10 @@ local function accumulate(dirs, waterq)
end
for i1=1, X*Y do
-- Find sources (nodes that have no donor)
if ndonors[i1] == 0 then
local i2 = i1
local dir = dirs[i2]
local w = waterq[i2]
-- Follow the water flow downstream: move 'i2' to the next node according to its flow direction
while dir > 0 do
if dir == 1 then
i2 = i2 + X
@ -417,12 +363,9 @@ local function accumulate(dirs, waterq)
elseif dir == 4 then
i2 = i2 - 1
end
-- Increment the water quantity of i2
w = w + waterq[i2]
waterq[i2] = w
-- Stop on an unresolved confluence (node with >1 donors) and decrease the number of remaining donors
-- When the ndonors of a confluence has decreased to 1, it means that its water quantity has already been incremented by its tributaries, so it can be resolved like a standard river section. However, do not decrease ndonors to zero to avoid considering it as a source.
if ndonors[i2] > 1 then
ndonors[i2] = ndonors[i2] - 1
break

View File

@ -1,4 +1,4 @@
-- twist.lua
-- bounds.lua
local function get_bounds(dirs, rivers)
local X, Y = dirs.X, dirs.Y