CHANGELOG.md

@@ -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

README.md

@@ -1,5 +1,5 @@
 # Map Generator with Rivers
-`mapgen_rivers v1.0.2` by Gaël de Sailly.
+`mapgen_rivers v1.0` 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.

compatibility.lua

@@ -1,24 +1,3 @@
--- Fix compatibility for settings-related changes
--- Only loaded if the versions of the mod and the world mismatch
-
-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 ""
 
@@ -42,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)
@@ -62,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

geometry.lua

@@ -0,0 +1,37 @@
+local function distance_to_segment(x1, y1, x2, y2, x, y)
+	-- get the distance between point (x,y) and segment (x1,y1)-(x2,y2)
+	local a = (x1-x2)^2 + (y1-y2)^2 -- square of distance
+	local b = (x1-x)^2 + (y1-y)^2
+	local c = (x2-x)^2 + (y2-y)^2
+	if a + b < c then
+		-- The closest point of the segment is the extremity 1
+		return math.sqrt(b)
+	elseif a + c < b then
+		-- The closest point of the segment is the extremity 2
+		return math.sqrt(c)
+	else
+		-- The closest point is on the segment
+		return math.abs(x1 * (y2-y) + x2 * (y-y1) + x * (y1-y2)) / math.sqrt(a)
+	end
+end
+
+local function transform_quadri(X, Y, x, y)
+	-- To index points in an irregular quadrilateral, giving x and y between 0 (one edge) and 1 (opposite edge)
+	-- X, Y 4-vectors giving the coordinates of the 4 vertices
+	-- x, y position to index.
+	local x1, x2, x3, x4 = unpack(X)
+	local y1, y2, y3, y4 = unpack(Y)
+
+	-- Compare distance to 2 opposite edges, they give the X coordinate
+	local d23 = distance_to_segment(x2,y2,x3,y3,x,y)
+	local d41 = distance_to_segment(x4,y4,x1,y1,x,y)
+	local xc = d41 / (d23+d41)
+
+	-- Same for the 2 other edges, they give the Y coordinate
+	local d12 = distance_to_segment(x1,y1,x2,y2,x,y)
+	local d34 = distance_to_segment(x3,y3,x4,y4,x,y)
+	local yc = d12 / (d12+d34)
+	return xc, yc
+end
+
+return transform_quadri

gridmanager.lua

