Adapt to multithreading, now it is working!

Also remove globally stored path for modpath and data.
Basically use less global variables for better adaptation to multithreading.

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,40 @@
+local sqrt, abs = math.sqrt, math.abs
+local unpk = unpack
+
+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
+
+return transform_quadri

gridio.lua

@@ -1,114 +0,0 @@
--- Input and output functions for grid maps
-
-local worldpath = mapgen_rivers.world_data_path
-
-local floor = math.floor
-local sbyte, schar = string.byte, string.char
-local unpk = unpack
-
--- Loading files
-local function load_full_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 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 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,
-}
-
-local function interactive_loader(filename, bytes, signed, size, converter)
-	local file = io.open(worldpath .. 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
-
-local load_methods = {
-	full        = load_full_map,
-	interactive = interactive_loader,
-}
-
-function mapgen_rivers.load_file(...)
-	local load_method = mapgen_rivers.settings.load_method
-	local load_func = load_methods[load_method]
-	if load_func then
-		return load_func(...)
-	else
-		minetest.log("error", ("[mapgen_rivers] Unknown load method %s"):format(load_method))
-	end
-end
-
--- Writing files
-function mapgen_rivers.write_file(filename, data, bytes)
-	local size = #data
-	local file = io.open(worldpath .. filename, 'wb')
-
-	local bytelist = {}
-	for j=1, bytes do
-		bytelist[j] = 0
-	end
-
-	for i=1, size do
-		local n = floor(data[i])
-		data[i] = n
-		for j=bytes, 2, -1 do
-			bytelist[j] = n % 256
-			n = floor(n / 256)
-		end
-		bytelist[1] = n % 256
-
-		file:write(schar(unpk(bytelist)))
-	end
-
-	file:close()
-end

gridmanager.lua

@@ -1,148 +0,0 @@
--- Manages grid loading, writing and generation
-
-local world_data = 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 grid_maps_list = {
-	dem = {bytes=2, signed=true},
-	lakes = {bytes=2, signed=true},
-	dirs = {bytes=1, signed=false},
-	rivers = {bytes=4, signed=false},
-
-	offset_x = {bytes=1, signed=true, conv=offset_conv},
-	offset_y = {bytes=1, signed=true, conv=offset_conv},
-}
-
-local function apply_grid_conversion(grid)
-	if grid.load_method ~= "full" then
-		minetest.log("warning", ("Could not apply data conversion for load method %s"):format(grid.load_method))
-		return false
-	end
-	if grid.conv_applied then
-		return true
-	end
-
-	local size = grid.size.x * grid.size.y
-	for mapname, params in pairs(grid_maps_list) do
-		local conv = params.conv
-		if conv then
-			local map = grid[mapname]
-			for i=1, size do
-				map[i] = conv(map[i])
-			end
-		end
-	end
-	grid.conv_applied = true
-
-	return true
-end
-
-function mapgen_rivers.try_load_grid(grid)
-	local load_method = mapgen_rivers.settings.load_method
-
-	-- First, check whether a grid is already loaded with the appropriate method
-	if mapgen_rivers.grid and mapgen_rivers.grid.load_method == load_method then
-		if not mapgen_rivers.grid.conv_applied then
-			apply_grid_conversion(mapgen_rivers.grid)
-		end
-		return true
-	-- Then, check the provided argument is a valid grid
-	elseif grid and grid.load_method == load_method then
-		if not mapgen_rivers.grid.conv_applied then
-			apply_grid_conversion(grid)
-		end
-		mapgen_rivers.grid = grid
-		on_grid_loaded_callback(grid)
-		return true
-	end
-
-	-- Fall back to loading the grid from the files
-	local sfile = io.open(world_data .. 'size', 'r')
-	if not sfile then
-		return false
-	end
-
-	local x, z = sfile:read('*n'), sfile:read('*n')
-	if not x or not z then
-		return false
-	end
-
-	if load_method == "full" then
-		minetest.log("action", '[mapgen_rivers] Loading full grid')
-	elseif load_method == "interactive" then
-		minetest.log("action", '[mapgen_rivers] Loading grid as interactive loaders')
-	end
-
-	grid = {
-		load_method = load_method,
-		size = {x=x, y=z},
-	}
-
-	for map, params in pairs(grid_maps_list) do
-		grid[map] = mapgen_rivers.load_file(map, params.bytes, params.signed, x*z, params.conv)
-	end
-	grid.conv_applied = true
-
-	mapgen_rivers.grid = grid
-	on_grid_loaded_callback(grid)
-
-	return true
-end
-
-function mapgen_rivers.generate_grid()
-	minetest.log("action", '[mapgen_rivers] Generating grid, this may take a while...')
-	local grid = {}
-
-	local blocksize = mapgen_rivers.settings.blocksize
-	local xsize = math.floor(mapgen_rivers.settings.map_x_size / blocksize)
-	local zsize = math.floor(mapgen_rivers.settings.map_z_size / blocksize)
-	grid.size = {x=xsize, y=zsize}
-	grid.conv_applied = false
-
-	if not mapgen_rivers.pregenerate then
-		minetest.log("error", "[mapgen_rivers] Pre-generation function is not available.")
-		return false
-	end
-
-	mapgen_rivers.pregenerate(grid)
-
-	return grid
-end
-
-function mapgen_rivers.write_grid(grid)
-	minetest.mkdir(world_data)
-
-	if grid.conv_applied then
-		minetest.log("error", '[mapgen_rivers] Could not write grid if data conversion is already done')
-		return false
-	end
-
-	for map, params in pairs(grid_maps_list) do
-		mapgen_rivers.write_file(map, grid[map], params.bytes)
-	end
-
-	local sfile = io.open(world_data .. 'size', "w")
-	sfile:write(grid.size.x..'\n'..grid.size.y)
-	sfile:close()
-
-	return true
-end
-

heightmap.lua

@@ -1,53 +1,17 @@
--- Transform polygon data into a heightmap
+local modpath = minetest.get_modpath(minetest.get_current_modname())
 
-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 sea_level = tonumber(minetest.get_mapgen_setting("water_level"))
+local riverbed_slope = tonumber(mapgen_rivers.settings:get("riverbed_slope")) * tonumber(mapgen_rivers.settings:get("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 floor, min, max = math.floor, math.min, math.max
 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)
 	local v0 = v01*xf + v00*(1-xf)
@@ -55,9 +19,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 = {}
@@ -82,57 +46,50 @@ function mapgen_rivers.make_heightmaps(minp, maxp)
 				-- = 0: on riverbank
 				-- > 0: inside river
 				local depth_factors = {
-					r_west  -    xf , -- West edge (1)
-					r_north -    zf , -- North edge (2)
-					r_east  - (1-xf), -- East edge (3)
-					r_south - (1-zf), -- South edge (4)
-					c_NW -    xf  -    zf , -- North-West corner (5)
-					c_NE - (1-xf) -    zf , -- North-East corner (6)
-					c_SE - (1-xf) - (1-zf), -- South-East corner (7)
-					c_SW -    xf  - (1-zf), -- South-West corner (8)
+					r_west - xf,
+					r_north - zf,
+					xf - r_east,
+					zf - r_south,
+					c_NW-xf-zf,
+					xf-zf-c_NE,
+					xf+zf-c_SE,
+					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.
-				-- If imax is still at 0, it means that we are not in a river.
-				local dpmax = 0
+				-- Find the maximal depth factor and determine to which river it belongs
+				local depth_factor_max = 0
 				local imax = 0
 				for i=1, 8 do
-					if depth_factors[i] > dpmax then
-						dpmax = depth_factors[i]
+					if depth_factors[i] >= depth_factor_max then
+						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)
-				local xfc, zfc
-				-- xfc:
-				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 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 x1 = 1-max(r_east-dpmax, c_NE-zf-dpmax, c_SE-(1-zf)-dpmax, 0) -- and bounded to 1 by eastern riverbank
-					if x0 >= x1 then
-						xfc = 0.5
-					else
-						xfc = (xf-x0) / (x1-x0)
-					end
-				elseif imax == 1 or imax == 5 or imax == 8 then -- river at the western side of the polygon
-					xfc = 0
-				else -- 3, 6, 7 : river at the eastern side of the polygon
-					xfc = 1
-				end
-
-				-- Same for zfc:
-				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
-					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
-					local z1 = 1-max(r_south-dpmax, c_SW-xf-dpmax, c_SE-(1-xf)-dpmax, 0) -- and bounded to 1 by southern riverbank
-					if z0 >= z1 then
-						zfc = 0.5
-					else
-						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
+				-- 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)
+				if imax == 0 then
+					local x0 = max(r_west, c_NW-zf, zf-c_SW)
+					local x1 = min(r_east, c_NE+zf, c_SE-zf)
+					local z0 = max(r_north, c_NW-xf, xf-c_NE)
+					local z1 = min(r_south, c_SW+xf, c_SE-xf)
+					xf = (xf-x0) / (x1-x0)
+					zf = (zf-z0) / (z1-z0)
+				elseif imax == 1 then
+					xf = 0
+				elseif imax == 2 then
+					zf = 0
+				elseif imax == 3 then
+					xf = 1
+				elseif imax == 4 then
+					zf = 1
+				elseif imax == 5 then
+					xf, zf = 0, 0
+				elseif imax == 6 then
+					xf, zf = 1, 0
+				elseif imax == 7 then
+					xf, zf = 1, 1
+				elseif imax == 8 then
+					xf, zf = 0, 1
 				end
 
 				-- Determine elevation by interpolation
@@ -142,12 +99,12 @@ function mapgen_rivers.make_heightmaps(minp, maxp)
 					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_z = xfc*(vdem[3]-vdem[2]) + (1-xfc)*(vdem[4]-vdem[1]) < 0
+				local slope_x = zf*(vdem[3]-vdem[4]) + (1-zf)*(vdem[2]-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
@@ -164,15 +121,15 @@ function mapgen_rivers.make_heightmaps(minp, maxp)
 				end
 				local lake_height = max(floor(poly.lake[lake_id]), terrain_height)
 
-				if imax > 0 and dpmax > 0 then
-					terrain_height = min(max(lake_height, sea_level) - floor(1+dpmax*riverbed_slope), terrain_height)
+				if imax > 0 and depth_factor_max > 0 then
+					terrain_height = min(max(lake_height, sea_level) - 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
-				lake_height_map[i] = out_elev
+				terrain_height_map[i] = MAP_BOTTOM
+				lake_height_map[i] = MAP_BOTTOM
 			end
 			i = i + 1
 		end
@@ -180,3 +137,5 @@ function mapgen_rivers.make_heightmaps(minp, maxp)
 
 	return terrain_height_map, lake_height_map
 end
+
+return heightmaps

init.lua

@@ -1,26 +1,43 @@
--- Main file, calls the other files and triggers main functions
-
 mapgen_rivers = {}
 
-local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
-mapgen_rivers.modpath = modpath
-mapgen_rivers.world_data_path = minetest.get_worldpath() .. '/river_data/'
+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
 
-dofile(modpath .. 'settings.lua')
-dofile(modpath .. 'gridmanager.lua')
-dofile(modpath .. 'gridio.lua')
-dofile(modpath .. 'polygons.lua')
-dofile(modpath .. 'heightmap.lua')
-dofile(modpath .. 'mapgen.lua')
+local modpath = minetest.get_modpath(minetest.get_current_modname())
+dofile(modpath .. '/settings.lua')
 
-minetest.register_on_mods_loaded(function()
-	local exist = mapgen_rivers.try_load_grid()
+local sfile = io.open(minetest.get_worldpath() .. '/river_data/size')
+if sfile then
+	sfile:close()
+else
+	dofile(modpath .. '/pregenerate.lua')
+	collectgarbage()
+end
 
-	if not exist then -- If grid does not exist yet, generate it
-		dofile(modpath .. 'pregenerate.lua')
-
-		local grid = mapgen_rivers.generate_grid()
-		mapgen_rivers.write_grid(grid)
-		mapgen_rivers.try_load_grid(grid) -- Reload if needed
+mapgen_rivers.use_mapgen_thread = minetest.settings:get_bool("mapgen_rivers_use_mapgen_thread")
+mapgen_rivers.thread = 'main'
+if mapgen_rivers.use_mapgen_thread then
+	if minetest.register_mapgen_dofile then
+		minetest.register_mapgen_dofile(modpath .. '/mapgen.lua')
+	else
+		minetest.log("warning", "[mapgen_rivers] Mapgen thread not available on this Minetest version.")
+		mapgen_rivers.use_mapgen_thread = false
 	end
-end)
+end
+
+if not mapgen_rivers.use_mapgen_thread then
+	dofile(modpath .. '/mapgen.lua')
+end
+
+-- Setup a metatable to load grid on request if not present
+local mt = {
+	__index = function(_, field)
+		if field == 'grid' then
+			dofile(modpath .. '/load_grid.lua')
+			return mapgen_rivers.grid
+		end
+	end,
+}
+setmetatable(mapgen_rivers, mt)

load_grid.lua

@@ -0,0 +1,107 @@
+local datapath = minetest.get_worldpath() .. "/river_data/"
+
+local floor = math.floor
+local sbyte = string.byte
+local unpk = unpack
+
+local load_map
+local use_interactive_loader
+if minetest.settings:has("mapgen_rivers_use_interactive_loader") then
+	use_interactive_loader = minetest.settings:get_bool("mapgen_rivers_use_interactive_loader")
+else
+	use_interactive_loader = not minetest.settings:get_bool("mapgen_rivers_load_all")
+end
+
+if use_interactive_loader then
+	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,
+	}
+
+	load_map = function(filename, bytes, signed, size, converter)
+		local file = io.open(datapath .. 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
+
+else
+
+	load_map = function(filename, bytes, signed, size, converter)
+		local file = io.open(datapath .. filename, 'rb')
+		local data = file:read('*all')
+		if #data < bytes*size then
+			data = minetest.decompress(data)
+		end
+
+		local map = {}
+
+		for i=1, size do
+			local i0 = (i-1)*bytes+1
+			local i1 = i*bytes
+			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
+end
+
+local sfile = io.open(datapath..'size', 'r')
+assert(sfile)
+local X, Z = tonumber(sfile:read('*l')), tonumber(sfile:read('*l'))
+sfile:close()
+
+local function offset_converter(o)
+	return (o + 0.5) * (1/256)
+end
+
+mapgen_rivers.grid = {
+	size = {x=X, y=Z},
+	dem = load_map('dem', 2, true, X*Z),
+	lakes = load_map('lakes', 2, true, X*Z),
+	dirs = load_map('dirs', 1, false, X*Z),
+	rivers = load_map('rivers', 4, false, X*Z),
+
+	offset_x = load_map('offset_x', 1, true, X*Z, offset_converter),
+	offset_y = load_map('offset_y', 1, true, X*Z, offset_converter);
+}

mapgen.lua

@@ -1,21 +1,40 @@
--- 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")
+-- Recreate mod table if we are in a separate environment
+if not minetest.global_exists("mapgen_rivers") then
+	mapgen_rivers = {use_mapgen_thread=true, thread='mapgen'}
+	mapgen_rivers.settings = Settings(minetest.get_worldpath() .. '/mapgen_rivers.conf')
 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
+if not mapgen_rivers.grid then
+	dofile(minetest.get_modpath(minetest.get_current_modname()) .. '/load_grid.lua')
+end
+
+local settings = mapgen_rivers.settings
+
+local sea_level = tonumber(minetest.get_mapgen_setting("water_level"))
+local elevation_chill = tonumber(settings:get('elevation_chill'))
+local use_distort = settings:get_bool('distort')
+local use_biomes = settings:get_bool('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 minetest.get_modpath("default") and not use_biomegen_mod
+
+local noiseparams = {
+	distort_x = settings:get_np_group('np_distort_x'),
+	distort_z = settings:get_np_group('np_distort_z'),
+	distort_amplitude = settings:get_np_group('np_distort_amplitude'),
+}
+
+if use_biomes then
+	noiseparams.heat = minetest.get_mapgen_setting_noiseparams("mg_biome_np_heat")
+	noiseparams.heat.offset = noiseparams.heat.offset + sea_level / elevation_chill
+	noiseparams.heat_blend = minetest.get_mapgen_setting_noiseparams("mg_biome_np_heat_blend")
+end
 
 if use_biomegen_mod then
 	biomegen.set_elevation_chill(elevation_chill)
 end
 
+local heightmaps = dofile(minetest.get_modpath(minetest.get_current_modname()) .. '/heightmap.lua')
+
 -- Linear interpolation
 local function interp(v00, v01, v11, v10, xf, zf)
 	local v0 = v01*xf + v00*(1-xf)
@@ -41,7 +60,24 @@ local sumtime = 0
 local sumtime2 = 0
 local ngen = 0
 
-function mapgen_rivers.make_chunk(minp, maxp, seed)
+local function init_mapgen(chulens)
+	mapsize = {
+		x = chulens.x,
+		y = chulens.y+1,
+		z = chulens.z,
+	}
+	if use_distort then
+		noise_x_obj = minetest.get_perlin_map(noiseparams.distort_x, mapsize)
+		noise_z_obj = minetest.get_perlin_map(noiseparams.distort_z, mapsize)
+		noise_distort_obj = minetest.get_perlin_map(noiseparams.distort_amplitude, chulens)
+	end
+	if use_biomes then
+		noise_heat_obj = minetest.get_perlin_map(noiseparams.heat, chulens)
+		noise_heat_blend_obj = minetest.get_perlin_map(noiseparams.heat_blend, chulens)
+	end
+end
+
+local function generate(vm, 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 = {
@@ -51,20 +87,7 @@ function mapgen_rivers.make_chunk(minp, maxp, seed)
 	}
 
 	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_mapgen(chulens)
 		init = true
 	end
 
@@ -109,9 +132,9 @@ function mapgen_rivers.make_chunk(minp, maxp, seed)
 		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)
+		terrain_map, lake_map = heightmaps(pminp, pmaxp)
 	else
-		terrain_map, lake_map = mapgen_rivers.make_heightmaps(minp, maxp)
+		terrain_map, lake_map = heightmaps(minp, maxp)
 	end
 
 	-- Check that there is at least one position that reaches min y
@@ -152,7 +175,7 @@ function mapgen_rivers.make_chunk(minp, maxp, seed)
 		c_ice = minetest.get_content_id("default:ice")
 	end
 
-	local vm, emin, emax = minetest.get_mapgen_object("voxelmanip")
+	local emin, emax = vm:get_emerged_area()
 	vm:get_data(data)
 
 	local a = VoxelArea:new({MinEdge = emin, MaxEdge = emax})
@@ -253,7 +276,9 @@ function mapgen_rivers.make_chunk(minp, maxp, seed)
 	vm:set_lighting({day = 0, night = 0})
 	vm:calc_lighting()
 	vm:update_liquids()
-	vm:write_to_map()
+	if mapgen_rivers.thread == "main" then
+		vm:write_to_map()
+	end
 
 	local t = os.clock()-t0
 	ngen = ngen + 1
@@ -262,9 +287,20 @@ function mapgen_rivers.make_chunk(minp, maxp, seed)
 	minetest.log("verbose", ("[mapgen_rivers] Done in %5.3f s"):format(t))
 end
 
-minetest.register_on_generated(mapgen_rivers.make_chunk)
+if mapgen_rivers.thread == "main" then
+	minetest.register_on_generated(function(minp, maxp, seed)
+		local vm = minetest.get_mapgen_object("voxelmanip")
+		generate(vm, minp, maxp, seed)
+	end)
+elseif mapgen_rivers.thread == "mapgen" then
+	minetest.register_on_generated(generate)
+end
 
 minetest.register_on_shutdown(function()
+	if ngen == 0 then
+		return
+	end
+
 	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))

polygons.lua

@@ -1,29 +1,26 @@
--- Fetch polygons from a given areas, and compute their properties
--- and find to which polygon every point belongs
-
-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.grid.size.x
+local Z = mapgen_rivers.grid.size.y
 
 local function index(x, z)
 	return z*X+x+1
 end
 
+local settings = mapgen_rivers.settings
+
+local blocksize = tonumber(settings:get('blocksize'))
+local min_catchment = tonumber(settings:get('min_catchment'))
+
 local map_offset = {x=0, z=0}
-mapgen_rivers.register_on_grid_loaded(function(grid)
-	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)
+if settings:get_bool('center') then
+	map_offset.x = blocksize*X/2
+	map_offset.z = blocksize*Z/2
+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
+min_catchment = min_catchment / (blocksize*blocksize)
+local wpower = settings:get('river_widening_power')
 local wfactor = 1/(2*blocksize * min_catchment^wpower)
 local function river_width(flow)
 	flow = abs(flow)
@@ -35,20 +32,20 @@ local function river_width(flow)
 end
 
 local noise_heat -- Need a large-scale noise here so no heat blend
-local elevation_chill = mapgen_rivers.settings.elevation_chill
+local elevation_chill = settings:get_bool('elevation_chill')
 local function get_temperature(x, y, z)
 	local pos = {x=x, y=z}
 	return noise_heat:get2d(pos) - y*elevation_chill
 end
 
-local glaciers = mapgen_rivers.settings.glaciers
-local glacier_factor = mapgen_rivers.settings.glacier_factor
+local glaciers = settings:get_bool('glaciers')
+local glacier_factor = tonumber(settings:get('glacier_factor'))
 
 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
@@ -165,15 +162,15 @@ function mapgen_rivers.make_polygons(minp, maxp)
 				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_south = (dirD==2 and riverD or 0) + (dirC==4 and riverC or 0)
+			local river_east = 1 - (dirB==1 and riverB or 0) - (dirC==3 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 +178,5 @@ function mapgen_rivers.make_polygons(minp, maxp)
 
 	return polygons
 end
+
+return make_polygons

pregenerate.lua

@@ -1,116 +1,131 @@
--- Generate the grid using terrainlib_lua
--- Only called on first mapgen, if there is no grid yet
+local modpath = minetest.get_modpath(minetest.get_current_modname())
 
-local EvolutionModel = dofile(mapgen_rivers.modpath .. '/terrainlib_lua/erosion.lua')
-local twist = dofile(mapgen_rivers.modpath .. '/terrainlib_lua/twist.lua')
+local EvolutionModel = dofile(modpath .. '/terrainlib_lua/erosion.lua')
+local twist = dofile(modpath .. '/terrainlib_lua/twist.lua')
 
-local blocksize = mapgen_rivers.settings.blocksize
-local tectonic_speed = mapgen_rivers.settings.tectonic_speed
+local blocksize = tonumber(mapgen_rivers.settings:get("blocksize"))
+local tectonic_speed = tonumber(mapgen_rivers.settings:get("tectonic_speed"))
 
-local np_base = table.copy(mapgen_rivers.noise_params.base)
-np_base.spread = vector.divide(np_base.spread, blocksize)
+local np_base = mapgen_rivers.settings:get_np_group("np_base")
+np_base.spread.x = np_base.spread.x / blocksize
+np_base.spread.y = np_base.spread.y / blocksize
+np_base.spread.z = np_base.spread.z / blocksize
 
-local evol_params = mapgen_rivers.settings.evol_params
+local evol_params = {
+	K = tonumber(mapgen_rivers.settings:get("river_erosion_coef")),
+	m = tonumber(mapgen_rivers.settings:get("river_erosion_power")),
+	d = tonumber(mapgen_rivers.settings:get("difusive_erosion")),
+}
 
-local time = mapgen_rivers.settings.evol_time
-local time_step = mapgen_rivers.settings.evol_time_step
+local time = tonumber(mapgen_rivers.settings:get("evol_time"))
+local time_step = tonumber(mapgen_rivers.settings:get("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
+-- Setup the model
+local size = {
+	x = tonumber(mapgen_rivers.settings:get("grid_x_size")),
+	y = tonumber(mapgen_rivers.settings:get("grid_z_size")),
+}
 
-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 seed = tonumber(minetest.get_mapgen_setting("seed"))
+np_base.seed = (np_base.seed or 0) + seed
 
-function mapgen_rivers.pregenerate(grid)
-	local size = grid.size
+local nobj_base = PerlinNoiseMap(np_base, {x=size.x, y=1, z=size.y})
 
-	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 dem = nobj_base:get_3d_map_flat({x=0, y=0, z=0})
+dem.X = size.x
+dem.Y = size.y
 
-	local seed = tonumber(minetest.get_mapgen_setting("seed"))
-	np_base.seed = (np_base.seed or 0) + seed
+local model = EvolutionModel(evol_params)
+model.dem = dem
+local ref_dem = model:define_isostasy(dem)
 
-	local nobj_base = PerlinNoiseMap(np_base, {x=size.x, y=1, z=size.y})
+local tectonic_step = tectonic_speed * time_step
+collectgarbage()
 
-	local dem = nobj_base:get_3d_map_flat({x=0, y=0, z=0})
-	dem.X = size.x
-	dem.Y = size.y
+-- Run the model
+for i=1, niter do
+	minetest.log("info", "[mapgen_rivers] Iteration " .. i .. " of " .. niter)
 
-	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()
-		model:erode(time_step)
-		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
+	model:diffuse(time_step)
 	model:flow()
-
-	local mfloor = math.floor
-	local mmin, mmax = math.min, math.max
-	local offset_x, offset_y = twist(model.dirs, model.rivers, 5)
-	for i=1, size.x*size.y do
-		offset_x[i] = mmin(mmax(offset_x[i]*256, -128), 127)
-		offset_y[i] = mmin(mmax(offset_y[i]*256, -128), 127)
+	model:erode(time_step)
+	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)
+		end
+		model:isostasy()
 	end
 
-	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
-
-	grid.load_method = "full"
-	grid.conv_applied = false
 	collectgarbage()
 end
+model:flow()
+
+local mfloor = math.floor
+local mmin, mmax = math.min, math.max
+local offset_x, offset_y = twist(model.dirs, model.rivers, 5)
+for i=1, size.x*size.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
+
+-- Write the results in the world directory
+local datapath = minetest.get_worldpath() .. "/river_data/"
+minetest.mkdir(datapath)
+
+local function write_map(filename, data, bytes)
+	local size = #data
+	local file = io.open(datapath .. filename, 'wb')
+
+	local bytelist = {}
+	for j=1, bytes do
+		bytelist[j] = 0
+	end
+
+    local unpk = unpack
+    local schar = string.char
+    local floor = math.floor
+	for i=1, size do
+		local n = floor(data[i])
+		data[i] = n
+		for j=bytes, 2, -1 do
+			bytelist[j] = n % 256
+			n = floor(n / 256)
+		end
+		bytelist[1] = n % 256
+
+		file:write(schar(unpk(bytelist)))
+	end
+
+	file:close()
+end
+
+write_map('dem', model.dem, 2)
+write_map('lakes', model.lakes, 2)
+write_map('dirs', model.dirs, 1)
+write_map('rivers', model.rivers, 4)
+write_map('offset_x', offset_x, 1)
+write_map('offset_y', offset_y, 1)
+local sfile = io.open(datapath .. 'size', "w")
+sfile:write(size.x..'\n'..size.y)
+sfile:close()
+
+local use_interactive_loader
+if minetest.settings:has("mapgen_rivers_use_interactive_loader") then
+	use_interactive_loader = minetest.settings:get_bool("mapgen_rivers_use_interactive_loader")
+else
+	use_interactive_loader = not minetest.settings:get_bool("mapgen_rivers_load_all")
+end
+
+if not use_interactive_loader then
+    mapgen_rivers.grid = {
+        size = size,
+	    dem = model.dem,
+	    lakes = model.lakes,
+	    dirs = model.dirs,
+	    rivers = model.rivers,
+	    offset_x = offset_x,
+	    offset_y = offset_y,
+    }
+end

settings.lua

@@ -1,12 +1,11 @@
--- Read global and per-world settings
-
 local mtsettings = minetest.settings
-local mgrsettings = Settings(minetest.get_worldpath() .. '/mapgen_rivers.conf')
+local settings = Settings(minetest.get_worldpath() .. '/mapgen_rivers.conf')
 
-mapgen_rivers.version = "1.0.2-dev1"
+mapgen_rivers.version = "1.0.2"
+mapgen_rivers.settings = settings
 
 local previous_version_mt = mtsettings:get("mapgen_rivers_version") or "0.0"
-local previous_version_mgr = mgrsettings:get("version") or "0.0"
+local previous_version_mgr = settings:get("version") or "0.0"
 
 if mapgen_rivers.version ~= previous_version_mt or mapgen_rivers.version ~= previous_version_mgr then
 	local compat_mt, compat_mgr = dofile(minetest.get_modpath(minetest.get_current_modname()) .. "/compatibility.lua")
@@ -14,124 +13,103 @@ if mapgen_rivers.version ~= previous_version_mt or mapgen_rivers.version ~= prev
 		compat_mt(mtsettings)
 	end
 	if mapgen_rivers.version ~= previous_version_mgr then
-		compat_mgr(mgrsettings)
+		compat_mgr(settings)
 	end
 end
 
 mtsettings:set("mapgen_rivers_version", mapgen_rivers.version)
-mgrsettings:set("version", mapgen_rivers.version)
+settings: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
+local modified = false
+local function verify_setting(name, dtype, default)
+	if settings:has(name) then
+		return
 	end
-end
 
-function mapgen_rivers.define_setting(name, dtype, default)
+	modified = true
+	local v = default
+	local mtname = 'mapgen_rivers_' .. name
+	local mthas = mtsettings:has(mtname)
 	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]
+		if mthas then
+			v = mtsettings:get(mtname)
+			if dtype == "number" and tonumber(v) == nil then
+				v = default
 			end
-			mgrsettings:set(name, v)
-		end
-		if dtype == "number" then
-			return tonumber(v)
-		else
-			return v
 		end
+		settings:set(name, v)
 	elseif dtype == "bool" then
-		local v = mgrsettings:get_bool(name)
-		if v == nil then
-			v = mtsettings:get_bool('mapgen_rivers_' .. name)
-			if v == nil then
-				v = defaults[name]
-			end
-			mgrsettings:set_bool(name, v)
+		if mthas then
+			v = mtsettings:get(mtname)
 		end
-		return v
+		settings:set_bool(name, v)
 	elseif dtype == "noise" then
-		local v = mgrsettings:get_np_group(name)
-		if v == nil then
-			v = mtsettings:get_np_group('mapgen_rivers_' .. name)
-			if v == nil then
-				v = defaults[name]
-			end
-			mgrsettings:set_np_group(name, v)
+		if mthas then
+			v = mtsettings:get_np_group(mtname)
 		end
-		clean_np(v)
-		return v
+		settings:set_np_group(name, v)
 	end
 end
 
-local def_setting = mapgen_rivers.define_setting
+verify_setting('center', 'bool', true)
+verify_setting('blocksize', 'number', 15)
+verify_setting('min_catchment', 'number', 3600)
+verify_setting('river_widening_power', 'number', 0.5)
+verify_setting('riverbed_slope', 'number', 0.4)
+verify_setting('distort', 'bool', true)
+verify_setting('biomes', 'bool', true)
+verify_setting('glaciers', 'bool', false)
+verify_setting('glacier_factor', 'number', 8)
+verify_setting('elevation_chill', 'number', 0.25)
+verify_setting('grid_x_size', 'number', 1000)
+verify_setting('grid_z_size', 'number', 1000)
+verify_setting('river_erosion_coef', 'number', 0.5)
+verify_setting('river_erosion_power', 'number', 0.4)
+verify_setting('diffusive_erosion', 'number', 0.5)
+verify_setting('compensation_radius', 'number', 50)
+verify_setting('tectonic_speed', 'number', 70)
+verify_setting('evol_time', 'number', 10)
+verify_setting('evol_time_step', 'number', 1)
 
-mapgen_rivers.settings = {
-	center = def_setting('center', 'bool'),
-	blocksize = def_setting('blocksize', 'number'),
-	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'),
+verify_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",
+})
+verify_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,
+})
+verify_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,
+})
+verify_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",
+})
 
