From c99b8338e0d2daefba1a19f4f481aaf45a3c87ac Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ga=C3=ABl=20C?= Date: Tue, 1 Jun 2021 19:07:09 +0200 Subject: [PATCH] Lua Terrainlib: added first Lua files for erosion and flow routing Tested, but not linked with the mod, yet. --- terrainlib_lua/erosion.lua | 218 ++++++++++++++++++ terrainlib_lua/gaussian.lua | 111 +++++++++ terrainlib_lua/rivermapper.lua | 402 +++++++++++++++++++++++++++++++++ 3 files changed, 731 insertions(+) create mode 100644 terrainlib_lua/erosion.lua create mode 100644 terrainlib_lua/gaussian.lua create mode 100644 terrainlib_lua/rivermapper.lua diff --git a/terrainlib_lua/erosion.lua b/terrainlib_lua/erosion.lua new file mode 100644 index 0000000..b7eeb67 --- /dev/null +++ b/terrainlib_lua/erosion.lua @@ -0,0 +1,218 @@ +-- erosion.lua + +local function erode(model, time) + --local tinsert = table.insert + local mmin, mmax = math.min, math.max + local dem = model.dem + local dirs = model.dirs + local lakes = model.lakes + local rivers = model.rivers + local sea_level = model.params.sea_level + local K = model.params.K + local m = model.params.m + local X, Y = dem.X, dem.Y + local scalars = type(K) == "number" and type(m) == "number" + + local erosion_time + if model.params.variable_erosion then + erosion_time = {} + else + erosion_time = model.erosion_time or {} + end + + if scalars then + for i=1, X*Y do + local etime = 1 / (K*rivers[i]^m) + erosion_time[i] = etime + lakes[i] = mmax(lakes[i], dem[i], sea_level) + end + else + for i=1, X*Y do + local etime = 1 / (K[i]*rivers[i]^m[i]) + erosion_time[i] = etime + lakes[i] = mmax(lakes[i], dem[i], sea_level) + end + end + + for i=1, X*Y do + local iw = i + local remaining = time + local new_elev + while true do + local inext = iw + local d = dirs[iw] + + if d == 0 then + new_elev = lakes[iw] + break + elseif d == 1 then + inext = iw+1 + elseif d == 2 then + inext = iw+X + elseif d == 3 then + inext = iw-1 + elseif d == 4 then + inext = iw-X + end + + local etime = erosion_time[iw] + if remaining <= etime then + local c = remaining / etime + new_elev = (1-c) * lakes[iw] + c * lakes[inext] + break + end + + remaining = remaining - etime + iw = inext + end + + dem[i] = mmin(dem[i], new_elev) + end +end + +local function diffuse(model, time) + 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 + 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 + 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 + -- TODO Test this +end + +local rivermapper, gaussian +if minetest then + rivermapper = dofile(minetest.get_modpath(minetest.get_current_modname()) .. '/terrainlib_lua/rivermapper.lua') + gaussian = dofile(minetest.get_modpath(minetest.get_current_modname()) .. '/terrainlib_lua/gaussian.lua') +else + rivermapper = dofile('rivermapper.lua') + gaussian = dofile('gaussian.lua') +end + +local function flow(model) + model.dirs, model.lakes = rivermapper.flow_routing(model.dem, model.dirs, model.lakes, 'semirandom') + model.rivers = rivermapper.accumulate(model.dirs, model.rivers) +end + +local function uplift(model, time) + local dem = model.dem + local X, Y = dem.X, dem.Y + local uplift_rate = model.params.uplift + if type(uplift_rate) == "number" then + local uplift_total = uplift_rate * time + for i=1, X*Y do + dem[i] = dem[i] + uplift_total + end + else + for i=1, X*Y do + dem[i] = dem[i] + uplift_rate[i]*time + end + end +end + +local function noise(model, time) + local random = math.random + local dem = model.dem + local noise_depth = model.params.noise * 2 * time + local X, Y = dem.X, dem.Y + for i=1, X*Y do + dem[i] = dem[i] + (random()-0.5) * noise_depth + end +end + +local function define_isostasy(model, ref, link) + 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 + + return ref2 +end + +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] + end + + gaussian.gaussian_blur_approx(temp, model.params.compensation_radius, 4) + for i=1, X*Y do + dem[i] = dem[i] + temp[i] + end +end + +local evol_model_mt = { + erode = erode, + diffuse = diffuse, + flow = flow, + uplift = uplift, + noise = noise, + isostasy = isostasy, + define_isostasy = define_isostasy, +} + +evol_model_mt.__index = evol_model_mt + +local defaults = { + K = 1, + m = 0.5, + d = 1, + variable_erosion = false, + sea_level = 0, + uplift = 10, + noise = 0.001, + compensation_radius = 50, +} + +local function EvolutionModel(params) + params = params or {} + local o = {params = params} + for k, v in pairs(defaults) do + if params[k] == nil then + params[k] = v + end + end + o.dem = params.dem + return setmetatable(o, evol_model_mt) +end + +return EvolutionModel diff --git a/terrainlib_lua/gaussian.lua b/terrainlib_lua/gaussian.lua new file mode 100644 index 0000000..f08399b --- /dev/null +++ b/terrainlib_lua/gaussian.lua @@ -0,0 +1,111 @@ +-- gaussian.lua + +local function get_box_size(sigma, n) + local v = sigma^2 / n + --print('v: '..v) + local r_ideal = ((12*v + 1) ^ 0.5 - 1) / 2 + --print('r_ideal: '..r_ideal) + local r_down = math.floor(r_ideal) + --print('r_down: '..r_down) + local r_up = math.ceil(r_ideal) + --print('r_up: '..r_up) + local v_down = ((2*r_down+1)^2 - 1) / 12 + --print('v_down: '..v_down) + local v_up = ((2*r_up+1)^2 - 1) / 12 + --print('v_up: '..v_up) + local m_ideal = (v - v_down) / (v_up - v_down) * n + --print('m_ideal: '..m_ideal) + local m = math.floor(m_ideal+0.5) + --print('m: '..m) + + local sizes = {} + for i=1, n do + sizes[i] = i<=m and 2*r_up+1 or 2*r_down+1 + end + --print('sizes: '..table.concat(sizes, ', ')) + + return sizes +end + +local function box_blur_1d(map, size, first, incr, len, map2) + local n = math.ceil(size/2) + first = first or 1 + incr = incr or 1 + len = len or math.floor((#map-first)/incr)+1 + local last = first + (len-1)*incr + + local nth = first+(n-1)*incr + local sum = 0 + for i=first, nth, incr do + if i == first then + sum = sum + map[i] + else + sum = sum + 2*map[i] + end + end + + local i_left = nth + local incr_left = -incr + local i_right = nth + local incr_right = incr + + map2 = map2 or {} + for i=first, last, incr do + map2[i] = sum / size + i_right = i_right + incr_right + sum = sum - map[i_left] + map[i_right] + i_left = i_left + incr_left + + if i_left == first then + incr_left = incr + end + if i_right == last then + incr_right = -incr + end + end + + return map2 +end + +local function box_blur_2d(map1, size, map2) + local X, Y = map1.X, map1.Y + map2 = map2 or {} + for y=1, Y do + box_blur_1d(map1, size, (y-1)*X+1, 1, X, map2) + end + for x=1, X do + box_blur_1d(map2, size, x, X, Y, map1) + end + + 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) + for i=1, n do + box_blur_2d(map, sizes[i], map2) + end + return map +end + +return { + get_box_size = get_box_size, + box_blur_1d = box_blur_1d, + box_blur_2d = box_blur_2d, + gaussian_blur_approx = gaussian_blur_approx, +} diff --git a/terrainlib_lua/rivermapper.lua b/terrainlib_lua/rivermapper.lua new file mode 100644 index 0000000..ecf0203 --- /dev/null +++ b/terrainlib_lua/rivermapper.lua @@ -0,0 +1,402 @@ +-- rivermapper.lua + +local function flow_local_semirandom(plist) + local sum = 0 + for i=1, #plist do + 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 + r = r - p + end + return 0 +end + +local flow_methods = { + semirandom = flow_local_semirandom, +} + +local function flow_routing(dem, dirs, lakes, method) + method = method or 'semirandom' + flow_local = flow_methods[method] or flow_local_semirandom + + dirs = dirs or {} + lakes = lakes or {} + + -- Localize for performance + --local tinsert = table.insert + local tremove = table.remove + local mmax = math.max + + local X, Y = dem.X, dem.Y + dirs.X = X + dirs.Y = Y + lakes.X = X + lakes.Y = Y + --print(X, Y) + local i = 1 + local dirs2 = {} + for i=1, X*Y do + dirs2[i] = 0 + end + + local singular = {} + for y=1, Y do + for x=1, X do + local zi = dem[i] + local plist = { + x1 and mmax(zi-dem[i-1], 0) or 0, -- Westward + y>1 and mmax(zi-dem[i-X], 0) or 0, -- Northward + } + + local d = flow_local(plist) + dirs[i] = d + if d == 0 then + singular[#singular+1] = i + elseif d == 1 then + dirs2[i+1] = dirs2[i+1] + 1 + elseif d == 2 then + dirs2[i+X] = dirs2[i+X] + 2 + elseif d == 3 then + dirs2[i-1] = dirs2[i-1] + 4 + elseif d == 4 then + dirs2[i-X] = dirs2[i-X] + 8 + end + i = i + 1 + end + end + + -- Compute basins and links + local nbasins = #singular + print(nbasins) + local basin_id = {} + local links = {} + local basin_links + + local function add_link(i1, i2, b1, isY) + local b2 + if i2 == 0 then + b2 = 0 + else + b2 = basin_id[i2] + if b2 == 0 then + return + end + end + 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 + l2.elev = elev + l2.i = mmax(i1,i2) + l2.is_y = isY + l2[1] = b2 + l2[2] = b1 + end + end + end + + for i=1, X*Y do + basin_id[i] = 0 + end + --for ib, s in ipairs(singular) do + for ib=1, nbasins do + --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] + + if d >= 8 then -- River coming from South + d = d - 8 + queue[#queue+1] = i+X + --tinsert(queue, i+X) + elseif i <= X*(Y-1) then + add_link(i, i+X, ib, true) + else + add_link(i, 0, ib, true) + end + + if d >= 4 then -- River coming from East + d = d - 4 + queue[#queue+1] = i+1 + --tinsert(queue, i+1) + elseif i%X > 0 then + add_link(i, i+1, ib, false) + else + add_link(i, 0, ib, false) + end + + if d >= 2 then -- River coming from North + d = d - 2 + queue[#queue+1] = i-X + --tinsert(queue, i-X) + elseif i > X then + add_link(i, i-X, ib, true) + else + add_link(i, 0, ib, true) + end + + if d >= 1 then -- River coming from West + queue[#queue+1] = i-1 + --tinsert(queue, i-1) + elseif i%X ~= 1 then + add_link(i, i-1, ib, false) + else + add_link(i, 0, ib, false) + end + end + end + dirs2 = nil + + links[0] = {} + local nlinks = {} + for i=0, nbasins do + nlinks[i] = 0 + end + + --for ib1, blinks in ipairs(links) do + for ib1=1, #links do + for ib2, link in pairs(links[ib1]) do + if ib2 < ib1 then + links[ib2][ib1] = link + nlinks[ib1] = nlinks[ib1] + 1 + nlinks[ib2] = nlinks[ib2] + 1 + end + end + end + + local lowlevel = {} + for i, n in pairs(nlinks) do + if n <= 8 then + lowlevel[i] = links[i] + end + end + + basin_graph = {} + for n=1, nbasins do + --print(n, nbasins) + local b1, lnk1 = next(lowlevel) + lowlevel[b1] = nil + + local b2 + local lowest = math.huge + local lnk1 = links[b1] + local i = 0 + --print('Scanning basin '..b1) + for bn, bdata in pairs(lnk1) do + --print('- Link '..bn) + i = i + 1 + if bdata.elev < lowest then + lowest = bdata.elev + b2 = bn + end + end + --print('Number of links: '..i..' vs '..nlinks[b1]) + + -- Add link to the graph + local bound = lnk1[b2] + local bb1, bb2 = bound[1], bound[2] + if not basin_graph[bb1] then + basin_graph[bb1] = {} + end + if not basin_graph[bb2] then + basin_graph[bb2] = {} + end + basin_graph[bb1][bb2] = bound + basin_graph[bb2][bb1] = bound + --if bb1 == 0 then + -- print(bb2) + --elseif bb2 == 0 then + -- print(bb1) + --end + + -- Merge basin b1 into b2 + --print("Merging "..b1.." into "..b2) + local lnk2 = links[b2] + -- First, remove the link between b1 and b2 + lnk1[b2] = nil + lnk2[b1] = nil + nlinks[b2] = nlinks[b2] - 1 + --print('Decreasing link count of '..b2..' ('..nlinks[b2]..')') + if nlinks[b2] == 8 then + --print('Added to lowlevel') + lowlevel[b2] = lnk2 + end + --print('Scanning neighbourg of '..b1..' to fix links') + -- Look for basin 1's neighbours, and add them to basin 2 if they have a lower pass + for bn, bdata in pairs(lnk1) do + --print('- Neighbour '..bn) + local lnkn = links[bn] + lnkn[b1] = nil + + if lnkn[b2] then + nlinks[bn] = nlinks[bn] - 1 + --print('Decreasing link count of '..bn..' ('..nlinks[bn]..')') + if nlinks[bn] == 8 then + --print('Added to lowlevel') + lowlevel[bn] = lnkn + end + else + nlinks[b2] = nlinks[b2] + 1 + --print('Increasing link count of '..b2..' ('..nlinks[b2]..')') + if nlinks[b2] == 9 then + --print('Removed from lowlevel') + lowlevel[b2] = nil + end + end + + if not lnkn[b2] or lnkn[b2].elev > bdata.elev then + --print(' - Redirecting link') + lnkn[b2] = bdata + lnk2[bn] = bdata + end + end + end + + local queue = {[0] = -math.huge} + local basin_lake = {} + for n=1, nbasins do + basin_lake[n] = 0 + end + local reverse = {3, 4, 1, 2, [0]=0} + for n=1, nbasins do + --print(n, nbasins) + local b1, elev1 = next(queue) + queue[b1] = nil + basin_lake[b1] = elev1 + --print('Scanning basin '..b1) + for b2, bound in pairs(basin_graph[b1]) do + --print('Flow '..b2..' into '..b1) + -- Make b2 flow into b1 + local i = bound.i + local dir = bound.is_y and 4 or 3 + --print(basin_id[i]) + if basin_id[i] ~= b2 then + dir = dir - 2 + if bound.is_y then + i = i - X + else + i = i - 1 + end + elseif b1 == 0 then + dir = 0 + end + --print(basin_id[i]) + --print('Reversing directions') + repeat + dir, dirs[i] = dirs[i], dir + if dir == 1 then + i = i + 1 + elseif dir == 2 then + i = i + X + elseif dir == 3 then + i = i - 1 + elseif dir == 4 then + i = i - X + end + dir = reverse[dir] + until dir == 0 + -- Add b2 into the queue + queue[b2] = mmax(elev1, bound.elev) + basin_graph[b2][b1] = nil + end + basin_graph[b1] = nil + end + + for i=1, X*Y do + lakes[i] = basin_lake[basin_id[i]] + end + + return dirs, lakes +end + +local function accumulate(dirs, waterq) + waterq = waterq or {} + local X, Y = dirs.X, dirs.Y + --local tinsert = table.insert + + local ndonors = {} + local waterq = {X=X, Y=Y} + for i=1, X*Y do + ndonors[i] = 0 + waterq[i] = 1 + end + + --for i1, dir in ipairs(dirs) do + for i1=1, X*Y do + local i2 + local dir = dirs[i1] + if dir == 1 then + i2 = i1+1 + elseif dir == 2 then + i2 = i1+X + elseif dir == 3 then + i2 = i1-1 + elseif dir == 4 then + i2 = i1-X + end + if i2 then + ndonors[i2] = ndonors[i2] + 1 + end + end + + for i1=1, X*Y do + --print(i1, ndonors[i1]) + if ndonors[i1] == 0 then + local i2 = i1 + local dir = dirs[i2] + local w = waterq[i2] + --print(dir) + while dir > 0 do + if dir == 1 then + i2 = i2 + 1 + elseif dir == 2 then + i2 = i2 + X + elseif dir == 3 then + i2 = i2 - 1 + elseif dir == 4 then + i2 = i2 - X + end + --print('Incrementing '..i2) + w = w + waterq[i2] + waterq[i2] = w + + if ndonors[i2] > 1 then + ndonors[i2] = ndonors[i2] - 1 + break + end + dir = dirs[i2] + end + end + end + + return waterq +end + +return { + flow_routing = flow_routing, + accumulate = accumulate, + flow_methods = flow_methods, +}