@@ -1,105 +0,0 @@
--- Manages grid loading, writing and generation
-
-local datapath = mapgen_rivers.world_data_path
-
-local registered_on_grid_loaded = {}
-function mapgen_rivers.register_on_grid_loaded(func)
-	if type(func) == "function" then
-		registered_on_grid_loaded[#registered_on_grid_loaded+1] = func
-	else
-		minetest.log("error", "[mapgen_rivers] register_on_grid_loaded can only register functions!")
-	end
-end
-
-local function on_grid_loaded_callback(grid)
-	for _, func in ipairs(registered_on_grid_loaded) do
-		func(grid)
-	end
-end
-
-local function offset_conv(o)
-	return (o + 0.5) * (1/256)
-end
-
-local floor = math.floor
-local sbyte, schar = string.byte, string.char
-local unpk = unpack
-
--- Loading files
-
--- Never load the full map during mapgen. Instead, create an empty lookup table
--- and read the file on-the-fly when an element is requested for the first time,
--- using __index metamethod.
-local loader_mt = {
-	__index = function(loader, i) -- Called when accessing a missing key
-		local file = loader.file
-		local bytes = loader.bytes
-		file:seek('set', (i-1)*bytes)
-		local strnum = file:read(bytes)
-
-		local n = sbyte(strnum, 1)
-		if loader.signed and n >= 128 then
-			n = n - 256
-		end
-
-		for j=2, bytes do
-			n = n*256 + sbyte(strnum, j)
-		end
-
-		if loader.conv then
-			n = loader.conv(n)
-		end
-		-- Cache key for next use
-		loader[i] = n
-		return n
-	end,
-}
-
-local function load_file(filename, bytes, signed, size, converter)
-	local file = io.open(datapath .. filename, 'rb')
-	if file then
-		converter = converter or false
-		return setmetatable({file=file, bytes=bytes, signed=signed, size=size, conv=converter}, loader_mt)
-	end
-end
-
-function mapgen_rivers.load_or_generate_grid()
-	-- First, check whether a grid is already loaded
-	if mapgen_rivers.grid then
-		return true
-	end
-
-	-- If not, try to load the grid from the files
-	local sfile = io.open(datapath .. 'size', 'r')
-	if not sfile then
-		dofile(mapgen_rivers.modpath .. "/pregenerate.lua")
-		collectgarbage()
-
-		sfile = io.open(datapath .. 'size', 'r')
-		if not sfile then
-			return false
-		end
-	end
-
-	local x, z = sfile:read('*n'), sfile:read('*n')
-	if not x or not z then
-		return false
-	end
-
-	minetest.log("action", '[mapgen_rivers] Loading grid')
-
-	local grid = {
-		size = {x=x, y=z},
-		dem = load_file('dem', 2, true, x*z),
-		lakes = load_file('lakes', 2, true, x*z),
-		dirs = load_file('dirs', 1, false, x*z),
-		rivers = load_file('rivers', 4, false, x*z),
-		offset_x = load_file('offset_x', 1, true, x*z, offset_conv),
-		offset_y = load_file('offset_y', 1, true, x*z, offset_conv),
-	}
-
-	mapgen_rivers.grid = grid
-	on_grid_loaded_callback(grid)
-
-	return true
-end

heightmap.lua

@@ -1,52 +1,12 @@
--- Transform polygon data into a heightmap
-
 local modpath = mapgen_rivers.modpath
 
+local make_polygons = dofile(modpath .. 'polygons.lua')
+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 out_elev = mapgen_rivers.settings.margin_elev
+local MAP_BOTTOM = -31000
-
--- Localize for performance
-local floor, min, max, sqrt, abs = math.floor, math.min, math.max, math.sqrt, math.abs
-local unpk = unpack
-
--- Geometrical helpers
-local function distance_to_segment(x1, y1, x2, y2, x, y)
-	-- get the distance between point (x,y) and segment (x1,y1)-(x2,y2)
-	local a = (x1-x2)^2 + (y1-y2)^2 -- square of distance
-	local b = (x1-x)^2 + (y1-y)^2
-	local c = (x2-x)^2 + (y2-y)^2
-	if a + b < c then
-		-- The closest point of the segment is the extremity 1
-		return sqrt(b)
-	elseif a + c < b then
-		-- The closest point of the segment is the extremity 2
-		return sqrt(c)
-	else
-		-- The closest point is on the segment
-		return abs(x1 * (y2-y) + x2 * (y-y1) + x * (y1-y2)) / sqrt(a)
-	end
-end
-
-local function transform_quadri(X, Y, x, y)
-	-- To index points in an irregular quadrilateral, giving x and y between 0 (one edge) and 1 (opposite edge)
-	-- X, Y 4-vectors giving the coordinates of the 4 vertices
-	-- x, y position to index.
-	local x1, x2, x3, x4 = unpk(X)
-	local y1, y2, y3, y4 = unpk(Y)
-
-	-- Compare distance to 2 opposite edges, they give the X coordinate
-	local d23 = distance_to_segment(x2,y2,x3,y3,x,y)
-	local d41 = distance_to_segment(x4,y4,x1,y1,x,y)
-	local xc = d41 / (d23+d41)
-
-	-- Same for the 2 other edges, they give the Y coordinate
-	local d12 = distance_to_segment(x1,y1,x2,y2,x,y)
-	local d34 = distance_to_segment(x3,y3,x4,y4,x,y)
-	local yc = d12 / (d12+d34)
-	return xc, yc
-end
 
 -- Linear interpolation
 local function interp(v00, v01, v11, v10, xf, zf)
@@ -55,9 +15,9 @@ local function interp(v00, v01, v11, v10, xf, zf)
 	return v1*zf + v0*(1-zf)
 end
 
-function mapgen_rivers.make_heightmaps(minp, maxp)
+local function heightmaps(minp, maxp)
 
-	local polygons = mapgen_rivers.make_polygons(minp, maxp)
+	local polygons = make_polygons(minp, maxp)
 	local incr = maxp.z-minp.z+1
 
 	local terrain_height_map = {}
@@ -70,11 +30,11 @@ function mapgen_rivers.make_heightmaps(minp, maxp)
 
 			if poly then
 				local xf, zf = transform_quadri(poly.x, poly.z, x, z)
-				local i00, i01, i11, i10 = unpk(poly.i)
+				local i00, i01, i11, i10 = unpack(poly.i)
 
 				-- Load river width on 4 edges and corners
-				local r_west, r_north, r_east, r_south = unpk(poly.rivers)
+				local r_west, r_north, r_east, r_south = unpack(poly.rivers)
-				local c_NW, c_NE, c_SE, c_SW = unpk(poly.river_corners)
+				local c_NW, c_NE, c_SE, c_SW = unpack(poly.river_corners)
 
 				-- Calculate the depth factor for each edge and corner.
 				-- Depth factor:
@@ -82,72 +42,65 @@ function mapgen_rivers.make_heightmaps(minp, maxp)
 				-- = 0: on riverbank
 				-- > 0: inside river
 				local depth_factors = {
-					r_west  -    xf , -- West edge (1)
+					r_west - xf,
-					r_north -    zf , -- North edge (2)
+					r_north - zf,
-					r_east  - (1-xf), -- East edge (3)
+					xf - r_east,
-					r_south - (1-zf), -- South edge (4)
+					zf - r_south,
-					c_NW -    xf  -    zf , -- North-West corner (5)
+					c_NW-xf-zf,
-					c_NE - (1-xf) -    zf , -- North-East corner (6)
+					xf-zf-c_NE,
-					c_SE - (1-xf) - (1-zf), -- South-East corner (7)
+					xf+zf-c_SE,
-					c_SW -    xf  - (1-zf), -- South-West corner (8)
+					zf-xf-c_SW,
 				}
 
-				-- Find the maximal depth factor, which determines to which of the 8 river sections (4 edges + 4 corners) the current point belongs.
+				-- Find the maximal depth factor and determine to which river it belongs
-				-- If imax is still at 0, it means that we are not in a river.
+				local depth_factor_max = 0
-				local dpmax = 0
 				local imax = 0
 				for i=1, 8 do
-					if depth_factors[i] > dpmax then
+					if depth_factors[i] >= depth_factor_max then
-						dpmax = depth_factors[i]
+						depth_factor_max = depth_factors[i]
 						imax = i
 					end
 				end
 
-				-- Transform the coordinates to have xfc and zfc = 0 or 1 in rivers (to avoid rivers having lateral slope and to accomodate the riverbanks smoothly)
+				-- Transform the coordinates to have xf and zf = 0 or 1 in rivers (to avoid rivers having lateral slope and to accomodate the surrounding smoothly)
-				local xfc, zfc
+				if imax == 0 then
-				-- xfc:
+					local x0 = math.max(r_west, c_NW-zf, zf-c_SW)
-				if imax == 0 or imax == 2 or imax == 4 then -- river segment does not constrain X coordinate, so accomodate xf in function of other river sections
+					local x1 = math.min(r_east, c_NE+zf, c_SE-zf)
-					local x0 =   max(r_west-dpmax, c_NW-zf-dpmax, c_SW-(1-zf)-dpmax, 0) -- new xf will be bounded to 0 by western riverbank
+					local z0 = math.max(r_north, c_NW-xf, xf-c_NE)
-					local x1 = 1-max(r_east-dpmax, c_NE-zf-dpmax, c_SE-(1-zf)-dpmax, 0) -- and bounded to 1 by eastern riverbank
+					local z1 = math.min(r_south, c_SW+xf, c_SE-xf)
-					if x0 >= x1 then
+					xf = (xf-x0) / (x1-x0)
-						xfc = 0.5
+					zf = (zf-z0) / (z1-z0)
-					else
+				elseif imax == 1 then
-						xfc = (xf-x0) / (x1-x0)
+					xf = 0
-					end
+				elseif imax == 2 then
-				elseif imax == 1 or imax == 5 or imax == 8 then -- river at the western side of the polygon
+					zf = 0
-					xfc = 0
+				elseif imax == 3 then
-				else -- 3, 6, 7 : river at the eastern side of the polygon
+					xf = 1
-					xfc = 1
+				elseif imax == 4 then
-				end
+					zf = 1
-
+				elseif imax == 5 then
-				-- Same for zfc:
+					xf, zf = 0, 0
-				if imax == 0 or imax == 1 or imax == 3 then -- river segment does not constrain Z coordinate, so accomodate zf in function of other river sections
+				elseif imax == 6 then
-					local z0 =   max(r_north-dpmax, c_NW-xf-dpmax, c_NE-(1-xf)-dpmax, 0) -- new zf will be bounded to 0 by northern riverbank
+					xf, zf = 1, 0
-					local z1 = 1-max(r_south-dpmax, c_SW-xf-dpmax, c_SE-(1-xf)-dpmax, 0) -- and bounded to 1 by southern riverbank
+				elseif imax == 7 then
-					if z0 >= z1 then
+					xf, zf = 1, 1
-						zfc = 0.5
+				elseif imax == 8 then
-					else
+					xf, zf = 0, 1
-						zfc = (zf-z0) / (z1-z0)
-					end
-				elseif imax == 2 or imax == 5 or imax == 6 then -- river at the northern side of the polygon
-					zfc = 0
-				else -- 4, 7, 8 : river at the southern side of the polygon
-					zfc = 1
 				end
 
 				-- Determine elevation by interpolation
 				local vdem = poly.dem
-				local terrain_height = floor(0.5+interp(
+				local terrain_height = math.floor(0.5+interp(
 					vdem[1],
 					vdem[2],
 					vdem[3],
 					vdem[4],
-					xfc, zfc
+					xf, zf
 				))
 
 				-- Spatial gradient of the interpolation
-				local slope_x = zfc*(vdem[3]-vdem[4]) + (1-zfc)*(vdem[2]-vdem[1]) < 0
+				local slope_x = zf*(vdem[3]-vdem[4]) + (1-zf)*(vdem[2]-vdem[1]) < 0
-				local slope_z = xfc*(vdem[3]-vdem[2]) + (1-xfc)*(vdem[4]-vdem[1]) < 0
+				local slope_z = xf*(vdem[3]-vdem[2]) + (1-xf)*(vdem[4]-vdem[1]) < 0
 				local lake_id = 0
 				if slope_x then
 					if slope_z then
@@ -162,17 +115,17 @@ function mapgen_rivers.make_heightmaps(minp, maxp)
 						lake_id = 1
 					end
 				end
-				local lake_height = max(floor(poly.lake[lake_id]), terrain_height)
+				local lake_height = math.max(math.floor(poly.lake[lake_id]), terrain_height)
 
-				if imax > 0 and dpmax > 0 then
+				if imax > 0 and depth_factor_max > 0 then
-					terrain_height = min(max(lake_height, sea_level) - floor(1+dpmax*riverbed_slope), terrain_height)
+					terrain_height = math.min(math.max(lake_height, sea_level) - math.floor(1+depth_factor_max*riverbed_slope), terrain_height)
 				end
 
 				terrain_height_map[i] = terrain_height
 				lake_height_map[i] = lake_height
 			else
-				terrain_height_map[i] = out_elev
+				terrain_height_map[i] = MAP_BOTTOM
-				lake_height_map[i] = out_elev
+				lake_height_map[i] = MAP_BOTTOM
 			end
 			i = i + 1
 		end
@@ -180,3 +133,5 @@ function mapgen_rivers.make_heightmaps(minp, maxp)
 
 	return terrain_height_map, lake_height_map
 end
+
+return heightmaps

init.lua

@@ -1,5 +1,3 @@
--- Main file, calls the other files and triggers main functions
-
 mapgen_rivers = {}
 
 local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
@@ -7,10 +5,241 @@ mapgen_rivers.modpath = modpath
 mapgen_rivers.world_data_path = minetest.get_worldpath() .. '/river_data/'
 
 dofile(modpath .. 'settings.lua')
-dofile(modpath .. 'gridmanager.lua')
-dofile(modpath .. 'polygons.lua')
-dofile(modpath .. 'heightmap.lua')
-dofile(modpath .. 'mapgen.lua')
-dofile(modpath .. 'spawn.lua')
 
-minetest.register_on_mods_loaded(mapgen_rivers.load_or_generate_grid)
+local sea_level = mapgen_rivers.settings.sea_level
+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 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')
+
+-- Linear interpolation
+local function interp(v00, v01, v11, v10, xf, zf)
+	local v0 = v01*xf + v00*(1-xf)
+	local v1 = v11*xf + v10*(1-xf)
+	return v1*zf + v0*(1-zf)
+end
+
+local data = {}
+
+local noise_x_obj, noise_z_obj, noise_distort_obj, noise_heat_obj, noise_heat_blend_obj
+local noise_x_map = {}
+local noise_z_map = {}
+local noise_distort_map = {}
+local noise_heat_map = {}
+local noise_heat_blend_map = {}
+local mapsize
+local init = false
+
+local sumtime = 0
+local sumtime2 = 0
+local ngen = 0
+
+local function generate(minp, maxp, seed)
+	print(("[mapgen_rivers] Generating from %s to %s"):format(minetest.pos_to_string(minp), minetest.pos_to_string(maxp)))
+
+	local chulens = {
+		x = maxp.x-minp.x+1,
+		y = maxp.y-minp.y+1,
+		z = maxp.z-minp.z+1,
+	}
+
+	if not init then
+		mapsize = {
+			x = chulens.x,
+			y = chulens.y+1,
+			z = chulens.z,
+		}
+		if use_distort then
+			noise_x_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.distort_x, mapsize)
+			noise_z_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.distort_z, mapsize)
+			noise_distort_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.distort_amplitude, chulens)
+		end
+		if use_biomes then
+			noise_heat_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.heat, chulens)
+			noise_heat_blend_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.heat_blend, chulens)
+		end
+		init = true
+	end
+
+	local t0 = os.clock()
+	local minp2d = {x=minp.x, y=minp.z}
+	if use_distort then
+		noise_x_obj:get_3d_map_flat(minp, noise_x_map)
+		noise_z_obj:get_3d_map_flat(minp, noise_z_map)
+		noise_distort_obj:get_2d_map_flat(minp2d, noise_distort_map)
+	end
+	if use_biomes then
+		noise_heat_obj:get_2d_map_flat(minp2d, noise_heat_map)
+		noise_heat_blend_obj:get_2d_map_flat(minp2d, noise_heat_blend_map)
+	end
+
+	local terrain_map, lake_map, incr, i_origin
+
+	if use_distort then
+		local xmin, xmax, zmin, zmax = minp.x, maxp.x, minp.z, maxp.z
+		local i = 0
+		local i2d = 0
+		for z=minp.z, maxp.z do
+			for y=minp.y, maxp.y+1 do
+				for x=minp.x, maxp.x do
+					i = i+1
+					i2d = i2d+1
+					local distort = noise_distort_map[i2d]
+					local xv = noise_x_map[i]*distort + x
+					if xv < xmin then xmin = xv end
+					if xv > xmax then xmax = xv end
+					noise_x_map[i] = xv
+					local zv = noise_z_map[i]*distort + z
+					if zv < zmin then zmin = zv end
+					if zv > zmax then zmax = zv end
+					noise_z_map[i] = zv
+				end
+				i2d = i2d-chulens.x
+			end
+		end
+
+		local pminp = {x=math.floor(xmin), z=math.floor(zmin)}
+		local pmaxp = {x=math.floor(xmax)+1, z=math.floor(zmax)+1}
+		incr = pmaxp.x-pminp.x+1
+		i_origin = 1 - pminp.z*incr - pminp.x
+		terrain_map, lake_map = heightmaps(pminp, pmaxp)
+	else
+		terrain_map, lake_map = heightmaps(minp, maxp)
+	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)
+
+	local a = VoxelArea:new({MinEdge = emin, MaxEdge = emax})
+	local ystride = a.ystride -- Tip : the ystride of a VoxelArea is the number to add to the array index to get the index of the position above. It's faster because it avoids to completely recalculate the index.
+
+	local nid = mapsize.x*(mapsize.y-1) + 1
+	local incrY = -mapsize.x
+	local incrX = 1 - mapsize.y*incrY
+	local incrZ = mapsize.x*mapsize.y - mapsize.x*incrX - mapsize.x*mapsize.y*incrY
+
+	local i2d = 1
+
+	for z = minp.z, maxp.z do
+		for x = minp.x, maxp.x do
+			local ivm = a:index(x, minp.y, z)
+			local ground_above = false
+			local temperature
+			if use_biomes then
+				temperature = noise_heat_map[i2d]+noise_heat_blend_map[i2d]
+			end
+			local terrain, lake
+			if not use_distort then
+				terrain = terrain_map[i2d]
+				lake = lake_map[i2d]
+			end
+
+			for y = maxp.y+1, minp.y, -1 do
+				if use_distort then
+					local xn = noise_x_map[nid]
+					local zn = noise_z_map[nid]
+					local x0 = math.floor(xn)
+					local z0 = math.floor(zn)
+
+					local i0 = i_origin + z0*incr + x0
+					local i1 = i0+1
+					local i2 = i1+incr
+					local i3 = i2-1
+
+					terrain = interp(terrain_map[i0], terrain_map[i1], terrain_map[i2], terrain_map[i3], xn-x0, zn-z0)
+					lake = math.min(lake_map[i0], lake_map[i1], lake_map[i2], lake_map[i3])
+				end
+
+				if y <= maxp.y then
+
+					local is_lake = lake > terrain
+					local ivm = a:index(x, y, z)
+					if y <= terrain then
+						if not use_biomes or y <= terrain-1 or ground_above then
+							data[ivm] = c_stone
+						elseif is_lake or y < sea_level then
+							data[ivm] = c_sand
+						else
+							local temperature_y = temperature - y*elevation_chill
+							if temperature_y >= 15 then
+								data[ivm] = c_lawn
+							elseif temperature_y >= 0 then
+								data[ivm] = c_dirtsnow
+							else
+								data[ivm] = c_snow
+							end
+						end
+					elseif y <= lake and lake > sea_level then
+						if not use_biomes or temperature - y*elevation_chill >= 0 then
+							data[ivm] = c_rwater
+						else
+							data[ivm] = c_ice
+						end
+					elseif y <= sea_level then
+						data[ivm] = c_water
+					end
+				end
+
+				ground_above = y <= terrain
+
+				ivm = ivm + ystride
+				if use_distort then
+					nid = nid + incrY
+				end
+			end
+
+			if use_distort then
+				nid = nid + incrX
+			end
+			i2d = i2d + 1
+		end
+
+		if use_distort then
+			nid = nid + incrZ
+		end
+	end
+
+	if use_biomegen_mod then
+		biomegen.generate_all(data, a, vm, minp, maxp, seed)
+	else
+		vm:set_data(data)
+		minetest.generate_ores(vm, minp, maxp)
+	end
+
+	vm:set_lighting({day = 0, night = 0})
+	vm:calc_lighting()
+	vm:update_liquids()
+	vm:write_to_map()
+	local t1 = os.clock()
+
+	local t = t1-t0
+	ngen = ngen + 1
+	sumtime = sumtime + t
+	sumtime2 = sumtime2 + t*t
+	print(("[mapgen_rivers] Done in %5.3f s"):format(t))
+end
+
+minetest.register_on_generated(generate)
+minetest.register_on_shutdown(function()
+	local avg = sumtime / ngen
+	local std = math.sqrt(sumtime2/ngen - avg*avg)
+	print(("[mapgen_rivers] Mapgen statistics:\n- Mapgen calls: %4d\n- Mean time: %5.3f s\n- Standard deviation: %5.3f s"):format(ngen, avg, std))
+end)