-	map_x_size = def_setting('map_x_size', 'number'),
-	map_z_size = def_setting('map_z_size', 'number'),
-	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'),
-		m = def_setting('river_erosion_power', 'number'),
-		d = def_setting('diffusive_erosion', 'number'),
-		compensation_radius = def_setting('compensation_radius', 'number'),
-	},
-	tectonic_speed = def_setting('tectonic_speed', 'number'),
-	evol_time = def_setting('evol_time', 'number'),
-	evol_time_step = def_setting('evol_time_step', 'number'),
-
-	load_all = mtsettings:get_bool('mapgen_rivers_load_all')
-}
-
-mapgen_rivers.settings.load_method =
-		mapgen_rivers.settings.load_all and "full" or "interactive"
-
-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()
+if modified then
+	settings:write()
 end
-
-minetest.register_on_mods_loaded(write_settings)
-minetest.register_on_shutdown(write_settings)

settings_default.json

@@ -1,71 +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,
-	"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"
-	}
-}

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.
-#    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
-
-#    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
+#    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
 
 #    Minimal catchment area for a river to be drawn, in square nodes
 #    Lower value means bigger river density

terrainlib_lua/erosion.lua

@@ -123,18 +123,14 @@ end
 
 local modpath = ""
 if minetest then
-	if minetest.global_exists('mapgen_rivers') then
-		modpath = mapgen_rivers.modpath .. "terrainlib_lua/"
-	else
-		modpath = minetest.get_modpath(minetest.get_current_modname()) .. "terrainlib_lua/"
-	end
+	modpath = minetest.get_modpath(minetest.get_current_modname()) .. "/terrainlib_lua/"
 end
 
 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
 

