terrainlib: More optimizations on flow routing

terrainlib_lua/rivermapper.lua

@@ -103,10 +103,9 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 		local elev = i2 == 0 and dem[i1] or mmax(dem[i1], dem[i2]) -- Elevation of the highest of the two sides of the link (or only i1 if b2 is map outside)
 		local l2 = basin_links[b2]
 		if not l2 then
-			l2 = {}
+			l2 = {b2, b1, elev=elev, i=mmax(i1,i2), is_y=isY}
 			basin_links[b2] = l2
-		end
-		if not l2.elev or l2.elev > elev then -- If this link is lower than the lowest registered link between these two basins, register it as the new lowest pass
+		elseif l2.elev > elev then -- If this link is lower than the lowest registered link between these two basins, register it as the new lowest pass
 			l2.elev = elev
 			l2.i = mmax(i1,i2)
 			l2.is_y = isY
@@ -210,26 +209,28 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 	-- Mareš' optimizations mainly consist in skipping elements that have over 8 links, until extra links are removed when other elements are merged.
 	-- Note that for this step we are only working on basins, not grid nodes.
 	local lowlevel = {}
-	for i, n in pairs(nlinks) do
-		if n <= 8 then
-			lowlevel[i] = links[i]
+	cur = 0
+	local ref = singular -- Reuse table
+	for i=0, nbasins do
+		if nlinks[i] <= 8 then
+			cur = cur + 1
+			lowlevel[cur] = i
+			ref[i] = cur
 		end
 	end
 
 	local basin_graph = {}
-	for n=1, nbasins do
+	while cur > 1 do
 		-- Iterate in lowlevel but its contents may change during the loop
-		-- 'next' called with only one argument always returns an element if table is not empty
-		local b1, lnk1 = next(lowlevel)
-		lowlevel[b1] = nil
+		local b1 = lowlevel[cur]
+		cur = cur - 1
+		local lnk1 = links[b1]
 
 		local b2
 		local lowest = math.huge
 		local lnk1 = links[b1]
-		local i = 0
 		-- Look for lowest link
 		for bn, bdata in pairs(lnk1) do
-			i = i + 1
 			if bdata.elev < lowest then
 				lowest = bdata.elev
 				b2 = bn
@@ -256,7 +257,9 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 		nlinks[b2] = nlinks[b2] - 1
 		-- When the number of links is changing, we need to check whether the basin can be added to / removed from 'lowlevel'
 		if nlinks[b2] == 8 then
-			lowlevel[b2] = lnk2
+			cur = cur + 1
+			lowlevel[cur] = b2
+			ref[b2] = cur
 		end
 		-- 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
@@ -266,12 +269,16 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 			if lnkn[b2] then -- If bassin bn is also linked to b2
 				nlinks[bn] = nlinks[bn] - 1 -- Then bassin bn is losing a link because it keeps only one link toward b1/b2 after the merge
 				if nlinks[bn] == 8 then
-					lowlevel[bn] = lnkn
+					cur = cur + 1
+					lowlevel[cur] = bn
+					ref[bn] = cur
 				end
 			else -- If bn was linked to b1 but not to b2
 				nlinks[b2] = nlinks[b2] + 1 -- Then b2 is gaining a link to bn because of the merge
 				if nlinks[b2] == 9 then
-					lowlevel[b2] = nil
+					lowlevel[ref[b2]] = lowlevel[cur]
+					ref[lowlevel[cur]] = ref[b2]
+					cur = cur - 1
 				end
 			end
 
@@ -289,15 +296,17 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 	-- To orient the basin graph, we will consider that the ultimate basin water should flow into is the map outside (basin #0). We will start from it and recursively walk upstream to the neighbouring basins, using only links that are in the minimal spanning tree. This gives the flow direction of the links, and thus, the outlet of every basin.
 	-- This will also give lake elevation, which is the highest link encountered between map outside and the given basin on the spanning tree.
 	-- And within each basin, we need to modify flow directions to connect the singular node to the outlet.
-	local queue = {[0] = -math.huge}
+	local queue = {0}
+	local queuevalues = {-math.huge}
+	cur = 1
 	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
-		local b1, elev1 = next(queue) -- Pop from queue
-		queue[b1] = nil
+	while cur > 0 do
+		local b1, elev1 = queue[cur], queuevalues[cur] -- Pop from queue
+		cur = cur - 1
 		basin_lake[b1] = elev1
 		-- Iterate through b1's neighbours (according to the spanning tree)
 		for b2, bound in pairs(basin_graph[b1]) do
@@ -336,7 +345,9 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 			until dir == 0 -- Stop when reaching the singular node
 
 			-- Add basin b2 into the queue, and keep the highest link elevation, that will define the elevation of the lake in b2
-			queue[b2] = mmax(elev1, bound.elev)
+			cur = cur + 1
+			queue[cur] = b2
+			queuevalues[cur] = mmax(elev1, bound.elev)
 			-- Remove b1 from b2's neighbours to avoid coming back to b1
 			basin_graph[b2][b1] = nil
 		end