view.py: use a different colormap for sea, to distinguish it from lakes

2025-07-24 19:30:24 +02:00 · 2021-09-07 12:00:47 +02:00
16 changed files with 202 additions and 303 deletions
--- a/CHANGELOG.md
+++ b/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
--- a/README.md
+++ b/README.md
@ -1,5 +1,5 @@
 # Map Generator with Rivers
-`mapgen_rivers v1.0.2` by Gaël de Sailly.
+`mapgen_rivers v1.0` by Gaël de Sailly.

 Semi-procedural map generator for Minetest 5.x. It aims to create realistic and nice-looking landscapes for the game, focused on river networks. It is based on algorithms modelling water flow and river erosion at a broad scale, similar to some used by researchers in Earth Sciences. It is taking some inspiration from [Fastscape](https://github.com/fastscape-lem/fastscape).

@ -13,7 +13,6 @@ It used to be composed of a Python script doing pre-generation, and a Lua mod re
 License: GNU LGPLv3.0

 Code: Gaël de Sailly
-
 Flow routing algorithm concept (in `terrainlib/rivermapper.lua`): Cordonnier, G., Bovy, B., & Braun, J. (2019). A versatile, linear complexity algorithm for flow routing in topographies with depressions. Earth Surface Dynamics, 7(2), 549-562.

 # Requirements
--- a/geometry.lua
+++ b/geometry.lua
@ -1,6 +1,3 @@
-local sqrt, abs = math.sqrt, math.abs
-local unpk = unpack
-
 local function distance_to_segment(x1, y1, x2, y2, x, y)
 	-- get the distance between point (x,y) and segment (x1,y1)-(x2,y2)
 	local a = (x1-x2)^2 + (y1-y2)^2 -- square of distance
@ -8,13 +5,13 @@ local function distance_to_segment(x1, y1, x2, y2, x, y)
 	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)
+		return math.sqrt(b)
 	elseif a + c < b then
 		-- The closest point of the segment is the extremity 2
-		return sqrt(c)
+		return math.sqrt(c)
 	else
 		-- The closest point is on the segment
-		return abs(x1 * (y2-y) + x2 * (y-y1) + x * (y1-y2)) / sqrt(a)
+		return math.abs(x1 * (y2-y) + x2 * (y-y1) + x * (y1-y2)) / math.sqrt(a)
 	end
 end

@ -22,8 +19,8 @@ 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)
+	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)
--- a/heightmap.lua
+++ b/heightmap.lua
@ -8,10 +8,6 @@ local riverbed_slope = mapgen_rivers.settings.riverbed_slope * mapgen_rivers.set

 local MAP_BOTTOM = -31000

-- Localize for performance
-local floor, min, max = math.floor, math.min, math.max
-local unpk = unpack
-
 -- Linear interpolation
 local function interp(v00, v01, v11, v10, xf, zf)
 	local v0 = v01*xf + v00*(1-xf)
@ -34,11 +30,11 @@ local function 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:
@ -68,10 +64,10 @@ local function heightmaps(minp, maxp)

 				-- Transform the coordinates to have xf and zf = 0 or 1 in rivers (to avoid rivers having lateral slope and to accomodate the surrounding smoothly)
 				if imax == 0 then
-					local x0 = max(r_west, c_NW-zf, zf-c_SW)
-					local x1 = min(r_east, c_NE+zf, c_SE-zf)
-					local z0 = max(r_north, c_NW-xf, xf-c_NE)
-					local z1 = min(r_south, c_SW+xf, c_SE-xf)
+					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
@ -94,7 +90,7 @@ local function heightmaps(minp, maxp)

 				-- 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],
@ -119,10 +115,10 @@ local function 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 depth_factor_max > 0 then
-					terrain_height = min(max(lake_height, sea_level) - floor(1+depth_factor_max*riverbed_slope), terrain_height)
+					terrain_height = math.min(math.max(lake_height, sea_level) - math.floor(1+depth_factor_max*riverbed_slope), terrain_height)
 				end

 				terrain_height_map[i] = terrain_height
--- a/init.lua
+++ b/init.lua
@ -6,7 +6,7 @@ mapgen_rivers.world_data_path = minetest.get_worldpath() .. '/river_data/'

 if minetest.get_mapgen_setting("mg_name") ~= "singlenode" then
 	minetest.set_mapgen_setting("mg_name", "singlenode", true)
-	minetest.log("warning", "[mapgen_rivers] Mapgen set to singlenode")
+	print("[mapgen_rivers] Mapgen set to singlenode")
 end

 dofile(modpath .. 'settings.lua')
