Add Grid and Flow implementations to the tests#221
Merged
Conversation
See TopoToolbox#55 These two structs are implementations of the GRIDobj and FLOWobj interfaces to be used in the tests. They should make adding and modifying tests much easier by hiding some of the details of creating and accessing data sets. Grid is a two dimensional array with an attached cellsize but no coordinate reference information. It provides methods for generating random data, blank datasets or loading through GDAL. It also provides arithmetic and comparison operations and TopoToolbox functions needed for flow routing. Flow is a topologically sorted edge list with source, target and fraction fields as well as the topological node ordering in stream. FlowRouteD8Carve runs D8 with morphological carving on the provided DEM and FlowAccumulation and FlowDrainageBasins use flow algebras. The random_dem test has been heavily modified to use the Grid and Flow structs where possible. Test functions take Grid or Flow objects, and the flow routing, gradient computation, flow accumulation, etc. has been converted to generate Grids or Flows. Some refactoring has been carried out along the way such as - Implement a simpler topological sorting test using a Flow. The previous one used flow directions and is a relic of the older flow routing interface. - Remove calls to flow_accumulation, flow_accumulation_edgelist and drainagebasins in favor of the flow algebra implementation of Flow (see TopoToolbox#212). The test_flow_accumulation_multimethod test is therefore removed as there is really only one supported flow accumulation method now. - Separate fillsinks/fillsinks_hybrid from the flow routing function. We previously had two versions of flow routing to be able to compare the performance of fillsinks and fillsinks_hybrid. The profiler now profiles the two fillsinks algorithms separately. There is also a Stream implementation on the way, but it is a little more challenging to ensure the memory allocation details are correct. Once that lands, the remaining random_dem tests of the stream network can also be refactored. GDAL is now required for the random_dem tests because Grid depends on it. It would be nice to figure out how to compile GridFromFile and GridToFile only if GDAL is available, so we can run random_dem without GDAL. Signed-off-by: William Kearney <william.kearney@uni-potsdam.de>
The GDAL requirement for random_dem means that it currently is not run on macos or Windows.
wkearn
added a commit
to wkearn/libtopotoolbox
that referenced
this pull request
Apr 9, 2026
This PR implements the D-infinity flow routing algorithm of Tarboton (1997). D-infinity routes flow to up to two downstream neighbors of each pixel, with flow proportioned according to the angle of the normal vector of 8 triangular facets surrounding the center pixel. This PR also introduces a major redesign of the flow routing algorithms to accommodate better multiple flow direction algorithms and to separate flow routing from flat resolution. This design is used in the D-infinity implementation and several helper functions are added. The existing D8 implementation has not yet been modified. The new flow routing procedure is: 1. Use a flow routing algorithm to determine flow directions over the DEM. These are stored in a bit array of type uint8_t where bit k is set in pixel (i, j) if there is an outgoing edge from that pixel in the kth direction. The directions are numbered starting from 1 to the right and moving counterclockwise. Any unresolved pixels have a direction of 0. These could be missing data in the DEM, boundary pixels or sinks. 2. Count the number of edges in the flow direction bitmap. We provide a function edgeset_count to do this, but it could also be done externally. It is the number of set bits in the direction bitmap. 3. Allocate a weight array of type float containing the counted number of edges. Fill it up according to the weighting scheme of the flow routing algorithm. These weights are packed into the array to minimize space usage. This is why the number of edges need to be counted prior to computing the weights. The flow_routing_dinf_weights function computes and stores them in the correct order. 4. If you don't want to resolve flats, skip to step 7. 5. Apply a flat resolution algorithm to create a new direction bitmap and weight array for the unresolved pixels. This process looks a lot like the regular flow routing: fill up the bitmap, count the edges, allocate and fill up the weights. 6. Count the total number of edges in the originally resolved direction bitmap and the bitmap of directions on flats. edgeset_count_merged is provided for this purpose. Allocate a weight array and merge the originally resolved directions and the resolved directions on flats using the merge_edges function. You have to do this after flat resolution because of the way the edges are sorted in the weight array. There might be a better way to do it that would allow you to skip this step and merge on the fly. 7. Now you have a bitmap of directions and float array of weights for all the edges you want in your Flow object. Use the edgeset_scan function to prepare an array of type ptrdiff_t with the offsets of each pixels edges in the weight array. This is done for you by edgeset_merge if you used that. edgeset_count, edgeset_scan and edgeset_merge all return the total number of edges in the edge set. 8. Allocate your Flow object arrays: - ptrdiff_t *stream: the topologically sorted list of nodes. Size is the number of pixels in your DEM - ptrdiff_t *source: the topologically sorted list of edge sources. Size is the number of edges - ptrdiff_t *target: the topologically sorted list of edge targets. Size is the number of edges. - float *weight: the topologically sorted list of edge weights. Size is the number of edges. 9. Use direction_tsort to traverse the flow network and topologically sort the edges. The node, source, target and weight arrays are filled up in the right order. The bitmap edge set provides a uniform interface for flow routing and flat resolution algorithms to exchange data. This could also be done with topologically sorted edge lists as stored in the Flow objects. However, it is harder to merge two topological sorts than it is to merge two edge sets. The topological sort has to be redone with another graph traversal, and edge lists don't let you compute traversals from an arbitrary pixel so easily. This design lets you interchange flow routing and flat resolution algorithms as desired. Each algorithm returns an edge set that can be merged and then topologically sorted. We can, for example, decouple the D8 flow routing from the least cost auxiliary topography carving that we currently do. And you could conceivably use the LCAT carving with Dinf if you wanted. The D-inf implementation consists of two functions in src/dinf.c. - flow_routing_dinf_directions computes the direction bitmap and an array of outgoing flow proportions for each pixel - flow_routing_dinf_weights computes the edge weights from the array of flow proportions. A new function flow_routing_tsort is added to src/flow_routing.c to topologically sort bitmap edge sets. This is intended to replace the topological sorting functionality of flow_routing_d8_carve. A shortest-path flat resolution algorithm is added to src/flow_routing.c with two functions - resolve_flats_shortest_path computes the direction bitmap using Dijkstra's algorithm to route flow along flats following the shortest distance to a sill - resolve_flats_shortest_path_weights computes the flow weights for the shortest-path flat resolution. These are all 1s, but this function is provided to maintain the direction/weights interface. The shortest path resolution has been added because it is what is described by Tarboton (1997). He uses an iterative algorithm, but it amounts to the same thing. This is not what is used by default in MATLAB TopoToolbox or, as far as I can tell, in Tarboton's TauDEM. Moreover this implementation has been done entirely from scratch based on the description in Tarboton (1997) without reliance on the D-inf implementation in MATLAB TopoToolbox or TauDEM. It probably does not behave exactly as in MATLAB TopoToolbox. The flat resolution approach is certainly different. MATLAB TopoToolbox's D-inf implementation also fails on some DEMs with NaNs such as the Taiwan example. I think the problem there is when there are embedded NaNs that are not contiguous with the border. libtopotoolbox handles embedded NaNs as missing data. Because of the way D-inf computes flow angles based on triangular facets, this can lead to some surprising results, such as the following situation: 8 9 7 5 5 6 NaN 1 NaN The center pixel should probably flow down, but because both facets containing the downward neighbor have NaNs, the slope of each facet is NaN, and the flow is unresolved. The flat resolution algorithm will probably flow to the neighbor of equal elevation to the left and from there to the bottom right. This actually occurs in the Kunashiri example DEM. It could be possible to fix this by special casing the situation where one of the two neighboring pixels for a facet is NaN. Functions for dealing with the edge sets have been added to src/helpers/edgeset.c. These could conceivably also be done in higher level languages, but making sure you get all the ordering right is a bit tricky. Test cases have been added to test/dinf.c. These tests use the Grid and Flow object interface added in TopoToolbox#221 and test on a variety of small test cases, including those from the Tarboton paper, the snapshot DEMs and randomly generated ones. Memory order tests are also conducted. I have also taken this opportunity to prototype a new testing harness that provides TAP output. Signed-off-by: William Kearney <william.kearney@uni-potsdam.de>
wkearn
added a commit
that referenced
this pull request
Apr 27, 2026
This PR implements the D-infinity flow routing algorithm of Tarboton (1997). D-infinity routes flow to up to two downstream neighbors of each pixel, with flow proportioned according to the angle of the normal vector of 8 triangular facets surrounding the center pixel. This PR also introduces a major redesign of the flow routing algorithms to accommodate better multiple flow direction algorithms and to separate flow routing from flat resolution. This design is used in the D-infinity implementation and several helper functions are added. The existing D8 implementation has not yet been modified. The new flow routing procedure is: 1. Use a flow routing algorithm to determine flow directions over the DEM. These are stored in a bit array of type uint8_t where bit k is set in pixel (i, j) if there is an outgoing edge from that pixel in the kth direction. The directions are numbered starting from 1 to the right and moving counterclockwise. Any unresolved pixels have a direction of 0. These could be missing data in the DEM, boundary pixels or sinks. 2. Count the number of edges in the flow direction bitmap. We provide a function edgeset_count to do this, but it could also be done externally. It is the number of set bits in the direction bitmap. 3. Allocate a weight array of type float containing the counted number of edges. Fill it up according to the weighting scheme of the flow routing algorithm. These weights are packed into the array to minimize space usage. This is why the number of edges need to be counted prior to computing the weights. The flow_routing_dinf_weights function computes and stores them in the correct order. 4. If you don't want to resolve flats, skip to step 7. 5. Apply a flat resolution algorithm to create a new direction bitmap and weight array for the unresolved pixels. This process looks a lot like the regular flow routing: fill up the bitmap, count the edges, allocate and fill up the weights. 6. Count the total number of edges in the originally resolved direction bitmap and the bitmap of directions on flats. edgeset_count_merged is provided for this purpose. Allocate a weight array and merge the originally resolved directions and the resolved directions on flats using the merge_edges function. You have to do this after flat resolution because of the way the edges are sorted in the weight array. There might be a better way to do it that would allow you to skip this step and merge on the fly. 7. Now you have a bitmap of directions and float array of weights for all the edges you want in your Flow object. Use the edgeset_scan function to prepare an array of type ptrdiff_t with the offsets of each pixels edges in the weight array. This is done for you by edgeset_merge if you used that. edgeset_count, edgeset_scan and edgeset_merge all return the total number of edges in the edge set. 8. Allocate your Flow object arrays: - ptrdiff_t *stream: the topologically sorted list of nodes. Size is the number of pixels in your DEM - ptrdiff_t *source: the topologically sorted list of edge sources. Size is the number of edges - ptrdiff_t *target: the topologically sorted list of edge targets. Size is the number of edges. - float *weight: the topologically sorted list of edge weights. Size is the number of edges. 9. Use direction_tsort to traverse the flow network and topologically sort the edges. The node, source, target and weight arrays are filled up in the right order. The bitmap edge set provides a uniform interface for flow routing and flat resolution algorithms to exchange data. This could also be done with topologically sorted edge lists as stored in the Flow objects. However, it is harder to merge two topological sorts than it is to merge two edge sets. The topological sort has to be redone with another graph traversal, and edge lists don't let you compute traversals from an arbitrary pixel so easily. This design lets you interchange flow routing and flat resolution algorithms as desired. Each algorithm returns an edge set that can be merged and then topologically sorted. We can, for example, decouple the D8 flow routing from the least cost auxiliary topography carving that we currently do. And you could conceivably use the LCAT carving with Dinf if you wanted. The D-inf implementation consists of two functions in src/dinf.c. - flow_routing_dinf_directions computes the direction bitmap and an array of outgoing flow proportions for each pixel - flow_routing_dinf_weights computes the edge weights from the array of flow proportions. A new function flow_routing_tsort is added to src/flow_routing.c to topologically sort bitmap edge sets. This is intended to replace the topological sorting functionality of flow_routing_d8_carve. A shortest-path flat resolution algorithm is added to src/flow_routing.c with two functions - resolve_flats_shortest_path computes the direction bitmap using Dijkstra's algorithm to route flow along flats following the shortest distance to a sill - resolve_flats_shortest_path_weights computes the flow weights for the shortest-path flat resolution. These are all 1s, but this function is provided to maintain the direction/weights interface. The shortest path resolution has been added because it is what is described by Tarboton (1997). He uses an iterative algorithm, but it amounts to the same thing. This is not what is used by default in MATLAB TopoToolbox or, as far as I can tell, in Tarboton's TauDEM. Moreover this implementation has been done entirely from scratch based on the description in Tarboton (1997) without reliance on the D-inf implementation in MATLAB TopoToolbox or TauDEM. It probably does not behave exactly as in MATLAB TopoToolbox. The flat resolution approach is certainly different. MATLAB TopoToolbox's D-inf implementation also fails on some DEMs with NaNs such as the Taiwan example. I think the problem there is when there are embedded NaNs that are not contiguous with the border. libtopotoolbox handles embedded NaNs as missing data. Because of the way D-inf computes flow angles based on triangular facets, this can lead to some surprising results, such as the following situation: 8 9 7 5 5 6 NaN 1 NaN The center pixel should probably flow down, but because both facets containing the downward neighbor have NaNs, the slope of each facet is NaN, and the flow is unresolved. The flat resolution algorithm will probably flow to the neighbor of equal elevation to the left and from there to the bottom right. This actually occurs in the Kunashiri example DEM. It could be possible to fix this by special casing the situation where one of the two neighboring pixels for a facet is NaN. Functions for dealing with the edge sets have been added to src/helpers/edgeset.c. These could conceivably also be done in higher level languages, but making sure you get all the ordering right is a bit tricky. Test cases have been added to test/dinf.c. These tests use the Grid and Flow object interface added in #221 and test on a variety of small test cases, including those from the Tarboton paper, the snapshot DEMs and randomly generated ones. Memory order tests are also conducted. I have also taken this opportunity to prototype a new testing harness that provides TAP output. Signed-off-by: William Kearney <william.kearney@uni-potsdam.de>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
See #55
These two structs are C implementations of the GRIDobj and FLOWobj interfaces to be used in the tests. They should make adding and modifying tests much easier by hiding some of the details of creating and accessing data sets. This is part of a larger project to clean up the tests.
Grid is a two dimensional array with an attached cellsize but no coordinate reference information. It provides methods for generating random data, blank datasets or loading through GDAL. It also provides arithmetic and comparison operations and TopoToolbox functions needed for flow routing.
Flow is a topologically sorted edge list with source, target and fraction fields as well as the topological node ordering in stream. FlowRouteD8Carve runs D8 with morphological carving on the provided DEM and FlowAccumulation and FlowDrainageBasins use flow algebras.
The random_dem test has been heavily modified to use the Grid and Flow structs where possible. Test functions take Grid or Flow objects, and the flow routing, gradient computation, flow accumulation, etc. has been converted to generate Grids or Flows.
Some refactoring has been carried out along the way such as
Implement a simpler topological sorting test using a Flow. The previous one used flow directions and is a relic of the older flow routing interface.
Remove calls to flow_accumulation, flow_accumulation_edgelist and drainagebasins in favor of the flow algebra implementation of Flow (see Deprecate
flow_accumulationanddrainagebasins#212). The test_flow_accumulation_multimethod test is therefore removed as there is really only one supported flow accumulation method now.Separate fillsinks/fillsinks_hybrid from the flow routing function. We previously had two versions of flow routing to be able to compare the performance of fillsinks and fillsinks_hybrid. The profiler now profiles the two fillsinks algorithms separately.
There is also a Stream implementation on the way, but it is a little more challenging to ensure the memory allocation details are correct. Once that lands, the remaining random_dem tests of the stream network can also be refactored.
GDAL is now required for the random_dem tests because Grid depends on it. It would be nice to figure out how to compile GridFromFile and GridToFile only if GDAL is available, so we can run random_dem without GDAL.