//erode: Finish the initial round of bugfixing, but I'm on the fence about it.

Specifically, I'm unsure about whether I'm happy with the effects of the 
algorithm.
Also, we convolve with a 3x3 gaussian kernel after erosion is complete - 
and we have verified that the erosion is having an positive effect at 
"roughening up" a terrain surface.
It seems like the initial blog post was correct: the algorithm does tend 
to make steep surfaces steeper.
It also appears that it's more effective on larger areas, and 'gentler' 
curves. THis might be because the surface normals are more conducive to 
making the snowballs roll.
Finally, 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.

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] matrix_size", matrix_size[0], matrix_size[1])
-	print("[convolve] border_size", border_size[0], border_size[1])
+	-- print("[convolve] border_size", border_size[0], border_size[1])
-	print("[convolve] heightmap_size: ", heightmap_size[0], heightmap_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] 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] 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

worldeditadditions/lib/erode/erode.lua

@@ -13,15 +13,17 @@ function worldeditadditions.erode.run(pos1, pos2, algorithm, params)
 	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("[erode.run] algorithm: "..algorithm..", params:");
-	print(worldeditadditions.map_stringify(params))
+	-- print(worldeditadditions.map_stringify(params))
 	worldeditadditions.print_2d(heightmap, heightmap_size[1])
 	
 	if algorithm == "snowballs" then
-		local success, msg = worldeditadditions.erode.snowballs(heightmap_eroded, heightmap_size, params)
+		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!"

worldeditadditions/lib/erode/snowballs.lua

@@ -1,49 +1,69 @@
 
--- Test command: //multi //fp set1 1312 5 5543 //fp set2 1336 18 5521 //erode//multi //fp set1 1312 5 5543 //fp set2 1336 18 5521 //erode
+-- 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 = 0, z = 0 }
+	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..")")
+	-- print("[snowball] startpos ("..pos.x..", "..pos.z.."), velocity: ("..velocity.x..", "..velocity.z..")")
 	
-	for i = 1, params.snowball_max_steps do
+	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[1]-1 or z >= heightmap[0]-1 then
+			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]..", "..heightmap_size[0]..")")
+			-- 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
+			return true, i
 		end
-		-- print("[snowball] now at ("..x..", "..z..") (bounds @ "..heightmap_size[1]..", "..heightmap_size[0]..")")
 		
-		if hi > heightmap_length then print("[snowball] out-of-bounds on the array, hi: "..hi..", heightmap_length: "..heightmap_length) return 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
 		
 		-- 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)
 		
-		local step_diff = step_deposit - step_erode
-		
 		-- 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_diff
+			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_diff
+		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
 
 --[[
@@ -52,26 +72,32 @@ Note that this *mutates* the given heightmap.
 @source https://jobtalle.com/simulating_hydraulic_erosion.html
 
 ]]--
-function worldeditadditions.erode.snowballs(heightmap, heightmap_size, params)
+function worldeditadditions.erode.snowballs(heightmap_initial, heightmap, heightmap_size, region_height, params_custom)
-	-- Apply the default settings
+	local params = {
-	worldeditadditions.table_apply({
+		rate_deposit = 0.03, -- 0.03
-		rate_deposit = 0.03,
+		rate_erosion = 0.04, -- 0.04
-		rate_erosion = 0.04,
 		friction = 0.07,
-		speed = 0.15,
+		speed = 1,
-		radius = 0.8,
+		max_steps = 80,
-		snowball_max_steps = 80,
+		velocity_hist_count = 3,
+		init_velocity = 0.25,
 		scale_iterations = 0.04,
-		drops_per_cell = 0.4,
+		maxdiff = 0.4,
-		snowball_count = 50000
+		count = 50000
-	}, params)
+	}
+	-- Apply the default settings
+	worldeditadditions.table_apply(params_custom, params)
+	
+	print("[erode/snowballs] params: ")
+	print(worldeditadditions.map_stringify(params))
 	
-	print("[erode/snowballs] params: "..worldeditadditions.map_stringify(params))
 	
 	local normals = worldeditadditions.calculate_normals(heightmap, heightmap_size)
 	
-	for i = 1, params.snowball_count do
+	local stats_steps = {}
-		snowball(
+	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),
@@ -79,14 +105,38 @@ function worldeditadditions.erode.snowballs(heightmap, heightmap_size, params)
 			},
 			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
 	
-	return true, params.snowball_count.." snowballs simulated"
+	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

worldeditadditions/utils/strings.lua

@@ -216,7 +216,7 @@ function worldeditadditions.parse_map(params_text)
 	
 	local last_key = nil
 	for i, part in ipairs(parts) do
-		if i % 2 == 1 then
+		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

worldeditadditions/utils/tables.lua

@@ -21,7 +21,10 @@ end
 -- @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

worldeditadditions/utils/vector.lua

@@ -8,6 +8,7 @@ end
 -- @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
 
@@ -18,6 +19,10 @@ end
 -- @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,

worldeditadditions_commands/commands/erode.lua

@@ -1,8 +1,8 @@
---  ██████  ██    ██ ███████ ██████  ██       █████  ██    ██
--- ██    ██ ██    ██ ██      ██   ██ ██      ██   ██  ██  ██
 -- ███████ ██████   ██████  ██████  ███████
--- ██    ██  ██  ██  ██      ██   ██ ██      ██   ██    ██
+-- ██      ██   ██ ██    ██ ██   ██ ██
---  ██████    ████   ███████ ██   ██ ███████ ██   ██    ██
+-- █████   ██████  ██    ██ ██   ██ █████
+-- ██      ██   ██ ██    ██ ██   ██ ██
+-- ███████ ██   ██  ██████  ██████  ███████
 worldedit.register_command("erode", {
 	params = "[<snowballs|...> [<key_1> [<vaue_1>]] [<key_2> [<value_2>]] ...]",
 	description = "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.",
@@ -35,6 +35,7 @@ worldedit.register_command("erode", {
 			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")