From cdba38d37dae68dfc85ef873cdf9a79305674894 Mon Sep 17 00:00:00 2001
From: Starbeamrainbowlabs <sbrl@starbeamrainbowlabs.com>
Date: Thu, 20 Aug 2020 01:53:26 +0100
Subject: [PATCH] Start implementing snowball erosion algorithm. There's still
 a long way to go - we're only just getting warmed up!

---
 worldeditadditions/init.lua                |  1 +
 worldeditadditions/lib/erode/snowballs.lua | 67 +++++++++++++++++++-
 worldeditadditions/utils.lua               | 26 +++++++-
 worldeditadditions/utils/nodes.lua         | 38 ------------
 worldeditadditions/utils/tables.lua        | 10 +++
 worldeditadditions/utils/terrain.lua       | 72 ++++++++++++++++++++++
 6 files changed, 173 insertions(+), 41 deletions(-)
 create mode 100644 worldeditadditions/utils/terrain.lua

diff --git a/worldeditadditions/init.lua b/worldeditadditions/init.lua
index 8be8c1a..571d68d 100644
--- a/worldeditadditions/init.lua
+++ b/worldeditadditions/init.lua
@@ -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")
diff --git a/worldeditadditions/lib/erode/snowballs.lua b/worldeditadditions/lib/erode/snowballs.lua
index e834f10..f9fc82e 100644
--- a/worldeditadditions/lib/erode/snowballs.lua
+++ b/worldeditadditions/lib/erode/snowballs.lua
@@ -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
diff --git a/worldeditadditions/utils.lua b/worldeditadditions/utils.lua
index 375c3b1..3ab8504 100644
--- a/worldeditadditions/utils.lua
+++ b/worldeditadditions/utils.lua
@@ -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
diff --git a/worldeditadditions/utils/nodes.lua b/worldeditadditions/utils/nodes.lua
index e39b7ae..fe7cf32 100644
--- a/worldeditadditions/utils/nodes.lua
+++ b/worldeditadditions/utils/nodes.lua
@@ -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
diff --git a/worldeditadditions/utils/tables.lua b/worldeditadditions/utils/tables.lua
index c41b47d..2c868b3 100644
--- a/worldeditadditions/utils/tables.lua
+++ b/worldeditadditions/utils/tables.lua
@@ -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
diff --git a/worldeditadditions/utils/terrain.lua b/worldeditadditions/utils/terrain.lua
new file mode 100644
index 0000000..2148f39
--- /dev/null
+++ b/worldeditadditions/utils/terrain.lua
@@ -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