@ -32,9 +32,6 @@ local function interp(v00, v01, v11, v10, xf, zf)
 	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
@ -51,7 +48,7 @@ local sumtime2 = 0
 local ngen = 0

 local function generate(minp, maxp, seed)
-	minetest.log("info", ("[mapgen_rivers] Generating from %s to %s"):format(minetest.pos_to_string(minp), minetest.pos_to_string(maxp)))
+	print(("[mapgen_rivers] Generating from %s to %s"):format(minetest.pos_to_string(minp), minetest.pos_to_string(maxp)))

 	local chulens = {
 		x = maxp.x-minp.x+1,
@ -114,8 +111,8 @@ local function generate(minp, maxp, seed)
 			end
 		end

-		local pminp = {x=floor(xmin), z=floor(zmin)}
-		local pmaxp = {x=floor(xmax)+1, z=floor(zmax)+1}
+		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)
@ -123,30 +120,6 @@ local function generate(minp, maxp, seed)
 		terrain_map, lake_map = 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("default:stone")
 	local c_dirt = minetest.get_content_id("default:dirt")
 	local c_lawn = minetest.get_content_id("default:dirt_with_grass")
@ -172,7 +145,7 @@ local function generate(minp, maxp, seed)

 	for z = minp.z, maxp.z do
 		for x = minp.x, maxp.x do
-			local ivm = a:index(x, maxp.y+1, z)
+			local ivm = a:index(x, minp.y, z)
 			local ground_above = false
 			local temperature
 			if use_biomes then
@ -188,8 +161,8 @@ local function generate(minp, maxp, seed)
 				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 x0 = math.floor(xn)
+					local z0 = math.floor(zn)

 					local i0 = i_origin + z0*incr + x0
 					local i1 = i0+1
@ -197,12 +170,13 @@ local function generate(minp, maxp, seed)
 					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])
+					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
@ -231,7 +205,7 @@ local function generate(minp, maxp, seed)

 				ground_above = y <= terrain

-				ivm = ivm - ystride
+				ivm = ivm + ystride
 				if use_distort then
 					nid = nid + incrY
 				end
@ -259,17 +233,18 @@ local function generate(minp, maxp, seed)
 	vm:calc_lighting()
 	vm:update_liquids()
 	vm:write_to_map()
+	local t1 = os.clock()

-	local t = os.clock()-t0
+	local t = t1-t0
 	ngen = ngen + 1
 	sumtime = sumtime + t
 	sumtime2 = sumtime2 + t*t
-	minetest.log("verbose", ("[mapgen_rivers] Done in %5.3f s"):format(t))
+	print(("[mapgen_rivers] Done in %5.3f s"):format(t))
 end

 minetest.register_on_generated(generate)
 minetest.register_on_shutdown(function()
 	local avg = sumtime / ngen
 	local std = math.sqrt(sumtime2/ngen - avg*avg)
-	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))
+	print(("[mapgen_rivers] Mapgen statistics:\n- Mapgen calls: %4d\n- Mean time: %5.3f s\n- Standard deviation: %5.3f s"):format(ngen, avg, std))
 end)
--- a/load.lua
+++ b/load.lua
@ -1,15 +1,12 @@
 local worldpath = mapgen_rivers.world_data_path

-local floor = math.floor
-local sbyte, schar = string.byte, string.char
-local unpk = unpack
-
 function mapgen_rivers.load_map(filename, bytes, signed, size, converter)
 	local file = io.open(worldpath .. filename, 'rb')
 	local data = file:read('*all')
 	if #data < bytes*size then
 		data = minetest.decompress(data)
 	end
+	local sbyte = string.byte

 	local map = {}

@ -38,6 +35,8 @@ function mapgen_rivers.load_map(filename, bytes, signed, size, converter)
 	return map
 end

