Use iterative finite differences for diffusion, instead of gaussian blur

Allows variable diffusion coefficients

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,41 +1,114 @@
 # Map Generator with Rivers
-`mapgen_rivers v1.0.2` by Gaël de Sailly.
+`mapgen_rivers v0.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).
+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).
 
 Its main particularity compared to conventional Minetest mapgens is that rivers that flow strictly downhill, and combine together to form wider rivers, until they reach the sea. Another notable feature is the possibility of large lakes above sea level.
 
-![Screenshot](https://content.minetest.net/uploads/fff09f2269.png)
+![Screenshot](https://user-images.githubusercontent.com/6905002/98825953-6289d980-2435-11eb-9e0b-704a95663ce0.png)
 
-It used to be composed of a Python script doing pre-generation, and a Lua mod reading the pre-generation output and generating the map. The code has been rewritten in full Lua for version 1.0 (July 2021), and is now usable out-of-the-box as any other Minetest mod.
-
-# Author and license
-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.
+**Important to know**: Unlike most other Minetest mods, it does not contain standalone Lua code, but does part of its processing with a separate Python program (included).
+- The Python part does pre-processing: it creates large-scale terrain data and applies landscape evolution algorithms, then outputs a grid of data in the mod's or world's folder. The grid is typically an array of 1000x1000 points of data, each of them representing a cell (by default 12x12 nodes). This pre-processing is long and should be run in advance.
+- The Lua part does actual map generation on Minetest. It reads grid data, upscales it (by a factor 12 by default), and adds small-scale features.
 
 # 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.
+
+Map pre-generation requires Python 3 with the following libraries installed:
+- `numpy`, widely used library for numerical calculations
+- `scipy`, a library for advanced data treatments, that is used here for Gaussian filtering
+- `noise`, implementing Perlin/Simplex noises
+
+Also, the following are optional (for map preview)
+- `matplotlib`, a famous library for graphical plotting
+- `colorcet` if you absolutely need better colormaps for preview :-)
+
+They are all commonly found on `pip` or `conda` Python distributions.
 
 # Installation
 This mod should be placed in the `mods/` directory of Minetest like any other mod.
 
 # Usage
-It is recommended to use it **only in new worlds, with `singlenode` mapgen**. On first start, it runs pre-generation to produce a grid, from which the map will be generated. This usually takes a few seconds, but depending on custom settings this can grow considerably longer.
+By default, the mod contains a demo 400x400 grid (so you can start the game directly), but it is recommended to run the pre-processing script to generate a new grid before world creation, if you can.
 
-By default, it only generates a 15k x 15k map, centered around the origin. To obtain a bigger map, you can increase grid size and/or block size in settings, but this can be more ressource-intensive (as the map has to be loaded in full at pre-generation).
+1. Run the script `generate.py` to generate a grid, preferentially from inside the mod's directory, but you can also run it directly in a Minetest world. See next paragraph for details about parameters.
+```
+./generate.py
+```
+2. Start Minetest, create a world with `singlenode` mapgen, enable `mapgen_rivers` mod, and launch the game. If you generated a grid in the world directory, it will copy it. If not, it will use the demo grid.
 
-## Settings
-Settings can be found in Minetest in the `Settings` tab, `All settings` -> `Mods` -> `mapgen_rivers`.
+## Parameters for `generate.py`
+For a basic use you do not need to append any argument:
+```
+./generate.py
+```
+By default this will produce a 1000x1000 grid and save it in `river_data/`. Expect a computing time of about 30 minutes.
 
-Most settings are world-specific and a copy is made in `mapgen_rivers.conf` in the world folder, during world first use, which means that further modification of global settings will not alter existing worlds.
+### Parameters and config files
+This pre-processing takes many parameters. Instead of asking all these parameters to the end user, they are grouped in `.conf` files for usability, but the script still allows to override individual settings.
+
+Generic usage:
+```
+./generate.py conf_file output_dir
+```
+
+- `conf_file`: Path to configuration file from which parameters should be read. If omitted, attempts to read in `terrain.conf`.
+- `output_dir`: Directory in which to save the grid data, defaults to `river_data/`. If it does not exist, it is created. If it already contains previous grid data, they are overwritten.
+
+#### Config files
+The mod currently includes 3 config files, providing different terrain styles:
+- `terrain_default.conf` generates the standard terrain, with highest elevations around 250 with sharp peaks, and otherwise hilly terrain.
+- `terrain_higher.conf` generates higher mountains (up to 400 nodes), and wider valleys.
+- `terrain_original.conf` provides a terrain similar to what was generated with the first release of `mapgen_rivers`.
+
+More work is needed to find better and more varied terrain styles.
+
+### Complete list of parameters
+Other parameters can be specified by `--parameter value`. Syntax `--parameter=value` is also supported.
+
+| Parameter     | Description | Example |
+|---------------|-------------|---------|
+| | **Generic parameters** |
+| `mapsize`     | Size of the grid, in number of cells per edge. Usually `1000`, so to have 1000x1000 cells, the grid will have 1001x1001 nodes. Note that the grid is upscaled 12x in the game (this ratio can be changed), so that a `mapsize` of 1000 will result in a 12000x12000 map by default. | `--mapsize 1000` |
+| `sea_level`   | Height of the sea; height below which a point is considered under water even if it is not in a closed depression. | `--sea_level 1` |
+| | **Noise parameters** |
+| `scale`       | Horizontal variation wavlength of the largest noise octave, in grid cells (equivalent to the `spread` of a `PerlinNoise`). | `--scale 400` |
+| `vscale`      | Elevation coefficient, determines the approximate height difference between deepest seas and highest mountains. | `--vscale 300` |
+| `offset`      | Offset of the noise, will determine mean elevation. | `--offset 0` |
+| `persistence` | Relative height of smaller noise octaves compared to bigger ones. | `--persistence 0.6` |
+| `lacunarity`  | Relative reduction of wavelength between octaves. If `lacunarity`×`persistence` is larger than 1 (usual case), smaller octaves result in higher slopes than larger ones. This case is interesting for rivers networks because slopes determine rivers position. | `--lacunarity 2` |
+| | **Landscape evolution parameters**|
+| `K`           | Abstract erosion constant. Increasing it will increase erosive intensity. | `--K 1` |
+| `m`           | Parameter representing the influence of river flux on erosion. For `m=0`, small and big rivers are equal contributors to erosion. For `m=1` the erosive capability is proportional to river flux (assumed to be catchment area). Usual values: `0.25`-`0.60`. Be careful, this parameter is *highly sensitive*. | `--m 0.35` |
+| `d`           | Diffusion coefficient acting on sea/lake floor. Usual values `0`-`1`. | `--d 0.2` |
+| `flex_radius` | Flexure radius. Wavelength over which loss/gain of mass is compensated by uplift/subsidence. This ensures that mountain ranges will not get eventually flattened by erosion, and that an equilibrium is reached. Geologically speaking, this implements [isostatic rebound](https://en.wikipedia.org/wiki/Isostasy). | `--flex_radius 20` |
+| `time`        | Simulated time of erosion modelling, in abstract units. | `--time 10` |
+| `niter`       | Number of iterations. Each iteration represents a time `time/niter`. | `--niter 10` |
+| `sea_level_variations` | Amplitude of sea level variations throughout the simulation (if any). | `--sea_level_variations 10` |
+| `sea_level_variations_time` | Characteristic time of variation for sea level, in the same units than `time`. Increasing it will result in slower variations between iterations. | `--sea_level_variations_time 1` |
+| `flow_method` | Algorithm used for local flow calculation. Possible values are `steepest` (every node flows toward the steepest neighbour when possible), and `semirandom` (default, flow direction is determined randomly between lower neighbours, with lowest ones having greater probability). | `--flow_method semirandom` |
+| | **Alternatives** |
+| `config`      | Another way to specify configuration file | `--config terrain_higher.conf` |
+| `output`      | Another way to specify output dir | `--output ~/.minetest/worlds/my_world/river_data` |
+
+### Example
+```
+./generate.py terrain_higher.conf --mapsize 700 --K 0.4 --m 0.5
+```
+Reads parameters in `terrain_higher.conf`, and will generate a 700x700 grid using custom values for `K` and `m`.
 
 ## Map preview
-The Python script `view_map.py` can display the full map. You need to have Python 3 installed, as well as the libraries `numpy`, `matplotlib`, and optionally `colorcet`. For `conda` users, an `environment.yml` file is provided.
+If you have `matplotlib` installed, `generate.py` will automatically show the grid aspect in real time during the erosion simulation.
 
-It can be run from command line by passing the world folder. Example:
+There is also a script to view a generated map afterwards: `view_map.py`. Its syntax is the following:
 ```
