From 9115ce2f71e0056902415aa805b13afcc3b46615 Mon Sep 17 00:00:00 2001 From: William Kearney Date: Mon, 27 Apr 2026 15:49:48 +0200 Subject: [PATCH] Implement D8 flow routing following new interface Resolves #223 This adds the D8 flow routing functions - `flow_routing_d8_directions` - `flow_routing_d8_weights` in src/d8.c that implement the new staged flow routing interface for D8. This function does not resolve flow over flats. The old flat resolution behavior of flow_routing_d8_carve has been migrated to src/lcat.c (for "Least Cost Auxiliary Topography"). This includes the two functions - `resolve_flats_lcat` - `resolve_flats_lcat_weights` that use the auxiliary topography, which must still be computed separately by `gwdt_computecosts` and `gwdt`. The new implementations do NOT preserve the behavior of the previous `flow_routing_d8_carve`. The order in which neighbors are checked has been changed to be consistent with the dinf implementation and with `flow_routing_tsort`. If the neighbor-checking order is changed to that of `flow_routing_d8_carve`, the results are identical to the previous ones. The tests introduced by #224 have been modified to use the new interface, and they should currently pass. Signed-off-by: William Kearney --- include/topotoolbox.h | 139 ++++++++++++++++++++++++++++++++ src/CMakeLists.txt | 2 + src/d8.c | 47 +++++++++++ src/lcat.c | 54 +++++++++++++ test/d8.c | 180 +++++++++++++++++++++++++++++++----------- 5 files changed, 376 insertions(+), 46 deletions(-) create mode 100644 src/d8.c create mode 100644 src/lcat.c diff --git a/include/topotoolbox.h b/include/topotoolbox.h index 901c50b..51bb22e 100644 --- a/include/topotoolbox.h +++ b/include/topotoolbox.h @@ -3146,4 +3146,143 @@ void resolve_flats_shortest_path(uint8_t *direction, uint8_t *resolved, TOPOTOOLBOX_API void resolve_flats_shortest_path_weights(float *weight, ptrdiff_t count); +/** + @brief Compute flow directions using the D8 method + + @details D8 chooses the steepest downstream neighbor of the 8 + pixels surrounding the center pixel. + + @param[out] direction A bitmap edge set + @parblock + + A uint8_t array of the same size as the DEM. Each of the 8 bits of + each pixel corresponds to an outgoing edge. + + @endparblock + + @param[in] dem The input DEM + @parblock + + A float array with size dims[0] x dims[1]. + + @endparblock + + @param[in] dims The dimensions of the DEM + @parblock + A pair of ptrdiff_t, fastest changing dimension first. + @endparblock + + @param[in] order The memory order of the underlying arrays + @parblock + 0 for column-major, 1 for row-major + @endparblock + */ +TOPOTOOLBOX_API +void flow_routing_d8_directions(uint8_t *direction, float *dem, + ptrdiff_t dims[2], int order); + +/** + @brief Compute and store the D8 edge weights + + @details Once the flow directions have been computed, the edge + weights need to be stored in an array in an order that will be used + by functions like flow_routing_tsort. The weights for D8 are all + equal to 1, so this function can be replaced in applications by + filling an appropriately sized array with ones. + + @param[out] weight The weight array + @parblock + + An array of floats that will be filled with the edge weights. This + should be preallocated with a size equal to the number of edges in + the direction bitmap. This size can be computed using + edgeset_count. + + @endparblock + + @param[in] count The number of edges in the weight array + */ +TOPOTOOLBOX_API +void flow_routing_d8_weights(float *weight, ptrdiff_t count); + +/** + @brief Resolve flats using least cost auxiliary topography carving. + + @details This flat resolution algorithm uses the auxiliary + topography in `aux` to carve a path through flat areas and sinks. + + See + + Schwanghart, W., Groom, G., Kuhn, N.J. and Heckrath, G. (2013), + Flow network derivation from a high resolution DEM in a low relief, + agrarian landscape. Earth Surf. Process. Landforms, 38: + 1576-1586. https://doi.org/10.1002/esp.3452 + + and references therein for more details. + + @param[out] direction The output direction bitmap + @param[in] resolved nonzero for pixels that have resolved directions + @parblock + + An array of uint8_t flagging unresolved pixels with a zero. Only + pixels with a zero in this array will be resolved by the shortest + path algorithm. + + @endparblock + + @param aux The auxiliary topography + @parblock + + An array of floats of size dims[0] x dims[1]. + + This should be generated using `gwdt_computecosts` and `gwdt`. + @endparblock + + @param[in] dem The input dem + @parblock + + This is used to make sure that flat pixels only flow to neighbors + of the same elevation. + + @endparblock + + @param[in] dims The dimensions of the DEM + @parblock + A pair of ptrdiff_t, fastest changing dimension first. + @endparblock + + @param[in] order The memory order of the underlying arrays + @parblock + 0 for column-major, 1 for row-major + @endparblock + */ +TOPOTOOLBOX_API +void resolve_flats_lcat(uint8_t *direction, uint8_t *resolved, float *aux, + float *dem, ptrdiff_t dims[2], int order); + +/** + @brief Compute the weights of the least cost auxiliary topography + flat resolution algorithm + + @details The weights are all 1. This function is provided for + consistency with other flat resolution implementations. + + @param[out] weight The array of edge weights + @parblock + + An array of floats preallocated with a size count. + + @endparblock + + @param[in] count The number of edges + @parblock + + This can be computed from the output direction bitmap of + resolve_flats_shortest_path using edgeset_count. + + @endparblock + */ +TOPOTOOLBOX_API +void resolve_flats_lcat_weights(float *weight, ptrdiff_t count); + #endif // TOPOTOOLBOX_H diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 55289aa..ca50d90 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -38,6 +38,8 @@ add_library(topotoolbox helpers/deque.h helpers/edgeset.c dinf.c + d8.c + lcat.c ) # Define the include directory diff --git a/src/d8.c b/src/d8.c new file mode 100644 index 0000000..4cd0993 --- /dev/null +++ b/src/d8.c @@ -0,0 +1,47 @@ +#define TOPOTOOLBOX_BUILD + +#include "topotoolbox.h" + +#define SQRT2 1.41421356237309504880f + +void flow_routing_d8_directions(uint8_t *direction, float *dem, + ptrdiff_t dims[2], int order) { + // Basic D8 flow routing: flow to the maximum downstream neighbor until you + // hit a flat or a sink. + + // These are the offsets for the appropriate neighbors. + ptrdiff_t e[2][8] = {{0, -1, -1, -1, 0, 1, 1, 1}, + {1, 1, 0, -1, -1, -1, 0, 1}}; + + float chamfer[8] = {1.0f, SQRT2, 1.0f, SQRT2, 1.0f, SQRT2, 1.0f, SQRT2}; + + for (ptrdiff_t j = 0; j < dims[1]; j++) { + for (ptrdiff_t i = 0; i < dims[0]; i++) { + float z = dem[j * dims[0] + i]; + float g = 0.0; + direction[j * dims[0] + i] = 0; + + for (ptrdiff_t n = 0; n < 8; n++) { + ptrdiff_t in = i + e[order & 1][n]; + ptrdiff_t jn = j + e[(order ^ 1) & 1][n]; + + if (in < 0 || in >= dims[0] || jn < 0 || jn >= dims[1]) { + continue; + } + + float gn = (z - dem[jn * dims[0] + in]) / chamfer[n]; + + if (gn > g) { + direction[j * dims[0] + i] = 1 << n; + g = gn; + } + } + } + } +} + +void flow_routing_d8_weights(float *weight, ptrdiff_t count) { + for (ptrdiff_t e = 0; e < count; e++) { + weight[e] = 1.0; + } +} diff --git a/src/lcat.c b/src/lcat.c new file mode 100644 index 0000000..98894b0 --- /dev/null +++ b/src/lcat.c @@ -0,0 +1,54 @@ +#define TOPOTOOLBOX_BUILD + +#include "topotoolbox.h" + +#define SQRT2 1.41421356237309504880f + +void resolve_flats_lcat(uint8_t *direction, uint8_t *resolved, float *aux, + float *dem, ptrdiff_t dims[2], int order) { + // Least cost auxiliary topography carving. The hard part is + // computing the auxiliary topography in aux, but you do that + // first. Otherwise this is D8 flow routing over unresolved pixels + // using the auxiliary topography. + ptrdiff_t e[2][8] = {{0, -1, -1, -1, 0, 1, 1, 1}, + {1, 1, 0, -1, -1, -1, 0, 1}}; + + float chamfer[8] = {1.0f, SQRT2, 1.0f, SQRT2, 1.0f, SQRT2, 1.0f, SQRT2}; + + for (ptrdiff_t j = 0; j < dims[1]; j++) { + for (ptrdiff_t i = 0; i < dims[0]; i++) { + if (resolved[j * dims[0] + i]) continue; + + float z = dem[j * dims[0] + i]; + float t = aux[j * dims[0] + i]; + float g = 0.0; + + for (ptrdiff_t n = 0; n < 8; n++) { + ptrdiff_t in = i + e[order & 1][n]; + ptrdiff_t jn = j + e[(order ^ 1) & 1][n]; + + if (in < 0 || in >= dims[0] || jn < 0 || jn >= dims[1]) { + continue; + } + + // Skip pixels that are greater than the current elevation + if (dem[jn * dims[0] + in] > z) { + continue; + } + + float gn = (t - aux[jn * dims[0] + in]) / chamfer[n]; + + if (gn > g) { + direction[j * dims[0] + i] = 1 << n; + g = gn; + } + } + } + } +} + +void resolve_flats_lcat_weights(float *weight, ptrdiff_t count) { + for (ptrdiff_t e = 0; e < count; e++) { + weight[e] = 1.0; + } +} diff --git a/test/d8.c b/test/d8.c index 448807a..4ee64b5 100644 --- a/test/d8.c +++ b/test/d8.c @@ -213,35 +213,94 @@ int test_sum_edge_weights(Flow fd) { return res; } -Flow d8_carve_flow(Grid dem, int order) { +int test_one_neighbor(Flow fd) { + // D8 networks should have a maximum of one downstream neighbor. + + int res = 1; + Grid neighbors = GridCreate(GridU8, NULL, fd.cellsize, fd.dims); + + for (ptrdiff_t e = 0; e < fd.count; e++) { + ptrdiff_t u = fd.source[e]; + + if (((uint8_t *)neighbors.data)[u] > 0) { + res = 0; + } + ((uint8_t *)neighbors.data)[u] += 1; + } + + GridFree(&neighbors); + return res; +} + +Flow d8_flow(Grid dem, int order) { Flow fd = {0}; + Grid direction = GridCreate(GridU8, NULL, dem.cellsize, dem.dims); + flow_routing_d8_directions(direction.data, dem.data, dem.dims, order); + + ptrdiff_t count = edgeset_count(direction.data, direction.dims); + float *weight = calloc(count, sizeof *weight); + flow_routing_d8_weights(weight, count); + Grid demf = GridFillsinks(dem); Grid flats = GridIdentifyFlats(dem); Grid costs = GridGWDTComputeCosts(flats, dem, demf); - Grid dist = GridGWDT(costs, flats); + Grid aux = GridGWDT(costs, flats); + + Grid rdirection = GridCreate(GridU8, NULL, dem.cellsize, dem.dims); + Grid resolved = GridCreate(GridU8, NULL, dem.cellsize, dem.dims); + + // Resolved pixels are either NaNs, or they have a flow direction + for (ptrdiff_t j = 0; j < dem.dims[1]; j++) { + for (ptrdiff_t i = 0; i < dem.dims[0]; i++) { + if (isnan(((float *)dem.data)[j * dem.dims[0] + i])) { + // If DEM is a NaN, treat the flow direction as resolved + ((uint8_t *)resolved.data)[j * dem.dims[0] + i] = 0xff; + } else if (isinf(((float *)aux.data)[j * dem.dims[0] + i])) { + // If aux is infinite, this is a sink/flat that has no sills. + ((uint8_t *)resolved.data)[j * dem.dims[0] + i] = 0xff; + } else { + ((uint8_t *)resolved.data)[j * dem.dims[0] + i] = + ((uint8_t *)direction.data)[j * dem.dims[0] + i] ? 0xff : 0; + } + } + } - Grid stream = GridCreate(GridIdx, NULL, dem.cellsize, dem.dims); - Grid direction = GridCreate(GridU8, NULL, dem.cellsize, dem.dims); + resolve_flats_lcat(rdirection.data, resolved.data, aux.data, demf.data, + dem.dims, order); - flow_routing_d8_carve(stream.data, direction.data, dem.data, dist.data, - flats.data, dem.dims, order); + ptrdiff_t rcount = edgeset_count(rdirection.data, rdirection.dims); + float *rweight = calloc(rcount, sizeof *rweight); - ptrdiff_t *source = calloc(GridElementCount(dem), sizeof *source); - ptrdiff_t *target = calloc(GridElementCount(dem), sizeof *target); - ptrdiff_t edge_count = flow_routing_d8_edgelist( - source, target, stream.data, direction.data, dem.dims, order); + resolve_flats_lcat_weights(rweight, rcount); - float *weight = calloc(edge_count, sizeof *weight); - for (ptrdiff_t e = 0; e < edge_count; e++) { - weight[e] = 1.0; - } + ptrdiff_t merged_count = + edgeset_count_merged(direction.data, rdirection.data, direction.dims); + + float *merged_weights = calloc(merged_count, sizeof *merged_weights); + Grid mergedscan = GridCreate(GridIdx, NULL, dem.cellsize, dem.dims); + + edgeset_merge(merged_weights, mergedscan.data, direction.data, weight, + rdirection.data, rweight, direction.dims); + + Grid stream = GridCreate(GridIdx, NULL, dem.cellsize, dem.dims); + ptrdiff_t *source = calloc(merged_count, sizeof *source); + ptrdiff_t *target = calloc(merged_count, sizeof *target); + float *sorted_weight = calloc(merged_count, sizeof *sorted_weight); + + Grid stack = GridCreate(GridIdx, NULL, dem.cellsize, dem.dims); + Grid stackdir = GridCreate(GridU8, NULL, dem.cellsize, dem.dims); + Grid visited = GridCreate(GridU8, NULL, dem.cellsize, dem.dims); + flow_routing_tsort(stream.data, source, target, sorted_weight, stack.data, + stackdir.data, direction.data, merged_weights, + mergedscan.data, visited.data, merged_count, dem.dims, + order); fd.stream = stream.data; fd.source = source; fd.target = target; - fd.fraction = weight; - fd.count = edge_count; + fd.fraction = sorted_weight; + fd.count = merged_count; fd.cellsize = dem.cellsize; fd.dims[0] = dem.dims[0]; fd.dims[1] = dem.dims[1]; @@ -249,7 +308,21 @@ Flow d8_carve_flow(Grid dem, int order) { GridFree(&demf); GridFree(&flats); GridFree(&costs); - GridFree(&dist); + GridFree(&aux); + GridFree(&resolved); + + GridFree(&visited); + GridFree(&stackdir); + GridFree(&stack); + + free(merged_weights); + GridFree(&mergedscan); + + free(rweight); + GridFree(&rdirection); + + free(weight); + GridFree(&direction); return fd; @@ -269,7 +342,7 @@ int main(int argc, char *argv[]) { ptrdiff_t dims[2] = {3, 3}; float data[9] = {2, 2, 2, 1, 1, 1, 0, 0, 0}; Grid dem = GridCreate(GridF32, data, 1.0, dims); - Flow fd = d8_carve_flow(dem, 0); + Flow fd = d8_flow(dem, 0); if (!test_node_once(fd)) test_fail; if (!test_tsort_node(fd)) test_fail; @@ -278,6 +351,7 @@ int main(int argc, char *argv[]) { if (!test_no_boundary_edges(fd)) test_fail; if (!test_outlets(fd)) test_fail; if (!test_sum_edge_weights(fd)) test_fail; + if (!test_one_neighbor(fd)) test_fail; FlowFree(&fd); } @@ -287,7 +361,7 @@ int main(int argc, char *argv[]) { float data[25] = {4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0}; Grid dem = GridCreate(GridF32, data, 1.0, dims); - Flow fd = d8_carve_flow(dem, 0); + Flow fd = d8_flow(dem, 0); if (!test_node_once(fd)) test_fail; if (!test_tsort_node(fd)) test_fail; @@ -296,6 +370,7 @@ int main(int argc, char *argv[]) { if (!test_no_boundary_edges(fd)) test_fail; if (!test_outlets(fd)) test_fail; if (!test_sum_edge_weights(fd)) test_fail; + if (!test_one_neighbor(fd)) test_fail; FlowFree(&fd); } @@ -304,7 +379,7 @@ int main(int argc, char *argv[]) { ptrdiff_t dims[2] = {3, 3}; float data[9] = {2, 3, 4, 1, 2, 3, 0, 1, 2}; Grid dem = GridCreate(GridF32, data, 1.0, dims); - Flow fd = d8_carve_flow(dem, 0); + Flow fd = d8_flow(dem, 0); if (!test_node_once(fd)) test_fail; if (!test_tsort_node(fd)) test_fail; @@ -313,6 +388,7 @@ int main(int argc, char *argv[]) { if (!test_no_boundary_edges(fd)) test_fail; if (!test_outlets(fd)) test_fail; if (!test_sum_edge_weights(fd)) test_fail; + if (!test_one_neighbor(fd)) test_fail; FlowFree(&fd); } @@ -321,7 +397,7 @@ int main(int argc, char *argv[]) { ptrdiff_t dims[2] = {3, 3}; float data[9] = {2, 2, 2, 2, 1, 2, 2, 2, 2}; Grid dem = GridCreate(GridF32, data, 1.0, dims); - Flow fd = d8_carve_flow(dem, 0); + Flow fd = d8_flow(dem, 0); if (!test_node_once(fd)) test_fail; if (!test_tsort_node(fd)) test_fail; @@ -330,12 +406,13 @@ int main(int argc, char *argv[]) { if (!test_no_boundary_edges(fd)) test_fail; if (!test_outlets(fd)) test_fail; if (!test_sum_edge_weights(fd)) test_fail; + if (!test_one_neighbor(fd)) test_fail; FlowFree(&fd); } test("Outward sloping cone") { - ptrdiff_t dims[2] = {16, 16}; + ptrdiff_t dims[2] = {17, 17}; float cellsize = 150.0 / (dims[0] - 1); Grid dem = GridCreate(GridF32, NULL, cellsize, dims); for (ptrdiff_t j = 0; j < dims[1]; j++) { @@ -346,7 +423,7 @@ int main(int argc, char *argv[]) { } } - Flow fd = d8_carve_flow(dem, 0); + Flow fd = d8_flow(dem, 0); if (!test_node_once(fd)) test_fail; if (!test_tsort_node(fd)) test_fail; @@ -355,13 +432,14 @@ int main(int argc, char *argv[]) { if (!test_no_boundary_edges(fd)) test_fail; if (!test_outlets(fd)) test_fail; if (!test_sum_edge_weights(fd)) test_fail; + if (!test_one_neighbor(fd)) test_fail; FlowFree(&fd); GridFree(&dem); } test("Outward sloping cone with flats") { - ptrdiff_t dims[2] = {16, 16}; + ptrdiff_t dims[2] = {17, 17}; float cellsize = 150.0 / (dims[0] - 1); Grid dem = GridCreate(GridF32, NULL, cellsize, dims); for (ptrdiff_t j = 0; j < dims[1]; j++) { @@ -373,7 +451,7 @@ int main(int argc, char *argv[]) { } } - Flow fd = d8_carve_flow(dem, 0); + Flow fd = d8_flow(dem, 0); if (!test_node_once(fd)) test_fail; if (!test_tsort_node(fd)) test_fail; @@ -382,6 +460,7 @@ int main(int argc, char *argv[]) { if (!test_no_boundary_edges(fd)) test_fail; if (!test_outlets(fd)) test_fail; if (!test_sum_edge_weights(fd)) test_fail; + if (!test_one_neighbor(fd)) test_fail; FlowFree(&fd); GridFree(&dem); @@ -399,7 +478,7 @@ int main(int argc, char *argv[]) { Grid dem = GridFromFile(path); free(path); Grid demf = GridFillsinks(dem); - Flow fd = d8_carve_flow(demf, 0); + Flow fd = d8_flow(demf, 0); if (!test_node_once(fd)) test_fail; if (!test_tsort_node(fd)) test_fail; @@ -408,6 +487,7 @@ int main(int argc, char *argv[]) { if (!test_no_boundary_edges(fd)) test_fail; if (!test_outlets(fd)) test_fail; if (!test_sum_edge_weights(fd)) test_fail; + if (!test_one_neighbor(fd)) test_fail; FlowFree(&fd); GridFree(&demf); @@ -425,7 +505,7 @@ int main(int argc, char *argv[]) { Grid dem = GridFromFile(path); free(path); Grid demf = GridFillsinks(dem); - Flow fd = d8_carve_flow(demf, 0); + Flow fd = d8_flow(demf, 0); if (!test_node_once(fd)) test_fail; if (!test_tsort_node(fd)) test_fail; @@ -434,6 +514,7 @@ int main(int argc, char *argv[]) { if (!test_no_boundary_edges(fd)) test_fail; if (!test_outlets(fd)) test_fail; if (!test_sum_edge_weights(fd)) test_fail; + if (!test_one_neighbor(fd)) test_fail; FlowFree(&fd); GridFree(&demf); @@ -451,7 +532,7 @@ int main(int argc, char *argv[]) { Grid dem = GridFromFile(path); free(path); Grid demf = GridFillsinks(dem); - Flow fd = d8_carve_flow(demf, 0); + Flow fd = d8_flow(demf, 0); if (!test_node_once(fd)) test_fail; if (!test_tsort_node(fd)) test_fail; @@ -460,6 +541,7 @@ int main(int argc, char *argv[]) { if (!test_no_boundary_edges(fd)) test_fail; if (!test_outlets(fd)) test_fail; if (!test_sum_edge_weights(fd)) test_fail; + if (!test_one_neighbor(fd)) test_fail; FlowFree(&fd); GridFree(&demf); @@ -477,7 +559,7 @@ int main(int argc, char *argv[]) { Grid dem = GridFromFile(path); free(path); Grid demf = GridFillsinks(dem); - Flow fd = d8_carve_flow(demf, 0); + Flow fd = d8_flow(demf, 0); if (!test_node_once(fd)) test_fail; if (!test_tsort_node(fd)) test_fail; @@ -486,6 +568,7 @@ int main(int argc, char *argv[]) { if (!test_no_boundary_edges(fd)) test_fail; if (!test_outlets(fd)) test_fail; if (!test_sum_edge_weights(fd)) test_fail; + if (!test_one_neighbor(fd)) test_fail; FlowFree(&fd); GridFree(&demf); @@ -503,7 +586,7 @@ int main(int argc, char *argv[]) { Grid dem = GridFromFile(path); free(path); Grid demf = GridFillsinks(dem); - Flow fd = d8_carve_flow(demf, 0); + Flow fd = d8_flow(demf, 0); if (!test_node_once(fd)) test_fail; if (!test_tsort_node(fd)) test_fail; @@ -512,6 +595,7 @@ int main(int argc, char *argv[]) { if (!test_no_boundary_edges(fd)) test_fail; if (!test_outlets(fd)) test_fail; if (!test_sum_edge_weights(fd)) test_fail; + if (!test_one_neighbor(fd)) test_fail; FlowFree(&fd); GridFree(&demf); @@ -529,7 +613,7 @@ int main(int argc, char *argv[]) { Grid dem = GridFromFile(path); free(path); Grid demf = GridFillsinks(dem); - Flow fd = d8_carve_flow(demf, 0); + Flow fd = d8_flow(demf, 0); if (!test_node_once(fd)) test_fail; if (!test_tsort_node(fd)) test_fail; @@ -538,6 +622,7 @@ int main(int argc, char *argv[]) { if (!test_no_boundary_edges(fd)) test_fail; if (!test_outlets(fd)) test_fail; if (!test_sum_edge_weights(fd)) test_fail; + if (!test_one_neighbor(fd)) test_fail; FlowFree(&fd); GridFree(&demf); @@ -556,7 +641,7 @@ int main(int argc, char *argv[]) { free(path); Grid demf = GridFillsinks(dem); - Flow fd = d8_carve_flow(demf, 0); + Flow fd = d8_flow(demf, 0); if (!test_node_once(fd)) test_fail; if (!test_tsort_node(fd)) test_fail; @@ -565,6 +650,7 @@ int main(int argc, char *argv[]) { if (!test_no_boundary_edges(fd)) test_fail; if (!test_outlets(fd)) test_fail; if (!test_sum_edge_weights(fd)) test_fail; + if (!test_one_neighbor(fd)) test_fail; FlowFree(&fd); GridFree(&demf); @@ -586,7 +672,7 @@ int main(int argc, char *argv[]) { } ((float *)dem.data)[2 * dims[0] + 2] = NAN; - Flow fd = d8_carve_flow(dem, 0); + Flow fd = d8_flow(dem, 0); if (!test_node_once(fd)) test_fail; if (!test_tsort_node(fd)) test_fail; @@ -595,6 +681,7 @@ int main(int argc, char *argv[]) { if (!test_no_boundary_edges(fd)) test_fail; if (!test_outlets(fd)) test_fail; if (!test_sum_edge_weights(fd)) test_fail; + if (!test_one_neighbor(fd)) test_fail; FlowFree(&fd); GridFree(&dem); @@ -619,8 +706,8 @@ int main(int argc, char *argv[]) { Grid cdemf = GridFillsinks(cdem); GridFree(&cdem); - Flow cfd = d8_carve_flow(cdemf, 1); - Flow ffd = d8_carve_flow(fdemf, 0); + Flow cfd = d8_flow(cdemf, 1); + Flow ffd = d8_flow(fdemf, 0); Grid cacc = FlowAccumulation(cfd); Grid facc = FlowAccumulation(ffd); @@ -631,11 +718,12 @@ int main(int argc, char *argv[]) { GridFree(&caccT); GridFree(&cacc); GridFree(&facc); - GridFree(&fdemf); - GridFree(&cdemf); - FlowFree(&ffd); FlowFree(&cfd); + FlowFree(&ffd); + + GridFree(&fdemf); + GridFree(&cdemf); } test("Memory order saddle (ties)") { @@ -647,8 +735,8 @@ int main(int argc, char *argv[]) { Grid fdem = GridCreate(GridF32, fdem_buffer, 1.0, fdims); Grid cdem = GridTranspose(fdem); - Flow cfd = d8_carve_flow(cdem, 1); - Flow ffd = d8_carve_flow(fdem, 0); + Flow cfd = d8_flow(cdem, 1); + Flow ffd = d8_flow(fdem, 0); Grid cacc = FlowAccumulation(cfd); Grid facc = FlowAccumulation(ffd); @@ -659,10 +747,11 @@ int main(int argc, char *argv[]) { GridFree(&caccT); GridFree(&cacc); GridFree(&facc); - GridFree(&cdem); - FlowFree(&ffd); FlowFree(&cfd); + FlowFree(&ffd); + + GridFree(&cdem); } test("Memory order random") { @@ -677,12 +766,11 @@ int main(int argc, char *argv[]) { Grid cdemf = GridFillsinks(cdem); GridFree(&cdem); - Flow cfd = d8_carve_flow(cdemf, 1); - Flow ffd = d8_carve_flow(fdemf, 0); + Flow cfd = d8_flow(cdem, 1); + Flow ffd = d8_flow(fdem, 0); Grid cacc = FlowAccumulation(cfd); Grid facc = FlowAccumulation(ffd); - Grid caccT = GridTranspose(cacc); if (!GridAllClose(caccT, facc, 1e-5)) test_fail; @@ -691,8 +779,8 @@ int main(int argc, char *argv[]) { GridFree(&cacc); GridFree(&facc); - FlowFree(&ffd); FlowFree(&cfd); + FlowFree(&ffd); GridFree(&fdemf); GridFree(&cdemf);