Use settings from minetest.conf

This prevents rounding errors and improves interpolation on nearly flat areas

.gitignore

@@ -1,7 +1,10 @@
 __pycache__/
 dem
 lakes
-links
-rivers
 size
+offset_x
+offset_y
+bounds_x
+bounds_y
+rivers
 unused/

README.md

@@ -5,6 +5,8 @@ Procedural map generator for Minetest 5.x. Still experimental and basic.
 
 Contains two distinct programs: Python scripts for pre-processing, and Lua scripts to generate the map on Minetest.
 
+![Screenshot](https://user-images.githubusercontent.com/6905002/79073532-7a567f00-7ce7-11ea-9791-8fb453f5175d.png)
+
 # Installation
 This mod should be placed in the `/mods` directory like any other Minetest mod.
 

bounds.py

@@ -0,0 +1,62 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+def make_bounds(dirs, rivers):
+    (Y, X) = dirs.shape
+    bounds_h = np.zeros((Y, X-1), dtype='i4')
+    bounds_v = np.zeros((Y-1, X), dtype='i4')
+
+    bounds_v += (rivers * (dirs==1))[:-1,:]
+    bounds_h += (rivers * (dirs==2))[:,:-1]
+    bounds_v -= (rivers * (dirs==3))[1:,:]
+    bounds_h -= (rivers * (dirs==4))[:,1:]
+
+    return bounds_h, bounds_v
+
+def get_fixed(dirs):
+    borders = np.zeros(dirs.shape, dtype='?')
+    borders[-1,:] |= dirs[-1,:]==1
+    borders[:,-1] |= dirs[:,-1]==2
+    borders[0,:] |= dirs[0,:]==3
+    borders[:,0] |= dirs[:,0]==4
+
+    donors = np.zeros(dirs.shape, dtype='?')
+    donors[1:,:] |= dirs[:-1,:]==1
+    donors[:,1:] |= dirs[:,:-1]==2
+    donors[:-1,:] |= dirs[1:,:]==3
+    donors[:,:-1] |= dirs[:,1:]==4
+    return borders | ~donors
+
+def twist(bounds_x, bounds_y, fixed, d=0.1, n=5):
+    moveable = ~fixed
+
+    (Y, X) = fixed.shape
+    offset_x = np.zeros((Y, X))
+    offset_y = np.zeros((Y, X))
+
+    for i in range(n):
+        force_long = np.abs(bounds_x) * (1+np.diff(offset_x, axis=1))
+        force_trans = np.abs(bounds_y) * np.diff(offset_x, axis=0)
+
+        force_x = np.zeros((Y, X))
+        force_x[:,:-1] = force_long
+        force_x[:,1:] -= force_long
+        force_x[:-1,:]+= force_trans
+        force_x[1:,:] -= force_trans
+
+        force_long = np.abs(bounds_y) * (1+np.diff(offset_y, axis=0))
+        force_trans = np.abs(bounds_x) * np.diff(offset_y, axis=1)
+
+        force_y = np.zeros((Y, X))
+        force_y[:-1,:] = force_long
+        force_y[1:,:] -= force_long
+        force_y[:,:-1]+= force_trans
+        force_y[:,1:] -= force_trans
+
+        length = np.hypot(force_x, force_y)
+        length[length==0] = 1
+        coeff = d / length * moveable
+        offset_x += force_x * coeff
+        offset_y += force_y * coeff
+
+    return offset_x, offset_y

geometry.lua

@@ -0,0 +1,45 @@
+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
+	local b = (x1-x)^2 + (y1-y)^2
+	local c = (x2-x)^2 + (y2-y)^2
+	if a + b < c then
+		return math.sqrt(b)
+	elseif a + c < b then
+		return math.sqrt(c)
+	else
+		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)
+	-- X, Y 4-vectors giving the coordinates of the 4 nodes
+	-- x, y position to index.
+	local x1, x2, x3, x4 = unpack(X)
+	local y1, y2, y3, y4 = unpack(Y)
+
+	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)
+
+	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
+
+local function area(X, Y) -- Signed area of polygon, in function of direction of rotation. Clockwise = positive.
+	local n = #X
+	local sum = X[1]*Y[n] - X[n]*Y[1]
+	for i=2, n do
+		sum = sum + X[i]*Y[i-1] - X[i-1]*Y[i]
+	end
+
+	return sum/2
+end
+
+return {
+	distance_to_segment = distance_to_segment,
+	transform_quadri = transform_quadri,
+	area = area,
+}

init.lua