+local sbyte = string.byte
+
 local loader_mt = {
 	__index = function(loader, i)
 		local file = loader.file
@ -76,6 +75,9 @@ end
 function mapgen_rivers.write_map(filename, data, bytes)
 	local size = #data
 	local file = io.open(worldpath .. filename, 'wb')
+	local mfloor = math.floor
+	local schar = string.char
+	local upack = unpack

 	local bytelist = {}
 	for j=1, bytes do
@ -83,15 +85,15 @@ function mapgen_rivers.write_map(filename, data, bytes)
 	end

 	for i=1, size do
-		local n = floor(data[i])
+		local n = mfloor(data[i])
 		data[i] = n
 		for j=bytes, 2, -1 do
 			bytelist[j] = n % 256
-			n = floor(n / 256)
+			n = mfloor(n / 256)
 		end
 		bytelist[1] = n % 256

-		file:write(schar(unpk(bytelist)))
+		file:write(schar(upack(bytelist)))
 	end

 	file:close()
--- a/noises.lua
+++ b/noises.lua
@ -73,7 +73,7 @@ for name, np in pairs(mapgen_rivers.noise_params) do
 	if lac > 1 then
 		local omax = math.floor(math.log(math.min(np.spread.x, np.spread.y, np.spread.z)) / math.log(lac))+1
 		if np.octaves > omax then
-			minetest.log("warning", "[mapgen_rivers] Noise " .. name .. ": 'octaves' reduced to " .. omax)
+			print("[mapgen_rivers] Noise " .. name .. ": 'octaves' reduced to " .. omax)
 			np.octaves = omax
 		end
 	end
--- a/polygons.lua
+++ b/polygons.lua
@ -33,7 +33,7 @@ if first_mapgen then
 	-- Generate a map!!
 	local pregenerate = dofile(mapgen_rivers.modpath .. '/pregenerate.lua')
 	minetest.register_on_mods_loaded(function()
-		minetest.log("action", '[mapgen_rivers] Generating grid, this may take a while...')
+		print('[mapgen_rivers] Generating grid')
 		pregenerate(load_all)

 		if load_all then
@ -58,9 +58,9 @@ if not (first_mapgen and load_all) then

 	minetest.register_on_mods_loaded(function()
 		if load_all then
-			minetest.log("action", '[mapgen_rivers] Loading full grid')
+			print('[mapgen_rivers] Loading full grid')
 		else
-			minetest.log("action", '[mapgen_rivers] Loading grid as interactive loaders')
+			print('[mapgen_rivers] Loading grid as interactive loaders')
 		end
 		local grid = mapgen_rivers.grid

@ -90,19 +90,16 @@ if mapgen_rivers.settings.center then
 	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)
 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
@ -141,8 +138,8 @@ local function 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)

 	-- Iterate over the polygons
 	for xp = xpmin, xpmax do
@ -168,8 +165,8 @@ local function make_polygons(minp, maxp)

 			local bounds = {} -- Will be a list of the intercepts of polygon edges for every Z position (scanline algorithm)
 			-- Calculate the min and max Z positions 
-			local zmin = max(floor(min(unpack(poly_z)))+1, minp.z)
-			local 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] = {}
@ -179,14 +176,14 @@ local function 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
@ -197,11 +194,11 @@ local function 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
@ -223,16 +220,16 @@ local function 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

--- a/pregenerate.lua
+++ b/pregenerate.lua
@ -33,7 +33,7 @@ local function pregenerate(keep_loaded)
 	local tectonic_step = tectonic_speed * time_step
 	collectgarbage()
 	for i=1, niter do
-		minetest.log("info", "[mapgen_rivers] Iteration " .. i .. " of " .. niter)
+		print("[mapgen_rivers] Iteration " .. i .. " of " .. niter)

 		model:diffuse(time_step)
 		model:flow()
--- a/settings.lua
+++ b/settings.lua
@ -1,7 +1,7 @@
 local mtsettings = minetest.settings
 local mgrsettings = Settings(minetest.get_worldpath() .. '/mapgen_rivers.conf')

-mapgen_rivers.version = "1.0.2"
+mapgen_rivers.version = "1.0"

 local previous_version_mt = mtsettings:get("mapgen_rivers_version") or "0.0"
 local previous_version_mgr = mgrsettings:get("version") or "0.0"
--- a/settingtypes.txt
+++ b/settingtypes.txt
@ -23,7 +23,7 @@ mapgen_rivers_min_catchment (Minimal catchment area) float 3600.0 100.0 1000000.

 #    Coefficient describing how rivers widen when merging.
 #    Riwer width is a power law W = a*D^p. D is river flow and p is this parameter.
-#    Higher value means that a river will grow more when receiving a tributary.
+#    Higher value means a river needs to receive more tributaries to grow in width.
 #    Note that a river can never exceed 2*blocksize.
 mapgen_rivers_river_widening_power (River widening power) float 0.5 0.0 1.0

--- a/terrainlib_lua/erosion.lua
+++ b/terrainlib_lua/erosion.lua
@ -1,13 +1,7 @@
 -- erosion.lua

