diff --git a/CHANGELOG.md b/CHANGELOG.md
index 721a74c..28ed75f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,7 @@ It's about time I started a changelog! This will serve from now on as the master
 
 ## v1.9 (unreleased)
  - Add `//many` for executing a command many times in a row
+ - Add **experimental** `//erode` command
 
 
 ## v1.8: The Quality of Life Update (17th July 2020)
diff --git a/README.md b/README.md
index c547237..887d07c 100644
--- a/README.md
+++ b/README.md
@@ -27,6 +27,7 @@ If you can dream of it, it probably belongs here!
  - [`//overlay <node_name_a> [<chance_a>] <node_name_b> [<chance_b>] [<node_name_N> [<chance_N>]] ...`](#overlay-node_name_a-chance_a-node_name_b-chance_b-node_name_n-chance_n-)
  - [`//layers [<node_name_1> [<layer_count_1>]] [<node_name_2> [<layer_count_2>]] ...`](#layers-node_name_1-layer_count_1-node_name_2-layer_count_2-)
  - [`//convolve <kernel> [<width>[,<height>]] [<sigma>]`](#convolve-kernel-widthheight-sigma)
+ - [`//erode [<snowballs|...> [<key_1> [<value_1>]] [<key_2> [<value_2>]] ...]`](#erode-snowballs-key_1-value_1-key_2-value_2-) **experimental**
 
 ### Statistics
  - [`//count`](#count)
@@ -53,7 +54,7 @@ Floods all connected nodes of the same type starting at _pos1_ with <replace_nod
 //floodfill glass 25
 ```
 
-### `//overlay <node_name_a> [<chance_a>] <node_name_b> [<chance_b>] [<node_name_N> [<chance_N>]] ...`
+### `//overlay <node_name_a> [<chance[<snowballs|...> [<key_1> [<vaue_1>]] [<key_2> [<value_2>]] ...]_a>] <node_name_b> [<chance_b>] [<node_name_N> [<chance_N>]] ...`
 Places `<node_name_a>` in the last contiguous air space encountered above the first non-air node. In other words, overlays all top-most nodes in the specified area with `<node_name_a>`. Optionally supports a mix of node names and chances, as `//mix` (WorldEdit) and `//replacemix` (WorldEditAdditions) does.
 
 Will also work in caves, as it scans columns of nodes from top to bottom, skipping every non-air node until it finds one - and only then will it start searching for a node to place the target node on top of.
@@ -256,6 +257,43 @@ The sigma value is only applicable to the `gaussian` kernel, and can be thought
 //convolve gaussian 5 0.2
 ```
 
+## `//erode [<snowballs|...> [<key_1> [<value_1>]] [<key_2> [<value_2>]] ...]`
+Runs an erosion algorithm over the defined region, optionally passing a number of key - value pairs representing parameters that are passed to the chosen algorithm. This command is **experimental**, as the author is currently on-the-fence about the effects it produces.
+
+Currently implemented algorithms:
+
+Algorithm	| Mode	| Description
+------------|-------|-------------------
+`snowballs`	| 2D	| The default. Simulates snowballs rolling across the terrain, eroding & depositing material. Then runs a 3x3 gaussian kernel over the result (i.e. like the `//conv` / `//smoothadv` command).
+
+Usage examples:
+
+```
+//erode
+//erode snowballs
+//erode snowballs count 25000
+```
+
+Each of the algorithms above have 1 or more parameters that they support. These are detailed below.
+
+### Parameters: snowballs
+
+Parameter			| Type		| Default Value		| Description
+--------------------|-----------|-------------------|------------------------
+rate_deposit		| `float`	| 0.03				| The rate at which snowballs will deposit material
+rate_erosion		| `float`	| 0.04				| The rate at which snowballs will erode material
+friction			| `float`	| 0.07				| More friction slows snowballs down more.
+speed				| `float`	| 1					| Speed multiplier to apply to snowballs at each step.
+max_steps			| `float`	| 80				| The maximum number of steps to simulate each snowball for.
+velocity_hist_count	| `float`	| 3					| The number of previous history values to average when detecting whether a snowball has stopped or not
+init_velocity		| `float`	| 0.25				| The maximum random initial velocity of a snowball for each component of the velocity vector.
+scale_iterations	| `float`	| 0.04				| How much to scale erosion by as time goes on. Higher values mean that any given snowball will erode more later on as more steps pass.
+maxdiff				| `float`	| 0.4				| The maximum difference in height (between 0 and 1) that is acceptable as a percentage of the defined region's height.
+count				| `float`	| 50000				| The number of snowballs to simulate.
+
+If you find any good combinations of these parameters, please [open an issue](https://github.com/sbrl/Minetest-WorldEditAdditions/issues/new) (or a PR!) and let me know! I'll include good combinations here.
+
+
 ### `//count`
 Counts all the nodes in the defined region and returns the result along with calculated percentages (note that if the chat window used a monospace font, the returned result would be a perfect table. If someone has a ~~hack~~ solution to make the columns line up neatly, please [open an issue](https://github.com/sbrl/Minetest-WorldEditAdditions/issues/new) :D)
 
diff --git a/worldeditadditions/init.lua b/worldeditadditions/init.lua
index 7758048..a7c10f8 100644
--- a/worldeditadditions/init.lua
+++ b/worldeditadditions/init.lua
@@ -7,13 +7,14 @@
 
 worldeditadditions = {}
 worldeditadditions.modpath = minetest.get_modpath("worldeditadditions")
+dofile(worldeditadditions.modpath.."/utils/vector.lua")
 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")
 dofile(worldeditadditions.modpath.."/lib/floodfill.lua")
 dofile(worldeditadditions.modpath.."/lib/overlay.lua")
 dofile(worldeditadditions.modpath.."/lib/layers.lua")
@@ -23,7 +24,8 @@ dofile(worldeditadditions.modpath.."/lib/walls.lua")
 dofile(worldeditadditions.modpath.."/lib/replacemix.lua")
 dofile(worldeditadditions.modpath.."/lib/maze2d.lua")
 dofile(worldeditadditions.modpath.."/lib/maze3d.lua")
-dofile(worldeditadditions.modpath.."/lib/conv/convolution.lua")
+dofile(worldeditadditions.modpath.."/lib/conv/conv.lua")
+dofile(worldeditadditions.modpath.."/lib/erode/erode.lua")
 
 dofile(worldeditadditions.modpath.."/lib/count.lua")
 
diff --git a/worldeditadditions/lib/conv/convolution.lua b/worldeditadditions/lib/conv/conv.lua
similarity index 63%
rename from worldeditadditions/lib/conv/convolution.lua
rename to worldeditadditions/lib/conv/conv.lua
index eff91ab..93b0dbc 100644
--- a/worldeditadditions/lib/conv/convolution.lua
+++ b/worldeditadditions/lib/conv/conv.lua
@@ -75,45 +75,11 @@ function worldeditadditions.convolve(pos1, pos2, kernel, kernel_size)
 	-- worldeditadditions.print_2d(heightmap, (pos2.z - pos1.z) + 1)
 	-- print("transformed")
 	-- worldeditadditions.print_2d(heightmap_conv, (pos2.z - pos1.z) + 1)
-	-- It seems to be convolving as intended, but something's probably getting lost in translation below
 	
-	for z = heightmap_size[0], 0, -1 do
-		for x = heightmap_size[1], 0, -1 do
-			local hi = z*heightmap_size[1] + x
-			
-			
-			local height_old = heightmap[hi]
-			local height_new = heightmap_conv[hi]
-			-- print("[conv/save] hi", hi, "height_old", heightmap[hi], "height_new", heightmap_conv[hi], "z", z, "x", x, "pos1.y", pos1.y)
-			
-			-- Lua doesn't have a continue statement :-/
-			if height_old == height_new then
-				-- noop
-			elseif height_new < height_old then
-				stats.removed = stats.removed + (height_old - height_new)
-				local y = height_new
-				while y < height_old do
-					local ci = area:index(pos1.x + x, pos1.y + y, pos1.z + z)
-					-- print("[conv/save] remove at y", y, "→", pos1.y + y, "current:", minetest.get_name_from_content_id(data[ci]))
-					data[ci] = node_id_air
-					y = y + 1
-				end
-			else -- height_new > height_old
-				-- We subtract one because the heightmap starts at 1 (i.e. 1 = 1 node in the column), but the selected region is inclusive
-				local node_id = data[area:index(pos1.x + x, pos1.y + (height_old - 1), pos1.z + z)]
-				-- print("[conv/save] filling with ", node_id, "→", minetest.get_name_from_content_id(node_id))
-				
-				stats.added = stats.added + (height_new - height_old)
-				local y = height_old
-				while y < height_new do
-					local ci = area:index(pos1.x + x, pos1.y + y, pos1.z + z)
-					-- print("[conv/save] add at y", y, "→", pos1.y + y, "current:", minetest.get_name_from_content_id(data[ci]))
-					data[ci] = node_id
-					y = y + 1
-				end
-			end
-		end
-	end
+	worldeditadditions.apply_heightmap_changes(
+		pos1, pos2, area, data,
+		heightmap, heightmap_conv, heightmap_size
+	)
 	
 	worldedit.manip_helpers.finish(manip, data)
 	
diff --git a/worldeditadditions/lib/conv/convolve.lua b/worldeditadditions/lib/conv/convolve.lua
index 3b1bb70..9404bf0 100644
--- a/worldeditadditions/lib/conv/convolve.lua
+++ b/worldeditadditions/lib/conv/convolve.lua
@@ -16,12 +16,12 @@ function worldeditadditions.conv.convolve(heightmap, heightmap_size, matrix, mat
 	local border_size = {}
 	border_size[0] = (matrix_size[0]-1) / 2		-- height
 	border_size[1] = (matrix_size[1]-1) / 2		-- width
-	print("[convolve] matrix_size", matrix_size[0], matrix_size[1])
-	print("[convolve] border_size", border_size[0], border_size[1])
-	print("[convolve] heightmap_size: ", heightmap_size[0], heightmap_size[1])
-	
-	print("[convolve] z: from", (heightmap_size[0]-border_size[0]) - 1, "to", border_size[0], "step", -1)
-	print("[convolve] x: from", (heightmap_size[1]-border_size[1]) - 1, "to", border_size[1], "step", -1)
+	-- print("[convolve] matrix_size", matrix_size[0], matrix_size[1])
+	-- print("[convolve] border_size", border_size[0], border_size[1])
+	-- print("[convolve] heightmap_size: ", heightmap_size[0], heightmap_size[1])
+	-- 
+	-- print("[convolve] z: from", (heightmap_size[0]-border_size[0]) - 1, "to", border_size[0], "step", -1)
+	-- print("[convolve] x: from", (heightmap_size[1]-border_size[1]) - 1, "to", border_size[1], "step", -1)
 	
 	-- Convolve over only the bit that allows us to use the full convolution matrix
 	for z = (heightmap_size[0]-border_size[0]) - 1, border_size[0], -1 do
diff --git a/worldeditadditions/lib/erode/erode.lua b/worldeditadditions/lib/erode/erode.lua
new file mode 100644
index 0000000..eb2a24c
--- /dev/null
+++ b/worldeditadditions/lib/erode/erode.lua
@@ -0,0 +1,42 @@
+worldeditadditions.erode = {}
+
+dofile(worldeditadditions.modpath.."/lib/erode/snowballs.lua")
+
+
+function worldeditadditions.erode.run(pos1, pos2, algorithm, params)
+	pos1, pos2 = worldedit.sort_pos(pos1, pos2)
+	
+	local manip, area = worldedit.manip_helpers.init(pos1, pos2)
+	local data = manip:get_data()
+	
+	local heightmap_size = {}
+	heightmap_size[0] = (pos2.z - pos1.z) + 1
+	heightmap_size[1] = (pos2.x - pos1.x) + 1
+	
+	local region_height = (pos2.y - pos1.y) + 1
+	
+	local heightmap = worldeditadditions.make_heightmap(pos1, pos2, manip, area, data)
+	local heightmap_eroded = worldeditadditions.shallowcopy(heightmap)
+	
+	-- print("[erode.run] algorithm: "..algorithm..", params:");
+	-- print(worldeditadditions.map_stringify(params))
+	worldeditadditions.print_2d(heightmap, heightmap_size[1])
+	
+	if algorithm == "snowballs" then
+		local success, msg = worldeditadditions.erode.snowballs(heightmap, heightmap_eroded, heightmap_size, region_height, params)
+		if not success then return success, msg end
+	else
+		return false, "Error: Unknown algorithm '"..algorithm.."'. Currently implemented algorithms: snowballs (2d; hydraulic-like). Ideas for algorithms to implement are welcome!"
+	end
+	
+	local success, stats = worldeditadditions.apply_heightmap_changes(
+		pos1, pos2, area, data,
+		heightmap, heightmap_eroded, heightmap_size
+	)
+	if not success then return success, stats end
+	worldedit.manip_helpers.finish(manip, data)
+	
+	print("[erode] stats")
+	print(worldeditadditions.map_stringify(stats))
+	return true, stats
+end
diff --git a/worldeditadditions/lib/erode/snowballs.lua b/worldeditadditions/lib/erode/snowballs.lua
new file mode 100644
index 0000000..dd3e369
--- /dev/null
+++ b/worldeditadditions/lib/erode/snowballs.lua
@@ -0,0 +1,143 @@
+
+-- Test command: //multi //fp set1 1313 6 5540 //fp set2 1338 17 5521 //erode snowballs
+
+local function snowball(heightmap, normalmap, heightmap_size, startpos, params)
+	local sediment = 0
+	local pos = { x = startpos.x, z = startpos.z }
+	local pos_prev = { x = pos.x, z = pos.z }
+	local velocity = {
+		x = (math.random() * 2 - 1) * params.init_velocity,
+		z = (math.random() * 2 - 1) * params.init_velocity
+	}
+	local heightmap_length = #heightmap
+	
+	-- print("[snowball] startpos ("..pos.x..", "..pos.z.."), velocity: ("..velocity.x..", "..velocity.z..")")
+	
+	local hist_velocity = {}
+	
+	for i = 1, params.max_steps do
+		local x = pos.x
+		local z = pos.z
+		local hi = math.floor(z+0.5)*heightmap_size[1] + math.floor(x+0.5)
+		-- Stop if we go out of bounds
+		if x < 0 or z < 0
+			or x >= heightmap_size[1]-1 or z >= heightmap_size[0]-1 then
+			-- print("[snowball] hit edge; stopping at ("..x..", "..z.."), (bounds @ "..(heightmap_size[1]-1)..", "..(heightmap_size[0]-1)..")", "x", x, "/", heightmap_size[1]-1, "z", z, "/", heightmap_size[0]-1)
+			return true, i
+		end
+		
+		if #hist_velocity > 0 and i > 5
+			and worldeditadditions.average(hist_velocity) < 0.03 then
+			-- print("[snowball] It looks like we've stopped")
+			return true, i
+		end
+		
+		if normalmap[hi].y == 1 then return true, i end
+		
+		if hi > heightmap_length then return false, "Out-of-bounds on the array, hi: "..hi..", heightmap_length: "..heightmap_length end
+		
+		-- NOTE: We need to decide whether we want to keep the precomputed normals as we have now, or whether we want to dynamically compute them at the some of request.
+		-- print("[snowball] sediment", sediment, "rate_deposit", params.rate_deposit, "normalmap[hi].z", normalmap[hi].z)
+		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(x)
+			and math.floor(pos_prev.z) ~= math.floor(z) then
+			heightmap[hi] = heightmap[hi] + (step_deposit - step_erode)
+		end
+		
+		velocity.x = params.friction * velocity.x + normalmap[hi].x * params.speed
+		velocity.z = params.friction * velocity.z + normalmap[hi].y * params.speed
+		
+		-- print("[snowball] now at ("..x..", "..z..") velocity "..worldeditadditions.vector.lengthsquared(velocity)..", sediment "..sediment)
+		local new_vel_sq = worldeditadditions.vector.lengthsquared(velocity)
+		if new_vel_sq > 1 then
+			-- print("[snowball] velocity squared over 1, normalising")
+			velocity = worldeditadditions.vector.normalize(velocity)
+		end
+		table.insert(hist_velocity, new_vel_sq)
+		if #hist_velocity > params.velocity_hist_count then table.remove(hist_velocity, 1) end
+		pos_prev.x = x
+		pos_prev.z = z
+		pos.x = pos.x + velocity.x
+		pos.z = pos.z + velocity.z
+		sediment = sediment + (step_erode - step_deposit) -- Needs to be erosion - deposit, which is the opposite to the above
+	end
+	return true, params.max_steps
+end
+
+--[[
+2D erosion algorithm based on snowballs
+Note that this *mutates* the given heightmap.
+@source https://jobtalle.com/simulating_hydraulic_erosion.html
+
+]]--
+function worldeditadditions.erode.snowballs(heightmap_initial, heightmap, heightmap_size, region_height, params_custom)
+	local params = {
+		rate_deposit = 0.03, -- 0.03
+		rate_erosion = 0.04, -- 0.04
+		friction = 0.07,
+		speed = 1,
+		max_steps = 80,
+		velocity_hist_count = 3,
+		init_velocity = 0.25,
+		scale_iterations = 0.04,
+		maxdiff = 0.4,
+		count = 50000
+	}
+	-- Apply the default settings
+	worldeditadditions.table_apply(params_custom, params)
+	
+	print("[erode/snowballs] params: ")
+	print(worldeditadditions.map_stringify(params))
+	
+	
+	local normals = worldeditadditions.calculate_normals(heightmap, heightmap_size)
+	
+	local stats_steps = {}
+	for i = 1, params.count do
+		-- print("[snowballs] starting snowball ", i)
+		local success, steps = snowball(
+			heightmap, normals, heightmap_size,
+			{
+				x = math.random() * (heightmap_size[1] - 1),
+				z = math.random() * (heightmap_size[0] - 1)
+			},
+			params
+		)
+		table.insert(stats_steps, steps)
+		
+		if not success then return false, "Error: Failed at snowball "..i..":"..steps end
+	end
+	
+	print("[snowballs] "..#stats_steps.." snowballs simulated, max "..params.max_steps.." steps, averaged ~"..worldeditadditions.average(stats_steps).."")
+	
+	-- Round everything to the nearest int, since you can't really have
+	-- something like .141592671 of a node
+	-- Note that we do this after *all* the erosion is complete
+	local clamp_limit = math.floor(region_height * params.maxdiff + 0.5)
+	for i,v in ipairs(heightmap) do
+		heightmap[i] = math.floor(heightmap[i] + 0.5)
+		if heightmap[i] < 0 then heightmap[i] = 0 end
+		-- Limit the distance to params.maxdiff% of the region height
+		if math.abs(heightmap_initial[i] - heightmap[i]) > region_height * params.maxdiff then
+			if heightmap_initial[i] - heightmap[i] > 0 then
+				heightmap[i] = heightmap_initial[i] - clamp_limit
+			else
+				heightmap[i] = heightmap_initial[i] + clamp_limit
+			end
+		end
+	end
+	
+	local success, matrix = worldeditadditions.get_conv_kernel("gaussian", 3, 3)
+	if not success then return success, matrix end
+	matrix_size = {} matrix_size[0] = 3 matrix_size[1] = 3
+	worldeditadditions.conv.convolve(
+		heightmap, heightmap_size,
+		matrix,
+		matrix_size
+	)
+	
+	return true, params.count.." snowballs simulated"
+end
diff --git a/worldeditadditions/utils.lua b/worldeditadditions/utils.lua
deleted file mode 100644
index 375c3b1..0000000
--- a/worldeditadditions/utils.lua
+++ /dev/null
@@ -1,15 +0,0 @@
-worldeditadditions.vector = {}
-
-function worldeditadditions.vector.tostring(v)
-	return "(" .. v.x ..", " .. v.y ..", " .. v.z ..")"
-end
-
-function worldeditadditions.vector.lengthsquared(v)
-	return v.x*v.x + v.y*v.y + v.z*v.z
-end
-
-function worldeditadditions.vector.floor(v)
-	v.x = math.floor(v.x)
-	v.y = math.floor(v.y)
-	v.z = math.floor(v.z)
-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/strings.lua b/worldeditadditions/utils/strings.lua
index 52a0b2b..1b2f805 100644
--- a/worldeditadditions/utils/strings.lua
+++ b/worldeditadditions/utils/strings.lua
@@ -86,6 +86,7 @@ end
 -- @param	tbl		number[]	The ZERO-indexed list of numbers
 -- @param	width	number		The width of 2D array.
 function worldeditadditions.print_2d(tbl, width)
+	print("==== count: "..#tbl..", width:"..width.." ====")
 	local display_width = 1
 	for _i,value in pairs(tbl) do
 		display_width = math.max(display_width, #tostring(value))
@@ -209,6 +210,32 @@ function worldeditadditions.parse_weighted_nodes(parts, as_list)
 	return true, result
 end
 
+function worldeditadditions.parse_map(params_text)
+	local result = {}
+	local parts = worldeditadditions.split(params_text, "%s+", false)
+	
+	local last_key = nil
+	for i, part in ipairs(parts) do
+		if i % 2 == 0 then -- Lua starts at 1 :-/
+			-- Try converting to a number to see if it works
+			local part_converted = tonumber(part)
+			if as_number == nil then part_converted = part end
+			result[last_key] = part
+		else
+			last_key = part
+		end
+	end
+	return true, result
+end
+
+function worldeditadditions.map_stringify(map)
+	local result = {}
+	for key, value in pairs(map) do
+		table.insert(result, key.."\t"..value)
+	end
+	return table.concat(result, "\n")
+end
+
 --- Converts a float milliseconds into a human-readable string.
 -- Ported from PHP human_time from Pepperminty Wiki: https://github.com/sbrl/Pepperminty-Wiki/blob/fa81f0d/core/05-functions.php#L82-L104
 -- @param	ms		float	The number of milliseconds to convert.
diff --git a/worldeditadditions/utils/tables.lua b/worldeditadditions/utils/tables.lua
index c41b47d..99e41ac 100644
--- a/worldeditadditions/utils/tables.lua
+++ b/worldeditadditions/utils/tables.lua
@@ -15,3 +15,16 @@ 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)
+	print("[table_apply] start")
+	for key, value in pairs(source) do
+		print("[table_apply] Applying", key, "=", value)
+		target[key] = value
+	end
+	print("[table_apply] end")
+end
diff --git a/worldeditadditions/utils/terrain.lua b/worldeditadditions/utils/terrain.lua
new file mode 100644
index 0000000..d4353b6
--- /dev/null
+++ b/worldeditadditions/utils/terrain.lua
@@ -0,0 +1,124 @@
+
+--- 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)
+	print("heightmap_size: "..heightmap_size[1].."x"..heightmap_size[0])
+	local result = {}
+	for z = heightmap_size[0]-1, 0, -1 do
+		for x = heightmap_size[1]-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[0]-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[1]-1 then right = heightmap[z*heightmap_size[1] + (x+1)] end
+			
+			-- print("[normals] UP	| index", (z-1)*heightmap_size[1] + x, "z", z, "z-1", z - 1, "up", up, "limit", 0)
+			-- print("[normals] DOWN	| index", (z+1)*heightmap_size[1] + x, "z", z, "z+1", z + 1, "down", down, "limit", heightmap_size[1]-1)
+			-- print("[normals] LEFT	| index", z*heightmap_size[1] + (x-1), "x", x, "x-1", x - 1, "left", left, "limit", 0)
+			-- print("[normals] RIGHT	| index", z*heightmap_size[1] + (x+1), "x", x, "x+1", x + 1, "right", right, "limit", heightmap_size[1]-1)
+			
+			result[hi] = worldeditadditions.vector.normalize({
+				x = left - right,
+				y = 2, -- Z & Y are flipped
+				z = down - up
+			})
+			
+			-- print("[normals] at "..hi.." ("..x..", "..z..") normal "..worldeditadditions.vector.tostring(result[hi]))
+		end
+	end
+	return result
+end
+
+function worldeditadditions.apply_heightmap_changes(pos1, pos2, area, data, heightmap_old, heightmap_new, heightmap_size)
+	local stats = { added = 0, removed = 0 }
+	local node_id_air = minetest.get_content_id("air")
+	
+	for z = heightmap_size[0], 0, -1 do
+		for x = heightmap_size[1], 0, -1 do
+			local hi = z*heightmap_size[1] + x
+			
+			local height_old = heightmap_old[hi]
+			local height_new = heightmap_new[hi]
+			-- print("[conv/save] hi", hi, "height_old", heightmap_old[hi], "height_new", heightmap_new[hi], "z", z, "x", x, "pos1.y", pos1.y)
+			
+			-- Lua doesn't have a continue statement :-/
+			if height_old == height_new then
+				-- noop
+			elseif height_new < height_old then
+				stats.removed = stats.removed + (height_old - height_new)
+				local y = height_new
+				while y < height_old do
+					local ci = area:index(pos1.x + x, pos1.y + y, pos1.z + z)
+					-- print("[conv/save] remove at y", y, "→", pos1.y + y, "current:", minetest.get_name_from_content_id(data[ci]))
+					data[ci] = node_id_air
+					y = y + 1
+				end
+			else -- height_new > height_old
+				-- We subtract one because the heightmap starts at 1 (i.e. 1 = 1 node in the column), but the selected region is inclusive
+				local node_id = data[area:index(pos1.x + x, pos1.y + (height_old - 1), pos1.z + z)]
+				-- print("[conv/save] filling with ", node_id, "→", minetest.get_name_from_content_id(node_id))
+				
+				stats.added = stats.added + (height_new - height_old)
+				local y = height_old
+				while y < height_new do
+					local ci = area:index(pos1.x + x, pos1.y + y, pos1.z + z)
+					-- print("[conv/save] add at y", y, "→", pos1.y + y, "current:", minetest.get_name_from_content_id(data[ci]))
+					data[ci] = node_id
+					y = y + 1
+				end
+			end
+		end
+	end
+	
+	return true, stats
+end
diff --git a/worldeditadditions/utils/vector.lua b/worldeditadditions/utils/vector.lua
new file mode 100644
index 0000000..b7c2357
--- /dev/null
+++ b/worldeditadditions/utils/vector.lua
@@ -0,0 +1,42 @@
+worldeditadditions.vector = {}
+
+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)
+	if not v.y then return v.x*v.x + v.z*v.z end
+	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.normalize(v)
+	local length = math.sqrt(worldeditadditions.vector.lengthsquared(v))
+	if not v.y then return {
+		x = v.x / length,
+		z = v.z / length
+	} end
+	return {
+		x = v.x / length,
+		y = v.y / length,
+		z = v.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)
+	-- 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_commands/commands/erode.lua b/worldeditadditions_commands/commands/erode.lua
new file mode 100644
index 0000000..e5c4ad7
--- /dev/null
+++ b/worldeditadditions_commands/commands/erode.lua
@@ -0,0 +1,44 @@
+-- ███████ ██████   ██████  ██████  ███████
+-- ██      ██   ██ ██    ██ ██   ██ ██
+-- █████   ██████  ██    ██ ██   ██ █████
+-- ██      ██   ██ ██    ██ ██   ██ ██
+-- ███████ ██   ██  ██████  ██████  ███████
+worldedit.register_command("erode", {
+	params = "[<snowballs|...> [<key_1> [<value_1>]] [<key_2> [<value_2>]] ...]",
+	description = "**experimental** Runs the specified erosion algorithm over the given defined region. This may occur in 2d or 3d. Currently implemented algorithms: snowballs (default;2d hydraulic-like). Also optionally takes an arbitrary set of key - value pairs representing parameters to pass to the algorithm. See the full documentation for details.",
+	privs = { worldedit = true },
+	require_pos = 2,
+	parse = function(params_text)
+		if not params_text or params_text == "" then
+			return true, "snowballs", {}
+		end
+		
+		if params_text:find("%s") == nil then
+			return true, params_text, {}
+		end
+		
+		local algorithm, params = params_text:match("([^%s]+)%s(.+)")
+		if algorithm == nil then
+			return false, "Failed to split params_text into 2 parts (this is probably a bug)"
+		end
+		
+		local success, map = worldeditadditions.parse_map(params)
+		if not success then return success, map end
+		return true, algorithm, map
+	end,
+	nodes_needed = function(name)
+		return worldedit.volume(worldedit.pos1[name], worldedit.pos2[name])
+	end,
+	func = function(name, algorithm, params)
+		local start_time = worldeditadditions.get_ms_time()
+		local success, stats = worldeditadditions.erode.run(
+			worldedit.pos1[name], worldedit.pos2[name],
+			algorithm, params
+		)
+		if not success then return success, stats end
+		local time_taken = worldeditadditions.get_ms_time() - start_time
+		
+		minetest.log("action", name .. " used //erode "..algorithm.." at " .. worldeditadditions.vector.tostring(worldedit.pos1[name]) .. ", adding " .. stats.added .. " nodes and removing " .. stats.removed .. " nodes in " .. time_taken .. "s")
+		return true, stats.added .. " nodes added and " .. stats.removed .. " nodes removed in " .. worldeditadditions.human_time(time_taken)
+	end
+})
diff --git a/worldeditadditions_commands/init.lua b/worldeditadditions_commands/init.lua
index 91c7366..229f460 100644
--- a/worldeditadditions_commands/init.lua
+++ b/worldeditadditions_commands/init.lua
@@ -30,6 +30,7 @@ dofile(we_c.modpath.."/commands/walls.lua")
 dofile(we_c.modpath.."/commands/maze.lua")
 dofile(we_c.modpath.."/commands/replacemix.lua")
 dofile(we_c.modpath.."/commands/convolve.lua")
+dofile(we_c.modpath.."/commands/erode.lua")
 
 dofile(we_c.modpath.."/commands/count.lua")