load.lua

@@ -0,0 +1,100 @@
+local worldpath = mapgen_rivers.world_data_path
+
+function mapgen_rivers.load_map(filename, bytes, signed, size, converter)
+	local file = io.open(worldpath .. filename, 'rb')
+	local data = file:read('*all')
+	if #data < bytes*size then
+		data = minetest.decompress(data)
+	end
+	local sbyte = string.byte
+
+	local map = {}
+
+	for i=1, size do
+		local i0 = (i-1)*bytes+1
+		local elements = {data:byte(i0, i1)}
+		local n = sbyte(data, i0)
+		if signed and n >= 128 then
+			n = n - 256
+		end
+
+		for j=1, bytes-1 do
+			n = n*256 + sbyte(data, i0+j)
+		end
+
+		map[i] = n
+	end
+	file:close()
+
+	if converter then
+		for i=1, size do
+			map[i] = converter(map[i])
+		end
+	end
+
+	return map
+end
+
+local sbyte = string.byte
+
+local loader_mt = {
+	__index = function(loader, i)
+		local file = loader.file
+		local bytes = loader.bytes
+		file:seek('set', (i-1)*bytes)
+		local strnum = file:read(bytes)
+
+		local n = sbyte(strnum, 1)
+		if loader.signed and n >= 128 then
+			n = n - 256
+		end
+
+		for j=2, bytes do
+			n = n*256 + sbyte(strnum, j)
+		end
+
+		if loader.conv then
+			n = loader.conv(n)
+		end
+		loader[i] = n
+		return n
+	end,
+}
+
+function mapgen_rivers.interactive_loader(filename, bytes, signed, size, converter)
+	local file = io.open(worldpath .. filename, 'rb')
+	if file then
+		minetest.register_on_shutdown(function()
+			file:close()
+		end)
+		converter = converter or false
+		return setmetatable({file=file, bytes=bytes, signed=signed, size=size, conv=converter}, loader_mt)
+	end
+end
+
+function mapgen_rivers.write_map(filename, data, bytes)
+	local size = #data
+	local file = io.open(worldpath .. filename, 'wb')
+	local mfloor = math.floor
+	local schar = string.char
+	local upack = unpack
+
+	local bytelist = {}
+	for j=1, bytes do
+		bytelist[j] = 0
+	end
+
+	for i=1, size do
+		local n = mfloor(data[i])
+		data[i] = n
+		for j=bytes, 2, -1 do
+			bytelist[j] = n % 256
+			n = mfloor(n / 256)
+		end
+		bytelist[1] = n % 256
+
+		file:write(schar(upack(bytelist)))
+	end
+
+	file:close()
+end

mapgen.lua