@@ -1,6 +1,7 @@
 local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
 local worldpath = minetest.get_worldpath() .. '/'
 local load_map = dofile(modpath .. 'load.lua')
+local geometry = dofile(modpath .. 'geometry.lua')
 
 local function copy_if_needed(filename)
 	local wfilename = worldpath..filename
@@ -23,45 +24,70 @@ local X = tonumber(sfile:read('*l'))
 local Z = tonumber(sfile:read('*l'))
 
 copy_if_needed('dem')
-local dem = load_map(worldpath..'dem', 2, true)
+local dem = load_map(worldpath..'dem', 2, true, X*Z)
 copy_if_needed('lakes')
-local lakes = load_map(worldpath..'lakes', 2, true)
+local lakes = load_map(worldpath..'lakes', 2, true, X*Z)
-copy_if_needed('links')
+copy_if_needed('bounds_x')
-local links = load_map(worldpath..'links', 1, false)
+local bounds_x = load_map(worldpath..'bounds_x', 4, false, (X-1)*Z)
-copy_if_needed('rivers')
+copy_if_needed('bounds_y')
-local rivers = load_map(worldpath..'rivers', 4, false)
+local bounds_z = load_map(worldpath..'bounds_y', 4, false, X*(Z-1))
+
+copy_if_needed('offset_x')
+local offset_x = load_map(worldpath..'offset_x', 1, true, X*Z)
+for k, v in ipairs(offset_x) do
+	offset_x[k] = (v+0.5)/256
+end
+
+copy_if_needed('offset_y')
+local offset_z = load_map(worldpath..'offset_y', 1, true, X*Z)
+for k, v in ipairs(offset_z) do
+	offset_z[k] = (v+0.5)/256
+end
+
 
 local function index(x, z)
 	return z*X+x+1
 end
 
-local function interp(v00, v01, v10, v11, xf, zf)
+local function get_point_location(x, z)
-	v0 = v01*xf + v00*(1-xf)
+	local i = index(x, z)
-	v1 = v11*xf + v10*(1-xf)
+	return x+offset_x[i], z+offset_z[i]
+end
+
+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 blocksize = 6
+local blocksize = 12
 local sea_level = 1
 local min_catchment = 25
+local max_catchment = 40000
+local riverbed_slope = 0.4
 
-local storage = minetest.get_mod_storage()
+local get_settings = dofile(modpath .. 'settings.lua')
-if storage:contains("blocksize") then
+
-	blocksize = storage:get_int("blocksize")
+blocksize = get_settings('blocksize', 'int', blocksize)
-else
+sea_level = get_settings('sea_level', 'int', sea_level)
-	storage:set_int("blocksize", blocksize)
+min_catchment = get_settings('min_catchment', 'float', min_catchment)
-end
+max_catchment = get_settings('max_catchment', 'float', max_catchment)
-if storage:contains("sea_level") then
+riverbed_slope = get_settings('riverbed_slope', 'float', riverbed_slope) * blocksize
-	sea_level = storage:get_int("sea_level")
+
-else
+-- Width coefficients: coefficients solving
-	storage:set_int("sea_level", sea_level)
+--   wfactor * min_catchment ^ wpower = 1/(2*blocksize)
-end
+--   wfactor * max_catchment ^ wpower = 1
-if storage:contains("min_catchment") then
+local wpower = math.log(2*blocksize)/math.log(max_catchment/min_catchment)
-	min_catchment = storage:get_float("min_catchment")
+local wfactor = 1 / max_catchment ^ wpower
-else
+local function river_width(flow)
-	storage:set_float("min_catchment", min_catchment)
+	flow = math.abs(flow)
+	if flow < min_catchment then
+		return 0
+	end
+
+	return math.min(wfactor * flow ^ wpower, 1)
 end
 
 local function generate(minp, maxp, seed)
@@ -77,69 +103,175 @@ local function generate(minp, maxp, seed)
 
 	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 chulens = maxp.z - minp.z + 1
 