-- This is the main file of terrainlib_lua. It registers the EvolutionModel object and some of the 
-
 local function erode(model, time)
-	-- Apply river erosion on the model
-	-- 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 tinsert = table.insert
 	local mmin, mmax = math.min, math.max
 	local dem = model.dem
 	local dirs = model.dirs
@ -28,9 +22,9 @@ local function erode(model, time)

 	if scalars then
 		for i=1, X*Y do
-			local etime = 1 / (K*rivers[i]^m) -- Inverse of erosion speed (Ex); time needed for the erosion wave to move through the river section.
+			local etime = 1 / (K*rivers[i]^m)
 			erosion_time[i] = etime
-			lakes[i] = mmax(lakes[i], dem[i], sea_level) -- Use lake/sea surface if higher than ground level, because rivers can not erode below.
+			lakes[i] = mmax(lakes[i], dem[i], sea_level)
 		end
 	else
 		for i=1, X*Y do
@ -45,12 +39,10 @@ local function erode(model, time)
 		local remaining = time
 		local new_elev
 		while true do
-			-- Explore downstream until we find the point 'iw' from which the erosion wave will reach 'i'
 			local inext = iw
 			local d = dirs[iw]

-			-- Follow the river downstream (move 'iw')
-			if d == 0 then -- If no flow direction, we reach the border of the map: set elevation to the latest node's elev and abort.
+			if d == 0 then
 				new_elev = lakes[iw]
 				break
 			elseif d == 1 then
@ -64,13 +56,13 @@ local function erode(model, time)
 			end

 			local etime = erosion_time[iw]
-			if remaining <= etime then -- We have found the node from which the erosion wave will take 'time' to arrive to 'i'.
+			if remaining <= etime then
 				local c = remaining / etime
-				new_elev = (1-c) * lakes[iw] + c * lakes[inext] -- Interpolate linearly between the two nodes
+				new_elev = (1-c) * lakes[iw] + c * lakes[inext]
 				break
 			end

-			remaining = remaining - etime -- If we still don't reach the target time, decrement time and move to next point.
+			remaining = remaining - etime
 			iw = inext
 		end

@ -79,55 +71,50 @@ local function erode(model, time)
 end

 local function diffuse(model, time)
-	-- Apply diffusion using finite differences methods
-	-- Adapted for small radiuses
-	local mmax = math.max
-	local dem = model.dem
-	local X, Y = dem.X, dem.Y
-	local d = model.params.d
-	-- 'd' is equal to 4 times the diffusion coefficient
-	local dmax = d
-	if type(d) == "table" then
-		dmax = -math.huge
-		for i=1, X*Y do
-			dmax = mmax(dmax, d[i])
-		end
-	end
+    local mmax = math.max
+    local dem = model.dem
+    local X, Y = dem.X, dem.Y
+    local d = model.params.d
+    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 diff = dmax * time
+    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
+    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
+                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
+        for i=1, X*Y do
+            dem[i] = dem[i] + temp[i]*ddiff
+        end
+    end
+    -- TODO Test this
 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
+    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")
@ -139,7 +126,6 @@ local function flow(model)
 end

 local function uplift(model, time)
-	-- Raises the terrain according to uplift rate (model.params.uplift)
 	local dem = model.dem
 	local X, Y = dem.X, dem.Y
 	local uplift_rate = model.params.uplift
@ -156,7 +142,6 @@ local function uplift(model, time)
 end

 local function noise(model, time)
-	-- Adds noise to the terrain according to noise depth (model.params.noise)
 	local random = math.random
 	local dem = model.dem
 	local noise_depth = model.params.noise * 2 * time
@ -166,43 +151,35 @@ local function noise(model, time)
 	end
 end

-- Isostasy
-- This is the geological phenomenon that makes the lithosphere "float" over the underlying layers.
-- One of the key implications is that when a very large mass is removed from the ground, the lithosphere reacts by moving upward. This compensation only occurs at large scale (as the lithosphere is not flexible enough for small scale adjustments) so the implementation is using a very large-window Gaussian blur of the elevation array.
-
-- This implementation is quite simplistic, it does not do a mass balance of the lithosphere as this would introduce too many parameters. Instead, it defines a reference equilibrium elevation, and the ground will react toward this elevation (at the scale of the gaussian window).
-- A change in reference isostasy during the run can also be used to simulate tectonic forcing, like making a new mountain range appear.
 local function define_isostasy(model, ref, link)