@@ -1,278 +0,0 @@
--- Mapgen loop and mapgen-related things
-
-if minetest.get_mapgen_setting("mg_name") ~= "singlenode" then
-	minetest.set_mapgen_setting("mg_name", "singlenode", true)
-	minetest.log("warning", "[mapgen_rivers] Mapgen set to singlenode")
-end
-
-local sea_level = mapgen_rivers.settings.sea_level
-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
-
-if use_biomegen_mod then
-	biomegen.set_elevation_chill(elevation_chill)
-end
-
--- Linear interpolation
-local function interp(v00, v01, v11, v10, xf, zf)
-	local v0 = v01*xf + v00*(1-xf)
-	local v1 = v11*xf + v10*(1-xf)
-	return v1*zf + v0*(1-zf)
-end
-
--- Localize for performance
-local floor, min = math.floor, math.min
-
-local data = {}
-
-local noise_x_obj, noise_z_obj, noise_distort_obj, noise_heat_obj, noise_heat_blend_obj
-local noise_x_map = {}
-local noise_z_map = {}
-local noise_distort_map = {}
-local noise_heat_map = {}
-local noise_heat_blend_map = {}
-local mapsize
-local init = false
-
-local sumtime = 0
-local sumtime2 = 0
-local ngen = 0
-
-function mapgen_rivers.make_chunk(minp, maxp, seed)
-	minetest.log("info", ("[mapgen_rivers] Generating from %s to %s"):format(minetest.pos_to_string(minp), minetest.pos_to_string(maxp)))
-
-	local chulens = {
-		x = maxp.x-minp.x+1,
-		y = maxp.y-minp.y+1,
-		z = maxp.z-minp.z+1,
-	}
-
-	if not init then
-		mapsize = {
-			x = chulens.x,
-			y = chulens.y+1,
-			z = chulens.z,
-		}
-		if use_distort then
-			noise_x_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.distort_x, mapsize)
-			noise_z_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.distort_z, mapsize)
-			noise_distort_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.distort_amplitude, chulens)
-		end
-		if use_biomes then
-			noise_heat_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.heat, chulens)
-			noise_heat_blend_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.heat_blend, chulens)
-		end
-		init = true
-	end
-
-	local t0 = os.clock()
-	local minp2d = {x=minp.x, y=minp.z}
-	if use_distort then
-		noise_x_obj:get_3d_map_flat(minp, noise_x_map)
-		noise_z_obj:get_3d_map_flat(minp, noise_z_map)
-		noise_distort_obj:get_2d_map_flat(minp2d, noise_distort_map)
-	end
-	if use_biomes then
-		noise_heat_obj:get_2d_map_flat(minp2d, noise_heat_map)
-		noise_heat_blend_obj:get_2d_map_flat(minp2d, noise_heat_blend_map)
-	end
-
-	local terrain_map, lake_map, incr, i_origin
-
-	if use_distort then
-		local xmin, xmax, zmin, zmax = minp.x, maxp.x, minp.z, maxp.z
-		local i = 0
-		local i2d = 0
-		for z=minp.z, maxp.z do
-			for y=minp.y, maxp.y+1 do
-				for x=minp.x, maxp.x do
-					i = i+1
-					i2d = i2d+1
-					local distort = noise_distort_map[i2d]
-					local xv = noise_x_map[i]*distort + x
-					if xv < xmin then xmin = xv end
-					if xv > xmax then xmax = xv end
-					noise_x_map[i] = xv
-					local zv = noise_z_map[i]*distort + z
-					if zv < zmin then zmin = zv end
-					if zv > zmax then zmax = zv end
-					noise_z_map[i] = zv
-				end
-				i2d = i2d-chulens.x
-			end
-			i2d = i2d+chulens.x
-		end
-
-		local pminp = {x=floor(xmin), z=floor(zmin)}
-		local pmaxp = {x=floor(xmax)+1, z=floor(zmax)+1}
-		incr = pmaxp.x-pminp.x+1
-		i_origin = 1 - pminp.z*incr - pminp.x
-		terrain_map, lake_map = mapgen_rivers.make_heightmaps(pminp, pmaxp)
-	else
-		terrain_map, lake_map = mapgen_rivers.make_heightmaps(minp, maxp)
-	end
-
-	-- Check that there is at least one position that reaches min y
-	if minp.y > sea_level then
-		local y0 = minp.y
-		local is_empty = true
-		for i=1, #terrain_map do
-			if terrain_map[i] >= y0 or lake_map[i] >= y0 then
-				is_empty = false
-				break
-			end
-		end
-
-		-- If not, skip chunk
-		if is_empty then
-			if use_biomegen_mod and biomegen.skip_chunk then
-				biomegen.skip_chunk(minp, maxp)
-			end
-
-			local t = os.clock() - t0
-			ngen = ngen + 1
-			sumtime = sumtime + t
-			sumtime2 = sumtime2 + t*t
-
-			minetest.log("verbose", "[mapgen_rivers] Skipping empty chunk (fully above ground level)")
-			minetest.log("verbose", ("[mapgen_rivers] Done in %5.3f s"):format(t))
-			return
-		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 vm, emin, emax = minetest.get_mapgen_object("voxelmanip")
-	vm:get_data(data)
-
-	local a = VoxelArea:new({MinEdge = emin, MaxEdge = emax})
-	local ystride = a.ystride -- Tip : the ystride of a VoxelArea is the number to add to the array index to get the index of the position above. It's faster because it avoids to completely recalculate the index.
-
-	local nid = mapsize.x*(mapsize.y-1) + 1
-	local incrY = -mapsize.x
-	local incrX = 1 - mapsize.y*incrY
-	local incrZ = mapsize.x*mapsize.y - mapsize.x*incrX - mapsize.x*mapsize.y*incrY
-
-	local i2d = 1
-
-	for z = minp.z, maxp.z do
-		for x = minp.x, maxp.x do
-			local ivm = a:index(x, maxp.y+1, z)
-			local ground_above = false
-			local temperature
-			if use_biomes then
-				temperature = noise_heat_map[i2d]+noise_heat_blend_map[i2d]
-			end
-			local terrain, lake
-			if not use_distort then
-				terrain = terrain_map[i2d]
-				lake = lake_map[i2d]
-			end
-
-			for y = maxp.y+1, minp.y, -1 do
-				if use_distort then
-					local xn = noise_x_map[nid]
-					local zn = noise_z_map[nid]
-					local x0 = floor(xn)
-					local z0 = floor(zn)
-
-					local i0 = i_origin + z0*incr + x0
-					local i1 = i0+1
-					local i2 = i1+incr
-					local i3 = i2-1
-
-					terrain = interp(terrain_map[i0], terrain_map[i1], terrain_map[i2], terrain_map[i3], xn-x0, zn-z0)
-					lake = min(lake_map[i0], lake_map[i1], lake_map[i2], lake_map[i3])
-				end
-
-				if y <= maxp.y then
-
-					local is_lake = lake > terrain
-					if y <= terrain then
-						if not use_biomes or y <= terrain-1 or ground_above then
-							data[ivm] = c_stone
-						elseif is_lake or y < sea_level then
-							data[ivm] = c_sand
-						else
-							local temperature_y = temperature - y*elevation_chill
-							if temperature_y >= 15 then
-								data[ivm] = c_lawn
-							elseif temperature_y >= 0 then
-								data[ivm] = c_dirtsnow
-							else
-								data[ivm] = c_snow
-							end
-						end
-					elseif y <= lake and lake > sea_level then
-						if not use_biomes or temperature - y*elevation_chill >= 0 then
-							data[ivm] = c_rwater
-						else
-							data[ivm] = c_ice
-						end
-					elseif y <= sea_level then
-						data[ivm] = c_water
-					end
-				end
-
-				ground_above = y <= terrain
-
-				ivm = ivm - ystride
-				if use_distort then
-					nid = nid + incrY
-				end
-			end
-
-			if use_distort then
-				nid = nid + incrX
-			end
-			i2d = i2d + 1
-		end
-
-		if use_distort then
-			nid = nid + incrZ
-		end
-	end
-
-	if use_biomegen_mod then
-		biomegen.generate_all(data, a, vm, minp, maxp, seed)
-	else
-		vm:set_data(data)
-		minetest.generate_ores(vm, minp, maxp)
-	end
-
-	vm:set_lighting({day = 0, night = 0})
-	vm:calc_lighting()
-	vm:update_liquids()
-	vm:write_to_map()
-
-	local t = os.clock()-t0
-	ngen = ngen + 1
-	sumtime = sumtime + t
-	sumtime2 = sumtime2 + t*t
-	minetest.log("verbose", ("[mapgen_rivers] Done in %5.3f s"):format(t))
-end
-
--- Enforce first position in mapgen callbacks
-table.insert(minetest.registered_on_generateds, 1, mapgen_rivers.make_chunk)
---minetest.register_on_generated(mapgen_rivers.make_chunk)
-
-minetest.register_on_shutdown(function()
-	local avg = sumtime / ngen
-	local std = math.sqrt(sumtime2/ngen - avg*avg)
-	minetest.log("action", ("[mapgen_rivers] Mapgen statistics:\n- Mapgen calls: %4d\n- Mean time: %5.3f s\n- Standard deviation: %5.3f s"):format(ngen, avg, std))
-end)

mod.conf

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

noises.lua

@@ -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
+			print("[mapgen_rivers] Noise " .. name .. ": 'octaves' reduced to " .. omax)
+			np.octaves = omax
+		end
+	end
+end

polygons.lua

@@ -1,37 +1,105 @@
--- Fetch polygons from a given areas, and compute their properties
+local modpath = mapgen_rivers.modpath
--- and find to which polygon every point belongs
+local mod_data_path = modpath .. 'river_data/'
+if not io.open(mod_data_path .. 'size', 'r') then
+	mod_data_path = modpath .. 'demo_data/'
+end
 
-local blocksize = mapgen_rivers.settings.blocksize
+local world_data_path = mapgen_rivers.world_data_path
-local X = math.floor(mapgen_rivers.settings.map_x_size / blocksize)
+minetest.mkdir(world_data_path)
-local Z = math.floor(mapgen_rivers.settings.map_z_size / blocksize)
+
+dofile(modpath .. 'load.lua')
+
+mapgen_rivers.grid = {}
+
+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)
+end
+
+local load_all = mapgen_rivers.settings.load_all
+
+-- Try to read file 'size'
+local sfile = io.open(world_data_path..'size', 'r')
+local first_mapgen = true
+if sfile then
+	X, Z = tonumber(sfile:read('*l')), tonumber(sfile:read('*l'))
+	sfile:close()
+	first_mapgen = false
+end
+
+if first_mapgen then
+	-- Generate a map!!
+	local pregenerate = dofile(mapgen_rivers.modpath .. '/pregenerate.lua')
+	minetest.register_on_mods_loaded(function()
+		print('[mapgen_rivers] Generating grid')
+		pregenerate(load_all)
+
+		if load_all then
+			local offset_x = mapgen_rivers.grid.offset_x
+			local offset_y = mapgen_rivers.grid.offset_y
+			for i=1, X*Z do
+				offset_x[i] = offset_converter(offset_x[i])
+				offset_y[i] = offset_converter(offset_y[i])
+			end
+		end
+	end)
+end
+
+-- if data not already loaded
+if not (first_mapgen and load_all) then
+	local load_map
+	if load_all then
+		load_map = mapgen_rivers.load_map
+	else
+		load_map = mapgen_rivers.interactive_loader
+	end
+
+	minetest.register_on_mods_loaded(function()
+		if load_all then
+			print('[mapgen_rivers] Loading full grid')
+		else
+			print('[mapgen_rivers] Loading grid as interactive loaders')
+		end
+		local grid = mapgen_rivers.grid
+
+		grid.dem = load_map('dem', 2, true, X*Z)
+		grid.lakes = load_map('lakes', 2, true, X*Z)
+		grid.dirs = load_map('dirs', 1, false, X*Z)
+		grid.rivers = load_map('rivers', 4, false, X*Z)
+
+		grid.offset_x = load_map('offset_x', 1, true, X*Z, offset_converter)
+		grid.offset_y = load_map('offset_y', 1, true, X*Z, offset_converter)
+	end)
+end
+
+mapgen_rivers.grid.size = {x=X, y=Z}
 
 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