-./view_map.py ~/.minetest/worlds/test_mg_rivers
+./view_map.py grid blocksize
+```
+
+- `grid` is the path to the grid directory to view. For example `river_data/`.
+- `blocksize` is the size at which 1 grid cell will be upscaled, in order to match game coordinates. If you use default settings, use `12`.
+
+Example:
+```
+./view_map.py river_data 12
 ```

compatibility.lua

@@ -1,79 +0,0 @@
--- 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 ""
-
-	local min_catchment = settings:get(prefix.."min_catchment")
-	if min_catchment then
-		min_catchment = tonumber(min_catchment)
-		local blocksize = tonumber(settings:get(prefix.."blocksize") or 15)
-		settings:set(prefix.."min_catchment", tonumber(min_catchment) * blocksize*blocksize)
-		local max_catchment = settings:get(prefix.."max_catchment")
-		if max_catchment then
-			max_catchment = tonumber(max_catchment)
-			local wpower = math.log(2*blocksize)/math.log(max_catchment/min_catchment)
-			settings:set(prefix.."river_widening_power", wpower)
-		end
-	end
-end
-
-local function fix_compatibility_minetest(settings)
-	local previous_version = settings:get("mapgen_rivers_version") or "0.0"
-
-	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)
-	local previous_version = settings:get("version") or "0.0"
-
-	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

demo_data/size

@@ -0,0 +1,2 @@
+401
+401

environment.yml

@@ -1,10 +0,0 @@
-name: mapgen_rivers
-
-channels:
-    - conda-forge
-
-dependencies:
-    - python
-    - matplotlib
-    - numpy
-    - colorcet

generate.py

@@ -0,0 +1,233 @@
+#!/usr/bin/env python3
+
+import numpy as np
+from noise import snoise2, snoise3
+import os
+import sys
+
+import terrainlib
+
+class noisemap:
+    def __init__(self, X, Y, scale=0.01, vscale=1.0, tscale=1.0, offset=0.0, log=None, xbase=None, ybase=None, **params):
+        # Determine noise offset randomly
+        if xbase is None:
+            xbase = np.random.randint(8192)-4096
+        if ybase is None:
+            ybase = np.random.randint(8192)-4096
+        self.xbase = xbase
+        self.ybase = ybase
+        self.X = X
+        self.Y = Y
+        self.scale = scale
+        if log:
+            vscale /= offset
+        self.vscale = vscale
+        self.tscale = tscale
+        self.offset = offset
+        self.log = log
+        self.params = params
+
+    def get2d(self):
+        n = np.zeros((self.X, self.Y))
+        for x in range(self.X):
+            for y in range(self.Y):
+                n[x,y] = snoise2(x/self.scale + self.xbase, y/self.scale + self.ybase, **self.params)
+
+        if self.log:
+            return np.exp(n*self.vscale) * self.offset
+        else:
+            return n*self.vscale + self.offset
+
+    def get3d(self, t=0):
+        t /= self.tscale
+        n = np.zeros((self.X, self.Y))
+        for x in range(self.X):
+            for y in range(self.Y):
+                n[x,y] = snoise3(x/self.scale + self.xbase, y/self.scale + self.ybase, t, **self.params)
+
+        if self.log:
+            return np.exp(n*self.vscale) * self.offset
+        else:
+            return n*self.vscale + self.offset
+
+### PARSE COMMAND-LINE ARGUMENTS
+argc = len(sys.argv)
+
+config_file = 'terrain_default.conf'
+output_dir = 'river_data'
+params_from_args = {}
+i = 1 # Index of arguments
+j = 1 # Number of 'orphan' arguments (the ones that are not preceded by '--something')
+while i < argc:
+    arg = sys.argv[i]
+    if arg[:2] == '--':
+        pname = arg[2:]
+        v = None
+        split = pname.split('=', maxsplit=1)
+        if len(split) == 2:
+            pname, v = split
+            i += 1
+        elif i+1 < argc:
+            v = sys.argv[i+1]
+            i += 2
+
+        if v is not None:
+            if pname == 'config':
+                config_file = v
+            elif pname == 'output':
+                output_dir = v
+            else:
+                params_from_args[pname] = v
+    else:
+        if j == 1:
+            config_file = arg
+        elif j == 2:
+            output_dir = arg
+        i += 1
+        j += 1
+
+print(config_file, output_dir)
+
+params = terrainlib.read_config_file(config_file)
+params.update(params_from_args) # Params given from args prevail against conf file
+
+### READ SETTINGS
+def get_setting(name, default):
+    if name in params:
+        return params[name]
+    return default
+
+mapsize = int(get_setting('mapsize', 1000))
+scale = float(get_setting('scale', 400.0))
+vscale = float(get_setting('vscale', 300.0))
+offset = float(get_setting('offset', 0.0))
+persistence = float(get_setting('persistence', 0.6))
+lacunarity = float(get_setting('lacunarity', 2.0))
+
+K = float(get_setting('K', 0.5))
+m = float(get_setting('m', 0.5))
+d = float(get_setting('d', 0.5))
+sea_level = float(get_setting('sea_level', 0.0))
+sea_level_variations = float(get_setting('sea_level_variations', 0.0))
+sea_level_variations_time = float(get_setting('sea_level_variations_time', 1.0))
+flex_radius = float(get_setting('flex_radius', 20.0))
+tectonics_time = float(get_setting('tectonics_time', 0.0))
+flow_method = get_setting('flow_method', 'semirandom')
+
+time = float(get_setting('time', 10.0))
+niter = int(get_setting('niter', 10))
+
+### MAKE INITIAL TOPOGRAPHY
+n = np.zeros((mapsize+1, mapsize+1))
+
+# Set noise parameters
+params = {
+    "offset" : offset,
+    "vscale" : vscale,
+    "scale" : scale,
+    "octaves" : int(np.ceil(np.log2(mapsize)))+1,
+    "persistence" : persistence,
+    "lacunarity" : lacunarity,
+}
+
+params_sealevel = {
+    "octaves" : 1,
+    "persistence" : 1,
+    "lacunarity" : 2,
+}
+
+params_K = {
+    "offset" : K,
+    "vscale" : K,
+    "scale" : 400,
+    "octaves" : 1,
+    "persistence" : 0.5,
+    "lacunarity" : 2,
+    "log" : True,
+}
+
+params_m = {
+    "offset" : m,
+    "vscale" : m*0.2,
+    "scale" : 400,
+    "octaves" : 1,
+    "persistence" : 0.5,
+    "lacunarity" : 2,
+    "log" : False,
+}
+
+if sea_level_variations != 0.0:
+    sea_ybase = np.random.randint(8192)-4096
+    sea_level_ref = snoise2(time * (1-1/niter) / sea_level_variations, sea_ybase, **params_sealevel) * sea_level_variations
+    params['offset'] -= (sea_level_ref + sea_level)
+
+if tectonics_time == 0.0:
+    n = noisemap(mapsize+1, mapsize+1, **params).get2d()
+else:
+    terrain_noisemap = noisemap(mapsize+1, mapsize+1, tscale=tectonics_time, **params)
+    n = terrain_noisemap.get3d()
+m_map = noisemap(mapsize+1, mapsize+1, **params_m).get2d()
+K_map = noisemap(mapsize+1, mapsize+1, **params_K).get2d()
+
+### COMPUTE LANDSCAPE EVOLUTION
+# Initialize landscape evolution model
+print('Initializing model')
+model = terrainlib.EvolutionModel(n, K=K_map, m=m_map, d=K_map, sea_level=sea_level, flex_radius=flex_radius, flow_method=flow_method)
+terrainlib.update(model.dem, model.lakes, t=5, sea_level=model.sea_level, title='Initializing...')
+
+dt = time/niter
+
+# Run the model's processes: the order in which the processes are run is arbitrary and could be changed.
+
+for i in range(niter):
+    disp_niter = 'Iteration {:d} of {:d}...'.format(i+1, niter)
+    if sea_level_variations != 0:
+        model.sea_level = snoise2((i*dt)/sea_level_variations_time, sea_ybase, **params_sealevel) * sea_level_variations - sea_level_ref
+    terrainlib.update(model.dem, model.lakes, sea_level=model.sea_level, title=disp_niter)
+    print(disp_niter)
+    print('Diffusion')
+    model.diffusion(dt)
+    print('Flow calculation')
+    model.calculate_flow()
+    terrainlib.update(model.dem, model.lakes, sea_level=model.sea_level, title=disp_niter)
+    print('Advection')
+    model.advection(dt)
+    if tectonics_time != 0.0:
+        print('Isostasy reference redefinition')
+        model.define_isostasy(terrain_noisemap.get3d((i+1)*dt))
+    print('Isostatic equilibration')
+    model.adjust_isostasy()
+
+print('Last flow calculation')
+model.calculate_flow()
+
+print('Done!')
+
+# Twist the grid
+bx, by = terrainlib.make_bounds(model.dirs, model.rivers)
+offset_x, offset_y = terrainlib.twist(bx, by, terrainlib.get_fixed(model.dirs))
+
+# Convert offset in 8-bits
+offset_x = np.clip(np.floor(offset_x * 256), -128, 127)
+offset_y = np.clip(np.floor(offset_y * 256), -128, 127)
+
+### SAVE OUTPUT
+if not os.path.isdir(output_dir):
+    os.mkdir(output_dir)
+os.chdir(output_dir)
+# Save the files
+terrainlib.save(model.dem, 'dem', dtype='>i2')
+terrainlib.save(model.lakes, 'lakes', dtype='>i2')
+terrainlib.save(offset_x, 'offset_x', dtype='i1')
+terrainlib.save(offset_y, 'offset_y', dtype='i1')
+
+terrainlib.save(model.dirs, 'dirs', dtype='u1')
+terrainlib.save(model.rivers, 'rivers', dtype='>u4')
+
+with open('size', 'w') as sfile:
+    sfile.write('{:d}\n{:d}'.format(mapsize+1, mapsize+1))
+
+terrainlib.stats(model.dem, model.lakes)
+print()
+print('Grid is ready for use!')
+terrainlib.plot(model.dem, model.lakes, title='Final grid, ready for use!')

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

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,52 +1,13 @@
--- 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 blocksize = mapgen_rivers.blocksize
+local sea_level = mapgen_rivers.sea_level
+local riverbed_slope = mapgen_rivers.riverbed_slope
 
-local out_elev = mapgen_rivers.settings.margin_elev
-
--- 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
+local MAP_BOTTOM = -31000
 
 -- Linear interpolation
 local function interp(v00, v01, v11, v10, xf, zf)
@@ -55,9 +16,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 +31,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 c_NW, c_NE, c_SE, c_SW = unpk(poly.river_corners)
+				local r_west, r_north, r_east, r_south = unpack(poly.rivers)
+				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 +43,65 @@ 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 = math.max(r_west, c_NW-zf, zf-c_SW)
+					local x1 = math.min(r_east, c_NE+zf, c_SE-zf)
+					local z0 = math.max(r_north, c_NW-xf, xf-c_NE)
+					local z1 = math.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
 				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_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
@@ -162,17 +116,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
-					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 = 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
-				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 +134,5 @@ function mapgen_rivers.make_heightmaps(minp, maxp)
 
 	return terrain_height_map, lake_height_map
 end
+
+return heightmaps

init.lua

@@ -1,26 +1,226 @@
--- 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/'
 
 dofile(modpath .. 'settings.lua')
-dofile(modpath .. 'gridmanager.lua')
-dofile(modpath .. 'gridio.lua')
-dofile(modpath .. 'polygons.lua')
-dofile(modpath .. 'heightmap.lua')
-dofile(modpath .. 'mapgen.lua')
 
-minetest.register_on_mods_loaded(function()
-	local exist = mapgen_rivers.try_load_grid()
+local blocksize = mapgen_rivers.blocksize
+local sea_level = mapgen_rivers.sea_level
+local riverbed_slope = mapgen_rivers.riverbed_slope
+local elevation_chill = mapgen_rivers.elevation_chill
+local use_distort = mapgen_rivers.distort
+local use_biomes = mapgen_rivers.biomes
+local use_biomegen_mod = use_biomes and minetest.global_exists('biomegen')
+use_biomes = use_biomes and not use_biomegen_mod
 
-	if not exist then -- If grid does not exist yet, generate it
-		dofile(modpath .. 'pregenerate.lua')
+if use_biomegen_mod then
+	biomegen.set_elevation_chill(elevation_chill)
+end
+dofile(modpath .. 'noises.lua')
 
-		local grid = mapgen_rivers.generate_grid()
-		mapgen_rivers.write_grid(grid)
-		mapgen_rivers.try_load_grid(grid) -- Reload if needed
+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 function generate(minp, maxp, seed)
+	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
-end)
+
+	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()
+end
+
+minetest.register_on_generated(generate)

load.lua

@@ -0,0 +1,31 @@
+local worldpath = minetest.get_worldpath() .. "/river_data/"
+
+local function load_map(filename, bytes, signed, size)
+	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, i1 = (i-1)*bytes+1, i*bytes
+		local elements = {data:byte(i0, i1)}
+		local n = elements[1]
+		if signed and n >= 128 then
+			n = n - 256
+		end
+
+		for j=2, bytes do
+			n = n*256 + elements[j]
+		end
+
+		map[i] = n
+	end
+	file:close()
+
+	return map
+end
+
+return load_map

mapgen.lua

@@ -1,271 +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
-		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
-			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
-
-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,37 @@
+mapgen_rivers.noise_params = {
+	distort_x = {
+		offset = 0,
+		scale = 1,
+		seed = -4574,
+		spread = {x=64, y=32, z=64},
+		octaves = 3,
+		persistence = 0.75,
+		lacunarity = 2,
+	},
+
+	distort_z = {
+		offset = 0,
+		scale = 1,
+		seed = -7940,
+		spread = {x=64, y=32, z=64},
+		octaves = 3,
+		persistence = 0.75,
+		lacunarity = 2,
+	},
+
+	distort_amplitude = {
+		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'),
+}
+
+mapgen_rivers.noise_params.heat.offset = mapgen_rivers.noise_params.heat.offset + mapgen_rivers.sea_level*mapgen_rivers.elevation_chill

polygons.lua

@@ -1,63 +1,102 @@
--- Fetch polygons from a given areas, and compute their properties
--- and find to which polygon every point belongs
+local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
+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 X = math.floor(mapgen_rivers.settings.map_x_size / blocksize)
-local Z = math.floor(mapgen_rivers.settings.map_z_size / blocksize)
+local world_data_path = minetest.get_worldpath() .. '/river_data/'
+minetest.mkdir(world_data_path)
 
+local load_map = dofile(modpath .. 'load.lua')
+
+local function copy_if_needed(filename)
+	local wfilename = world_data_path..filename
+	local wfile = io.open(wfilename, 'rb')
+	if wfile then
+		wfile:close()
+		return
+	end
+	local mfilename = mod_data_path..filename
+	local mfile = io.open(mfilename, 'rb')
+	local wfile = io.open(wfilename, 'wb')
+	wfile:write(mfile:read("*all"))
+	mfile:close()
+	wfile:close()
+end
+
+copy_if_needed('size')
+local sfile = io.open(world_data_path..'size', 'r')
+local X = tonumber(sfile:read('*l'))
+local Z = tonumber(sfile:read('*l'))
+sfile:close()
+
+copy_if_needed('dem')
+local dem = load_map('dem', 2, true, X*Z)
+copy_if_needed('lakes')
+local lakes = load_map('lakes', 2, true, X*Z)
+copy_if_needed('dirs')
+local dirs = load_map('dirs', 1, false, X*Z)
+copy_if_needed('rivers')
+local rivers = load_map('rivers', 4, false, X*Z)
+
+copy_if_needed('offset_x')
+local offset_x = load_map('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('offset_y', 1, true, X*Z)
+for k, v in ipairs(offset_z) do
+	offset_z[k] = (v+0.5)/256
+end
+
+-- To index a flat array representing a 2D map
 local function index(x, z)
 	return z*X+x+1
 end
 
+local blocksize = mapgen_rivers.blocksize
+local min_catchment = mapgen_rivers.min_catchment
+local max_catchment = mapgen_rivers.max_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 mapgen_rivers.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
-local wfactor = 1/(2*blocksize * min_catchment^wpower)
+-- Width coefficients: coefficients solving
+--   wfactor * min_catchment ^ wpower = 1/(2*blocksize)
+--   wfactor * max_catchment ^ wpower = 1
+local wpower = math.log(2*blocksize)/math.log(max_catchment/min_catchment)
+local wfactor = 1 / max_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
-local elevation_chill = mapgen_rivers.settings.elevation_chill
+local elevation_chill = mapgen_rivers.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 = mapgen_rivers.glaciers
+local glacier_factor = mapgen_rivers.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 grid = mapgen_rivers.grid
-	local dem = grid.dem
-	local lakes = grid.lakes
-	local dirs = grid.dirs
-	local rivers = grid.rivers
-
-	local offset_x = grid.offset_x
-	local offset_z = grid.offset_y
+local function make_polygons(minp, maxp)
+	print("Generating polygon map")
+	print(minp.x, maxp.x, minp.z, maxp.z)
 
 	if not init then
 		if glaciers then
@@ -70,8 +109,9 @@ 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 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 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 = 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)
+	print(xpmin, xpmax, zpmin, zpmax)
 
 	-- Iterate over the polygons
 	for xp = xpmin, xpmax do
@@ -93,12 +133,15 @@ function mapgen_rivers.make_polygons(minp, maxp)
 				(offset_z[iC]+zp+1) * blocksize - map_offset.z,
 				(offset_z[iD]+zp+1) * blocksize - map_offset.z,
 			}
+			if xp==xpmin and zp==zpmin then
+				print(xp, zp, poly_x[1], poly_z[1])
+			end
 			local polygon = {x=poly_x, z=poly_z, i={iA, iB, iC, iD}}
 
 			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 zmax = min(floor(max(unpack(poly_z))), maxp.z)
+			local zmin = math.max(math.floor(math.min(unpack(poly_z)))+1, minp.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 +151,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 lzmax = floor(max(z1, z2))
+				local lzmin = math.floor(math.min(z1, z2))+1
+				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 +169,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 xmax = min(floor(zlist[l*2]), maxp.x)
+					local xmin = math.max(math.floor(zlist[l*2-1])+1, minp.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 +195,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_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 +224,5 @@ function mapgen_rivers.make_polygons(minp, maxp)
 
 	return polygons
 end
+
+return make_polygons

pregenerate.lua

@@ -1,116 +0,0 @@
--- Generate the grid using terrainlib_lua
--- Only called on first mapgen, if there is no grid yet
-
-local EvolutionModel = dofile(mapgen_rivers.modpath .. '/terrainlib_lua/erosion.lua')
-local twist = dofile(mapgen_rivers.modpath .. '/terrainlib_lua/twist.lua')
-
-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
-
-local time = mapgen_rivers.settings.evol_time
-local time_step = mapgen_rivers.settings.evol_time_step
-local niter = math.ceil(time/time_step)
-time_step = time / niter
-
-local use_margin = mapgen_rivers.settings.margin
-local margin_width = mapgen_rivers.settings.margin_width / blocksize
-local margin_elev = mapgen_rivers.settings.margin_elev
-
-local function margin(dem, width, elev)
-	local X, Y = dem.X, dem.Y
-	for i=1, width do
-		local c1 = ((i-1)/width) ^ 0.5
-		local c2 = (1-c1) * elev
-		local index = (i-1)*X + 1
-		for x=1, X do
-			dem[index] = dem[index] * c1 + c2
-			index = index + 1
-		end
-		index = i
-		for y=1, Y do
-			dem[index] = dem[index] * c1 + c2
-			index = index + X
-		end
-		index = X*(Y-i) + 1
-		for x=1, X do
-			dem[index] = dem[index] * c1 + c2
-			index = index + 1
-		end
-		index = X-i + 1
-		for y=1, Y do
-			dem[index] = dem[index] * c1 + c2
-			index = index + X
-		end
-	end
-end
-
-function mapgen_rivers.pregenerate(grid)
-	local size = grid.size
-
-	if size.x * size.y > 4e6 then
-		minetest.log("warning", "[mapgen_rivers] You are going to generate a very large grid (>4M nodes). If you experience problems, you should increase blocksize or reduce map size.")
-	end
-
-	local seed = tonumber(minetest.get_mapgen_setting("seed"))
-	np_base.seed = (np_base.seed or 0) + seed
-
-	local nobj_base = PerlinNoiseMap(np_base, {x=size.x, y=1, z=size.y})
-
-	local dem = nobj_base:get_3d_map_flat({x=0, y=0, z=0})
-	dem.X = size.x
-	dem.Y = size.y
-
-	if use_margin then
-		margin(dem, margin_width, margin_elev)
-	end
-
-	local model = EvolutionModel(evol_params)
-	model.dem = dem
-	local ref_dem = model:define_isostasy(dem)
-
-	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: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
-
-	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

readconfig.py

@@ -1,42 +0,0 @@
-def read_conf_file(filename):
-    f = open(filename, 'r')
-    return read_conf(f)
-
-def read_conf(f, end_tag=None):
-    conf = {}
-    while True:
-        line = f.readline()
-        if len(line) == 0:
-            return conf
-        line = line.strip()
-        if line == end_tag:
-            return conf
-        if len(line) == 0 or line[0] == '#':
-            continue
-
-        eqpos = line.find('=')
-        if eqpos < 0:
-            continue
-            
-        name, value = line[:eqpos].rstrip(), line[eqpos+1:].lstrip()
-        if value == '{':
-            # Group
-            conf[name] = read_conf(f, end_tag='}')
-
-        elif value == '"""':
-            # Multiline
-            conf[value] = read_multiline(f)
-
-        else:
-            conf[name] = value
-
-def read_multiline(f):
-    mline = ''
-    while True:
-        line = f.readline()
-        if len(line) == 0:
-            return mline
-        line = line.strip()
-        if line == '"""':
-            return mline
-        mline += line + '\n'

settings.lua

@@ -1,137 +1,58 @@
--- Read global and per-world settings
+local storage = minetest.get_mod_storage()
+local settings = minetest.settings
 
-local mtsettings = minetest.settings
-local mgrsettings = Settings(minetest.get_worldpath() .. '/mapgen_rivers.conf')
-
-mapgen_rivers.version = "1.0.2-dev1"
-
-local previous_version_mt = mtsettings:get("mapgen_rivers_version") or "0.0"
-local previous_version_mgr = mgrsettings: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")
-	if mapgen_rivers.version ~= previous_version_mt then
-		compat_mt(mtsettings)
-	end
-	if mapgen_rivers.version ~= previous_version_mgr then
-		compat_mgr(mgrsettings)
-	end
-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
+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)
+		elseif dtype == "bool" then
+			return storage:get_string(key) == 'true'
 		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]
-			end
-			mgrsettings:set(name, v)
-		end
-		if dtype == "number" then
-			return tonumber(v)
+	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)
 		else
-			return v
-		end
-	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]
+			storage:set_string(key, conf_val)
+			if dtype == "bool" then
+				conf_val = conf_val == 'true'
 			end
-			mgrsettings:set_bool(name, v)
 		end
-		return 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)
+
+		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)
+		elseif dtype == "bool" then
+			storage:set_string(key, tostring(default))
 		end
-		clean_np(v)
-		return v
+
+		return default
 	end
 end
 
-local def_setting = mapgen_rivers.define_setting
-
-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'),
-
-	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()
-end
-
-minetest.register_on_mods_loaded(write_settings)
-minetest.register_on_shutdown(write_settings)
+mapgen_rivers.center = get_settings('center', 'bool', false)
+mapgen_rivers.blocksize = get_settings('blocksize', 'int', 12)
+mapgen_rivers.sea_level = get_settings('sea_level', 'int', 1)
+mapgen_rivers.min_catchment = get_settings('min_catchment', 'float', 25)
+mapgen_rivers.max_catchment = get_settings('max_catchment', 'float', 40000)
+mapgen_rivers.riverbed_slope = get_settings('riverbed_slope', 'float', 0.4) * mapgen_rivers.blocksize
+mapgen_rivers.distort = get_settings('distort', 'bool', true)
+mapgen_rivers.biomes = get_settings('biomes', 'bool', true)
+mapgen_rivers.glaciers = get_settings('glaciers', 'bool', false)
+mapgen_rivers.glacier_factor = get_settings('glacier_factor', 'float', 8)
+mapgen_rivers.elevation_chill = get_settings('elevation_chill', 'float', 0.25)

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

@@ -1,47 +1,26 @@
 # File containing all settings for 'mapgen_rivers' mod.
 
 #    Whether the map should be centered at x=0, z=0.
-mapgen_rivers_center (Center map) bool true
+mapgen_rivers_center (Center map) bool false
 
-#    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
-mapgen_rivers_blocksize (Block size) float 15.0 2.0 100.0
+#    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 15000.
+mapgen_rivers_blocksize (Block size) float 12.0 2.0 40.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
+#    Sea level used by mapgen_rivers
+mapgen_rivers_sea_level (Sea level) int 1
 
-#    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
-
-#    Minimal catchment area for a river to be drawn, in square nodes
+#    Minimal catchment area for a river to be drawn, in grid cells
+#    (1 cell = blocksize x blocksize).
 #    Lower value means bigger river density
-mapgen_rivers_min_catchment (Minimal catchment area) float 3600.0 100.0 1000000.0
+mapgen_rivers_min_catchment (Minimal catchment area) float 25.0 1.0 1000.0
 
-#    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.
-#    Note that a river can never exceed 2*blocksize.
-mapgen_rivers_river_widening_power (River widening power) float 0.5 0.0 1.0
+#    Catchment area in grid cells (1 grid cell = blocksize x blocksize)
+#    at which rivers reach their maximal width of 2*blocksize.
+#    Higher value means a river needs to receive more tributaries to grow in width.
+mapgen_rivers_max_catchment (Maximal catchment area) float 40000.0 1000.0 10000000.0
 
 #    Lateral slope of the riverbed.
 #    Higher value means deeper rivers.
@@ -70,64 +49,4 @@ 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
-mapgen_rivers_evol_time (Landscape evolution time) float 10.0 0.0 100.0
-
-#    Model time steps in arbitrary units
-#    Smaller values will result in more time steps to be necessary to
-#    complete the simulation, taking more time.
-mapgen_rivers_evol_time_step (Landscape evolution time step) float 1.0 0.0 50.0
-
-#    To adjust river erosion proportionnally.
-#    This type of erosion acts by deepening the valleys.
-mapgen_rivers_river_erosion_coef (River erosion coefficient) float 0.5 0.0 10.0
-
-#    Represents how much river erosion depends on river flow (catchment area).
-#    Catchment area is elevated to this power.
-#    Extreme cases: 0.0 -> All rivers have the same erosive capabilities
-#    1.0 -> Erosion is proportional to river flow
-#    Reasonable values are generally between 0.4 and 0.7.
-#
-#    This parameter is extremely sensitive, and changes may require to adjust
-#    'river_erosion_coef' as well.
-mapgen_rivers_river_erosion_power (River erosion power) float 0.4 0.0 1.0
-
-#    Intensity of diffusive erosion.
-#    Smoothes peaks and valleys, and tends to prevent sharp cliffs from forming.
-mapgen_rivers_diffusive_erosion (Diffusive erosion) float 0.5 0.0 10.0
-
-#    Radius of compensation for isostatic/tectonic processes
-#    Tectonic uplift forces will have a diffuse effect over this radius
-mapgen_rivers_compensation_radius (Tectonic compensation radius) float 50 1.0 1000.0
-
-#    Speed of evolution of tectonic conditions between steps
-#    Higher values means tectonics will be very different from one step to the other,
-#    resulting in geologically unstable and more varied landforms (plateau, gorge, lake...)
-mapgen_rivers_tectonic_speed (Tectonic speed) float 70 0 10000
-
-[Noises]
-
-#    Y level of terrain at a very large scale. Only used during pre-generation.
-#    X and Z axes correspond to map's X and Z directions, and Y axis is time.
-#    Successive XZ slices of this noise represent successive tectonic states.
-mapgen_rivers_np_base (Terrain base noise) noise_params_3d 0, 300, (2048, 2048, 2048), 2469, 8, 0.6, 2.0, eased
-
-#    This noise will shear the terrain on the X axis,
-#    to break the regularity of the river grid.
-mapgen_rivers_np_distort_x (X-axis distorsion noise) noise_params_3d 0, 1, (64, 32, 64), -4574, 3, 0.75, 2.0
-
-#    This noise will shear the terrain on the Z axis,
-#    to break the regularity of the river grid.
-mapgen_rivers_np_distort_z (Z-axis distorsion noise) noise_params_3d 0, 1, (64, 32, 64), -7940, 3, 0.75, 2.0
-
-#    Amplitude of the distorsion.
-#    Too small values may leave the grid pattern apparent,
-#    and too high values could make the terrain insanely twisted.
-mapgen_rivers_np_distort_amplitude (Distorsion amplitude noise) noise_params_2d 0, 10, (1024, 1024, 1024), 676, 5, 0.5, 2.0, absvalue
+# Noises: to be added. For now they are hardcoded.

terrain_default.conf

@@ -0,0 +1,17 @@
+mapsize = 1000
+scale = 400
+vscale = 300
+offset = 0
+persistence = 0.6
+lacunarity = 2.0
+
+K = 0.5
+m = 0.5
+d = 0.5
+sea_level = 0
+sea_level_variations = 8
+sea_level_variations_time = 2
+flex_radius = 20
+
+time = 10
+niter = 10

terrain_higher.conf

@@ -0,0 +1,17 @@
+mapsize = 1000
+scale = 400
+vscale = 600
+offset = 0
+persistence = 0.65
+lacunarity = 2.0
+
+K = 0.5
+m = 0.45
+d = 0.55
+sea_level = 0
+sea_level_variations = 12
+sea_level_variations_time = 2
+flex_radius = 50
+
+time = 10
+niter = 10

terrain_original.conf

@@ -0,0 +1,16 @@
+mapsize = 1000
+scale = 400
+vscale = 300
+offset = 0
+persistence = 0.6
+lacunarity = 2.0
+flow_method = steepest
+
+K = 1
+m = 0.35
+d = 0
+sea_level = 0
+flex_radius = 20
+
+time = 10
+niter = 10

terrainlib/__init__.py

@@ -0,0 +1,7 @@
+# Load packages and provide easy access to important functions
+
+from .settings import read_config_file
+from .erosion import EvolutionModel
+from .save import save
+from .bounds import make_bounds, twist, get_fixed
+from .view import stats, update, plot

terrainlib/bounds.py

@@ -0,0 +1,74 @@
+import numpy as np
+
+def make_bounds(dirs, rivers):
+    """
+    Give an array of all horizontal and vertical bounds
+    """
+
+    (Y, X) = dirs.shape
+    bounds_h = np.zeros((Y, X-1), dtype=rivers.dtype)
+    bounds_v = np.zeros((Y-1, X), dtype=rivers.dtype)
+
+    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):
+    """
+    Give the list of points that should not be twisted
+    """
+
+    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):
+    """
+    Twist the grid (define an offset for every node). Model river bounds as if they were elastics.
+    Smoothes preferentially big rivers.
+    """
+
+    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 # Normalize, take into account the direction only
+        offset_x += force_x * coeff
+        offset_y += force_y * coeff
+
+    return offset_x, offset_y

terrainlib/erosion.py

@@ -0,0 +1,121 @@
+import numpy as np
+import scipy.ndimage as im
+import scipy.signal as si
+from .rivermapper import flow
+
+def advection(dem, dirs, rivers, time, K=1, m=0.5, sea_level=0):
+    """
+    Simulate erosion by rivers.
+    This models erosion as an upstream advection of elevations ("erosion waves").
+    Advection speed depends on water flux and parameters:
+
+    v = K * flux^m
+    """
+
+    adv_time = 1 / (K*rivers**m) # For every pixel, calculate the time an "erosion wave" will need to cross it.
+    dem = np.maximum(dem, sea_level)
+    dem_new = np.zeros(dem.shape)
+
+    for y in range(dirs.shape[0]):
+        for x in range(dirs.shape[1]):
+            # Elevations propagate upstream, so for every pixel we seek the downstream pixel whose erosion wave just reached the current pixel.
+            # This means summing the advection times downstream until we reach the erosion time.
+            x0, y0 = x, y
+            x1, y1 = x, y
+            remaining = time
+            while True:
+                # Move one pixel downstream
+                flow_dir = dirs[y0,x0]
+                if flow_dir == 0:
+                    remaining = 0
+                    break
+                elif flow_dir == 1:
+                    y1 += 1
+                elif flow_dir == 2:
+                    x1 += 1
+                elif flow_dir == 3:
+                    y1 -= 1
+                elif flow_dir == 4:
+                    x1 -= 1
+
+                if remaining <= adv_time[y0,x0]: # Time is over, we found it.
+                    break
+                remaining -= adv_time[y0,x0]
+                x0, y0 = x1, y1
+
+            c = remaining / adv_time[y0,x0]
+            dem_new[y,x] = c*dem[y1,x1] + (1-c)*dem[y0,x0] # If between 2 pixels, perform linear interpolation.
+
+    return dem_new
+
+second_derivative_matrix = np.array([
+    [0., 0.25, 0.],
+    [0.25,-1., 0.25],
+    [0., 0.25, 0.],
+])
+
+diff_max = 1.0
+
+def diffusion(dem, time, d=1.0):
+    if isinstance(d, np.ndarray):
+        dmax = d.max()
+    else:
+        dmax = d
+    diff = time * dmax
+    print(diff)
+    niter = int(diff//diff_max) + 1
+    ddiff = d * (time / niter)
+
+    #print('{:d} iterations'.format(niter))
+    for i in range(niter):
+        dem[1:-1,1:-1] += si.convolve2d(dem, second_derivative_matrix, mode='valid') * ddiff
+        #print('iteration {:d}'.format(i+1))
+
+    return dem
+    #return im.gaussian_filter(dem, radius, mode='reflect') # Diffusive erosion is a simple Gaussian blur
+
+class EvolutionModel:
+    def __init__(self, dem, K=1, m=0.5, d=1, sea_level=0, flow=False, flex_radius=100, flow_method='semirandom'):
+        self.dem = dem
+        #self.bedrock = dem
+        self.K = K
+        self.m = m
+        if isinstance(d, np.ndarray):
+            self.d = d[1:-1,1:-1]
+        else:
+            self.d = d
+        self.sea_level = sea_level
+        self.flex_radius = flex_radius
+        self.define_isostasy()
+        self.flow_method = flow_method
+        #set_flow_method(flow_method)
+        if flow:
+            self.calculate_flow()
+        else:
+            self.lakes = dem
+            self.dirs = np.zeros(dem.shape, dtype=int)
+            self.rivers = np.zeros(dem.shape, dtype=int)
+            self.flow_uptodate = False
+
+    def calculate_flow(self):
+        self.dirs, self.lakes, self.rivers = flow(self.dem, method=self.flow_method)
+        self.flow_uptodate = True
+
+    def advection(self, time):
+        dem = advection(np.maximum(self.dem, self.lakes), self.dirs, self.rivers, time, K=self.K, m=self.m, sea_level=self.sea_level)
+        self.dem = np.minimum(dem, self.dem)
+        self.flow_uptodate = False
+
+    def diffusion(self, time):
+        self.dem = diffusion(self.dem, time, d=self.d)
+        self.flow_uptodate = False
+
+    def define_isostasy(self, dem=None):
+        if dem is None:
+            dem = self.dem
+        self.ref_isostasy = im.gaussian_filter(dem, self.flex_radius, mode='reflect') # Define a blurred version of the DEM that will be considered as the reference isostatic elevation.
+
+    def adjust_isostasy(self, rate=1):
+        isostasy = im.gaussian_filter(self.dem, self.flex_radius, mode='reflect') # Calculate blurred DEM
+        correction = (self.ref_isostasy - isostasy) * rate # Compare it with the reference isostasy
+        self.dem = self.dem + correction # Adjust

terrainlib/rivermapper.py

@@ -0,0 +1,278 @@
+import numpy as np
+import numpy.random as npr
+from collections import defaultdict
+
+# This file provide functions to construct the river tree from an elevation model.
+# Based on a research paper:
+#   | Cordonnier, G., Bovy, B., and Braun, J.:
+#   | A versatile, linear complexity algorithm for flow routing in topographies with depressions,
+#   | Earth Surf. Dynam., 7, 549–562, https://doi.org/10.5194/esurf-7-549-2019, 2019.
+# 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 Boruvka algorithm.
+# Only flow_local and accumulate_flow are custom algorithms.
+
+# Define two different method for local flow routing
+def flow_local_steepest(plist):
+    vmax = 0.0
+    imax = 0.0
+    for i, p in enumerate(plist):
+        if p > vmax:
+            vmax = p
+            imax = i
+    if vmax > 0.0:
+        return imax+1
+    return 0
+
+def flow_local_semirandom(plist):
+    """
+    Determines a flow direction based on denivellation for every neighbouring node.
+    Denivellation must be positive for downward and zero for flat or upward:
+    dz = max(zref-z, 0)
+    """
+    psum = sum(plist)
+    if psum == 0:
+        return 0
+    r = npr.random() * psum
+    for i, p in enumerate(plist):
+        if r < p:
+            return i+1
+        r -= p
+
+flow_local_methods = {
+    'steepest' : flow_local_steepest,
+    'semirandom' : flow_local_semirandom,
+}
+
+def flow(dem, method='semirandom'):
+    if method in flow_local_methods:
+        flow_local = flow_local_methods[method]
+    else:
+        raise KeyError('Flow method \'{}\' does not exist'.format(method))
+
+    # Flow locally
+    dirs1 = np.zeros(dem.shape, dtype=int)
+    dirs2 = np.zeros(dem.shape, dtype=int)
+    (X, Y) = dem.shape
+    Xmax, Ymax = X-1, Y-1
+    singular = []
+    for x in range(X):
+        z0 = z1 = z2 = dem[x,0]
+        for y in range(Y):
+            z0 = z1
+            z1 = z2
+            if y < Ymax:
+                z2 = dem[x, y+1]
+
+            plist = [
+                max(z1-dem[x+1,y],0) if x<Xmax else 0, # 1: x -> x+1
+                max(z1-z2,0),                          # 2: y -> y+1
+                max(z1-dem[x-1,y],0) if x>0    else 0, # 3: x -> x-1
+                max(z1-z0,0),                          # 4: y -> y-1
+            ]
+            
+            pdir = flow_local(plist)
+            dirs2[x,y] = pdir
+            if pdir == 0:
+                singular.append((x,y))
+            elif pdir == 1:
+                dirs1[x+1,y] += 1
+            elif pdir == 2:
+                dirs1[x,y+1] += 2
+            elif pdir == 3:
+                dirs1[x-1,y] += 4
+            elif pdir == 4:
+                dirs1[x,y-1] += 8
+
+    # Compute basins
+    basin_id = np.zeros(dem.shape, dtype=int)
+    stack = []
+
+    for i, s in enumerate(singular):
+        queue = [s]
+        while queue:
+            x, y = queue.pop()
+            basin_id[x,y] = i
+            d = int(dirs1[x,y])
+
+            if d & 1:
+                queue.append((x-1,y))
+            if d & 2:
+                queue.append((x,y-1))
+            if d & 4:
+                queue.append((x+1,y))
+            if d & 8:
+                queue.append((x,y+1))
+
+    del dirs1
+
+    # Link basins
+    nsing = len(singular)
+    links = {}
+    def add_link(b0, b1, elev, bound):
+        b = (min(b0,b1),max(b0,b1))
+        if b not in links or links[b][0] > elev:
+            links[b] = (elev, bound)
+
+    for x in range(X):
+        b0 = basin_id[x,0]
+        add_link(-1, b0, dem[x,0], (True, x, 0))
+        for y in range(1,Y):
+            b1 = basin_id[x,y]
+            if b0 != b1:
+                add_link(b0, b1, max(dem[x,y-1],dem[x,y]), (True, x, y))
+            b0 = b1
+        add_link(-1, b1, dem[x,Ymax], (True, x, Y))
+    for y in range(Y):
+        b0 = basin_id[0,y]
+        add_link(-1, b0, dem[0,y], (False, 0, y))
+        for x in range(1,X):
+            b1 = basin_id[x,y]
+            if b0 != b1:
+                add_link(b0, b1, max(dem[x-1,y],dem[x,y]), (False, x, y))
+            b0 = b1
+        add_link(-1, b1, dem[Xmax,y], (False, X, y))
+
+    # Computing basin tree
+    graph = planar_boruvka(links)
+
+    basin_links = defaultdict(dict)
+    for elev, b1, b2, bound in graph:
+        basin_links[b1][b2] = basin_links[b2][b1] = (elev, bound)
+    basins = np.zeros(nsing+1)
+    stack = [(-1, float('-inf'))]
+
+    # Applying basin flowing
+    dir_reverse = (0, 3, 4, 1, 2)
+    while stack:
+        b1, elev1 = stack.pop()
+        basins[b1] = elev1
+
+        for b2, (elev2, bound) in basin_links[b1].items():
+            stack.append((b2, max(elev1, elev2)))
+
+            # Reverse flow direction in b2 (TODO)
+            isY, x, y = bound
+            backward = True # Whether water will escape the basin in +X/+Y direction
+            if not (x < X and y < Y and basin_id[x,y] == b2):
+                if isY:
+                    y -= 1
+                else:
+                    x -= 1
+                backward = False
+            d = 2*backward + isY + 1
+            while d > 0:
+                d, dirs2[x,y] = dirs2[x,y], d
+                if d == 1:
+                    x += 1
+                elif d == 2:
+                    y += 1
+                elif d == 3:
+                    x -= 1
+                elif d == 4:
+                    y -= 1
+                d = dir_reverse[d]
+
+            del basin_links[b2][b1]
+        del basin_links[b1]
+
+    # Calculating water quantity
+    dirs2[-1,:][dirs2[-1,:]==1] = 0
+    dirs2[:,-1][dirs2[:,-1]==2] = 0
+    dirs2[0,:][dirs2[0,:]==3] = 0
+    dirs2[:,0][dirs2[:,0]==4] = 0
+
+    waterq = accumulate_flow(dirs2)
+
+    return dirs2, basins[basin_id], waterq
+
+def accumulate_flow(dirs):
+    ndonors = np.zeros(dirs.shape, dtype=int)
+    ndonors[1:,:] += dirs[:-1,:] == 1
+    ndonors[:,1:] += dirs[:,:-1] == 2
+    ndonors[:-1,:] += dirs[1:,:] == 3
+    ndonors[:,:-1] += dirs[:,1:] == 4
+    waterq = np.ones(dirs.shape, dtype=int)
+
+    (X, Y) = dirs.shape
+    rangeX = range(X)
+    rangeY = range(Y)
+    for x in rangeX:
+        for y in rangeY:
+            if ndonors[x,y] > 0:
+                continue
+            xw, yw = x, y
+            w = waterq[xw, yw]
+            while 1:
+                d = dirs[xw, yw]
+                if d <= 0:
+                    break
+                elif d == 1:
+                    xw += 1
+                elif d == 2:
+                    yw += 1
+                elif d == 3:
+                    xw -= 1
+                elif d == 4:
+                    yw -= 1
+
+                w += waterq[xw, yw]
+                waterq[xw, yw] = w
+
+                if ndonors[xw, yw] > 1:
+                    ndonors[xw, yw] -= 1
+                    break
+
+    return waterq
+
+def planar_boruvka(links):
+    # Compute basin tree
+
+    basin_list = defaultdict(dict)
+
+    for (b1, b2), (elev, bound) in links.items():
+        basin_list[b1][b2] = basin_list[b2][b1] = (elev, b1, b2, bound)
+
+    threshold = 8
+    lowlevel = {}
+    for k, v in basin_list.items():
+        if len(v) <= threshold:
+            lowlevel[k] = v
+
+    basin_graph = []
+    n = len(basin_list)
+    while n > 1:
+        (b1, lnk1) = lowlevel.popitem()
+        b2 = min(lnk1, key=lnk1.get)
+        lnk2 = basin_list[b2]
+
+        # Add link to the graph
+        basin_graph.append(lnk1[b2])
+
+        # Union : merge basin 1 into basin 2
+        # First, delete the direct link
+        del lnk1[b2]
+        del lnk2[b1]
+
+        # Look for basin 1's neighbours, and add them to basin 2 if they have a lower pass
+        for k, v in lnk1.items():
+            bk = basin_list[k]
+            if k in lnk2 and lnk2[k] < v:
+                del bk[b1]
+            else:
+                lnk2[k] = v
+                bk[b2] = bk.pop(b1)
+
+            if k not in lowlevel and len(bk) <= threshold:
+                lowlevel[k] = bk
+
+        if b2 in lowlevel:
+            if len(lnk2) > threshold:
+                del lowlevel[b2]
+        elif len(lnk2) <= threshold:
+            lowlevel[b2] = lnk2
+        del lnk1
+
+        n -= 1
+
+    return basin_graph

terrainlib/save.py

@@ -0,0 +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(bin_data)

terrainlib/settings.py

@@ -0,0 +1,16 @@
+import os.path
+
+def read_config_file(fname):
+    settings = {}
+
+    if not os.path.isfile(fname):
+        return settings
+
+    with open(fname, 'r') as f:
+        for line in f:
+            slist = line.split('=', 1)
+            if len(slist) >= 2:
+                prefix, suffix = slist
+                settings[prefix.strip()] = suffix.strip()
+
+    return settings

terrainlib/view.py

@@ -19,7 +19,7 @@ except ImportError: # No module matplotlib
     has_matplotlib = False
 
 if has_matplotlib:
-    def view_map(dem, lakes, scale=1, center=False, sea_level=0.0, title=None):
+    def view_map(dem, lakes, scale=1, sea_level=0.0, title=None):
         lakes_sea = np.maximum(lakes, sea_level)
         water = np.maximum(lakes_sea - dem, 0)
         max_elev = dem.max()
@@ -31,10 +31,7 @@ if has_matplotlib:
         rgb = ls.shade(dem, cmap=cmap1, vert_exag=1/scale, blend_mode='soft', norm=norm_ground)
 
         (X, Y) = dem.shape
-        if center:
-            extent = (-(Y+1)*scale/2, (Y-1)*scale/2, -(X+1)*scale/2, (X-1)*scale/2)
-        else:
-            extent = (-0.5*scale, (Y-0.5)*scale, -0.5*scale, (X-0.5)*scale)
+        extent = (0, Y*scale, 0, X*scale)
         plt.imshow(np.flipud(rgb), extent=extent, interpolation='antialiased')
         alpha = (water > 0).astype('u1')
         plt.imshow(np.flipud(water), alpha=np.flipud(alpha), cmap=cmap2, extent=extent, vmin=0, vmax=max_depth, interpolation='antialiased')

terrainlib_lua/erosion.lua

@@ -1,244 +0,0 @@
--- erosion.lua
-
--- This is the main file of terrainlib_lua. It registers the EvolutionModel object and some of the 
-
-local function erode(model, time)
-	-- Apply river erosion on the model
-	-- Erosion model is based on the simplified version of the stream-power law Ey = K×A^m×S
-	-- where Ey is the vertical erosion speed, A catchment area of the river, S slope along the river, m and K local constants.
-	-- It is equivalent to considering a horizontal erosion wave travelling at Ex = K×A^m, and this latter approach allows much greather time steps so it is used here.
-	-- For each point, instead of moving upstream and see what point the erosion wave would reach, we move downstream and see from which point the erosion wave would reach the given point, then we can set the elevation.
-	local mmin, mmax = math.min, math.max
-	local dem = model.dem
-	local dirs = model.dirs
-	local lakes = model.lakes
-	local rivers = model.rivers
-	local sea_level = model.params.sea_level
-	local K = model.params.K
-	local m = model.params.m
-	local X, Y = dem.X, dem.Y
-	local scalars = type(K) == "number" and type(m) == "number"
-
-	local erosion_time
-	if model.params.variable_erosion then
-		erosion_time = {}
-	else
-		erosion_time = model.erosion_time or {}
-	end
-
-	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.
-			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.
-		end
-	else
-		for i=1, X*Y do
-			local etime = 1 / (K[i]*rivers[i]^m[i])
-			erosion_time[i] = etime
-			lakes[i] = mmax(lakes[i], dem[i], sea_level)
-		end
-	end
-
-	for i=1, X*Y do
-		local iw = i
-		local remaining = time
-		local new_elev
-		while true do
-			-- Explore downstream until we find the point 'iw' from which the erosion wave will reach 'i'
-			local inext = iw
-			local d = dirs[iw]
-
-			-- Follow the river downstream (move 'iw')
-			if d == 0 then -- If no flow direction, we reach the border of the map: set elevation to the latest node's elev and abort.
-				new_elev = lakes[iw]
-				break
-			elseif d == 1 then
-				inext = iw+X
-			elseif d == 2 then
-				inext = iw+1
-			elseif d == 3 then
-				inext = iw-X
-			elseif d == 4 then
-				inext = iw-1
-			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'.
-				local c = remaining / etime
-				new_elev = (1-c) * lakes[iw] + c * lakes[inext] -- Interpolate linearly between the two nodes
-				break
-			end
-
-			remaining = remaining - etime -- If we still don't reach the target time, decrement time and move to next point.
-			iw = inext
-		end
-
-		dem[i] = mmin(dem[i], new_elev)
-	end
-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
-		for i=1, X*Y do
-			dmax = mmax(dmax, d[i])
-		end
-	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
-
-	local temp = {}
-	for n=1, niter do
-		local i = 1
-		for y=1, Y do
-			local iN = (y==1) and 0 or -X
-			local iS = (y==Y) and 0 or X
-			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
-		end
-
-		for i=1, X*Y do
-			dem[i] = dem[i] + temp[i]*ddiff
-		end
-	end
-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
-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.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
-	if type(uplift_rate) == "number" then
-		local uplift_total = uplift_rate * time
-		for i=1, X*Y do
-			dem[i] = dem[i] + uplift_total
-		end
-	else
-		for i=1, X*Y do
-			dem[i] = dem[i] + uplift_rate[i]*time
-		end
-	end
-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
-	local X, Y = dem.X, dem.Y
-	for i=1, X*Y do
-		dem[i] = dem[i] + (random()-0.5) * noise_depth
-	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
-		model.isostasy_ref = ref
-		return
-	end
-
-	local X, Y = ref.X, ref.Y
-	local ref2 = model.isostasy_ref or {X=X, Y=Y}
-	model.isostasy_ref = ref2
-	for i=1, X*Y do
-		ref2[i] = ref[i]
-	end
-
-	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
-	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
-	end
-end
-
-local evol_model_mt = {
-	erode = erode,
-	diffuse = diffuse,
-	flow = flow,
-	uplift = uplift,
-	noise = noise,
-	isostasy = isostasy,
-	define_isostasy = define_isostasy,
-}
-
-evol_model_mt.__index = evol_model_mt
-
-local defaults = {
-	K = 1,
-	m = 0.5,
-	d = 1,
-	variable_erosion = false,
-	sea_level = 0,
-	uplift = 10,
-	noise = 0.001,
-	compensation_radius = 50,
-}
-
-local function EvolutionModel(params)
-	params = params or {}
-	local o = {params = params}
-	for k, v in pairs(defaults) do
-		if params[k] == nil then
-			params[k] = v
-		end
-	end
-	o.dem = params.dem
-	return setmetatable(o, evol_model_mt)
-end
-
-return EvolutionModel

terrainlib_lua/gaussian.lua

@@ -1,88 +0,0 @@
--- gaussian.lua
-
-local function get_box_size(sigma, n)
-	local v = sigma^2 / n
-	local r_ideal = ((12*v + 1) ^ 0.5 - 1) / 2
-	local r_down = math.floor(r_ideal)
-	local r_up = math.ceil(r_ideal)
-	local v_down = ((2*r_down+1)^2 - 1) / 12
-	local v_up = ((2*r_up+1)^2 - 1) / 12
-	local m_ideal = (v - v_down) / (v_up - v_down) * n
-	local m = math.floor(m_ideal+0.5)
-
-	local sizes = {}
-	for i=1, n do
-		sizes[i] = i<=m and 2*r_up+1 or 2*r_down+1
-	end
-
-	return sizes
-end
-
-local function box_blur_1d(map, size, first, incr, len, map2)
-	local n = math.ceil(size/2)
-	first = first or 1
-	incr = incr or 1
-	len = len or math.floor((#map-first)/incr)+1
-	local last = first + (len-1)*incr
-
-	local nth = first+(n-1)*incr
-	local sum = 0
-	for i=first, nth, incr do
-		if i == first then
-			sum = sum + map[i]
-		else
-			sum = sum + 2*map[i]
-		end
-	end
-
-	local i_left = nth
-	local incr_left = -incr
-	local i_right = nth
-	local incr_right = incr
-
-	map2 = map2 or {}
-	for i=first, last, incr do
-		map2[i] = sum / size
-		i_right = i_right + incr_right
-		sum = sum - map[i_left] + map[i_right]
-		i_left = i_left + incr_left
-
-		if i_left == first then
-			incr_left = incr
-		end
-		if i_right == last then
-			incr_right = -incr
-		end
-	end
-
-	return map2
-end
-
-local function box_blur_2d(map1, size, map2)
-	local X, Y = map1.X, map1.Y
-	map2 = map2 or {}
-	for y=1, Y do
-		box_blur_1d(map1, size, (y-1)*X+1, 1, X, map2)
-	end
-	for x=1, X do
-		box_blur_1d(map2, size, x, X, Y, map1)
-	end
-
-	return map1
-end
-
-local function gaussian_blur_approx(map, sigma, n, map2)
-	map2 = map2 or {}
-	local sizes = get_box_size(sigma, n)
-	for i=1, n do
-		box_blur_2d(map, sizes[i], map2)
-	end
-	return map
-end
-
-return {
-	get_box_size = get_box_size,
-	box_blur_1d = box_blur_1d,
-	box_blur_2d = box_blur_2d,
-	gaussian_blur_approx = gaussian_blur_approx,
-}

terrainlib_lua/rivermapper.lua

@@ -1,510 +0,0 @@
--- rivermapper.lua
-
--- This file provide functions to construct the river tree from an elevation model.
--- Based on a research paper:
---
---     Cordonnier, G., Bovy, B., and Braun, J.:
---     A versatile, linear complexity algorithm for flow routing in topographies with depressions,
---     Earth Surf. Dynam., 7, 549–562, https://doi.org/10.5194/esurf-7-549-2019, 2019.
---
--- 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.
-
-
--- 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 tremove = table.remove
-	local mmax = math.max
-	local mrand = math.random
-
-	local X, Y = dem.X, dem.Y
-	dirs.X = X
-	dirs.Y = Y
-	lakes.X = X
-	lakes.Y = Y
-	local i = 1
-	local dirs2 = {}
-	for i=1, X*Y do
-		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.
-			-- 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
-
-			-- '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
-				singular[#singular+1] = i
-			elseif d == 1 then
-				dirs2[i+X] = dirs2[i+X] + 1
-			elseif d == 2 then
-				dirs2[i+1] = dirs2[i+1] + 2
-			elseif d == 3 then
-				dirs2[i-X] = dirs2[i-X] + 4
-			elseif d == 4 then
-				dirs2[i-1] = dirs2[i-1] + 8
-			end
-			i = i + 1
-		end
-	end
-
-	--------------------------------------
-	-- STEP 2: Compute basins and links --
-	--------------------------------------
-	-- Now water can flow until it reaches a singular node (which is in most cases the bottom of a depression)
-	-- We will calculate the drainage basin of every singular node (all the nodes from which the water will flow in this singular node, directly or indirectly), make an adjacency list of basins, and find the lowest pass between each pair of adjacent basins (they are potential lake outlets)
-	local nbasins = #singular
-	local basin_id = {}
-	local links = {}
-	local basin_links
-
-	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
-			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
-			-- 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
-			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
-
-			-----------
-			-- SOUTH --
-			-----------
-			if d >= 4 then -- River coming from the South
-				d = d - 4
-				cur = cur + 1
-				singular[cur] = 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
-			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
-			end
-
-			----------
-			-- WEST --
-			----------
-			if d >= 2 then -- River coming from the West
-				d = d - 2
-				cur = cur + 1
-				singular[cur] = 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
-			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
-			end
-
-			-----------
-			-- NORTH --
-			-----------
-			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
-					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]
-				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
-	dirs2 = nil
-
-	links[0] = {}
-	local nlinks = {}
-	for i=0, nbasins do
-		nlinks[i] = 0
-	end
-
-	-- Iterate through pairs of adjacent basins, and make the links reciprocal
-	for ib1=1, #links do
-		for ib2, link in pairs(links[ib1]) do
-			if ib2 < ib1 then
-				links[ib2][ib1] = link
-				nlinks[ib1] = nlinks[ib1] + 1
-				nlinks[ib2] = nlinks[ib2] + 1
-			end
-		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
-	local ref = singular -- Reuse table
-	for i=0, nbasins do
-		if nlinks[i] <= 8 then
-			cur = cur + 1
-			lowlevel[cur] = i
-			ref[i] = cur
-		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
-		-- 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
-		for bn, bdata in pairs(lnk1) do
-			if bdata.elev < lowest then
-				lowest = bdata.elev
-				b2 = bn
-			end
-		end
-
-		-- 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
-		basin_graph[bb2][bb1] = bound
-
-		-- Merge basin b1 into b2
-		local lnk2 = links[b2]
-		-- First, remove the link between b1 and b2
-		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[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
-				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
-				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
-				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
-				lnkn[b2] = bdata
-				lnk2[bn] = bdata
-			end
-		end
-	end
-
-	--------------------------------------------------------------
-	-- STEP 4: Orient basin graph, and grid nodes inside basins --
-	--------------------------------------------------------------
-	-- We will finally solve those freaking singular nodes.
-	-- To orient the basin graph, we will consider that the ultimate basin water should flow into is the map outside (basin #0). We will start from it and recursively walk upstream to the neighbouring basins, using only links that are in the minimal spanning tree. This gives the flow direction of the links, and thus, the outlet of every basin.
-	-- This will also give lake elevation, which is the highest link encountered between map outside and the given basin on the spanning tree.
-	-- And within each basin, we need to modify flow directions to connect the singular node to the outlet.
-	local queue = {0}
-	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
-		local b1, elev1 = queue[cur], queuevalues[cur] -- Pop from queue
-		cur = cur - 1
-		basin_lake[b1] = elev1
-		-- Iterate through b1's neighbours (according to the spanning tree)
-		for b2, bound in pairs(basin_graph[b1]) do
-			-- Make b2 flow into b1
-			local i = bound.i -- Get the coordinate of the link (which is the basin's outlet)
-			local dir = bound.is_y and 3 or 4 -- And get the direction (S/E/N/W)
-			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
-					i = i - 1
-				end
-			elseif b1 == 0 then
-				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
-					i = i + 1
-				elseif dir == 3 then
-					i = i - X
-				elseif dir == 4 then
-					i = i - 1
-				end
-				-- Reverse the flow direction for the next node, which will flow into i
-				dir = reverse[dir]
-			until dir == 0 -- Stop when reaching the singular node
-
-			-- Add basin b2 into the queue, and keep the highest link elevation, that will define the elevation of the lake in b2
-			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
-	-- 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 ndonors = {}
-	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
-		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
-				elseif dir == 2 then
-					i2 = i2 + 1
-				elseif dir == 3 then
-					i2 = i2 - X
-				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
-				end
-				dir = dirs[i2]
-			end
-		end
-	end
-
-	return waterq
-end
-
-return {
-	flow_routing = flow_routing,
-	accumulate = accumulate,
-}

terrainlib_lua/twist.lua

@@ -1,102 +0,0 @@
--- twist.lua
-
-local function get_bounds(dirs, rivers)
-	local X, Y = dirs.X, dirs.Y
-	local bounds_x = {X=X, Y=Y}
-	local bounds_y = {X=X, Y=Y}
-	for i=1, X*Y do
-		bounds_x[i] = 0
-		bounds_y[i] = 0
-	end
-
-	for i=1, X*Y do
-		local dir = dirs[i]
-		local river = rivers[i]
-		if dir == 1 then -- South (+Y)
-			bounds_y[i] = river
-		elseif dir == 2 then -- East (+X)
-			bounds_x[i] = river
-		elseif dir == 3 then -- North (-Y)
-			bounds_y[i-X] = river
-		elseif dir == 4 then -- West (-X)
-			bounds_x[i-1] = river
-		end
-	end
-
-	return bounds_x, bounds_y
-end
-
-local function twist(dirs, rivers, n)
-	n = n or 5
-	local X, Y = dirs.X, dirs.Y
-	local bounds_x, bounds_y = get_bounds(dirs, rivers)
-	local dn = 0.5 / n
-
-	local offset_x = {X=X, Y=Y}
-	local offset_y = {X=X, Y=Y}
-	local offset_x_alt = {X=X, Y=Y}
-	local offset_y_alt = {X=X, Y=Y}
-	for i=1, X*Y do
-		offset_x[i] = 0
-		offset_y[i] = 0
-	end
-
-	for nn=1, n do
-		local i = 1
-		for y=1, Y do
-			for x=1, X do
-				local ox, oy = offset_x[i], offset_y[i]
-				if dirs[i] ~= 0 and rivers[i] > 1 then
-					local sum_fx = 0
-					local sum_fy = 0
-					local sum_w = 0
-					local b
-					if x < X then
-						b = bounds_x[i]
-						sum_fx = sum_fx + b*(offset_x[i+1]+1)
-						sum_fy = sum_fy + b*offset_y[i+1]
-						sum_w = sum_w + b
-					end
-					if y < Y then
-						b = bounds_y[i]
-						sum_fx = sum_fx + b*offset_x[i+X]
-						sum_fy = sum_fy + b*(offset_y[i+X]+1)
-						sum_w = sum_w + b
-					end
-					if x > 1 then
-						b = bounds_x[i-1]
-						sum_fx = sum_fx + b*(offset_x[i-1]-1)
-						sum_fy = sum_fy + b*offset_y[i-1]
-						sum_w = sum_w + b
-					end
-					if y > 1 then
-						b = bounds_y[i-X]
-						sum_fx = sum_fx + b*offset_x[i-X]
-						sum_fy = sum_fy + b*(offset_y[i-X]-1)
-						sum_w = sum_w + b
-					end
-					local fx, fy = sum_fx/sum_w - ox, sum_fy/sum_w - oy
-					local fd = (fx*fx+fy*fy) ^ 0.5
-					if fd > dn then
-						local c = dn/fd
-						fx, fy = fx*c, fy*c
-					end
-
-					offset_x_alt[i] = ox+fx
-					offset_y_alt[i] = oy+fy
-				else
-					offset_x_alt[i] = ox
-					offset_y_alt[i] = oy
-				end
-
-				i = i + 1
-			end
-		end
-		offset_x, offset_x_alt = offset_x_alt, offset_x
-		offset_y, offset_y_alt = offset_y_alt, offset_y
-	end
-
-	return offset_x, offset_y
-end
-
-return twist

view_map.py

@@ -5,20 +5,13 @@ import zlib
 import sys
 import os
 
-from view import stats, plot
-from readconfig import read_conf_file
+from terrainlib import stats, plot
 
-os.chdir(sys.argv[1])
-conf = read_conf_file('mapgen_rivers.conf')
-if 'center' in conf:
-    center = conf['center'] == 'true'
-else:
-    center = True
-
-if 'blocksize' in conf:
-    blocksize = float(conf['blocksize'])
-else:
-    blocksize = 15.0
+scale = 1
+if len(sys.argv) > 1:
+    os.chdir(sys.argv[1])
+if len(sys.argv) > 2:
+    scale = int(sys.argv[2])
 
 def load_map(name, dtype, shape):
     dtype = np.dtype(dtype)
@@ -28,10 +21,9 @@ def load_map(name, dtype, shape):
             data = zlib.decompress(data)
         return np.frombuffer(data, dtype=dtype).reshape(shape)
 
-shape = np.loadtxt('river_data/size', dtype='u4')
-shape = (shape[1], shape[0])
-dem = load_map('river_data/dem', '>i2', shape)
-lakes = load_map('river_data/lakes', '>i2', shape)
+shape = np.loadtxt('size', dtype='u4')
+dem = load_map('dem', '>i2', shape)
+lakes = load_map('lakes', '>i2', shape)
 
-stats(dem, lakes, scale=blocksize)
-plot(dem, lakes, scale=blocksize, center=center)
+stats(dem, lakes, scale=scale)
+plot(dem, lakes, scale=scale)