-	ref = ref or model.dem
-	if link then
-		model.isostasy_ref = ref
-		return
-	end
+    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
+    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
+    return ref2
 end

-- Apply isostasy
 local function isostasy(model)
 	local dem = model.dem
 	local X, Y = dem.X, dem.Y
 	local temp = {X=X, Y=Y}
 	local ref = model.isostasy_ref
 	for i=1, X*Y do
-		temp[i] = ref[i] - dem[i] -- Compute the difference between the ground level and the target level
+		temp[i] = ref[i] - dem[i]
 	end

-	-- Blur the difference map using Gaussian blur
 	gaussian.gaussian_blur_approx(temp, model.params.compensation_radius, 4)
 	for i=1, X*Y do
-		dem[i] = dem[i] + temp[i] -- Apply the difference
+		dem[i] = dem[i] + temp[i]
 	end
 end

--- a/terrainlib_lua/gaussian.lua
+++ b/terrainlib_lua/gaussian.lua
@ -71,6 +71,20 @@ local function box_blur_2d(map1, size, map2)
 	return map1
 end

+--[[local function gaussian_blur(map, std, tail)
+	local exp = math.exp
+
+	local kernel_mid = math.ceil(std*tail) + 1
+	local kernel_size = kernel_mid * 2 - 1
+	local kernel = {}
+	local cst1 = 1/(std*(2*math.pi)^0.5)
+	local cst2 = -1/(2*std^2)
+	for i=1, kernel_size do
+		kernel[i] = cst1 * exp((i-kernel_mid)^2 * cst2)
+	end
+
+	]]
+
 local function gaussian_blur_approx(map, sigma, n, map2)
 	map2 = map2 or {}
 	local sizes = get_box_size(sigma, n)
--- a/terrainlib_lua/rivermapper.lua
+++ b/terrainlib_lua/rivermapper.lua
@ -9,26 +9,24 @@
 --
 -- Big thanks to them for releasing this paper under a free license ! :)

-- The algorithm here makes use of most of the paper's concepts, including the Planar Borůvka algorithm.
+-- The algorithm here makes use of most of the paper's concepts, including the Planar Boruvka algorithm.
 -- Only flow_local and accumulate_flow are custom algorithms.

-
 local function flow_local_semirandom(plist)
-	-- Determines how water should flow at 1 node scale.
-	-- The straightforward approach would be "Water will flow to the lowest of the 4 neighbours", but here water flows to one of the lower neighbours, chosen randomly, but probability depends on height difference.
-	-- This makes rivers better follow the curvature of the topography at large scale, and be less biased by pure N/E/S/W directions.
-	-- 'plist': array of downward height differences (0 if upward)
 	local sum = 0
 	for i=1, #plist do
-		sum = sum + plist[i] -- Sum of probabilities
+		sum = sum + plist[i]
 	end
-
+	--for _, p in ipairs(plist) do
+		--sum = sum + p
+	--end
 	if sum == 0 then
 		return 0
 	end
 	local r = math.random() * sum
 	for i=1, #plist do
 		local p = plist[i]
+	--for i, p in ipairs(plist) do
 		if r < p then
 			return i
 		end
@ -37,14 +35,11 @@ local function flow_local_semirandom(plist)
 	return 0
 end

-- Maybe implement more flow methods in the future?
 local flow_methods = {
 	semirandom = flow_local_semirandom,
 }

-- Applies all steps of the flow routing, to calculate flow direction for every node, and lake surface elevation.
-- It's quite a hard piece of code, but we will go step by step and explain what's going on, so stay with me and... let's goooooooo!
-local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are optional tables to reuse for memory optimization, they may contain any data.
+local function flow_routing(dem, dirs, lakes, method)
 	method = method or 'semirandom'
 	local flow_local = flow_methods[method] or flow_local_semirandom

@ -52,6 +47,7 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
 	lakes = lakes or {}

 	-- Localize for performance
+	--local tinsert = table.insert
 	local tremove = table.remove
 	local mmax = math.max

@ -66,15 +62,11 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
 		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]