+
 local map_offset = {x=0, z=0}
-mapgen_rivers.register_on_grid_loaded(function(grid)
+if mapgen_rivers.settings.center then
-	X = grid.size.x
-	Z = grid.size.y
-	if mapgen_rivers.settings.center then
 	map_offset.x = blocksize*X/2
 	map_offset.z = blocksize*Z/2
-	end
+end
-end)
-
--- Localize for performance
-local floor, ceil, min, max, abs = math.floor, math.ceil, math.min, math.max, math.abs
 
 local min_catchment = mapgen_rivers.settings.min_catchment / (blocksize*blocksize)
 local wpower = mapgen_rivers.settings.river_widening_power
 local wfactor = 1/(2*blocksize * min_catchment^wpower)
 local function river_width(flow)
-	flow = abs(flow)
+	flow = math.abs(flow)
 	if flow < min_catchment then
 		return 0
 	end
 
-	return min(wfactor * flow ^ wpower, 1)
+	return math.min(wfactor * flow ^ wpower, 1)
 end
 
 local noise_heat -- Need a large-scale noise here so no heat blend
@@ -48,7 +116,7 @@ local init = false
 
 -- On map generation, determine into which polygon every point (in 2D) will fall.
 -- Also store polygon-specific data
-function mapgen_rivers.make_polygons(minp, maxp)
+local function make_polygons(minp, maxp)
 
 	local grid = mapgen_rivers.grid
 	local dem = grid.dem
@@ -70,8 +138,8 @@ function mapgen_rivers.make_polygons(minp, maxp)
 
 	local polygons = {}
 	-- Determine the minimum and maximum coordinates of the polygons that could be on the chunk, knowing that they have an average size of 'blocksize' and a maximal offset of 0.5 blocksize.
-	local xpmin, xpmax = max(floor((minp.x+map_offset.x)/blocksize - 0.5), 0), min(ceil((maxp.x+map_offset.x)/blocksize + 0.5), X-2)
+	local xpmin, xpmax = math.max(math.floor((minp.x+map_offset.x)/blocksize - 0.5), 0), math.min(math.ceil((maxp.x+map_offset.x)/blocksize + 0.5), X-2)
-	local zpmin, zpmax = max(floor((minp.z+map_offset.z)/blocksize - 0.5), 0), min(ceil((maxp.z+map_offset.z)/blocksize + 0.5), Z-2)
+	local zpmin, zpmax = math.max(math.floor((minp.z+map_offset.z)/blocksize - 0.5), 0), math.min(math.ceil((maxp.z+map_offset.z)/blocksize + 0.5), Z-2)
 
 	-- Iterate over the polygons
 	for xp = xpmin, xpmax do
@@ -97,8 +165,8 @@ function mapgen_rivers.make_polygons(minp, maxp)
 
 			local bounds = {} -- Will be a list of the intercepts of polygon edges for every Z position (scanline algorithm)
 			-- Calculate the min and max Z positions 
-			local zmin = max(floor(min(unpack(poly_z)))+1, minp.z)
+			local zmin = math.max(math.floor(math.min(unpack(poly_z)))+1, minp.z)
-			local zmax = min(floor(max(unpack(poly_z))), maxp.z)
+			local zmax = math.min(math.floor(math.max(unpack(poly_z))), maxp.z)
 			-- And initialize the arrays
 			for z=zmin, zmax do
 				bounds[z] = {}
@@ -108,14 +176,14 @@ function mapgen_rivers.make_polygons(minp, maxp)
 			for i2=1, 4 do -- Loop on 4 edges
 				local z1, z2 = poly_z[i1], poly_z[i2]
 				-- Calculate the integer Z positions over which this edge spans
-				local lzmin = floor(min(z1, z2))+1
+				local lzmin = math.floor(math.min(z1, z2))+1
-				local lzmax = floor(max(z1, z2))
+				local lzmax = math.floor(math.max(z1, z2))
 				if lzmin <= lzmax then -- If there is at least one position in it
 					local x1, x2 = poly_x[i1], poly_x[i2]
 					-- Calculate coefficient of the equation defining the edge: X=aZ+b
 					local a = (x1-x2) / (z1-z2)
 					local b = (x1 - a*z1)
-					for z=max(lzmin, minp.z), min(lzmax, maxp.z) do
+					for z=math.max(lzmin, minp.z), math.min(lzmax, maxp.z) do
 						-- For every Z position involved, add the intercepted X position in the table
 						table.insert(bounds[z], a*z+b)
 					end
@@ -126,11 +194,11 @@ function mapgen_rivers.make_polygons(minp, maxp)
 				-- Now sort the bounds list
 				local zlist = bounds[z]
 				table.sort(zlist)
