diff --git a/init.lua b/init.lua
index 957fc64..1e21aa4 100644
--- a/init.lua
+++ b/init.lua
@@ -45,55 +45,63 @@ 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
+				-- Load river width on 4 edges and corners
 				local r_west, r_north, r_east, r_south = unpack(poly.rivers)
-				if xf > r_east then
-					is_river = true
-					depth_factor = xf-r_east
-					xf = 1
-				elseif xf < r_west then
-					is_river = true
-					depth_factor = r_west-xf
-					xf = 0
-				end
-				if zf > r_south then
-					is_river = true
-					depth_factor = zf-r_south
-					zf = 1
-				elseif zf < r_north then
-					is_river = true
-					depth_factor = r_north-zf
-					zf = 0
-				end
+				local c_NW, c_NE, c_SE, c_SW = unpack(poly.river_corners)
 
-				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
-						depth_factor = c_NW-xf-zf
-						xf, zf = 0, 0
-					elseif 1-xf+zf < c_NE then
-						is_river = true
-						depth_factor = c_NE-1+xf-zf
-						xf, zf = 1, 0
-					elseif 2-xf-zf < c_SE then
-						is_river = true
-						depth_factor = c_SE-2+xf+zf
-						xf, zf = 1, 1
-					elseif xf+1-zf < c_SW then
-						is_river = true
-						depth_factor = c_SW-xf-1+zf
-						xf, zf = 0, 1
+				-- Calculate the depth factor for each edge and corner.
+				-- Depth factor:
+				-- < 0: outside river
+				-- = 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,
+					c_NE+xf-zf-1,
+					c_SE+xf+zf-2,
+					c_SW-xf+zf-1,
+				}
+
+				-- 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] >= depth_factor_max then
+						depth_factor_max = depth_factors[i]
+						imax = i
 					end
 				end
 
-				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)
+				-- 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, c_SW+zf-1)
+					local x1 = math.min(r_east, 1-c_NE+zf, 2-c_SE-zf)
+					local z0 = math.max(r_north, c_NW-xf, c_NE+xf-1)
+					local z1 = math.min(r_south, 1-c_SW+xf, 2-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 = math.floor(0.5+interp(
 					vdem[1],
@@ -104,8 +112,8 @@ local function generate(minp, maxp, seed)
 				))
 
 				local lake_height = math.max(math.floor(poly.lake), terrain_height)
-				if is_river then
-					terrain_height = math.min(math.max(lake_height, sea_level) - math.floor(1+depth_factor*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
 				local is_lake = lake_height > terrain_height
 				local ivm = a:index(x, minp.y-1, z)
diff --git a/polygons.lua b/polygons.lua
index 4824327..afa9f3f 100644
--- a/polygons.lua
+++ b/polygons.lua
@@ -153,17 +153,6 @@ local function make_polygons(minp, maxp)
 			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)
 
-			-- Only if opposite rivers overlap (should be rare)
-			if river_west > river_east then
-				local mean = (river_west + river_east) / 2
-				river_west = mean
-				river_east = mean
-			end
-			if river_north > river_south then
-				local mean = (river_north + river_south) / 2
-				river_north = mean
-				river_south = mean
-			end
 			polygon.rivers = {river_west, river_north, river_east, river_south}
 		end
 	end