Start implementing snowball erosion algorithm.

There's still a long way to go - we're only just getting warmed up!
This commit is contained in:
Starbeamrainbowlabs 2020-08-20 01:53:26 +01:00
parent 20ccf321ce
commit cdba38d37d
Signed by: sbrl
GPG Key ID: 1BE5172E637709C2
6 changed files with 173 additions and 41 deletions

View File

@ -11,6 +11,7 @@ dofile(worldeditadditions.modpath.."/utils/strings.lua")
dofile(worldeditadditions.modpath.."/utils/numbers.lua")
dofile(worldeditadditions.modpath.."/utils/nodes.lua")
dofile(worldeditadditions.modpath.."/utils/tables.lua")
dofile(worldeditadditions.modpath.."/utils/terrain.lua")
dofile(worldeditadditions.modpath.."/utils/raycast_adv.lua") -- For the farwand
dofile(worldeditadditions.modpath.."/utils.lua")

View File

@ -4,6 +4,71 @@ Note that this *mutates* the given heightmap.
@source https://jobtalle.com/simulating_hydraulic_erosion.html
]]--
function worldeditadditions.erode.snowballs(heightmap, heightmap_size, erosion_params)
function worldeditadditions.erode.snowballs(heightmap, heightmap_size, params)
-- Apply the default settings
worldeditadditions.table_apply({
rate_deposit = 0.03,
rate_erosion = 0.04,
friction = 0.07,
speed = 0.15,
radius = 0.8,
snowball_max_steps = 80,
scale_iterations = 0.04,
drops_per_cell = 0.4,
snowball_count = 50000
}, params)
local normals = worldeditadditions.calculate_normals(heightmap, heightmap_size)
for i = 1, params.snowball_count do
snowball(
heightmap, normals, heightmap_size,
{ x = math.random() }
)
end
end
local function snowball(heightmap, normalmap, heightmap_size, startpos, params)
local offset = { -- Random jitter - apparently helps to avoid snowballs from entrenching too much
x = (math.random() * 2 - 1) * params.radius,
z = (math.random() * 2 - 1) * params.radius
}
local sediment = 0
local pos = { x = startpos.x, z = startpos.z }
local pos_prev = { x = pos.x, z = pos.z }
local velocity = { x = 0, z = 0 }
local heightmap_length = #heightmap
for i = 1, params.snowball_max_steps do
local hi = math.floor(pos.z+offset.z+0.5)*heightmap_size[1] + math.floor(pos.x+offset.x+0.5)
if hi > heightmap_length then break end
-- Stop if we go out of bounds
if offset.x < 0 or offset.z < 0
or offset.x >= heightmap[1] or offset.z >= heightmap[0] then
break
end
local step_deposit = sediment * params.rate_deposit * normalmap[hi].z
local step_erode = params.rate_erosion * (1 - normalmap[hi].z) * math.min(1, i*params.scale_iterations)
-- Erode / Deposit, but only if we are on a different node than we were in the last step
if math.floor(pos_prev.x) ~= math.floor(pos.x)
and math.floor(pos_prev.z) ~= math.floor(pos.z) then
heightmap[hi] = heightmap[hi] + (deposit - erosion)
end
velocity.x = params.friction * velocity.x + normalmap[hi].x * params.speed
velocity.z = params.friction * velocity.z + normalmap[hi].y * params.speed
pos_prev.x = pos.x
pos_prev.z = pos.z
pos.x = pos.x + velocity.x
pos.z = pos.z + velocity.z
end
-- Round everything to the nearest int, since you can't really have
-- something like .141592671 of a node
for i,v in ipairs(heightmap) do
heightmap[i] = math.floor(heightmap[i] + 0.5)
end
end

View File

@ -4,12 +4,34 @@ function worldeditadditions.vector.tostring(v)
return "(" .. v.x ..", " .. v.y ..", " .. v.z ..")"
end
-- Calculates the length squared of the given vector.
-- @param v Vector The vector to operate on
-- @return number The length of the given vector squared
function worldeditadditions.vector.lengthsquared(v)
return v.x*v.x + v.y*v.y + v.z*v.z
end
--- Normalises the given vector such that its length is 1.
-- Also known as calculating the unit vector.
-- This method does *not* mutate.
-- @param v Vector The vector to calculate from.
-- @return Vector A new normalised vector.
function worldeditadditions.vector.normalise(v)
local length = math.sqrt(worldeditadditions.lengthsquared(v))
return {
x = x / length,
y = y / length,
z = z / length
}
end
--- Rounds the values in a vector down.
-- Warning: This MUTATES the given vector!
-- @param v Vector The vector to operate on
function worldeditadditions.vector.floor(v)
v.x = math.floor(v.x)
v.y = math.floor(v.y)
v.z = math.floor(v.z)
-- Some vectors are 2d, but on the x / z axes
if v.y then v.y = math.floor(v.y) end
-- Some vectors are 2d
if v.z then v.z = math.floor(v.z) end
end

View File