-				local c = floor(#zlist/2)
+				local c = math.floor(#zlist/2)
 				for l=1, c do
 					-- Take pairs of X coordinates: all positions between them belong to the polygon.
-					local xmin = max(floor(zlist[l*2-1])+1, minp.x)
+					local xmin = math.max(math.floor(zlist[l*2-1])+1, minp.x)
-					local xmax = min(floor(zlist[l*2]), maxp.x)
+					local xmax = math.min(math.floor(zlist[l*2]), maxp.x)
 					local i = (z-minp.z) * chulens + (xmin-minp.x) + 1
 					for x=xmin, xmax do
 						-- Fill the map at these places
@@ -152,28 +220,28 @@ function mapgen_rivers.make_polygons(minp, maxp)
 			local riverD = river_width(rivers[iD])
 			if glaciers then -- Widen the river
 				if get_temperature(poly_x[1], poly_dem[1], poly_z[1]) < 0 then
-					riverA = min(riverA*glacier_factor, 1)
+					riverA = math.min(riverA*glacier_factor, 1)
 				end
 				if get_temperature(poly_x[2], poly_dem[2], poly_z[2]) < 0 then
-					riverB = min(riverB*glacier_factor, 1)
+					riverB = math.min(riverB*glacier_factor, 1)
 				end
 				if get_temperature(poly_x[3], poly_dem[3], poly_z[3]) < 0 then
-					riverC = min(riverC*glacier_factor, 1)
+					riverC = math.min(riverC*glacier_factor, 1)
 				end
 				if get_temperature(poly_x[4], poly_dem[4], poly_z[4]) < 0 then
-					riverD = min(riverD*glacier_factor, 1)
+					riverD = math.min(riverD*glacier_factor, 1)
 				end
 			end
 
-			polygon.river_corners = {riverA, riverB, riverC, riverD}
+			polygon.river_corners = {riverA, 1-riverB, 2-riverC, 1-riverD}
 
 			-- Flow directions
 			local dirA, dirB, dirC, dirD = dirs[iA], dirs[iB], dirs[iC], dirs[iD]
 			-- Determine the river flux on the edges, by testing dirs values
 			local river_west = (dirA==1 and riverA or 0) + (dirD==3 and riverD or 0)
 			local river_north = (dirA==2 and riverA or 0) + (dirB==4 and riverB or 0)
-			local river_east = (dirB==1 and riverB or 0) + (dirC==3 and riverC or 0)
+			local river_east = 1 - (dirB==1 and riverB or 0) - (dirC==3 and riverC or 0)
-			local river_south = (dirD==2 and riverD or 0) + (dirC==4 and riverC or 0)
+			local river_south = 1 - (dirD==2 and riverD or 0) - (dirC==4 and riverC or 0)
 
 			polygon.rivers = {river_west, river_north, river_east, river_south}
 		end
@@ -181,3 +249,5 @@ function mapgen_rivers.make_polygons(minp, maxp)
 
 	return polygons
 end
+
+return make_polygons

pregenerate.lua

@@ -1,8 +1,3 @@
--- Generate the grid using terrainlib_lua
--- Only called on first mapgen, if there is no grid yet
-
--- Constants
-
 local EvolutionModel = dofile(mapgen_rivers.modpath .. '/terrainlib_lua/erosion.lua')
 local twist = dofile(mapgen_rivers.modpath .. '/terrainlib_lua/twist.lua')
 
@@ -10,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
 
@@ -19,70 +13,27 @@ 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 function pregenerate(keep_loaded)
-local margin_width = mapgen_rivers.settings.margin_width / blocksize
+	local grid = mapgen_rivers.grid
-local margin_elev = mapgen_rivers.settings.margin_elev
+	local size = grid.size
 
-local X = math.floor(mapgen_rivers.settings.map_x_size / blocksize)
+	local seed = tonumber(minetest.get_mapgen_setting("seed"))
-local Y = math.floor(mapgen_rivers.settings.map_z_size / blocksize)
+	np_base.seed = (np_base.seed or 0) + seed
 
-local function margin(dem, width, elev)
+	local nobj_base = PerlinNoiseMap(np_base, {x=size.x, y=1, z=size.y})
-	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
 
--- Generate grid
+	local dem = nobj_base:get_3d_map_flat({x=0, y=0, z=0})
+	dem.X = size.x
+	dem.Y = size.y
 
-minetest.log("action", '[mapgen_rivers] Generating grid, this may take a while...')
+	local model = EvolutionModel(evol_params)
+	model.dem = dem
+	local ref_dem = model:define_isostasy(dem)
 
-if X*Y > 4e6 then
+	local tectonic_step = tectonic_speed * time_step
-	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.")
+	collectgarbage()
-end
+	for i=1, niter do
-
+		print("[mapgen_rivers] Iteration " .. i .. " of " .. niter)
-local seed = tonumber(minetest.get_mapgen_setting("seed"):sub(-10))
-np_base.seed = (np_base.seed or 0) + seed
-
-local nobj_base = PerlinNoiseMap(np_base, {x=X, y=1, z=Y})
-
-local dem = nobj_base:get_3d_map_flat({x=0, y=0, z=0})
-dem.X = X
-dem.Y = 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)
-
-local tectonic_step = tectonic_speed * time_step
-collectgarbage()
-for i=1, niter do
-	minetest.log("info", "[mapgen_rivers] Iteration " .. i .. " of " .. niter)
 
 		model:diffuse(time_step)
 		model:flow()
@@ -90,61 +41,41 @@ for i=1, niter do
 		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
 
 		collectgarbage()
-end
+	end
-model:flow()
+	model:flow()
 
-local mfloor = math.floor
+	local mfloor = math.floor
-local mmin, mmax = math.min, math.max
+	local mmin, mmax = math.min, math.max
-local unpk, schar = unpack, string.char
+	local offset_x, offset_y = twist(model.dirs, model.rivers, 5)
-local offset_x, offset_y = twist(model.dirs, model.rivers, 5)
+	for i=1, size.x*size.y do
-for i=1, X*Y do
 		offset_x[i] = mmin(mmax(offset_x[i]*256, -128), 127)
 		offset_y[i] = mmin(mmax(offset_y[i]*256, -128), 127)
+	end
+
+	mapgen_rivers.write_map('dem', model.dem, 2)
+	mapgen_rivers.write_map('lakes', model.lakes, 2)
+	mapgen_rivers.write_map('dirs', model.dirs, 1)
+	mapgen_rivers.write_map('rivers', model.rivers, 4)
+	mapgen_rivers.write_map('offset_x', offset_x, 1)
+	mapgen_rivers.write_map('offset_y', offset_y, 1)
+	local sfile = io.open(mapgen_rivers.world_data_path .. 'size', "w")
+	sfile:write(size.x..'\n'..size.y)
+	sfile:close()
+
+	if keep_loaded then
+		grid.dem = model.dem
+		grid.lakes = model.lakes
+		grid.dirs = model.dirs
+		grid.rivers = model.rivers
+		grid.offset_x = offset_x
+		grid.offset_y = offset_y
+	end
+	collectgarbage()
 end
 
--- Write data
+return pregenerate
-
-local datapath = mapgen_rivers.world_data_path
-minetest.mkdir(datapath)
-
-local sfile = io.open(datapath .. 'size', "w")
-sfile:write(X..'\n'..Y)
-sfile:close()
-
-local function write_file(filename, data, bytes)
-	local file = io.open(datapath .. filename, 'wb')
-
-	local bytelist = {}
-	for j=1, bytes do
-		bytelist[j] = 0
-	end
-
-	for i=1, #data do
-		local n = mfloor(data[i])
-		data[i] = n
-		for j=bytes, 2, -1 do
-			bytelist[j] = n % 256
-			n = mfloor(n / 256)
-		end
-		bytelist[1] = n % 256
-
-		file:write(schar(unpk(bytelist)))
-	end
-
-	file:close()
-end
-
-write_file('dem', model.dem, 2)
-write_file('lakes', model.lakes, 2)
-write_file('dirs', model.dirs, 1)
-write_file('rivers', model.rivers, 4)
-write_file('offset_x', offset_x, 1)
-write_file('offset_y', offset_y, 1)

settings.lua

@@ -1,9 +1,7 @@
--- Read global and per-world settings
-
 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"
 
 local previous_version_mt = mtsettings:get("mapgen_rivers_version") or "0.0"
 local previous_version_mgr = mgrsettings:get("version") or "0.0"
@@ -21,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
@@ -61,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
@@ -71,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
@@ -83,47 +60,33 @@ end
 local def_setting = mapgen_rivers.define_setting
 
 mapgen_rivers.settings = {
-	center = def_setting('center', 'bool'),
+	center = def_setting('center', 'bool', true),
-	blocksize = def_setting('blocksize', 'number'),
+	blocksize = def_setting('blocksize', 'number', 15),
 	sea_level = tonumber(minetest.get_mapgen_setting('water_level')),
-	min_catchment = def_setting('min_catchment', 'number'),
+	min_catchment = def_setting('min_catchment', 'number', 3600),
-	river_widening_power = def_setting('river_widening_power', 'number'),
+	river_widening_power = def_setting('river_widening_power', 'number', 0.5),
-	riverbed_slope = def_setting('riverbed_slope', 'number'),
+	riverbed_slope = def_setting('riverbed_slope', 'number', 0.4),
-	distort = def_setting('distort', 'bool'),
+	distort = def_setting('distort', 'bool', true),
-	biomes = def_setting('biomes', 'bool'),
+	biomes = def_setting('biomes', 'bool', true),
-	glaciers = def_setting('glaciers', 'bool'),
+	glaciers = def_setting('glaciers', 'bool', false),
-	glacier_factor = def_setting('glacier_factor', 'number'),
+	glacier_factor = def_setting('glacier_factor', 'number', 8),
-	elevation_chill = def_setting('elevation_chill', 'number'),
+	elevation_chill = def_setting('elevation_chill', 'number', 0.25),
 
-	map_x_size = def_setting('map_x_size', 'number'),
+	grid_x_size = def_setting('grid_x_size', 'number', 1000),
-	map_z_size = def_setting('map_z_size', 'number'),
+	grid_z_size = def_setting('grid_z_size', 'number', 1000),
-	margin = def_setting('margin', 'bool'),
-	margin_width = def_setting('margin_width', 'number'),
-	margin_elev = def_setting('margin_elev', 'number'),
 	evol_params = {
-		K = def_setting('river_erosion_coef', 'number'),
+		K = def_setting('river_erosion_coef', 'number', 0.5),
-		m = def_setting('river_erosion_power', 'number'),
+		m = def_setting('river_erosion_power', 'number', 0.4),
-		d = def_setting('diffusive_erosion', 'number'),
+		d = def_setting('diffusive_erosion', 'number', 0.5),
-		compensation_radius = def_setting('compensation_radius', 'number'),
+		compensation_radius = def_setting('compensation_radius', 'number', 50),
 	},
-	tectonic_speed = def_setting('tectonic_speed', 'number'),
+	tectonic_speed = def_setting('tectonic_speed', 'number', 70),
-	evol_time = def_setting('evol_time', 'number'),
+	evol_time = def_setting('evol_time', 'number', 10),
-	evol_time_step = def_setting('evol_time_step', 'number'),
+	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

settings_default.json

@@ -1,70 +0,0 @@
-{
-	"version": "1.0.2-dev1",
-
-	"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,
-	"margin": true,
-	"margin_width": 2000,
-	"margin_elev": -200,
-	"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,
-
-	"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"
-	}
-}

settingtypes.txt

@@ -3,35 +3,19 @@
 #    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.
+#    Represents horizontal map scale. Every cell of the grid will be upscaled to
-#    A lower value will result in more detailed terrain and finer computation
+#    a square of this size.
-#    of rivers, but will be slower to generate and use more resources.
+#    For example if the grid size is 1000x1000 and block size is 12,
-#
+#    the actual size of the map will be 12000.
-#    WARNING: Excessively low values may cause crashes at pre-generation, due to
-#    memory issues
 mapgen_rivers_blocksize (Block size) float 15.0 2.0 100.0
 
-#    X size of the map being generated
+#    X size of the grid being generated
-#
+#    Actual size of the map is grid_x_size * blocksize
-#    X size of the river grid will be this/blocksize
+mapgen_rivers_grid_x_size (Grid X size) int 1000 50 5000
-#    When increasing, it is recommended to increase blocksize too
-mapgen_rivers_map_x_size (Map X size) int 15000 500 66000
 
-#    Z size of the map being generated
+#    Z size of the grid being generated
-#
+#    Actual size of the map is grid_z_size * blocksize
-#    Z size of the river grid will be this/blocksize
+mapgen_rivers_grid_z_size (Grid Z size) int 1000 50 5000
-#    When increasing, it is recommended to increase blocksize too
-mapgen_rivers_map_z_size (Map Z size) int 15000 500 66000
-
-#    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
@@ -39,7 +23,7 @@ mapgen_rivers_min_catchment (Minimal catchment area) float 3600.0 100.0 1000000.
 
 #    Coefficient describing how rivers widen when merging.
 #    Riwer width is a power law W = a*D^p. D is river flow and p is this parameter.
-#    Higher value means that a river will grow more when receiving a tributary.
+#    Higher value means a river needs to receive more tributaries to grow in width.
 #    Note that a river can never exceed 2*blocksize.
 mapgen_rivers_river_widening_power (River widening power) float 0.5 0.0 1.0
 
@@ -70,6 +54,11 @@ mapgen_rivers_glacier_widening_factor (Glacier widening factor) float 8.0 1.0 20
 #    This results in mountains being more covered by snow.
 mapgen_rivers_elevation_chill (Elevation chill) float 0.25 0.0 5.0
 
+#    If enabled, loads all grid data in memory at init time.
+#    If disabled, data will be loaded on request and cached in memory.
+#    It's recommended to disable it for very large maps (> 2000 grid nodes or so)
+mapgen_rivers_load_all (Load all data in memory) bool false
+
 [Landscape evolution parameters]
 
 #    Modelled landscape evolution time, in arbitrary units

spawn.lua

@@ -1,119 +0,0 @@
-local np_distort_x = mapgen_rivers.noise_params.distort_x
-local np_distort_z = mapgen_rivers.noise_params.distort_z
-local np_distort_amplitude = mapgen_rivers.noise_params.distort_amplitude
-
-local nobj_distort_x, nobj_distort_z, nobj_distort_amplitude
-
-local sea_level = mapgen_rivers.settings.sea_level
-local distort = mapgen_rivers.settings.distort
-
--- Linear interpolation
-local function interp(v00, v01, v11, v10, xf, zf)
-	local v0 = v01*xf + v00*(1-xf)
-	local v1 = v11*xf + v10*(1-xf)
-	return v1*zf + v0*(1-zf)
-end
-
-local function estimate_spawn_level(pos, use_distort)
-	local x, z = pos.x, pos.z
-	if distort and use_distort then
-		nobj_distort_x = nobj_distort_x or minetest.get_perlin(np_distort_x)
-		nobj_distort_z = nobj_distort_z or minetest.get_perlin(np_distort_z)
-		nobj_distort_amplitude = nobj_distort_amplitude or minetest.get_perlin(np_distort_amplitude)
-
-		local amplitude = nobj_distort_amplitude:get_2d({x=pos.x, y=pos.z})
-		x = x + nobj_distort_x:get_3d(pos)*amplitude
-		z = z + nobj_distort_z:get_3d(pos)*amplitude
-	end
-
-	local terrain, lakes = mapgen_rivers.make_heightmaps(
-		{x=math.floor(x),   z=math.floor(z)  },
-		{x=math.floor(x)+1, z=math.floor(z)+1}
-	)
-
-	local ex, ez = x % 1, z % 1
-	--local h = terrain[1]*(1-ex)*(1-ez) + terrain[2]*ex*(1-ez) + terrain[3]*(1-ex)*ez + terrain[4]*ex*ez
-	local h = interp(terrain[1], terrain[2], terrain[4], terrain[3], ex, ez)
-	local lake = math.min(lakes[1], lakes[2], lakes[3], lakes[4])
-
-	if h < lake or h <= sea_level then
-		return false, h
-	end
-
-	return true, h
-end
-
-local function get_spawn_level(x, z)
-	local pos = {x=x, z=z}
-	local suitable, y = estimate_spawn_level(pos, false)
-	if not suitable then
-		return
-	end
-
-	if not distort then
-		return math.floor(y) + 1
-	end
-
-	local low_bound = -math.huge
-	local high_bound = math.huge
-	local suitable_high = false
-
-	repeat
-		pos.y = math.max(math.min(math.floor(y+0.5), high_bound-1), low_bound+1)
-		suitable, y = estimate_spawn_level(pos, true)
-		if y > pos.y then
-			low_bound = pos.y
-		else
-			high_bound = pos.y
-			suitable_high = suitable
-		end
-	until high_bound - low_bound <= 1
-
-	if not suitable_high then
-		return
-	end
-
-	return high_bound + 1
-end
-minetest.get_spawn_level = get_spawn_level
-
-local rmax = 2000
-local function find_spawn_point(seed)
-	local level = get_spawn_level(0, 0)
-	if level then
-		return {x=0, y=level, z=0}
-	end
-
-	local pr = PcgRandom(seed or os.time())
-	local incr = 16
-	local r0 = 0
-
-	while r0 < rmax do
-		local r1 = r0 + incr
-		local r = pr:next(r0*r0+1, r1*r1) ^ 0.5
-		local a = pr:next() / 2147483648 * math.pi
-
-		local x = math.floor(math.cos(a) * r + 0.5)
-		local z = math.floor(math.sin(a) * r + 0.5)
-
-		level = get_spawn_level(x, z)
-		if level then
-			return {x=x, y=level, z=z}
-		end
-		r0 = r1
-	end
-end
-
-local function spawn_player(player)
-	if minetest.settings:get("static_spawnpoint") then
-		return
-	end
-
-	local spawn_point = find_spawn_point()
-	if spawn_point then
-		player:set_pos(spawn_point)
-	end
-end
-
-minetest.register_on_newplayer(spawn_player)
-minetest.register_on_respawnplayer(spawn_player)

terrainlib_lua/erosion.lua

@@ -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
+	--local tinsert = table.insert
-	-- 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 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 d == 0 then -- If no flow direction, we reach the border of the map: set elevation to the latest node's elev and abort.
 				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 = ""
@@ -134,12 +121,11 @@ local rivermapper = dofile(modpath .. "rivermapper.lua")
 local gaussian = dofile(modpath .. "gaussian.lua")
 
 local function flow(model)
-	model.dirs, model.lakes = rivermapper.flow_routing(model.dem, model.dirs, model.lakes)
+	model.dirs, model.lakes = rivermapper.flow_routing(model.dem, model.dirs, model.lakes, 'semirandom')
 	model.rivers = rivermapper.accumulate(model.dirs, model.rivers)
 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
 

terrainlib_lua/gaussian.lua

@@ -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)

