Lua Terrainlib: added first Lua files for erosion and flow routing

Tested, but not linked with the mod, yet.
This commit is contained in:
Gaël C 2021-06-01 19:07:09 +02:00
parent c175f2bbf7
commit c99b8338e0
3 changed files with 731 additions and 0 deletions

218
terrainlib_lua/erosion.lua Normal file
View File

@ -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

111
terrainlib_lua/gaussian.lua Normal file
View File

@ -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,
}

View File

@ -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 = {
x<X and mmax(zi-dem[i+1], 0) or 0, -- Eastward
y<Y and mmax(zi-dem[i+X], 0) or 0, -- Southward
x>1 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,
}