+	local polygons = {}
+	local xpmin, xpmax = math.max(math.floor(minp.x/blocksize - 0.5), 0), math.min(math.ceil(maxp.x/blocksize), X-2)
+	local zpmin, zpmax = math.max(math.floor(minp.z/blocksize - 0.5), 0), math.min(math.ceil(maxp.z/blocksize), Z-2)
+
+	for xp = xpmin, xpmax do
+		for zp=zpmin, zpmax do
+			local iA = index(xp, zp)
+			local iB = index(xp+1, zp)
+			local iC = index(xp+1, zp+1)
+			local iD = index(xp, zp+1)
+			local poly_x = {offset_x[iA]+xp, offset_x[iB]+xp+1, offset_x[iC]+xp+1, offset_x[iD]+xp}
+			local poly_z = {offset_z[iA]+zp, offset_z[iB]+zp, offset_z[iC]+zp+1, offset_z[iD]+zp+1}
+			local polygon = {x=poly_x, z=poly_z, i={iA, iB, iC, iD}}
+
+			local bounds = {}
+			local xmin = math.max(math.floor(blocksize*math.min(unpack(poly_x)))+1, minp.x)
+			local xmax = math.min(math.floor(blocksize*math.max(unpack(poly_x))), maxp.x)
+			for x=xmin, xmax do
+				bounds[x] = {}
+			end
+
+			local i1 = 4
+			for i2=1, 4 do -- Loop on 4 edges
+				local x1, x2 = poly_x[i1], poly_x[i2]
+				local lxmin = math.floor(blocksize*math.min(x1, x2))+1
+				local lxmax = math.floor(blocksize*math.max(x1, x2))
+				if lxmin <= lxmax then
+					local z1, z2 = poly_z[i1], poly_z[i2]
+					local a = (z1-z2) / (x1-x2)
+					local b = blocksize*(z1 - a*x1)
+					for x=math.max(lxmin, minp.x), math.min(lxmax, maxp.x) do
+						table.insert(bounds[x], a*x+b)
+					end
+				end
+				i1 = i2
+			end
+			for x=xmin, xmax do
+				local xlist = bounds[x]
+				table.sort(xlist)
+				local c = math.floor(#xlist/2)
+				for l=1, c do
+					local zmin = math.max(math.floor(xlist[l*2-1])+1, minp.z)
+					local zmax = math.min(math.floor(xlist[l*2]), maxp.z)
+					local i = (x-minp.x) * chulens + (zmin-minp.z) + 1
+					for z=zmin, zmax do
+						polygons[i] = polygon
+						i = i + 1
+					end
+				end
+			end
+
+			polygon.dem = {dem[iA], dem[iB], dem[iC], dem[iD]}
+			polygon.lake = math.min(lakes[iA], lakes[iB], lakes[iC], lakes[iD])
+
+			local river_west = river_width(bounds_z[iA])
+			local river_north = river_width(bounds_x[iA-zp])
+			local river_east = 1-river_width(bounds_z[iB])
+			local river_south = 1-river_width(bounds_x[iD-zp-1])
+			if river_west > river_east then
+				local mean = (river_west + river_east) / 2
+				river_west = mean
+				river_east = mean
+			end
+			if river_north > river_south then
+				local mean = (river_north + river_south) / 2
+				river_north = mean
+				river_south = mean
+			end
+			polygon.rivers = {river_west, river_north, river_east, river_south}
+
+			local around = {0,0,0,0,0,0,0,0}
+			if zp > 0 then
+				around[1] = river_width(bounds_z[iA-X])
+				around[2] = river_width(bounds_z[iB-X])
+			end
+			if xp < X-2 then
+				around[3] = river_width(bounds_x[iB-zp])
+				around[4] = river_width(bounds_x[iC-zp-1])
+			end
+			if zp < Z-2 then
+				around[5] = river_width(bounds_z[iC])
+				around[6] = river_width(bounds_z[iD])
+			end
+			if xp > 0 then
+				around[7] = river_width(bounds_x[iD-zp-2])
+				around[8] = river_width(bounds_x[iA-zp-1])
+			end
+			polygon.river_corners = {math.max(around[8], around[1]), math.max(around[2], around[3]), math.max(around[4], around[5]), math.max(around[6], around[7])}
+		end
+	end
+
+	local i = 1
 	for x = minp.x, maxp.x do
 		for z = minp.z, maxp.z do
-			local xb = x/blocksize
+			local poly = polygons[i]
-			local zb = z/blocksize
+			if poly then
+				local xf, zf = geometry.transform_quadri(poly.x, poly.z, x/blocksize, z/blocksize)
+				local i00, i01, i11, i10 = unpack(poly.i)
 
-			if xb >= 0 and xb < X-1 and zb >= 0 and zb < Z-1 then
+				local is_river = false
-				local x0 = math.floor(xb)
+				local depth_factor = 0
-				local x1 = x0+1
+				local r_west, r_north, r_east, r_south = unpack(poly.rivers)
-				local z0 = math.floor(zb)
+				if xf >= r_east then
-				local z1 = z0+1
+					is_river = true
+					depth_factor = xf-r_east
+					xf = 1
+				elseif xf <= r_west then
+					is_river = true
+					depth_factor = r_west-xf
+					xf = 0
+				end
+				if zf >= r_south then
+					is_river = true
+					depth_factor = zf-r_south
+					zf = 1
+				elseif zf <= r_north then
+					is_river = true
+					depth_factor = r_north-zf
+					zf = 0
+				end
 
-				local xf = xb - x0
+				if not is_river then
-				local zf = zb - z0
+					local c_NW, c_NE, c_SE, c_SW = unpack(poly.river_corners)
+					if xf+zf <= c_NW then
+						is_river = true
+						depth_factor = c_NW-xf-zf
+						xf, zf = 0, 0
+					elseif 1-xf+zf <= c_NE then
+						is_river = true
+						depth_factor = c_NE-1+xf-zf
+						xf, zf = 1, 0
+					elseif 2-xf-zf <= c_SE then
+						is_river = true
+						depth_factor = c_SE-2+xf+zf
+						xf, zf = 1, 1
+					elseif xf+1-zf <= c_SW then
+						is_river = true
+						depth_factor = c_SW-xf-1+zf
+						xf, zf = 0, 1
+					end
+				end
 
-				local i00 = index(x0,z0)
+				if not is_river then
-				local i01 = index(x1,z0)
+					xf = (xf-r_west) / (r_east-r_west)
-				local i10 = index(x0,z1)
+					zf = (zf-r_north) / (r_south-r_north)
-				local i11 = index(x1,z1)
+				end
 
-				local terrain_height = math.floor(interp(
+				local vdem = poly.dem
-					dem[i00],
+				local terrain_height = math.floor(0.5+interp(
-					dem[i01],
+					vdem[1],
-					dem[i10],
+					vdem[2],
-					dem[i11],
+					vdem[3],
+					vdem[4],
 					xf, zf
 				))
 
-				local lake_height = math.floor(math.min(
+				local lake_height = math.max(math.floor(poly.lake), terrain_height)
-					lakes[i00],
+				if is_river then
-					lakes[i01],
+					terrain_height = math.min(math.max(lake_height, sea_level) - math.floor(1+depth_factor*riverbed_slope), terrain_height)
-					lakes[i10],
+				end
-					lakes[i11]
-				))
-
 				local is_lake = lake_height > terrain_height
-
-				local is_river = false
-				if xf == 0 then
-					if links[i00] == 1 and rivers[i00] >= min_catchment then
-						is_river = true
-					elseif links[i10] == 3 and rivers[i10] >= min_catchment then
-						is_river = true
-					end
-				end
-
-				if zf == 0 then
-					if links[i00] == 2 and rivers[i00] >= min_catchment then
-						is_river = true
-					elseif links[i01] == 4 and rivers[i01] >= min_catchment then
-						is_river = true
-					end
-				end
-				
 				local ivm = a:index(x, minp.y-1, z)
-
 				if terrain_height >= minp.y then
 					for y=minp.y, math.min(maxp.y, terrain_height) do
 						if y == terrain_height then
 							if is_lake or y <= sea_level then
 								data[ivm] = c_sand
-							elseif is_river then
-								data[ivm] = c_rwater
 							else
 								data[ivm] = c_lawn
 							end
@@ -164,6 +296,7 @@ local function generate(minp, maxp, seed)
 					end
 				end
 			end
+			i = i + 1
 		end
 	end
 

load.lua

@@ -1,11 +1,12 @@
-local function load_map(filename, bytes, signed)
+local function load_map(filename, bytes, signed, size)
 	local file = io.open(filename, 'r')
 	local data = file:read('*all')
+	if #data < bytes*size then
+		data = minetest.decompress(data)
+	end
 
 	local map = {}
 
-	local size = math.floor(#data/bytes)
-
 	for i=1, size do
 		local i0, i1 = (i-1)*bytes+1, i*bytes
 		local elements = {data:byte(i0, i1)}

rivermapper.py

@@ -18,7 +18,7 @@ neighbours_dirs = np.array([
 
 neighbours_pattern = neighbours_dirs > 0
 
-def flow_dirs_lakes(dem, random=0.0625):
+def flow_dirs_lakes(dem, random=0):
     (Y, X) = dem.shape
 
     dem_margin = np.zeros((Y+2, X+2))

save.py

@@ -1,8 +1,13 @@
 import numpy as np
+import zlib
 
 def save(data, fname, dtype=None):
     if dtype is not None:
         data = data.astype(dtype)
 
+    bin_data = data.tobytes()
+    bin_data_comp = zlib.compress(bin_data, 9)
+    if len(bin_data_comp) < len(bin_data):
+        bin_data = bin_data_comp
     with open(fname, 'wb') as f:
-        f.write(data.tobytes())
+        f.write(bin_data)

settings.lua

@@ -0,0 +1,41 @@
+local storage = minetest.get_mod_storage()
+local settings = minetest.settings
+
+local function get_settings(key, dtype, default)
+	if storage:contains(key) then
+		if dtype == "string" then
+			return storage:get_string(key)
+		elseif dtype == "int" then
+			return storage:get_int(key)
+		elseif dtype == "float" then
+			return storage:get_float(key)
+		end
+	end
+
+	local conf_val = settings:get('mapgen_rivers_' .. key)
+	if conf_val then
+		if dtype == "int" then
+			conf_val = tonumber(conf_val)
+			storage:set_int(key, conf_val)
+		elseif dtype == "float" then
+			conf_val = tonumber(conf_val)
+			storage:set_float(key, conf_val)
+		elseif dtype == "string" then
+			storage:set_string(key, conf_val)
+		end
+
+		return conf_val
+	else
+		if dtype == "int" then
+			storage:set_int(key, default)
+		elseif dtype == "float" then
+			storage:set_float(key, default)
+		elseif dtype == "string" then
+			storage:set_string(key, default)
+		end
+
+		return default
+	end
+end
+
+return get_settings

terrain_rivers.py

@@ -4,6 +4,7 @@ import numpy as np
 import noise
 from save import save
 from erosion import EvolutionModel
+import bounds
 import os
 import sys
 
@@ -22,7 +23,7 @@ n = np.zeros((mapsize, mapsize))
 #micronoise_depth = 0.05
 
 params = {
-    "octaves" : 8,
+    "octaves" : int(np.log2(mapsize)),
     "persistence" : 0.5,
     "lacunarity" : 2.,
 }
@@ -67,9 +68,19 @@ model.calculate_flow()
 
 print('Done')
 
+bx, by = bounds.make_bounds(model.dirs, model.rivers)
+ox, oy = bounds.twist(bx, by, bounds.get_fixed(model.dirs))
+
+offset_x = np.clip(np.floor(ox * 256), -128, 127)
+offset_y = np.clip(np.floor(oy * 256), -128, 127)
+
 save(model.dem, 'dem', dtype='>i2')
 save(model.lakes, 'lakes', dtype='>i2')
-save(model.dirs, 'links', dtype='u1')
+save(np.abs(bx), 'bounds_x', dtype='>i4')
+save(np.abs(by), 'bounds_y', dtype='>i4')
+save(offset_x, 'offset_x', dtype='i1')
+save(offset_y, 'offset_y', dtype='i1')
+
 save(model.rivers, 'rivers', dtype='>u4')
 
 with open('size', 'w') as sfile:

view_map.py

@@ -0,0 +1,39 @@
+#!/usr/bin/env python3
+
+import numpy as np
+import zlib
+import matplotlib.pyplot as plt
+
+def load_map(name, dtype, shape):
+    dtype = np.dtype(dtype)
+    with open(name, 'rb') as f:
+        data = f.read()
+    if len(data) < shape[0]*shape[1]*dtype.itemsize:
+        data = zlib.decompress(data)
+    return np.frombuffer(data, dtype=dtype).reshape(shape)
+
+shape = np.loadtxt('size', dtype='u4')
+n = shape[0] * shape[1]
+dem = load_map('dem', '>i2', shape)
+lakes = load_map('lakes', '>i2', shape)
+rivers = load_map('rivers', '>u4', shape)
+
+plt.subplot(1,3,1)
+plt.pcolormesh(dem, cmap='viridis')
+plt.gca().set_aspect('equal', 'box')
+#plt.colorbar(orientation='horizontal')
+plt.title('Raw elevation')
+
+plt.subplot(1,3,2)
+plt.pcolormesh(lakes, cmap='viridis')
+plt.gca().set_aspect('equal', 'box')
+#plt.colorbar(orientation='horizontal')
+plt.title('Lake surface elevation')
+
+plt.subplot(1,3,3)
+plt.pcolormesh(np.log(rivers), vmin=0, vmax=np.log(n/25), cmap='Blues')
+plt.gca().set_aspect('equal', 'box')
+#plt.colorbar(orientation='horizontal')
+plt.title('Rivers discharge')
+
+plt.show()