terrainlib_lua/rivermapper.lua

@@ -9,19 +9,47 @@
 --
 -- 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)
+	local sum = 0
+	for i=1, #plist do
+		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
+		r = r - p
+	end
+	return 0
+end
+
+local flow_methods = {
+	semirandom = flow_local_semirandom,
+}
+
+local function flow_routing(dem, dirs, lakes, method)
+	method = method or 'semirandom'
+	local flow_local = flow_methods[method] or 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) -- 'dirs' and 'lakes' are optional tables to reuse for memory optimization, they may contain any data.
 	dirs = dirs or {}
 	lakes = lakes or {}
 
 	-- Localize for performance
+	--local tinsert = table.insert
 	local tremove = table.remove
 	local mmax = math.max
-	local mrand = math.random
 
 	local X, Y = dem.X, dem.Y
 	dirs.X = X
@@ -34,41 +62,20 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 		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]
-			-- Determine how water should flow at 1 node scale.
+			local plist = {
-			-- 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, with probability depending on height difference.
+				y<Y and mmax(zi-dem[i+X], 0) or 0, -- Southward
-			-- This makes rivers better follow the curvature of the topography at large scale, and be less biased by pure N/E/S/W directions.
+				x<X and mmax(zi-dem[i+1], 0) or 0, -- Eastward
-			local pSouth = y<Y and mmax(zi-dem[i+X], 0) or 0
+				y>1 and mmax(zi-dem[i-X], 0) or 0, -- Northward
-			local pEast = x<X and mmax(zi-dem[i+1], 0) or 0
+				x>1 and mmax(zi-dem[i-1], 0) or 0, -- Westward
-			local pNorth = y>1 and mmax(zi-dem[i-X], 0) or 0
+			}
-			local pWest = x>1 and mmax(zi-dem[i-1], 0) or 0
 
-			local d = 0
+			local d = flow_local(plist)
-			local sum = pSouth + pEast + pNorth + pWest
-			local r = mrand() * sum
-			if sum > 0 then
-				if r < pSouth then
-					d = 1
-				elseif r-pSouth < pEast then
-					d = 2
-				elseif r-pSouth-pEast < pNorth then
-					d = 3
-				else
-					d = 4
-				end
-			end
-
-			-- '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
@@ -83,180 +90,91 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 		end
 	end
 
-	--------------------------------------
+	-- Compute basins and links
-	-- 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)
 	local nbasins = #singular
 	local basin_id = {}
 	local links = {}
 	local basin_links
 
+	local function add_link(i1, i2, b1, isY)
+		local b2
+		if i2 == 0 then
+			b2 = 0
+		else
+			b2 = basin_id[i2]
+			if b2 == 0 then
+				return
+			end
+		end
+		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
+				l2.elev = elev
+				l2.i = mmax(i1,i2)
+				l2.is_y = isY
+				l2[1] = b2
+				l2[2] = b1
+			end
+		end
+	end
+
 	for i=1, X*Y do
 		basin_id[i] = 0
 	end