-			local plist = { -- Get the height difference of the 4 neighbours (and 0 if uphill)
+			local plist = {
 				y<Y and mmax(zi-dem[i+X], 0) or 0, -- Southward
 				x<X and mmax(zi-dem[i+1], 0) or 0, -- Eastward
 				y>1 and mmax(zi-dem[i-X], 0) or 0, -- Northward
@ -82,10 +74,8 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
 			}

 			local d = flow_local(plist)
-			-- 'dirs': Direction toward which water flow
-			-- 'dirs2': Directions from which water comes
 			dirs[i] = d
-			if d == 0 then -- If water can't flow from this node, add it to the list of singular nodes that will be resolved later
+			if d == 0 then
 				singular[#singular+1] = i
 			elseif d == 1 then
 				dirs2[i+X] = dirs2[i+X] + 1
@ -100,39 +90,30 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
 		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)
+	-- Compute basins and links
 	local nbasins = #singular
 	local basin_id = {}
 	local links = {}
 	local basin_links

-	-- Function to analyse a link between two nodes
 	local function add_link(i1, i2, b1, isY)
-		-- i1, i2: coordinates of two nodes
-		-- b1: basin that contains i1
-		-- isY: whether the link is in Y direction
 		local b2
-		-- Note that basin number #0 represents the outside of the map; or if the coordinate is inside the map, means that the basin number is uninitialized.
-		if i2 == 0 then -- If outside the map
+		if i2 == 0 then
 			b2 = 0
 		else
 			b2 = basin_id[i2]
-			if b2 == 0 then -- If basin of i2 is not already computed, skip
+			if b2 == 0 then
 				return
 			end
 		end
-		if b2 ~= b1 then -- If these two nodes don't belong to the same basin, we have found a link between two adjacent basins
-			local elev = i2 == 0 and dem[i1] or mmax(dem[i1], dem[i2]) -- Elevation of the highest of the two sides of the link (or only i1 if b2 is map outside)
+		if b2 ~= b1 then
+			local elev = i2 == 0 and dem[i1] or mmax(dem[i1], dem[i2])
 			local l2 = basin_links[b2]
 			if not l2 then
 				l2 = {}
 				basin_links[b2] = l2
 			end
-			if not l2.elev or l2.elev > elev then -- If this link is lower than the lowest registered link between these two basins, register it as the new lowest pass
+			if not l2.elev or l2.elev > elev then
 				l2.elev = elev
 				l2.i = mmax(i1,i2)
 				l2.is_y = isY
@ -145,48 +126,51 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
 	for i=1, X*Y do
 		basin_id[i] = 0
 	end
-
+	--for ib, s in ipairs(singular) do
 	for ib=1, nbasins do
-		-- Here we will recursively search upstream from the singular node to determine its drainage basin
-		local queue = {singular[ib]} -- Start with the singular node, then this queue will be filled with water donors neighbours
+		--local s = singular[ib]
+		local queue = {singular[ib]}
 		basin_links = {}
 		links[#links+1] = basin_links
+		--tinsert(links, basin_links)
 		while #queue > 0 do
 			local i = tremove(queue)
 			basin_id[i] = ib
-			local d = dirs2[i] -- Get the directions water is coming from
+			local d = dirs2[i]

-			-- Iterate through the 4 directions
-			if d >= 8 then -- River coming from the East
+			if d >= 8 then -- River coming from East
 				d = d - 8
 				queue[#queue+1] = i+1
-			-- If no river is coming from the East, we might be at the limit of two basins, thus we need to test adjacency.
+				--tinsert(queue, i+X)
 			elseif i%X > 0 then
 				add_link(i, i+1, ib, false)
-			else -- If the eastern neighbour is outside the map
+			else
 				add_link(i, 0, ib, false)
 			end

-			if d >= 4 then -- River coming from the South
+			if d >= 4 then -- River coming from South
 				d = d - 4
 				queue[#queue+1] = i+X
+				--tinsert(queue, i+1)
 			elseif i <= X*(Y-1) then
 				add_link(i, i+X, ib, true)
 			else
 				add_link(i, 0, ib, true)
 			end

-			if d >= 2 then -- River coming from the West
+			if d >= 2 then -- River coming from West
 				d = d - 2
 				queue[#queue+1] = i-1
+				--tinsert(queue, i-X)
 			elseif i%X ~= 1 then
 				add_link(i, i-1, ib, false)
 			else
 				add_link(i, 0, ib, false)
 			end

-			if d >= 1 then -- River coming from the North
+			if d >= 1 then -- River coming from North
 				queue[#queue+1] = i-X
+				--tinsert(queue, i-1)
 			elseif i > X then
 				add_link(i, i-X, ib, true)
 			else
@ -202,7 +186,7 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
 		nlinks[i] = 0
 	end

-	-- Iterate through pairs of adjacent basins, and make the links reciprocal
+	--for ib1, blinks in ipairs(links) do
 	for ib1=1, #links do
 		for ib2, link in pairs(links[ib1]) do
 			if ib2 < ib1 then
@ -213,15 +197,6 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
 		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 = {}
 	for i, n in pairs(nlinks) do
 		if n <= 8 then
@ -231,8 +206,6 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are

 	local basin_graph = {}
 	for n=1, nbasins do
-		-- Iterate in lowlevel but its contents may change during the loop
-		-- 'next' called with only one argument always returns an element if table is not empty
 		local b1, lnk1 = next(lowlevel)
 		lowlevel[b1] = nil

@ -240,7 +213,6 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
 		local lowest = math.huge
 		local lnk1 = links[b1]
 		local i = 0
-		-- Look for lowest link
 		for bn, bdata in pairs(lnk1) do
 			i = i + 1
 			if bdata.elev < lowest then
@ -249,7 +221,7 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
 			end
 		end

-		-- Add link to the graph, in both directions
+		-- Add link to the graph
 		local bound = lnk1[b2]
 		local bb1, bb2 = bound[1], bound[2]
 		if not basin_graph[bb1] then
@ -267,7 +239,6 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
 		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
 			lowlevel[b2] = lnk2
 		end
@ -276,32 +247,25 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
 			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 lnkn[b2] then
+				nlinks[bn] = nlinks[bn] - 1
 				if nlinks[bn] == 8 then
 					lowlevel[bn] = lnkn
 				end
-			else -- If bn was linked to b1 but not to b2
-				nlinks[b2] = nlinks[b2] + 1 -- Then b2 is gaining a link to bn because of the merge
+			else
+				nlinks[b2] = nlinks[b2] + 1
 				if nlinks[b2] == 9 then
 					lowlevel[b2] = nil
 				end
 			end

-			if not lnkn[b2] or lnkn[b2].elev > bdata.elev then -- If the link b1-bn will become the new lowest link between b2 and bn, redirect the link to b2
+			if not lnkn[b2] or lnkn[b2].elev > bdata.elev then
 				lnkn[b2] = bdata
 				lnk2[bn] = bdata
 			end
 		end
 	end

-	--------------------------------------------------------------
-	-- 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] = -math.huge}
 	local basin_lake = {}
 	for n=1, nbasins do
@ -309,17 +273,15 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
 	end
 	local reverse = {3, 4, 1, 2, [0]=0}
 	for n=1, nbasins do
-		local b1, elev1 = next(queue) -- Pop from queue
+		local b1, elev1 = next(queue)
 		queue[b1] = nil
 		basin_lake[b1] = elev1
-		-- Iterate through b1's neighbours (according to the spanning tree)
 		for b2, bound in pairs(basin_graph[b1]) do
 			-- 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)
+			local i = bound.i
+			local dir = bound.is_y and 3 or 4
 			if basin_id[i] ~= b2 then
 				dir = dir - 2
-				-- Coordinate 'i' refers to the side of the link with the highest X/Y position. In case it is in the wrong basin, take the other side by decrementing by one row/column.
 				if bound.is_y then
 					i = i - X
 				else
@ -329,12 +291,8 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
 				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
@ -344,36 +302,26 @@ local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are
 				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
+			until dir == 0
+			-- Add b2 into the queue
 			queue[b2] = 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.
 	waterq = waterq or {}
 	local X, Y = dirs.X, dirs.Y
+	--local tinsert = table.insert

 	local ndonors = {}
 	local waterq = {X=X, Y=Y}
@ -382,7 +330,7 @@ local function accumulate(dirs, waterq)
 		waterq[i] = 1
 	end

-	-- Calculate the number of direct donors
+	--for i1, dir in ipairs(dirs) do
 	for i1=1, X*Y do
 		local i2
 		local dir = dirs[i1]
@ -401,12 +349,10 @@ local function accumulate(dirs, waterq)
 	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
@ -417,12 +363,9 @@ local function accumulate(dirs, waterq)
 				elseif dir == 4 then
 					i2 = i2 - 1
 				end
-				-- Increment the water quantity of i2
 				w = w + waterq[i2]
 				waterq[i2] = w

-				-- Stop on an unresolved confluence (node with >1 donors) and decrease the number of remaining donors
-				-- When the ndonors of a confluence has decreased to 1, it means that its water quantity has already been incremented by its tributaries, so it can be resolved like a standard river section. However, do not decrease ndonors to zero to avoid considering it as a source.
 				if ndonors[i2] > 1 then
 					ndonors[i2] = ndonors[i2] - 1
 					break
--- a/terrainlib_lua/twist.lua
+++ b/terrainlib_lua/twist.lua
@ -1,4 +1,4 @@
-- twist.lua
+-- bounds.lua

 local function get_bounds(dirs, rivers)
 	local X, Y = dirs.X, dirs.Y
--- a/view.py
+++ b/view.py
@ -11,10 +11,12 @@ try:
        import colorcet as cc
        cmap1 = cc.cm.CET_L11
        cmap2 = cc.cm.CET_L12
+        cmap3 = cc.cm.CET_L6.reversed()
    except ImportError: # No module colorcet
        import matplotlib.cm as cm
        cmap1 = cm.summer
-        cmap2 = cm.Blues
+        cmap2 = cm.ocean.reversed()
+        cmap3 = cm.Blues
 except ImportError: # No module matplotlib
    has_matplotlib = False

@ -24,10 +26,12 @@ if has_matplotlib:
        water = np.maximum(lakes_sea - dem, 0)
        max_elev = dem.max()
        max_depth = water.max()
+        max_lake_depth = lakes.max()

        ls = mcl.LightSource(azdeg=315, altdeg=45)
        norm_ground = plt.Normalize(vmin=sea_level, vmax=max_elev)
        norm_sea = plt.Normalize(vmin=0, vmax=max_depth)
+        norm_lake = plt.Normalize(vmin=0, vmax=max_lake_depth)
        rgb = ls.shade(dem, cmap=cmap1, vert_exag=1/scale, blend_mode='soft', norm=norm_ground)

        (X, Y) = dem.shape
@ -37,13 +41,23 @@ if has_matplotlib:
            extent = (-0.5*scale, (Y-0.5)*scale, -0.5*scale, (X-0.5)*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')
+        lakes_alpha = ((lakes_sea - np.maximum(dem,sea_level)) > 0).astype('u1')
+        # plt.imshow(np.flipud(water), alpha=np.flipud(alpha), cmap=cmap2, extent=extent, vmin=0, vmax=max_depth, interpolation='antialiased')
+        plt.imshow(np.flipud(water), alpha=np.flipud(alpha), cmap=cmap3, extent=extent, vmin=0, vmax=max_depth, interpolation='antialiased')
+        plt.imshow(np.flipud(water), alpha=np.flipud(lakes_alpha), cmap=cmap2, extent=extent, vmin=0, vmax=max_depth, interpolation='antialiased')

        sm1 = plt.cm.ScalarMappable(cmap=cmap1, norm=norm_ground)
        plt.colorbar(sm1).set_label('Elevation')

-        sm2 = plt.cm.ScalarMappable(cmap=cmap2, norm=norm_sea)
-        plt.colorbar(sm2).set_label('Water depth')
+        sm2 = plt.cm.ScalarMappable(cmap=cmap2, norm=norm_lake)
+        cb2 = plt.colorbar(sm2)
+        cb2.ax.invert_yaxis()
+        cb2.set_label('Lake Depth')
+
+        sm3 = plt.cm.ScalarMappable(cmap=cmap3, norm=norm_sea)
+        cb3 = plt.colorbar(sm3)
+        cb3.ax.invert_yaxis()
+        cb3.set_label('Ocean Depth')

        plt.xlabel('X')
        plt.ylabel('Z')
@ -84,9 +98,10 @@ def stats(dem, lakes, scale=1):
    lake_surface = lake.sum()

    print('---   General    ---')
-    print('Grid size:    {:5d}x{:5d}'.format(dem.shape[0], dem.shape[1]))
+    print('Grid size (dem):     {:5d}x{:5d}'.format(dem.shape[0], dem.shape[1]))
+    print('Grid size (lakes):   {:5d}x{:5d}'.format(lakes.shape[0], lakes.shape[1]))
    if scale > 1:
-        print('Map size:     {:5d}x{:5d}'.format(int(dem.shape[0]*scale), int(dem.shape[1]*scale)))
+        print('Map size:            {:5d}x{:5d}'.format(int(dem.shape[0]*scale), int(dem.shape[1]*scale)))
    print()
    print('---   Surfaces   ---')
    print('Continents:        {:6.2%}'.format(continent_surface/surface))
@ -100,3 +115,5 @@ def stats(dem, lakes, scale=1):
    print('Mean continent elev: {:4.0f}'.format((dem*continent).sum()/continent_surface))
    print('Lowest elevation:    {:4.0f}'.format(dem.min()))
    print('Highest elevation:   {:4.0f}'.format(dem.max()))
+    print()
+