terrainlib_lua/rivermapper.lua

@@ -10,18 +10,50 @@
 -- 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.
+-- 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
+	end
+
+	if sum == 0 then
+		return 0
+	end
+	local r = math.random() * sum
+	for i=1, #plist do
+		local p = plist[i]
+		if r < p then
+			return i
+		end
+		r = r - p
+	end
+	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) -- 'dirs' and 'lakes' are optional tables to reuse for memory optimization, they may contain any data.
+local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are optional tables to reuse for memory optimization, they may contain any data.
+	method = method or 'semirandom'
+	local flow_local = flow_methods[method] or flow_local_semirandom
+
 	dirs = dirs or {}
 	lakes = lakes or {}
 
 	-- Localize for performance
 	local tremove = table.remove
 	local mmax = math.max
-	local mrand = math.random
 
 	local X, Y = dem.X, dem.Y
 	dirs.X = X
@@ -42,29 +74,14 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 	for y=1, Y do
 		for x=1, X do
 			local zi = dem[i]
-			-- Determine 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, with probability depending 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.
-			local pSouth = y<Y and mmax(zi-dem[i+X], 0) or 0
-			local pEast = x<X and mmax(zi-dem[i+1], 0) or 0
-			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 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
+			local plist = { -- Get the height difference of the 4 neighbours (and 0 if uphill)
+				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
+				x>1 and mmax(zi-dem[i-1], 0) or 0, -- Westward
+			}
 