-
+	--for ib, s in ipairs(singular) do
-	local cur = nbasins
+	for ib=1, nbasins do
-	local ib = 0
+		--local s = singular[ib]
-	while cur > 0 do
+		local queue = {singular[ib]}
-			local i = singular[cur]
-			cur = cur - 1
-			if dirs[i] == 0 then
 		basin_links = {}
 		links[#links+1] = basin_links
-				ib = ib + 1
+		--tinsert(links, basin_links)
-			end
+		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 East
-			-- Loop is unrolled on purpose, for performance (critical part!)
-			----------
-			-- EAST --
-			----------
-			if d >= 8 then -- River coming from the East
 				d = d - 8
-				cur = cur + 1
+				queue[#queue+1] = i+1
-				singular[cur] = i+1
+				--tinsert(queue, i+X)
-			-- If no river is coming from the East, we might be at the limit of two basins, thus we need to test adjacency.
 			elseif i%X > 0 then
-				if basin_id[i+1] ~= ib and basin_id[i+1] ~= 0 then
+				add_link(i, i+1, ib, false)
-					local b2 = basin_id[i+1]
+			else
-					local elev = mmax(dem[i], dem[i+1]) -- Elevation of the highest of the two sides of the link (or only i1 if b2 is map outside)
+				add_link(i, 0, ib, false)
-					local l2 = basin_links[b2]
-					if not l2 then
-						l2 = {b2, ib, elev=elev, i=i+1, is_y=false}
-						basin_links[b2] = l2 -- Potential non-linear complexity here
-					elseif 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
-						l2.elev = elev
-						l2.i = i+1
-						l2.is_y = false
-						l2[1] = b2
-						l2[2] = ib
-					end
-				end
-			else -- If the eastern neighbour is outside the map
-				local l2 = basin_links[0]
-				if not l2 then
-					l2 = {0, ib, elev=dem[i], i=i, is_y=false}
-					basin_links[0] = l2
-				elseif l2.elev > dem[i] then
-					l2.elev = dem[i]
-					l2.i = i
-					l2.is_y = false
-					l2[1] = 0
-					l2[2] = ib
-				end
 			end
 
-			-----------
+			if d >= 4 then -- River coming from South
-			-- SOUTH --
-			-----------
-			if d >= 4 then -- River coming from the South
 				d = d - 4
-				cur = cur + 1
+				queue[#queue+1] = i+X
-				singular[cur] = i+X
+				--tinsert(queue, i+1)
 			elseif i <= X*(Y-1) then
-				if basin_id[i+X] ~= ib and basin_id[i+X] ~= 0 then
+				add_link(i, i+X, ib, true)
-					local b2 = basin_id[i+X]
-					local elev = mmax(dem[i], dem[i+X])
-					local l2 = basin_links[b2]
-					if not l2 then
-						l2 = {b2, ib, elev=elev, i=i+X, is_y=true}
-						basin_links[b2] = l2
-					elseif l2.elev > elev then
-						l2.elev = elev
-						l2.i = i+X
-						l2.is_y = true
-						l2[1] = b2
-						l2[2] = ib
-					end
-				end
 			else
-				local l2 = basin_links[0]
+				add_link(i, 0, ib, true)
-				if not l2 then
-					l2 = {0, ib, elev=dem[i], i=i, is_y=true}
-					basin_links[0] = l2
-				elseif l2.elev > dem[i] then
-					l2.elev = dem[i]
-					l2.i = i
-					l2.is_y = true
-					l2[1] = 0
-					l2[2] = ib
-				end
 			end
 
-			----------
+			if d >= 2 then -- River coming from West
-			-- WEST --
-			----------
-			if d >= 2 then -- River coming from the West
 				d = d - 2
-				cur = cur + 1
+				queue[#queue+1] = i-1
-				singular[cur] = i-1
+				--tinsert(queue, i-X)
 			elseif i%X ~= 1 then
-				if basin_id[i-1] ~= ib and basin_id[i-1] ~= 0 then
+				add_link(i, i-1, ib, false)
-					local b2 = basin_id[i-1]
-					local elev = mmax(dem[i], dem[i-1])
-					local l2 = basin_links[b2]
-					if not l2 then
-						l2 = {b2, ib, elev=elev, i=i, is_y=false}
-						basin_links[b2] = l2
-					elseif l2.elev > elev then
-						l2.elev = elev
-						l2.i = i
-						l2.is_y = false
-						l2[1] = b2
-						l2[2] = ib
-					end
-				end
 			else
-				local l2 = basin_links[0]
+				add_link(i, 0, ib, false)
-				if not l2 then
-					l2 = {0, ib, elev=dem[i], i=i, is_y=false}
-					basin_links[0] = l2
-				elseif l2.elev > dem[i] then
-					l2.elev = dem[i]
-					l2.i = i
-					l2.is_y = false
-					l2[1] = 0
-					l2[2] = ib
-				end
 			end
 
-			-----------
+			if d >= 1 then -- River coming from North
-			-- NORTH --
+				queue[#queue+1] = i-X
-			-----------
+				--tinsert(queue, i-1)
-			if d >= 1 then -- River coming from the North
-				cur = cur + 1
-				singular[cur] = i-X
 			elseif i > X then
-				if basin_id[i-X] ~= ib and basin_id[i-X] ~= 0 then
+				add_link(i, i-X, ib, true)
-					local b2 = basin_id[i-X]
-					local elev = mmax(dem[i], dem[i-X])
-					local l2 = basin_links[b2]
-					if not l2 then
-						l2 = {b2, ib, elev=elev, i=i, is_y=true}
-						basin_links[b2] = l2
-					elseif l2.elev > elev then
-						l2.elev = elev
-						l2.i = i
-						l2.is_y = true
-						l2[1] = b2
-						l2[2] = ib
-					end
-				end
 			else
-				local l2 = basin_links[0]
+				add_link(i, 0, ib, true)
-				if not l2 then
-					l2 = {0, ib, elev=dem[i], i=i, is_y=true}
-					basin_links[0] = l2
-				elseif l2.elev > dem[i] then
-					l2.elev = dem[i]
-					l2.i = i
-					l2.is_y = true
-					l2[1] = 0
-					l2[2] = ib
 			end
 		end
 	end
@@ -268,7 +186,7 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 		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
@@ -279,51 +197,40 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 		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 = {}
-	cur = 0
+	for i, n in pairs(nlinks) do
-	local ref = singular -- Reuse table
+		if n <= 8 then
-	for i=0, nbasins do
+			lowlevel[i] = links[i]
-		if nlinks[i] <= 8 then
-			cur = cur + 1
-			lowlevel[cur] = i
-			ref[i] = cur
 		end
 	end
 
 	local basin_graph = {}
-	for i=0, nbasins do
+	for n=1, nbasins do
-		basin_graph[i] = {} -- Initialize (to ensure subtables don't go in the hash part)
+		local b1, lnk1 = next(lowlevel)
-	end
+		lowlevel[b1] = nil
-	for i=1, nbasins do
-		-- Iterate in lowlevel but its contents may change during the loop
-		local b1 = lowlevel[cur]
-		cur = cur - 1
-		local lnk1 = links[b1]
 
 		local b2
 		local lowest = math.huge
 		local lnk1 = links[b1]
-		-- Look for lowest link
+		local i = 0
 		for bn, bdata in pairs(lnk1) do
+			i = i + 1
 			if bdata.elev < lowest then
 				lowest = bdata.elev
 				b2 = bn
 			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]
-		basin_graph[bb1][bb2] = bound -- Potential non-linear complexity here
+		if not basin_graph[bb1] then
+			basin_graph[bb1] = {}
+		end
+		if not basin_graph[bb2] then
+			basin_graph[bb2] = {}
+		end
+		basin_graph[bb1][bb2] = bound
 		basin_graph[bb2][bb1] = bound
 
 		-- Merge basin b1 into b2
@@ -332,67 +239,49 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 		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
-			cur = cur + 1
+			lowlevel[b2] = lnk2
-			lowlevel[cur] = b2
-			ref[b2] = cur
 		end
 		-- Look for basin 1's neighbours, and add them to basin 2 if they have a lower pass
 		for bn, bdata in pairs(lnk1) do
 			local lnkn = links[bn]
 			lnkn[b1] = nil
 
-			if lnkn[b2] then -- If bassin bn is also linked to b2
+			if lnkn[b2] then
-				nlinks[bn] = nlinks[bn] - 1 -- Then bassin bn is losing a link because it keeps only one link toward b1/b2 after the merge
+				nlinks[bn] = nlinks[bn] - 1
 				if nlinks[bn] == 8 then
-					cur = cur + 1
+					lowlevel[bn] = lnkn
-					lowlevel[cur] = bn
-					ref[bn] = cur
 				end
-			else -- If bn was linked to b1 but not to b2
+			else
-				nlinks[b2] = nlinks[b2] + 1 -- Then b2 is gaining a link to bn because of the merge
+				nlinks[b2] = nlinks[b2] + 1
 				if nlinks[b2] == 9 then
-					lowlevel[ref[b2]] = lowlevel[cur]
+					lowlevel[b2] = nil
-					ref[lowlevel[cur]] = ref[b2]
-					cur = cur - 1
 				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
 
-	--------------------------------------------------------------
+	local queue = {[0] = -math.huge}
-	-- 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}
-	local queuevalues = {-math.huge}
-	cur = 1
 	local basin_lake = {}
 	for n=1, nbasins do
 		basin_lake[n] = 0
 	end
 	local reverse = {3, 4, 1, 2, [0]=0}
-	while cur > 0 do
+	for n=1, nbasins do
-		local b1, elev1 = queue[cur], queuevalues[cur] -- Pop from queue
+		local b1, elev1 = next(queue)
-		cur = cur - 1
+		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 i = bound.i
-			local dir = bound.is_y and 3 or 4 -- And get the direction (S/E/N/W)
+			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
@@ -402,12 +291,8 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 				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
@@ -417,65 +302,57 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 				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
+			until dir == 0
-
+			-- Add b2 into the queue
-			-- Add basin b2 into the queue, and keep the highest link elevation, that will define the elevation of the lake in b2
+			queue[b2] = mmax(elev1, bound.elev)
-			cur = cur + 1
-			queue[cur] = b2
-			queuevalues[cur] = 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
+	waterq = waterq or {}
-	-- 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.
 	local X, Y = dirs.X, dirs.Y
-	waterq = waterq or {X=X, Y=Y}
+	--local tinsert = table.insert
 
 	local ndonors = {}
+	local waterq = {X=X, Y=Y}
 	for i=1, X*Y do
 		ndonors[i] = 0
 		waterq[i] = 1
 	end
 
-	-- Calculate the number of direct donors
+	--for i1, dir in ipairs(dirs) do
-	for i=1, X*Y do
+	for i1=1, X*Y do
-		if dirs[i] == 1 then
+		local i2
-			ndonors[i+X] = ndonors[i+X] + 1
+		local dir = dirs[i1]
-		elseif dirs[i] == 2 then
+		if dir == 1 then
-			ndonors[i+1] = ndonors[i+1] + 1
+			i2 = i1+X
-		elseif dirs[i] == 3 then
+		elseif dir == 2 then
-			ndonors[i-X] = ndonors[i-X] + 1
+			i2 = i1+1
-		elseif dirs[i] == 4 then
+		elseif dir == 3 then
-			ndonors[i-1] = ndonors[i-1] + 1
+			i2 = i1-X
+		elseif dir == 4 then
+			i2 = i1-1
+		end
+		if i2 then
+			ndonors[i2] = ndonors[i2] + 1
 		end
 	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
@@ -486,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
@@ -507,4 +381,5 @@ end
 return {
 	flow_routing = flow_routing,
 	accumulate = accumulate,
+	flow_methods = flow_methods,
 }

terrainlib_lua/twist.lua

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