Fix river shape in confluences (less sharp riverbeds when a small rivers joins a big one)

Also cleaned and commented the code

heightmap.lua

@@ -46,50 +46,57 @@ local function heightmaps(minp, maxp)
 				-- = 0: on riverbank
 				-- > 0: inside river
 				local depth_factors = {
-					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,
+					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)
 				}
 
-				-- Find the maximal depth factor and determine to which river it belongs
-				local depth_factor_max = 0
+				-- 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
 				local imax = 0
 				for i=1, 8 do
-					if depth_factors[i] > depth_factor_max then
-						depth_factor_max = depth_factors[i]
+					if depth_factors[i] > dpmax then
+						dpmax = depth_factors[i]
 						imax = i
 					end
 				end
 
-				-- Transform the coordinates to have xf and zf = 0 or 1 in rivers (to avoid rivers having lateral slope and to accomodate the surrounding smoothly)
-				if imax == 0 then
-					local x0 = max(r_west, c_NW-zf, zf-c_SW)
-					local x1 = min(r_east, c_NE+zf, c_SE-zf)
-					local z0 = max(r_north, c_NW-xf, xf-c_NE)
-					local z1 = min(r_south, c_SW+xf, c_SE-xf)
-					xf = (xf-x0) / (x1-x0)
-					zf = (zf-z0) / (z1-z0)
-				elseif imax == 1 then
-					xf = 0
-				elseif imax == 2 then
-					zf = 0
-				elseif imax == 3 then
-					xf = 1
-				elseif imax == 4 then
-					zf = 1
-				elseif imax == 5 then
-					xf, zf = 0, 0
-				elseif imax == 6 then
-					xf, zf = 1, 0
-				elseif imax == 7 then
-					xf, zf = 1, 1
-				elseif imax == 8 then
-					xf, zf = 0, 1
+				-- 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
 				end
 
 				-- Determine elevation by interpolation
@@ -99,12 +106,12 @@ local function heightmaps(minp, maxp)
 					vdem[2],
 					vdem[3],
 					vdem[4],
-					xf, zf
+					xfc, zfc
 				))
 
 				-- Spatial gradient of the interpolation
-				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 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 lake_id = 0
 				if slope_x then
 					if slope_z then
@@ -121,8 +128,8 @@ local function heightmaps(minp, maxp)
 				end
 				local lake_height = max(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)
+				if imax > 0 and dpmax > 0 then
+					terrain_height = min(max(lake_height, sea_level) - floor(1+dpmax*riverbed_slope), terrain_height)
 				end
 
 				terrain_height_map[i] = terrain_height

polygons.lua

@@ -236,15 +236,15 @@ local function make_polygons(minp, maxp)
 				end
 			end
 
-			polygon.river_corners = {riverA, 1-riverB, 2-riverC, 1-riverD}
+			polygon.river_corners = {riverA, riverB, riverC, 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 = 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)
+			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)
 
 			polygon.rivers = {river_west, river_north, river_east, river_south}
 		end