+			local d = flow_local(plist)
 			-- 'dirs': Direction toward which water flow
 			-- 'dirs2': Directions from which water comes
 			dirs[i] = d
@@ -93,172 +110,89 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 	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
+			b2 = 0
+		else
+			b2 = basin_id[i2]
+			if b2 == 0 then -- If basin of i2 is not already computed, skip
+				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)
+			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
+				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
 
-	local cur = nbasins
-	local ib = 0
-	while cur > 0 do
-			local i = singular[cur]
-			cur = cur - 1
-			if dirs[i] == 0 then
-				basin_links = {}
-				links[#links+1] = basin_links
-				ib = ib + 1
-			end
+	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
+		basin_links = {}
+		links[#links+1] = 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
 
 			-- Iterate through the 4 directions
-			-- 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
-				singular[cur] = i+1
+				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.
 			elseif i%X > 0 then
-				if basin_id[i+1] ~= ib and basin_id[i+1] ~= 0 then
-					local b2 = basin_id[i+1]
-					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)
-					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
+				add_link(i, i+1, ib, false)
 			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
+				add_link(i, 0, ib, false)
 			end
 
-			-----------
-			-- SOUTH --
-			-----------
 			if d >= 4 then -- River coming from the South
 				d = d - 4
-				cur = cur + 1
-				singular[cur] = i+X
+				queue[#queue+1] = i+X
 			elseif i <= X*(Y-1) then
-				if basin_id[i+X] ~= ib and basin_id[i+X] ~= 0 then
-					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
+				add_link(i, i+X, ib, true)
 			else
-				local l2 = basin_links[0]
-				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
+				add_link(i, 0, ib, true)
 			end
 
-			----------
-			-- WEST --
-			----------
 			if d >= 2 then -- River coming from the West
 				d = d - 2
-				cur = cur + 1
-				singular[cur] = i-1
+				queue[#queue+1] = i-1
 			elseif i%X ~= 1 then
-				if basin_id[i-1] ~= ib and basin_id[i-1] ~= 0 then
-					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
+				add_link(i, i-1, ib, false)
 			else
-				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
+				add_link(i, 0, ib, false)
 			end
 
-			-----------
-			-- NORTH --
-			-----------
 			if d >= 1 then -- River coming from the North
-				cur = cur + 1
-				singular[cur] = i-X
+				queue[#queue+1] = i-X
 			elseif i > X then
-				if basin_id[i-X] ~= ib and basin_id[i-X] ~= 0 then
-					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
+				add_link(i, i-X, ib, true)
 			else
-				local l2 = basin_links[0]
-				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
+				add_link(i, 0, ib, true)
 			end
+		end
 	end
 	dirs2 = nil
 
@@ -289,31 +223,26 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 	-- 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
-	local ref = singular -- Reuse table
-	for i=0, nbasins do
-		if nlinks[i] <= 8 then
-			cur = cur + 1
-			lowlevel[cur] = i
-			ref[i] = cur
+	for i, n in pairs(nlinks) do
+		if n <= 8 then
+			lowlevel[i] = links[i]
 		end
 	end
 
 	local basin_graph = {}
-	for i=0, nbasins do
-		basin_graph[i] = {} -- Initialize (to ensure subtables don't go in the hash part)
-	end
-	for i=1, nbasins do
+	for n=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]
+		-- 'next' called with only one argument always returns an element if table is not empty
+		local b1, lnk1 = next(lowlevel)
+		lowlevel[b1] = nil
 
 		local b2
 		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
 				lowest = bdata.elev
 				b2 = bn
@@ -323,7 +252,13 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 		-- Add link to the graph, in both directions
 		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
@@ -334,9 +269,7 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 		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[cur] = b2
-			ref[b2] = cur
+			lowlevel[b2] = lnk2
 		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
@@ -346,16 +279,12 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 			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 nlinks[bn] == 8 then
-					cur = cur + 1
-					lowlevel[cur] = bn
-					ref[bn] = cur
+					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
 				if nlinks[b2] == 9 then
-					lowlevel[ref[b2]] = lowlevel[cur]
-					ref[lowlevel[cur]] = ref[b2]
-					cur = cur - 1
+					lowlevel[b2] = nil
 				end
 			end
 
@@ -373,17 +302,15 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 	-- 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 queue = {[0] = -math.huge}
 	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
-		local b1, elev1 = queue[cur], queuevalues[cur] -- Pop from queue
-		cur = cur - 1
+	for n=1, nbasins do
+		local b1, elev1 = next(queue) -- Pop from 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
@@ -422,9 +349,7 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 			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
-			cur = cur + 1
-			queue[cur] = b2
-			queuevalues[cur] = mmax(elev1, bound.elev)
+			queue[b2] = mmax(elev1, bound.elev)
 			-- Remove b1 from b2's neighbours to avoid coming back to b1
 			basin_graph[b2][b1] = nil
 		end
@@ -447,25 +372,31 @@ local function accumulate(dirs, waterq)
 	-- 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
-	waterq = waterq or {X=X, Y=Y}
 
 	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 i=1, X*Y do
-		if dirs[i] == 1 then
-			ndonors[i+X] = ndonors[i+X] + 1
-		elseif dirs[i] == 2 then
-			ndonors[i+1] = ndonors[i+1] + 1
-		elseif dirs[i] == 3 then
-			ndonors[i-X] = ndonors[i-X] + 1
-		elseif dirs[i] == 4 then
-			ndonors[i-1] = ndonors[i-1] + 1
+	for i1=1, X*Y do
+		local i2
+		local dir = dirs[i1]
+		if dir == 1 then
+			i2 = i1+X
+		elseif dir == 2 then
+			i2 = i1+1
+		elseif dir == 3 then
+			i2 = i1-X
+		elseif dir == 4 then
+			i2 = i1-1
+		end
+		if i2 then
+			ndonors[i2] = ndonors[i2] + 1
 		end
 	end
 
@@ -507,4 +438,5 @@ end
 return {
 	flow_routing = flow_routing,
 	accumulate = accumulate,
+	flow_methods = flow_methods,
 }