Comment and clarify

README.md

@@ -1,17 +1,17 @@
 mapgen_rivers
 =============
 
-Procedural map generator for Minetest 5.x. Still experimental and basic.
+Procedural map generator for Minetest 5.x. Focused on river networks, and features valley erosion and lakes.
 
 Contains two distinct programs: Python scripts for pre-processing, and Lua scripts to generate the map on Minetest.
 
-![Screenshot](https://user-images.githubusercontent.com/6905002/79073532-7a567f00-7ce7-11ea-9791-8fb453f5175d.png)
+![Screenshot](https://user-images.githubusercontent.com/6905002/79541028-687b3000-8089-11ea-9209-c23c15d75383.png)
 
 # Installation
 This mod should be placed in the `/mods` directory like any other Minetest mod.
 
 The Python part relies on external libraries that you need to install:
-- `numpy`, a widely used library for numerical calculations
+- `numpy` and `scipy`, widely used libraries for numerical calculations
 - `noise`, doing Perlin/Simplex noises
 - optionally, `matplotlib` (for map preview)
 
@@ -23,7 +23,7 @@ Run the script `terrain_rivers.py` via command line. You can optionally append t
 ```
 ./terrain_rivers.py 1000
 ```
-For a default 400x400 map, it should take between 1 and 2 minutes. It will generate 5 files directly in the mod folder, containing the map data (1.4 MB for the default size).
+For a default 401x401 map, it should take between 1 and 2 minutes. It will generate 5 files directly in the mod folder, containing the map data.
 
 ## Map generation
 Just create a Minetest world with `singlenode` mapgen, enable this mod and start the world. The data files are immediately copied in the world folder so you can re-generate them afterwards, it won't affect the old worlds.

bounds.py

@@ -1,7 +1,10 @@
 import numpy as np
-import matplotlib.pyplot as plt
 
 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='i4')
     bounds_v = np.zeros((Y-1, X), dtype='i4')
@@ -14,6 +17,10 @@ def make_bounds(dirs, rivers):
     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
@@ -28,6 +35,11 @@ def get_fixed(dirs):
     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
@@ -55,7 +67,7 @@ def twist(bounds_x, bounds_y, fixed, d=0.1, n=5):
 
         length = np.hypot(force_x, force_y)
         length[length==0] = 1
-        coeff = d / length * moveable
+        coeff = d / length * moveable # Normalize, take into account the direction only
         offset_x += force_x * coeff
         offset_y += force_y * coeff
 

erosion.py

@@ -3,22 +3,33 @@ import scipy.ndimage as im
 import rivermapper as rm
 
 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
+    """
+
     dirs = dirs.copy()
     dirs[0,:] = 0
     dirs[-1,:] = 0
     dirs[:,0] = 0
     dirs[:,-1] = 0
 
-    adv_time = 1 / (K*rivers**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
@@ -32,19 +43,19 @@ def advection(dem, dirs, rivers, time, K=1, m=0.5, sea_level=0):
                 elif flow_dir == 4:
                     x1 -= 1
 
-                if remaining <= adv_time[y0,x0]:
+                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]
+            dem_new[y,x] = c*dem[y1,x1] + (1-c)*dem[y0,x0] # If between 2 pixels, perform linear interpolation.
 
     return np.minimum(dem, dem_new)
 
 def diffusion(dem, time, d=1):
     radius = d * time**.5
-    return im.gaussian_filter(dem, radius, mode='reflect')
+    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):
@@ -78,9 +89,9 @@ class EvolutionModel:
         self.flow_uptodate = False
 
     def define_isostasy(self):
-        self.ref_isostasy = im.gaussian_filter(self.dem, self.flex_radius, mode='reflect')
+        self.ref_isostasy = im.gaussian_filter(self.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')
+        isostasy = im.gaussian_filter(self.dem, self.flex_radius, mode='reflect') # Calculate blurred DEM
-        correction = (self.ref_isostasy - isostasy) * rate
+        correction = (self.ref_isostasy - isostasy) * rate # Compare it with the reference isostasy
-        self.dem = self.dem + correction
+        self.dem = self.dem + correction # Adjust

geometry.lua

@@ -1,27 +1,33 @@
 local function distance_to_segment(x1, y1, x2, y2, x, y)
 	-- get the distance between point (x,y) and segment (x1,y1)-(x2,y2)
-	local a = (x1-x2)^2 + (y1-y2)^2
+	local 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)
-	-- X, Y 4-vectors giving the coordinates of the 4 nodes
+	-- 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)

init.lua

@@ -12,6 +12,7 @@ local make_polygons = dofile(modpath .. 'polygons.lua')
 
 local transform_quadri = dofile(modpath .. 'geometry.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)
@@ -44,6 +45,7 @@ local function generate(minp, maxp, seed)
 				local xf, zf = transform_quadri(poly.x, poly.z, x/blocksize, z/blocksize)
 				local i00, i01, i11, i10 = unpack(poly.i)
 
+				-- Test the 4 edges to see whether we are in a river or not
 				local is_river = false
 				local depth_factor = 0
 				local r_west, r_north, r_east, r_south = unpack(poly.rivers)
@@ -66,7 +68,7 @@ local function generate(minp, maxp, seed)
 					zf = 0
 				end
 
-				if not is_river then
+				if not is_river then -- Test corners also
 					local c_NW, c_NE, c_SE, c_SW = unpack(poly.river_corners)
 					if xf+zf < c_NW then
 						is_river = true
@@ -87,7 +89,7 @@ local function generate(minp, maxp, seed)
 					end
 				end
 
-				if not is_river then
+				if not is_river then -- Redefine indicesto have 0/1 on the riverbanks (avoids ugly edge cuts, at least for small rivers)
 					xf = (xf-r_west) / (r_east-r_west)
 					zf = (zf-r_north) / (r_south-r_north)
 				end

polygons.lua

@@ -44,6 +44,7 @@ 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
@@ -66,26 +67,33 @@ local function river_width(flow)
 	return math.min(wfactor * flow ^ wpower, 1)
 end
 
+-- On map generation, determine into which polygon every point (in 2D) will fall.
+-- Also store polygon-specific data
 local function make_polygons(minp, maxp)
 	local chulens = maxp.z - minp.z + 1
 
 	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 = math.max(math.floor(minp.x/blocksize - 0.5), 0), math.min(math.ceil(maxp.x/blocksize), X-2)
 	local zpmin, zpmax = math.max(math.floor(minp.z/blocksize - 0.5), 0), math.min(math.ceil(maxp.z/blocksize), Z-2)
 
+	-- Iterate over the polygons
 	for xp = xpmin, xpmax do
 		for zp=zpmin, zpmax do
 			local iA = index(xp, zp)
 			local iB = index(xp+1, zp)
 			local iC = index(xp+1, zp+1)
 			local iD = index(xp, zp+1)
+			-- Extract the vertices of the polygon
 			local poly_x = {offset_x[iA]+xp, offset_x[iB]+xp+1, offset_x[iC]+xp+1, offset_x[iD]+xp}
 			local poly_z = {offset_z[iA]+zp, offset_z[iB]+zp, offset_z[iC]+zp+1, offset_z[iD]+zp+1}
 			local polygon = {x=poly_x, z=poly_z, i={iA, iB, iC, iD}}
 
-			local bounds = {}
+			local bounds = {} -- Will be a list of the intercepts of polygon edges for every X position (scanline algorithm)
+			-- Calculate the min and max X positions 
 			local xmin = math.max(math.floor(blocksize*math.min(unpack(poly_x)))+1, minp.x)
 			local xmax = math.min(math.floor(blocksize*math.max(unpack(poly_x))), maxp.x)
+			-- And initialize the arrays
 			for x=xmin, xmax do
 				bounds[x] = {}
 			end
@@ -93,27 +101,33 @@ local function make_polygons(minp, maxp)
 			local i1 = 4
 			for i2=1, 4 do -- Loop on 4 edges
 				local x1, x2 = poly_x[i1], poly_x[i2]
+				-- Calculate the integer X positions over which this edge spans
 				local lxmin = math.floor(blocksize*math.min(x1, x2))+1
 				local lxmax = math.floor(blocksize*math.max(x1, x2))
-				if lxmin <= lxmax then
+				if lxmin <= lxmax then -- If there is at least one position in it
 					local z1, z2 = poly_z[i1], poly_z[i2]
+					-- Calculate coefficient of the equation defining the edge: Z=aX+b
 					local a = (z1-z2) / (x1-x2)
 					local b = blocksize*(z1 - a*x1)
 					for x=math.max(lxmin, minp.x), math.min(lxmax, maxp.x) do
+						-- For every X position involved, add the intercepted Z position in the table
 						table.insert(bounds[x], a*x+b)
 					end
 				end
 				i1 = i2
 			end
 			for x=xmin, xmax do
+				-- Now sort the bounds list
 				local xlist = bounds[x]
 				table.sort(xlist)
 				local c = math.floor(#xlist/2)
 				for l=1, c do
+					-- Take pairs of Z coordinates: all positions between them belong to the polygon.
 					local zmin = math.max(math.floor(xlist[l*2-1])+1, minp.z)
 					local zmax = math.min(math.floor(xlist[l*2]), maxp.z)
 					local i = (x-minp.x) * chulens + (zmin-minp.z) + 1
 					for z=zmin, zmax do
+						-- Fill the map at these places
 						polygons[i] = polygon
 						i = i + 1
 					end
@@ -123,10 +137,13 @@ local function make_polygons(minp, maxp)
 			polygon.dem = {dem[iA], dem[iB], dem[iC], dem[iD]}
 			polygon.lake = math.min(lakes[iA], lakes[iB], lakes[iC], lakes[iD])
 
+			-- Now, rivers.
+			-- Start by finding the river width (if any) for the polygon's 4 edges.
 			local river_west = river_width(bounds_z[iA])
 			local river_north = river_width(bounds_x[iA-zp])
 			local river_east = 1-river_width(bounds_z[iB])
 			local river_south = 1-river_width(bounds_x[iD-zp-1])
+			-- Only if opposite rivers overlap (should be rare)
 			if river_west > river_east then
 				local mean = (river_west + river_east) / 2
 				river_west = mean
@@ -139,6 +156,7 @@ local function make_polygons(minp, maxp)
 			end
 			polygon.rivers = {river_west, river_north, river_east, river_south}
 
+			-- Look for river corners
 			local around = {0,0,0,0,0,0,0,0}
 			if zp > 0 then
 				around[1] = river_width(bounds_z[iA-X])

rivermapper.py

@@ -19,9 +19,14 @@ neighbours_dirs = np.array([
 neighbours_pattern = neighbours_dirs > 0
 
 def flow_dirs_lakes(dem, random=0):
+    """
+    Calculates flow direction in D4 (4 choices) for every pixel of the DEM
+    Also returns an array of lake elevation
+    """
+
     (Y, X) = dem.shape
 
-    dem_margin = np.zeros((Y+2, X+2))
+    dem_margin = np.zeros((Y+2, X+2)) # We need a margin of one pixel at every edge, to prevent crashes when scanning the neighbour pixels on the borders
     dem_margin[1:-1,1:-1] = dem
     if random > 0:
         dem_margin += np.random.random(dem_margin.shape) * random
@@ -41,10 +46,11 @@ def flow_dirs_lakes(dem, random=0):
         dem_east = dem_margin[y,X]
         borders.append((dem_east, dem_east, y, X))
 
+    # Make a binary heap
     heapq.heapify(borders)
 
     dirs = np.zeros((Y+2, X+2), dtype='u1')
-    dirs[-2:,:] = 1
+    dirs[-2:,:] = 1 # Border pixels flow outside the map
     dirs[:,-2:] = 2
     dirs[ :2,:] = 3
     dirs[:, :2] = 4
@@ -56,21 +62,26 @@ def flow_dirs_lakes(dem, random=0):
         heapq.heappush(borders, (alt, altmax, y, x))
 
     while len(borders) > 0:
-        (alt, altmax, y, x) = heapq.heappop(borders)
+        (alt, altmax, y, x) = heapq.heappop(borders) # Take the lowest pixel in the queue
         neighbours = dirs[y-1:y+2, x-1:x+2]
-        empty_neighbours = (neighbours == 0) * neighbours_pattern
+        empty_neighbours = (neighbours == 0) * neighbours_pattern # Find the neighbours whose flow direction is not yet defined
-        neighbours += empty_neighbours * neighbours_dirs
+        neighbours += empty_neighbours * neighbours_dirs # They flow into the pixel being studied
 
-        lake = max(alt, altmax)
+        lake = max(alt, altmax) # Set lake elevation to the maximal height of the downstream section.
         lakes[y-1,x-1] = lake
 
         coords = np.transpose(empty_neighbours.nonzero())
-        for (dy,dx) in coords-1:
+        for (dy,dx) in coords-1: # Add these neighbours into the queue
             add_point(y+dy, x+dx, lake)
 
     return dirs[1:-1,1:-1], lakes
 
 def accumulate(dirs, dem=None):
+    """
+    Calculates the quantity of water that accumulates at every pixel,
+    following flow directions.
+    """
+
     (Y, X) = dirs.shape
     dirs_margin = np.zeros((Y+2,X+2))-1
     dirs_margin[1:-1,1:-1] = dirs
@@ -79,13 +90,13 @@ def accumulate(dirs, dem=None):
     def calculate_quantity(y, x):
         if quantity[y,x] > 0:
             return quantity[y,x]
-        q = 1
+        q = 1 # Consider that every pixel contains a water quantity of 1 by default.
         neighbours = dirs_margin[y:y+3, x:x+3]
-        donors = neighbours == neighbours_dirs
+        donors = neighbours == neighbours_dirs # Identify neighbours that give their water to the pixel being studied
 
         coords = np.transpose(donors.nonzero())
         for (dy,dx) in coords-1:
-            q += calculate_quantity(y+dy, x+dx)
+            q += calculate_quantity(y+dy, x+dx) # Add water quantity of the donors pixels (this triggers calculation for these pixels, recursively)
         quantity[y, x] = q
         return q
 
@@ -96,5 +107,9 @@ def accumulate(dirs, dem=None):
     return quantity
 
 def flow(dem):
+    """
+    Calculates flow directions and water quantity
+    """
+
     dirs, lakes = flow_dirs_lakes(dem)
     return dirs, lakes, accumulate(dirs, dem)

terrain_rivers.py

@@ -20,28 +20,29 @@ else:
 scale = (mapsize-1) / 2
 n = np.zeros((mapsize, mapsize))
 
-#micronoise_depth = 0.05
+# Set noise parameters
-
 params = {
     "octaves" : int(np.ceil(np.log2(mapsize-1)))+1,
     "persistence" : 0.5,
     "lacunarity" : 2.,
 }
 
+# Determine noise offset randomly
 xbase = np.random.randint(65536)
 ybase = np.random.randint(65536)
 
+# Generate the noise
 for x in range(mapsize):
     for y in range(mapsize):
         n[x,y] = noise.snoise2(x/scale + xbase, y/scale + ybase, **params)
 
-#micronoise = np.random.rand(mapsize, mapsize)
-#nn = np.exp(n*2) + micronoise*micronoise_depth
 nn = n*mapsize/5 + mapsize/20
 
+# Initialize landscape evolution model
 print('Initializing model')
 model = EvolutionModel(nn, K=1, m=0.35, d=1, sea_level=0)
 
+# Run the model's processes: the order in which the processes are run is arbitrary and could be changed.
 print('Flow calculation 1')
 model.calculate_flow()
 
@@ -68,12 +69,15 @@ model.calculate_flow()
 
 print('Done')
 
+# Twist the grid
 bx, by = bounds.make_bounds(model.dirs, model.rivers)
 ox, oy = bounds.twist(bx, by, bounds.get_fixed(model.dirs))
 
+# Convert offset in 8-bits
 offset_x = np.clip(np.floor(ox * 256), -128, 127)
 offset_y = np.clip(np.floor(oy * 256), -128, 127)
 
+# Save the files
 save(model.dem, 'dem', dtype='>i2')
 save(model.lakes, 'lakes', dtype='>i2')
 save(np.abs(bx), 'bounds_x', dtype='>i4')
@@ -86,6 +90,7 @@ save(model.rivers, 'rivers', dtype='>u4')
 with open('size', 'w') as sfile:
     sfile.write('{:d}\n{:d}'.format(mapsize, mapsize))
 
+# Display the map if matplotlib is found
 try:
     import matplotlib.pyplot as plt
 
@@ -111,7 +116,7 @@ try:
     plt.pcolormesh(model.rivers, vmin=0, vmax=mapsize**2/25, cmap='Blues')
     plt.gca().set_aspect('equal', 'box')
     #plt.colorbar(orientation='horizontal')
-    plt.title('Rivers discharge')
+    plt.title('Rivers flux')
 
     plt.show()
 except:

view_map.py

@@ -34,6 +34,6 @@ plt.subplot(1,3,3)
 plt.pcolormesh(np.log(rivers), vmin=0, vmax=np.log(n/25), cmap='Blues')
 plt.gca().set_aspect('equal', 'box')
 #plt.colorbar(orientation='horizontal')
-plt.title('Rivers discharge')
+plt.title('Rivers flux')
 
 plt.show()