@ -78,41 +78,3 @@ function worldeditadditions.is_liquidlike(id)
-- If it's not none, then it has to be a liquid as the only other values are source and flowing
return true
end
--- Given a manip object and associates, generates a 2D x/z heightmap.
-- Note that pos1 and pos2 should have already been pushed through
-- worldedit.sort_pos(pos1, pos2) before passing them to this function.
-- @param pos1 Vector Position 1 of the region to operate on
-- @param pos2 Vector Position 2 of the region to operate on
-- @param manip VoxelManip The VoxelManip object.
-- @param area area The associated area object.
-- @param data table The associated data object.
-- @return table The ZERO-indexed heightmap data (as 1 single flat array).
function worldeditadditions.make_heightmap(pos1, pos2, manip, area, data)
-- z y x (in reverse for little-endian machines) is the preferred loop order, but that isn't really possible here
local heightmap = {}
local hi = 0
local changes = { updated = 0, skipped_columns = 0 }
for z = pos1.z, pos2.z, 1 do
for x = pos1.x, pos2.x, 1 do
local found_node = false
-- Scan each column top to bottom
for y = pos2.y+1, pos1.y, -1 do
local i = area:index(x, y, z)
if not worldeditadditions.is_airlike(data[i]) then
-- It's the first non-airlike node in this column
-- Start heightmap values from 1 (i.e. there's at least 1 node in the column)
heightmap[hi] = (y - pos1.y) + 1
found_node = true
break
end
end
if not found_node then heightmap[hi] = -1 end
hi = hi + 1
end
end
return heightmap
end

View File

@ -15,3 +15,13 @@ function worldeditadditions.shallowcopy(orig)
end
return copy
end
--- SHALLOWLY applies the values in source to overwrite the equivalent keys in target.
-- Warning: This function mutates target!
-- @param source table The source to take values from
-- @param target table The target to write values to
function worldeditadditions.table_apply(source, target)
for key, value in pairs(source) do
target[key] = value
end
end

View File

@ -0,0 +1,72 @@
--- Given a manip object and associates, generates a 2D x/z heightmap.
-- Note that pos1 and pos2 should have already been pushed through
-- worldedit.sort_pos(pos1, pos2) before passing them to this function.
-- @param pos1 Vector Position 1 of the region to operate on
-- @param pos2 Vector Position 2 of the region to operate on
-- @param manip VoxelManip The VoxelManip object.
-- @param area area The associated area object.
-- @param data table The associated data object.
-- @return table The ZERO-indexed heightmap data (as 1 single flat array).
function worldeditadditions.make_heightmap(pos1, pos2, manip, area, data)
-- z y x (in reverse for little-endian machines) is the preferred loop order, but that isn't really possible here
local heightmap = {}
local hi = 0
local changes = { updated = 0, skipped_columns = 0 }
for z = pos1.z, pos2.z, 1 do
for x = pos1.x, pos2.x, 1 do
local found_node = false
-- Scan each column top to bottom
for y = pos2.y+1, pos1.y, -1 do
local i = area:index(x, y, z)
if not worldeditadditions.is_airlike(data[i]) then
-- It's the first non-airlike node in this column
-- Start heightmap values from 1 (i.e. there's at least 1 node in the column)
heightmap[hi] = (y - pos1.y) + 1
found_node = true
break
end
end
if not found_node then heightmap[hi] = -1 end
hi = hi + 1
end
end
return heightmap
end
--- Calculates a normal map for the given heightmap.
-- Caution: This method (like worldeditadditions.make_heightmap) works on
-- X AND Z, and NOT x and y. This means that the resulting 3d normal vectors
-- will have the z and y values swapped.
-- @param heightmap table A ZERO indexed flat heightmap. See worldeditadditions.make_heightmap().
-- @param heightmap_size int[] The size of the heightmap in the form [ z, x ]
-- @return Vector[] The calculated normal map, in the same form as the input heightmap. Each element of the array is a 3D Vector (i.e. { x, y, z }) representing a normal.
function worldeditadditions.calculate_normals(heightmap, heightmap_size)
local result = {}
for z = heightmap_size[0], 0, -1 do
for x = heightmap_size[1], 0, -1 do
-- Algorithm ref https://stackoverflow.com/a/13983431/1460422
-- Also ref Vector.mjs, which I implemented myself (available upon request)
local hi = z*heightmap_size[1] + x
-- Default to this pixel's height
local up = heightmap[hi]
local down = heightmap[hi]
local left = heightmap[hi]
local right = heightmap[hi]
if z - 1 > 0 then up = heightmap[(z-1)*heightmap_size[1] + x] end
if z + 1 < heightmap_size[1] then down = heightmap[(z+1)*heightmap_size[1] + x] end
if x - 1 > 0 then left = heightmap[z*heightmap_size[1] + (x-1)] end
if x + 1 < heightmap_size[0] then right = heightmap[z*heightmap_size[1] + (x+1)] end
result[hi] = worldeditadditions.vector.normalize({
x = left - right,
y = 2, -- Z & Y are flipped
z = down - up
})
end
end
return result
end