mapgen_rivers/heightmap.lua

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 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)
	local v1 = v11*xf + v10*(1-xf)
	return v1*zf + v0*(1-zf)
end

local function heightmaps(minp, maxp)

	local polygons = make_polygons(minp, maxp)
	local incr = maxp.z-minp.z+1

	local terrain_height_map = {}
	local lake_height_map = {}

	local i = 1
	for z=minp.z, maxp.z do
		for x=minp.x, maxp.x do
			local poly = polygons[i]

			if poly then
				local xf, zf = transform_quadri(poly.x, poly.z, x, z)
				local i00, i01, i11, i10 = unpk(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)

				-- 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 , -- 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, 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] > dpmax then
						dpmax = 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
				end

				-- Determine elevation by interpolation
				local vdem = poly.dem
				local terrain_height = floor(0.5+interp(
					vdem[1],
					vdem[2],
					vdem[3],
					vdem[4],
					xfc, zfc
				))

				-- 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 lake_id = 0
				if slope_x then
					if slope_z then
						lake_id = 3
					else
						lake_id = 2
					end
				else
					if slope_z then
						lake_id = 4
					else
						lake_id = 1
					end
				end
				local lake_height = max(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)
				end

				terrain_height_map[i] = terrain_height
				lake_height_map[i] = lake_height
			else
				terrain_height_map[i] = MAP_BOTTOM
				lake_height_map[i] = MAP_BOTTOM
			end
			i = i + 1
		end
	end

	return terrain_height_map, lake_height_map
end

return heightmaps