From 20527d2d108ea918abf1e016d66ab6ad7cc25216 Mon Sep 17 00:00:00 2001 From: William Kearney Date: Wed, 1 Apr 2026 15:58:16 +0200 Subject: [PATCH 1/3] Add Grid and Flow implementations to the tests See #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 #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 --- test/CMakeLists.txt | 23 +- test/flow.c | 143 ++++++++++ test/flow.h | 31 +++ test/grid.c | 636 ++++++++++++++++++++++++++++++++++++++++++++ test/grid.h | 64 +++++ test/random_dem.cpp | 521 ++++++++++++++++-------------------- 6 files changed, 1123 insertions(+), 295 deletions(-) create mode 100644 test/flow.c create mode 100644 test/flow.h create mode 100644 test/grid.c create mode 100644 test/grid.h diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index f30fb424..16963d75 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -31,16 +31,20 @@ set_tests_properties(versioninfo PROPERTIES ENVIRONMENT_MODIFICATION # TEST : random_dem # # Runs libtopotoolbox functions on randomly generated DEMs. -add_executable(random_dem random_dem.cpp utils.c utils.h profiler.h) -if(TT_SANITIZE AND NOT MSVC) - # Set up the AddressSanitizer flags. - target_compile_options(random_dem PRIVATE "$<$:-fsanitize=address>") - target_link_options(random_dem PRIVATE "$<$:-fsanitize=address>") +include(FindGDAL) + +if (GDAL_FOUND) + add_executable(random_dem random_dem.cpp utils.c utils.h profiler.h grid.c flow.c) + if(TT_SANITIZE AND NOT MSVC) + # Set up the AddressSanitizer flags. + target_compile_options(random_dem PRIVATE "$<$:-fsanitize=address>") + target_link_options(random_dem PRIVATE "$<$:-fsanitize=address>") + endif() + target_link_libraries(random_dem PRIVATE topotoolbox GDAL::GDAL) + add_test(NAME random_dem COMMAND random_dem) + set_tests_properties(random_dem PROPERTIES ENVIRONMENT_MODIFICATION + "PATH=path_list_prepend:$<$:$>") endif() -target_link_libraries(random_dem PRIVATE topotoolbox) -add_test(NAME random_dem COMMAND random_dem) -set_tests_properties(random_dem PROPERTIES ENVIRONMENT_MODIFICATION - "PATH=path_list_prepend:$<$:$>") # TEST : excesstopography # @@ -136,7 +140,6 @@ else() message(STATUS "Could not find snapshot data. Snapshot tests will not be run") endif() -include(FindGDAL) if(GDAL_FOUND AND TT_SNAPSHOT_DIR) add_executable(snapshot snapshot.cpp utils.c utils.h profiler.h) if(TT_SANITIZE AND NOT MSVC) diff --git a/test/flow.c b/test/flow.c new file mode 100644 index 00000000..7ea50c3f --- /dev/null +++ b/test/flow.c @@ -0,0 +1,143 @@ +#include "flow.h" + +#include + +#include "../include/topotoolbox.h" +#include "grid.h" + +void FlowFree(Flow *f) { + free(f->source); + free(f->target); + free(f->stream); + free(f->fraction); + + f->source = NULL; + f->target = NULL; + f->stream = NULL; + f->fraction = NULL; + f->dims[0] = 0; + f->dims[1] = 0; + f->count = 0; + f->cellsize = 0.0; +} + +Flow FlowInit(Grid dem) { + Flow fd = {0}; + fd.cellsize = dem.cellsize; + fd.dims[0] = dem.dims[0]; + fd.dims[1] = dem.dims[1]; + + fd.stream = calloc(GridElementCount(dem), sizeof *fd.stream); + fd.source = calloc(GridElementCount(dem), sizeof *fd.source); + fd.target = calloc(GridElementCount(dem), sizeof *fd.target); + fd.fraction = calloc(GridElementCount(dem), sizeof *fd.fraction); + + return fd; +} + +Flow FlowRouteD8Carve(Grid dem) { + Flow fd = FlowInit(dem); + + Grid demf = GridFillsinks(dem); + Grid flats = GridIdentifyFlats(demf); + Grid costs = GridGWDTComputeCosts(flats, dem, demf); + Grid dist = GridGWDT(costs, flats); + + if (demf.data && flats.data && costs.data && dist.data) { + Grid direction = GridCreate(GridU8, NULL, dem.cellsize, dem.dims); + if (direction.data && fd.stream) { + flow_routing_d8_carve(fd.stream, direction.data, demf.data, dist.data, + flats.data, dem.dims, 1); + + if (fd.source && fd.target && fd.fraction) { + fd.count = flow_routing_d8_edgelist(fd.source, fd.target, fd.stream, + direction.data, dem.dims, 1); + for (ptrdiff_t e = 0; e < fd.count; e++) { + fd.fraction[e] = 1.0; + } + } + } + GridFree(&direction); + } + + GridFree(&demf); + GridFree(&flats); + GridFree(&costs); + GridFree(&dist); + return fd; +} + +Grid FlowAccumulation(Flow f) { + Grid a = GridFill(1.0, f.cellsize, f.dims); + + if (a.data && f.count > 0) { + traverse_down_f32_add_mul(a.data, f.fraction, f.source, f.target, f.count); + } + + return a; +} + +Grid Outlets(Flow f) { + Grid outlets = GridCreate(GridU8, NULL, f.cellsize, f.dims); + if (outlets.data && f.count > 0) { + uint8_t *outletdata = outlets.data; + + Grid outdegree = GridCreate(GridU8, NULL, f.cellsize, f.dims); + if (outdegree.data) { + uint8_t *outdegreedata = outdegree.data; + + ptrdiff_t dims[2] = {f.dims[0], f.dims[1]}; + + edgelist_degree(outletdata, outdegreedata, f.source, f.target, + GridElementCount(outlets), f.count); + + for (ptrdiff_t j = 0; j < outlets.dims[1]; j++) { + for (ptrdiff_t i = 0; i < outlets.dims[0]; i++) { + outletdata[j * dims[0] + i] = (outletdata[j * dims[0] + i] > 0) & + (outdegreedata[j * dims[0] + i] == 0); + } + } + + GridFree(&outdegree); + } + } + return outlets; +} + +Grid FlowDrainageBasins(Flow f) { + Grid basins = GridCreate(GridU32, NULL, f.cellsize, f.dims); + + if (f.count > 0 && basins.data) { + uint32_t *basinsdata = basins.data; + + uint32_t *weights = calloc(f.count, sizeof *weights); + + if (weights) { + for (ptrdiff_t e = 0; e < f.count; e++) { + weights[e] = 0xffffffff; + } + + Grid outlets = Outlets(f); + if (outlets.data) { + uint8_t *outdata = outlets.data; + + uint32_t basincount = 1; + for (ptrdiff_t j = 0; j < basins.dims[1]; j++) { + for (ptrdiff_t i = 0; i < basins.dims[0]; i++) { + if (outdata[j * basins.dims[0] + i]) { + basinsdata[j * basins.dims[0] + i] = basincount; + basincount++; + } + } + } + + traverse_up_u32_or_and(basinsdata, weights, f.source, f.target, + f.count); + + GridFree(&outlets); + } + free(weights); + } + } + return basins; +} diff --git a/test/flow.h b/test/flow.h new file mode 100644 index 00000000..2d7219ef --- /dev/null +++ b/test/flow.h @@ -0,0 +1,31 @@ +#ifndef TT3_TEST_FLOW_H +#define TT3_TEST_FLOW_H + +#include + +#include "grid.h" + +typedef struct Flow { + ptrdiff_t *source; + ptrdiff_t *target; + ptrdiff_t *stream; + float *fraction; + ptrdiff_t count; + ptrdiff_t dims[2]; + float cellsize; +} Flow; + +// Destructor +void FlowFree(Flow *f); + +// Creation +Flow FlowInit(Grid g); + +// Flow routing +Flow FlowRouteD8Carve(Grid g); + +// TopoToolbox +Grid FlowAccumulation(Flow f); +Grid FlowDrainageBasins(Flow f); + +#endif diff --git a/test/grid.c b/test/grid.c new file mode 100644 index 00000000..6ffb780a --- /dev/null +++ b/test/grid.c @@ -0,0 +1,636 @@ +#include "grid.h" + +#include +#include +#include +#include +#include + +#include "../include/topotoolbox.h" +#include "utils.h" + +void GridFree(Grid *g) { + free(g->data); + g->data = NULL; + g->dims[0] = 0; + g->dims[0] = 0; + g->cellsize = 0.0; +} + +Grid GridCreate(GridDataType dtype, void *buffer, float cellsize, + ptrdiff_t dims[2]) { + Grid res = {0}; + + res.cellsize = cellsize; + res.dtype = dtype; + res.dims[0] = dims[0]; + res.dims[1] = dims[1]; + + if (res.dims[0] > 0 && (res.dims[1] > PTRDIFF_MAX / res.dims[0])) { + return res; + } + + if (buffer == NULL) { + buffer = calloc(GridElementCount(res), GridElementSize(res)); + } + res.data = buffer; + + return res; +} + +Grid GridFromFile(const char *path) { + Grid res = {0}; + + const GDALAccess eAccess = GA_ReadOnly; + + GDALDatasetH hDataset = GDALOpen(path, eAccess); + + if (hDataset == NULL) { + return res; + } + + GDALRasterBandH hBand = GDALGetRasterBand(hDataset, 1); + + ptrdiff_t nx = GDALGetRasterBandXSize(hBand); + ptrdiff_t ny = GDALGetRasterBandYSize(hBand); + + res.dims[0] = nx; + res.dims[1] = ny; + + if (nx > 0 && (ny > PTRDIFF_MAX / nx)) { + return res; + } + + double adfTransform[6]; + GDALGetGeoTransform(hDataset, adfTransform); + + if (fabs(fabs(adfTransform[1]) - fabs(adfTransform[5])) > 1e-6) { + return res; + } + + res.cellsize = adfTransform[1]; + res.dtype = GridF32; + + float *buffer = calloc(GridElementCount(res), GridElementSize(res)); + + if (buffer == NULL) { + return res; + } + + CPLErr err = GDALRasterIO(hBand, GF_Read, 0, 0, nx, ny, buffer, nx, ny, + GDT_Float32, 0, 0); + if (err) { + return res; + } + + res.data = buffer; + + GDALClose(hDataset); + + return res; +} + +GDALDataType GridGDALDataType(Grid g) { + switch (g.dtype) { + case GridF32: + return GDT_Float32; + case GridU8: + return GDT_Byte; + case GridU32: + return GDT_UInt32; + case GridI32: + return GDT_Int32; + case GridIdx: + return GDT_Int64; + default: + return GDT_Unknown; + } +} + +int GridToFile(Grid g, const char *dst) { + GDALDriverH driver = GDALGetDriverByName("GTiff"); + if (driver == NULL) return -1; + + GDALDataType dtype = GridGDALDataType(g); + GDALDatasetH ds = + GDALCreate(driver, dst, g.dims[0], g.dims[1], 1, dtype, NULL); + + if (ds == NULL) return -1; + + double gt[6] = {0, g.cellsize, 0, 0, 0, -g.cellsize}; + if (GDALSetGeoTransform(ds, gt)) return -1; + + GDALRasterBandH band = GDALGetRasterBand(ds, 1); + + if (band == NULL) return -1; + + CPLErr err = GDALRasterIO(band, GF_Write, 0, 0, g.dims[0], g.dims[1], g.data, + g.dims[0], g.dims[1], dtype, 0, 0); + if (err) return -1; + + GDALClose(ds); + + return 0; +} + +Grid GridRandom(uint64_t seed, float cellsize, ptrdiff_t dims[2]) { + Grid res = {0}; + + ptrdiff_t nx = dims[0]; + ptrdiff_t ny = dims[1]; + + res.dims[0] = nx; + res.dims[1] = ny; + + if (nx > 0 && (ny > PTRDIFF_MAX / nx)) { + return res; + } + + res.dtype = GridF32; + + float *buffer = calloc(GridElementCount(res), GridElementSize(res)); + + if (buffer == NULL) { + return res; + } + + uint32_t c = seed >> 32; + uint32_t d = seed & ((1UL << 32) - 1); + for (ptrdiff_t j = 0; j < dims[1]; j++) { + for (ptrdiff_t i = 0; i < dims[0]; i++) { + buffer[j * dims[0] + i] = 100.0f * pcg4d(i, j, c, d); + } + } + + res.data = buffer; + + return res; +} + +size_t GridElementSize(Grid g) { + switch (g.dtype) { + case GridF32: + return sizeof(float); + case GridU8: + return sizeof(uint8_t); + case GridU32: + return sizeof(uint32_t); + case GridI32: + return sizeof(int32_t); + case GridIdx: + return sizeof(ptrdiff_t); + default: + return 0; + } +} + +ptrdiff_t GridElementCount(Grid g) { return g.dims[0] * g.dims[1]; } + +size_t GridSize(Grid g) { return GridElementSize(g) * GridElementCount(g); } + +Grid GridCopy(Grid g) { + Grid res = {0}; + + res.dims[0] = g.dims[0]; + res.dims[1] = g.dims[1]; + + res.cellsize = g.cellsize; + + void *buffer = calloc(GridElementCount(g), GridElementSize(g)); + + if (buffer == NULL) { + return res; + } + + res.data = buffer; + res.dtype = g.dtype; + + memcpy(res.data, g.data, GridSize(g)); + + return res; +} + +Grid GridZeros(float cellsize, ptrdiff_t dims[2]) { + Grid res = {0}; + + res.cellsize = cellsize; + res.dtype = GridF32; + res.dims[0] = dims[0]; + res.dims[1] = dims[1]; + + if (res.dims[0] > 0 && (res.dims[1] > PTRDIFF_MAX / res.dims[0])) { + return res; + } + + void *buffer = calloc(GridElementCount(res), GridElementSize(res)); + + res.data = buffer; + + return res; +} + +Grid GridZerosLike(Grid g) { + Grid res = {0}; + + res.dims[0] = g.dims[0]; + res.dims[1] = g.dims[1]; + + res.cellsize = g.cellsize; + + res.dtype = g.dtype; + + res.data = calloc(GridElementCount(g), GridElementSize(g)); + + return res; +} + +Grid GridFill(float value, float cellsize, ptrdiff_t dims[2]) { + Grid res = {0}; + + res.cellsize = cellsize; + res.dtype = GridF32; + res.dims[0] = dims[0]; + res.dims[1] = dims[1]; + + if (res.dims[0] > 0 && (res.dims[1] > PTRDIFF_MAX / res.dims[0])) { + return res; + } + + float *buffer = calloc(GridElementCount(res), GridElementSize(res)); + + if (buffer == NULL) return res; + + for (ptrdiff_t j = 0; j < dims[1]; j++) { + for (ptrdiff_t i = 0; i < dims[0]; i++) { + buffer[j * dims[0] + i] = value; + } + } + + res.data = buffer; + + return res; +} + +Grid GridTranspose(Grid g) { + Grid res = {0}; + res.cellsize = g.cellsize; + res.dims[0] = g.dims[1]; + res.dims[1] = g.dims[0]; + res.dtype = g.dtype; + + float *buffer = calloc(GridElementCount(res), GridElementSize(res)); + + if (buffer == NULL) return res; + + for (ptrdiff_t j = 0; j < g.dims[1]; j++) { + for (ptrdiff_t i = 0; i < g.dims[0]; i++) { + buffer[i * res.dims[0] + j] = ((float *)g.data)[j * g.dims[0] + i]; + } + } + + res.data = buffer; + + return res; +} + +Grid GridSubtract(Grid g1, Grid g2) { + assert(g1.dims[0] == g2.dims[0]); + assert(g1.dims[1] == g2.dims[1]); + + // Arithmetic operations are only defined for float Grids + assert(g1.dtype == g2.dtype); + assert(g1.dtype == GridF32); + + Grid res = GridZerosLike(g1); + + if (g1.data && g2.data && res.data) { + float *r = res.data; + float *d1 = g1.data; + float *d2 = g2.data; + + for (ptrdiff_t j = 0; j < g1.dims[1]; j++) { + for (ptrdiff_t i = 0; i < g1.dims[0]; i++) { + r[j * g1.dims[0] + i] = d1[j * g1.dims[0] + i] - d2[j * g1.dims[0] + i]; + } + } + } + + return res; +} + +float GridMin(Grid g) { + assert(g.dtype == GridF32); + float *data = (float *)g.data; + + float min = INFINITY; + if (data) { + for (ptrdiff_t j = 0; j < g.dims[1]; j++) { + for (ptrdiff_t i = 0; i < g.dims[0]; i++) { + min = fminf(min, data[j * g.dims[0] + i]); + } + } + } + return min; +} + +float GridMax(Grid g) { + assert(g.dtype == GridF32); + float *data = (float *)g.data; + + float max = -INFINITY; + if (data) { + for (ptrdiff_t j = 0; j < g.dims[1]; j++) { + for (ptrdiff_t i = 0; i < g.dims[0]; i++) { + max = fmaxf(max, data[j * g.dims[0] + i]); + } + } + } + return max; +} + +int GridAllNonnegative(Grid g) { + assert(g.dtype == GridF32); + float *data = (float *)g.data; + + if (!data) return 1; + + for (ptrdiff_t j = 0; j < g.dims[1]; j++) { + for (ptrdiff_t i = 0; i < g.dims[0]; i++) { + if (data[j * g.dims[0] + i] < 0) { + return 0; + } + } + } + return 1; +} + +int GridAll(Grid g) { + assert(g.dtype == GridF32); + float *data = (float *)g.data; + + if (!data) return 0; + + for (ptrdiff_t j = 0; j < g.dims[1]; j++) { + for (ptrdiff_t i = 0; i < g.dims[0]; i++) { + if (data[j * g.dims[0] + i] == 0.0) { + return 0; + } + } + } + return 1; +} + +int GridAny(Grid g) { + assert(g.dtype == GridF32); + float *data = (float *)g.data; + + if (!data) return 0; + + for (ptrdiff_t j = 0; j < g.dims[1]; j++) { + for (ptrdiff_t i = 0; i < g.dims[0]; i++) { + if (data[j * g.dims[0] + i] != 0.0) { + return 1; + } + } + } + return 0; +} + +int GridEq(Grid g1, Grid g2) { + assert(g1.dtype == GridF32); + assert(g2.dtype == GridF32); + + if (g1.dims[0] != g2.dims[0] || g1.dims[1] != g2.dims[1]) { + return 1; + } + + float *data1 = g1.data; + float *data2 = g2.data; + + if (!data1 || !data2) return 0; + + for (ptrdiff_t j = 0; j < g1.dims[1]; j++) { + for (ptrdiff_t i = 0; i < g2.dims[0]; i++) { + if (data1[j * g1.dims[0] + i] != data2[j * g1.dims[0] + i]) { + return 0; + } + } + } + return 1; +} + +int GridAllClose(Grid g1, Grid g2, double rtol) { + assert(g1.dtype == GridF32); + assert(g2.dtype == GridF32); + + if (g1.dims[0] != g2.dims[0] || g1.dims[1] != g2.dims[1]) { + return 1; + } + + float *data1 = g1.data; + float *data2 = g2.data; + + if (!data1 || !data2) return 0; + + for (ptrdiff_t j = 0; j < g1.dims[1]; j++) { + for (ptrdiff_t i = 0; i < g2.dims[0]; i++) { + float z1 = data1[j * g1.dims[0] + i]; + float z2 = data2[j * g1.dims[0] + i]; + if (fabsf(z1 - z2) > fabsf(z2) * rtol) { + return 0; + } + } + } + + return 1; +} + +Grid GridInvert(Grid g) { + Grid res = GridCopy(g); + + assert(res.dtype == GridF32); + + if (res.data) { + float *data = (float *)res.data; + for (ptrdiff_t j = 0; j < g.dims[1]; j++) { + for (ptrdiff_t i = 0; i < g.dims[0]; i++) { + data[j * g.dims[0] + i] *= -1; + } + } + } + return res; +} + +Grid GridNot(Grid g) { + Grid res = GridCopy(g); + + assert(res.dtype == GridF32); + + if (res.data) { + float *data = (float *)res.data; + for (ptrdiff_t j = 0; j < g.dims[1]; j++) { + for (ptrdiff_t i = 0; i < g.dims[0]; i++) { + data[j * g.dims[0] + i] = 1 - data[j * g.dims[0] + i]; + } + } + } + return res; +} + +Grid GridIsNan(Grid g) { + Grid res = GridCopy(g); + + if (res.data) { + float *data = res.data; + for (ptrdiff_t j = 0; j < g.dims[1]; j++) { + for (ptrdiff_t i = 0; i < g.dims[0]; i++) { + data[j * g.dims[0] + i] = isnan(data[j * g.dims[0] + i]); + } + } + } + + return res; +} + +Grid GridFillsinks(Grid g) { + assert(g.dtype == GridF32); + float *data = g.data; + + Grid g2 = GridZerosLike(g); + + if (g2.data) { + Grid bc = GridCreate(GridU8, NULL, g.cellsize, g.dims); + + if (bc.data) { + for (ptrdiff_t j = 0; j < g.dims[1]; j++) { + for (ptrdiff_t i = 0; i < g.dims[0]; i++) { + if (i == 0 || i == g.dims[0] - 1 || j == 0 || j == g.dims[1] - 1) { + ((uint8_t *)bc.data)[j * g.dims[0] + i] = 1; + } + if (isnan(data[j * g.dims[0] + i])) { + ((uint8_t *)bc.data)[j * g.dims[0] + i] = 1; + data[j * g.dims[0] + i] = -INFINITY; + } + } + } + + Grid queue = GridCreate(GridIdx, NULL, g.cellsize, g.dims); + if (queue.data) { + fillsinks_hybrid(g2.data, queue.data, data, bc.data, g.dims); + + // Restore NaNs + for (ptrdiff_t j = 0; j < g.dims[1]; j++) { + for (ptrdiff_t i = 0; i < g.dims[0]; i++) { + if (data[j * g.dims[0] + i] == -INFINITY) { + data[j * g.dims[0] + i] = NAN; + ((float *)g2.data)[j * g.dims[0] + i] = NAN; + } + } + } + + GridFree(&queue); + } + GridFree(&bc); + } + } + return g2; +} + +Grid GridLabelSinks(Grid g) { + assert(g.dtype == GridF32); + + Grid res = GridZerosLike(g); + + if (!res.data || !g.data) return res; + + float *r = res.data; + float *data = g.data; + + ptrdiff_t j_offset[8] = {-1, -1, -1, 0, 1, 1, 1, 0}; + ptrdiff_t i_offset[8] = {1, 0, -1, -1, -1, 0, 1, 1}; + + for (ptrdiff_t j = 1; j < g.dims[1] - 1; j++) { + for (ptrdiff_t i = 1; i < g.dims[0] - 1; i++) { + float z = data[i + g.dims[0] * j]; + ptrdiff_t up_neighbor_count = 0; + for (int32_t neighbor = 0; neighbor < 8; neighbor++) { + ptrdiff_t neighbor_i = i + i_offset[neighbor]; + ptrdiff_t neighbor_j = j + j_offset[neighbor]; + + if (data[g.dims[0] * neighbor_j + neighbor_i] > z) { + up_neighbor_count++; + } + } + r[j * g.dims[0] + i] = up_neighbor_count == 8 ? 1.0 : 0.0; + } + } + + return res; +} + +Grid GridIdentifyFlats(Grid g) { + assert(g.dtype == GridF32); + + Grid res = GridCreate(GridI32, NULL, g.cellsize, g.dims); + + if (!g.data || !res.data) return res; + + identifyflats(res.data, g.data, g.dims); + + return res; +} + +Grid GridGWDTComputeCosts(Grid flats, Grid dem, Grid demf) { + assert(flats.dtype == GridI32); + assert(dem.dtype == GridF32); + assert(dem.dtype == GridF32); + + Grid costs = GridZerosLike(dem); + if (!flats.data || !dem.data || !costs.data) return costs; + + Grid conncomps = GridCreate(GridIdx, NULL, dem.cellsize, dem.dims); + if (!conncomps.data) return costs; + + gwdt_computecosts(costs.data, conncomps.data, flats.data, dem.data, demf.data, + dem.dims); + + GridFree(&conncomps); + + return costs; +} + +Grid GridGWDT(Grid costs, Grid flats) { + assert(costs.dtype == GridF32); + assert(flats.dtype == GridI32); + + Grid dist = GridZerosLike(costs); + if (dist.data && costs.data && flats.data) { + Grid prev = GridCreate(GridIdx, NULL, dist.cellsize, dist.dims); + if (prev.data) { + Grid heap = GridCreate(GridIdx, NULL, dist.cellsize, dist.dims); + + if (heap.data) { + Grid back = GridCreate(GridIdx, NULL, dist.cellsize, dist.dims); + + if (back.data) { + gwdt(dist.data, prev.data, costs.data, flats.data, heap.data, + back.data, costs.dims); + + GridFree(&back); + } + GridFree(&heap); + } + GridFree(&prev); + } + } + + return dist; +} + +Grid GridGradient8(Grid dem, int use_mp) { + Grid res = GridZerosLike(dem); + if (res.data) { + gradient8(res.data, dem.data, dem.cellsize, use_mp, dem.dims); + } + return res; +} diff --git a/test/grid.h b/test/grid.h new file mode 100644 index 00000000..a2bcaf48 --- /dev/null +++ b/test/grid.h @@ -0,0 +1,64 @@ +#ifndef TT3_TEST_GRID_H +#define TT3_TEST_GRID_H + +#include +#include + +typedef enum GridDataType { + GridF32, + GridU8, + GridU32, + GridI32, + GridIdx +} GridDataType; + +typedef struct Grid { + void *data; + ptrdiff_t dims[2]; + float cellsize; + GridDataType dtype; +} Grid; + +void GridFree(Grid *g); + +// Creation +Grid GridCreate(GridDataType, void *buffer, float cellsize, ptrdiff_t dims[2]); +Grid GridFromFile(const char *path); +Grid GridRandom(uint64_t seed, float cellsize, ptrdiff_t dims[2]); +Grid GridCopy(Grid g); +Grid GridZeros(float cellsize, ptrdiff_t dims[2]); +Grid GridZerosLike(Grid g); +Grid GridFill(float value, float cellsize, ptrdiff_t dims[2]); +Grid GridTranspose(Grid g); + +// Output +int GridToFile(Grid g, const char *dst); + +// Details +size_t GridElementSize(Grid g); +ptrdiff_t GridElementCount(Grid g); +size_t GridSize(Grid g); + +// Arithmetic +Grid GridSubtract(Grid g1, Grid g2); +int GridAll(Grid g); +int GridAny(Grid g); +int GridAllNonnegative(Grid g); +int GridEq(Grid g1, Grid g2); +int GridAllClose(Grid g1, Grid g2, double rtol); +Grid GridInvert(Grid g); +Grid GridNot(Grid g); +Grid GridIsNan(Grid g); + +float GridMin(Grid g); +float GridMax(Grid g); + +// TopoToolbox +Grid GridFillsinks(Grid g); +Grid GridLabelSinks(Grid g); +Grid GridIdentifyFlats(Grid g); +Grid GridGWDTComputeCosts(Grid flats, Grid dem, Grid demf); +Grid GridGWDT(Grid costs, Grid flats); +Grid GridGradient8(Grid dem, int use_mp); + +#endif diff --git a/test/random_dem.cpp b/test/random_dem.cpp index 568c502c..af278e5c 100644 --- a/test/random_dem.cpp +++ b/test/random_dem.cpp @@ -10,8 +10,8 @@ #include #include -// Include topotoolbox.h in its own namespace to help prevent naming -// conflicts in the global scope. +// Include topotoolbox in its own namespace to help prevent +// naming conflicts in the global scope. namespace tt { extern "C" { #include "topotoolbox.h" @@ -19,6 +19,8 @@ extern "C" { } // namespace tt extern "C" { +#include "flow.h" +#include "grid.h" #include "utils.h" } @@ -32,13 +34,10 @@ Profiler prof; // Global Profiler that our test functions can access Each pixel of the filled DEM should be greater than or equal to the corresponding pixel in the original DEM. */ -int32_t test_fillsinks_ge(float *original_dem, float *filled_dem, - ptrdiff_t dims[2]) { - for (ptrdiff_t j = 0; j < dims[1]; j++) { - for (ptrdiff_t i = 0; i < dims[0]; i++) { - assert(filled_dem[i + dims[0] * j] >= original_dem[i + dims[0] * j]); - } - } +int32_t test_fillsinks_ge(Grid original, Grid filled) { + Grid sub = GridSubtract(filled, original); + assert(GridAllNonnegative(sub)); + GridFree(&sub); return 0; } @@ -46,30 +45,10 @@ int32_t test_fillsinks_ge(float *original_dem, float *filled_dem, No pixel in the filled DEM should be completely surrounded by pixels higher than it. */ -int32_t test_fillsinks_filled(float *filled_dem, ptrdiff_t dims[2]) { - ptrdiff_t j_offset[8] = {-1, -1, -1, 0, 1, 1, 1, 0}; - ptrdiff_t i_offset[8] = {1, 0, -1, -1, -1, 0, 1, 1}; - - for (ptrdiff_t j = 0; j < dims[1]; j++) { - for (ptrdiff_t i = 0; i < dims[0]; i++) { - float z = filled_dem[i + dims[0] * j]; - ptrdiff_t up_neighbor_count = 0; - for (int32_t neighbor = 0; neighbor < 8; neighbor++) { - ptrdiff_t neighbor_i = i + i_offset[neighbor]; - ptrdiff_t neighbor_j = j + j_offset[neighbor]; - - if (neighbor_i < 0 || neighbor_i >= dims[0] || neighbor_j < 0 || - neighbor_j >= dims[1]) { - continue; - } - - if (filled_dem[neighbor_i + dims[0] * neighbor_j] > z) { - up_neighbor_count++; - } - } - assert(up_neighbor_count != 8); - } - } +int32_t test_fillsinks_filled(Grid filled) { + Grid sinks = GridLabelSinks(filled); + assert(GridAny(sinks) == 0); + GridFree(&sinks); return 0; } @@ -78,15 +57,19 @@ int32_t test_fillsinks_filled(float *filled_dem, ptrdiff_t dims[2]) { 8 higher neighbors should be labeled a flat. Likewise, every flat should have no lower neighbors and fewer than 8 higher neighbors. */ -int32_t test_identifyflats_flats(int32_t *flats, float *dem, - ptrdiff_t dims[2]) { +int32_t test_identifyflats_flats(Grid flats, Grid dem) { ptrdiff_t col_offset[8] = {-1, -1, -1, 0, 1, 1, 1, 0}; ptrdiff_t row_offset[8] = {1, 0, -1, -1, -1, 0, 1, 1}; + int32_t *flat_data = (int32_t *)flats.data; + float *dem_data = (float *)dem.data; + + ptrdiff_t dims[2] = {dem.dims[0], dem.dims[1]}; + for (ptrdiff_t j = 0; j < dims[1]; j++) { for (ptrdiff_t i = 0; i < dims[0]; i++) { - float z = dem[i + dims[0] * j]; - int32_t flat = flats[i + dims[0] * j]; + float z = dem_data[i + dims[0] * j]; + int32_t flat = flat_data[i + dims[0] * j]; int32_t up_neighbor_count = 0; int32_t down_neighbor_count = 0; @@ -99,7 +82,7 @@ int32_t test_identifyflats_flats(int32_t *flats, float *dem, continue; } - float neighbor_height = dem[neighbor_i + dims[0] * neighbor_j]; + float neighbor_height = dem_data[neighbor_i + dims[0] * neighbor_j]; if (neighbor_height > z) { up_neighbor_count++; @@ -122,15 +105,19 @@ int32_t test_identifyflats_flats(int32_t *flats, float *dem, Every pixel that has a lower neighbor, borders a flat, and has the same elevation as a flat that it touches should be labeled a sill. */ -int32_t test_identifyflats_sills(int32_t *flats, float *dem, - ptrdiff_t dims[2]) { +int32_t test_identifyflats_sills(Grid flats, Grid dem) { ptrdiff_t j_offset[8] = {-1, -1, -1, 0, 1, 1, 1, 0}; ptrdiff_t i_offset[8] = {1, 0, -1, -1, -1, 0, 1, 1}; + int32_t *flat_data = (int32_t *)flats.data; + float *dem_data = (float *)dem.data; + + ptrdiff_t dims[2] = {dem.dims[0], dem.dims[1]}; + for (ptrdiff_t j = 0; j < dims[1]; j++) { for (ptrdiff_t i = 0; i < dims[0]; i++) { - float z = dem[i + j * dims[0]]; - int32_t flat = flats[i + j * dims[0]]; + float z = dem_data[i + j * dims[0]]; + int32_t flat = flat_data[i + j * dims[0]]; int32_t down_neighbor_count = 0; int32_t equal_neighbor_flats = 0; @@ -147,8 +134,8 @@ int32_t test_identifyflats_sills(int32_t *flats, float *dem, continue; } - float neighbor_height = dem[neighbor_i + dims[0] * neighbor_j]; - int32_t neighbor_flat = flats[neighbor_i + dims[0] * neighbor_j]; + float neighbor_height = dem_data[neighbor_i + dims[0] * neighbor_j]; + int32_t neighbor_flat = flat_data[neighbor_i + dims[0] * neighbor_j]; if (neighbor_height < z) { down_neighbor_count++; @@ -170,15 +157,19 @@ int32_t test_identifyflats_sills(int32_t *flats, float *dem, Every pixel that is a flat and borders a sill of the same height should be labeled a presill */ -int32_t test_identifyflats_presills(int32_t *flats, float *dem, - ptrdiff_t dims[2]) { +int32_t test_identifyflats_presills(Grid flats, Grid dem) { ptrdiff_t j_offset[8] = {-1, -1, -1, 0, 1, 1, 1, 0}; ptrdiff_t i_offset[8] = {1, 0, -1, -1, -1, 0, 1, 1}; + int32_t *flat_data = (int32_t *)flats.data; + float *dem_data = (float *)dem.data; + + ptrdiff_t dims[2] = {dem.dims[0], dem.dims[1]}; + for (ptrdiff_t j = 0; j < dims[1]; j++) { for (ptrdiff_t i = 0; i < dims[0]; i++) { - float z = dem[i + j * dims[0]]; - int32_t flat = flats[i + j * dims[0]]; + float z = dem_data[i + j * dims[0]]; + int32_t flat = flat_data[i + j * dims[0]]; int32_t equal_neighbor_sills = 0; @@ -191,8 +182,8 @@ int32_t test_identifyflats_presills(int32_t *flats, float *dem, continue; } - float neighbor_height = dem[neighbor_i + dims[0] * neighbor_j]; - int32_t neighbor_flat = flats[neighbor_i + dims[0] * neighbor_j]; + float neighbor_height = dem_data[neighbor_i + dims[0] * neighbor_j]; + int32_t neighbor_flat = flat_data[neighbor_i + dims[0] * neighbor_j]; if (((neighbor_flat & 2) > 0) && (neighbor_height == z)) { equal_neighbor_sills++; @@ -209,11 +200,11 @@ int32_t test_identifyflats_presills(int32_t *flats, float *dem, /* Costs should be zero on nonflats and positive on flats. */ -int32_t test_gwdt_costs(float *costs, int32_t *flats, ptrdiff_t dims[2]) { - for (ptrdiff_t j = 0; j < dims[1]; j++) { - for (ptrdiff_t i = 0; i < dims[0]; i++) { - float cost = costs[i + dims[0] * j]; - int32_t flat = flats[i + dims[0] * j]; +int32_t test_gwdt_costs(Grid costs, Grid flats) { + for (ptrdiff_t j = 0; j < flats.dims[1]; j++) { + for (ptrdiff_t i = 0; i < flats.dims[0]; i++) { + float cost = ((float *)costs.data)[i + costs.dims[0] * j]; + int32_t flat = ((int32_t *)flats.data)[i + flats.dims[0] * j]; assert((cost >= 0) && (((flat & 1) > 0) == (cost > 0))); } @@ -226,15 +217,15 @@ int32_t test_gwdt_costs(float *costs, int32_t *flats, ptrdiff_t dims[2]) { Flats neighboring other flats should all have the same nonzero connected component label */ -int32_t test_gwdt_conncomps(ptrdiff_t *conncomps, int32_t *flats, - ptrdiff_t dims[2]) { +int32_t test_gwdt_conncomps(Grid conncomps, Grid flats) { ptrdiff_t j_offset[8] = {-1, -1, -1, 0, 1, 1, 1, 0}; ptrdiff_t i_offset[8] = {1, 0, -1, -1, -1, 0, 1, 1}; - for (ptrdiff_t j = 0; j < dims[1]; j++) { - for (ptrdiff_t i = 0; i < dims[0]; i++) { - ptrdiff_t label = conncomps[i + dims[0] * j]; - int32_t flat = flats[i + dims[0] * j]; + for (ptrdiff_t j = 0; j < flats.dims[1]; j++) { + for (ptrdiff_t i = 0; i < flats.dims[0]; i++) { + ptrdiff_t label = + ((ptrdiff_t *)conncomps.data)[i + conncomps.dims[0] * j]; + int32_t flat = ((int32_t *)flats.data)[i + flats.dims[0] * j]; if (!(flat & 1)) { continue; @@ -244,13 +235,16 @@ int32_t test_gwdt_conncomps(ptrdiff_t *conncomps, int32_t *flats, ptrdiff_t neighbor_i = i + i_offset[neighbor]; ptrdiff_t neighbor_j = j + j_offset[neighbor]; - if (neighbor_i < 0 || neighbor_i >= dims[0] || neighbor_j < 0 || - neighbor_j >= dims[1]) { + if (neighbor_i < 0 || neighbor_i >= flats.dims[0] || neighbor_j < 0 || + neighbor_j >= flats.dims[1]) { continue; } - int32_t neighbor_flat = flats[neighbor_i + dims[0] * neighbor_j]; - ptrdiff_t neighbor_label = conncomps[neighbor_i + dims[0] * neighbor_j]; + int32_t neighbor_flat = + ((int32_t *)flats.data)[neighbor_i + flats.dims[0] * neighbor_j]; + ptrdiff_t neighbor_label = + ((ptrdiff_t *) + conncomps.data)[neighbor_i + conncomps.dims[0] * neighbor_j]; if (neighbor_flat & 1) { assert(neighbor_label == label); } @@ -269,12 +263,13 @@ int32_t test_gwdt_conncomps(ptrdiff_t *conncomps, int32_t *flats, 4. equal to the distance to the parent pixel plus the geodesic distance between the parent and the current pixel. */ -int32_t test_gwdt(float *dist, ptrdiff_t *prev, float *costs, int32_t *flats, - ptrdiff_t dims[2]) { +int32_t test_gwdt(Grid dist, Grid prev, Grid costs, Grid flats) { + ptrdiff_t dims[2] = {dist.dims[0], dist.dims[1]}; + for (ptrdiff_t j = 0; j < dims[1]; j++) { for (ptrdiff_t i = 0; i < dims[0]; i++) { - float d = dist[i + j * dims[0]]; - int32_t flat = flats[i + j * dims[0]]; + float d = ((float *)dist.data)[j * dims[0] + i]; + int32_t flat = ((int32_t *)flats.data)[j * dims[0] + i]; assert(d >= 0); @@ -283,13 +278,15 @@ int32_t test_gwdt(float *dist, ptrdiff_t *prev, float *costs, int32_t *flats, } else if (flat & 1) { assert(d > 1.0); - ptrdiff_t parent = prev[i + j * dims[0]]; + ptrdiff_t parent = ((ptrdiff_t *)prev.data)[j * dims[0] + i]; ptrdiff_t parent_j = parent / dims[0]; ptrdiff_t parent_i = parent % dims[0]; float chamfer = (parent_j != j) && (parent_i != i) ? SQRT2f : 1.0f; - float proposed_dist = - dist[parent] + - chamfer * (costs[i + j * dims[0]] + costs[parent]) / 2; + float proposed_dist = ((float *)dist.data)[parent] + + chamfer * + (((float *)costs.data)[j * dims[0] + i] + + ((float *)costs.data)[parent]) / + 2; assert(proposed_dist == d); } } @@ -301,13 +298,18 @@ int32_t test_gwdt(float *dist, ptrdiff_t *prev, float *costs, int32_t *flats, For each cell, the saved gradient should be greater than or equal to the signed gradient to all 8 neighboring cells. */ -int32_t test_gradient8(float *gradient, float *dem, float cellsize, - ptrdiff_t dims[2]) { +int32_t test_gradient8(Grid gradient, Grid dem) { ptrdiff_t i_offset[8] = {0, 1, 1, 1, 0, -1, -1, -1}; ptrdiff_t j_offset[8] = {1, 1, 0, -1, -1, -1, 0, 1}; + + ptrdiff_t dims[2] = {dem.dims[0], dem.dims[1]}; + + float *gradient_data = (float *)gradient.data; + float *dem_data = (float *)dem.data; + for (ptrdiff_t j = 0; j < dims[1]; j++) { for (ptrdiff_t i = 0; i < dims[0]; i++) { - float max_gradient = gradient[j * dims[0] + i]; + float max_gradient = gradient_data[j * dims[0] + i]; // Iterate over all 8 neighboring cells for (int k = 0; k < 8; k++) { @@ -321,11 +323,11 @@ int32_t test_gradient8(float *gradient, float *dem, float cellsize, float vertical_dist; float local_gradient; - horizontal_dist = cellsize; + horizontal_dist = dem.cellsize; horizontal_dist *= neighbor_i != i && neighbor_j != j ? SQRT2 : 1.0f; - vertical_dist = - dem[j * dims[0] + i] - dem[neighbor_j * dims[0] + neighbor_i]; + vertical_dist = dem_data[j * dims[0] + i] - + dem_data[neighbor_j * dims[0] + neighbor_i]; local_gradient = vertical_dist / horizontal_dist; assert(max_gradient >= local_gradient); @@ -338,38 +340,35 @@ int32_t test_gradient8(float *gradient, float *dem, float cellsize, /* Computing the gradient8 using OpenMP should yield the same result as without. */ -int32_t test_gradient8_mp(float *gradient, float *gradient_mp, - ptrdiff_t dims[2]) { - for (ptrdiff_t j = 0; j < dims[1]; j++) { - for (ptrdiff_t i = 0; i < dims[0]; i++) { - ptrdiff_t position = j * dims[0] + i; - assert(gradient[position] == gradient_mp[position]); - } - } +int32_t test_gradient8_mp(Grid gradient, Grid gradient_mp) { + assert(GridEq(gradient, gradient_mp)); return 0; } /* Flow direction should point downstream or across flats */ -int32_t test_routeflowd8_direction(uint8_t *direction, float *filled_dem, - float *dist, int32_t *flats, - ptrdiff_t dims[2]) { +int32_t test_routeflowd8_direction(Grid direction, Grid filled_dem) { // 1<<5 1<<6 1<<7 // 1<<4 0 1<<0 // 1<<3 1<<2 1<<1 ptrdiff_t i_offset[8] = {0, 1, 1, 1, 0, -1, -1, -1}; ptrdiff_t j_offset[8] = {1, 1, 0, -1, -1, -1, 0, 1}; + ptrdiff_t dims[2] = {direction.dims[0], direction.dims[1]}; + + uint8_t *direction_data = (uint8_t *)direction.data; + float *filled_data = (float *)filled_dem.data; + for (ptrdiff_t j = 0; j < dims[1]; j++) { for (ptrdiff_t i = 0; i < dims[0]; i++) { - uint8_t flowdir = direction[j * dims[0] + i]; + uint8_t flowdir = direction_data[j * dims[0] + i]; if (flowdir == 0) { // Pixel is a sink continue; } - float z = filled_dem[j * dims[0] + i]; + float z = filled_data[j * dims[0] + i]; uint8_t v = flowdir; uint8_t r = 0; while (v >>= 1) { @@ -387,71 +386,41 @@ int32_t test_routeflowd8_direction(uint8_t *direction, float *filled_dem, // Neighbor elevation should be less than or equal to the // current elevation - assert(filled_dem[neighbor_j * dims[0] + neighbor_i] <= z); + assert(filled_data[neighbor_j * dims[0] + neighbor_i] <= z); } } return 0; } -/* - source should be a topological sort of the - graph defined by direction. - */ -int32_t test_routeflowd8_tsort(uint8_t *marks, ptrdiff_t *source, - uint8_t *direction, ptrdiff_t dims[2]) { - // 1<<5 1<<6 1<<7 - // 1<<4 0 1<<0 - // 1<<3 1<<2 1<<1 - ptrdiff_t i_offset[8] = {0, 1, 1, 1, 0, -1, -1, -1}; - ptrdiff_t j_offset[8] = {1, 1, 0, -1, -1, -1, 0, 1}; +int32_t test_tsort(Flow fd) { + // source and target arrays should be topologically sorted. All + // incoming edges to a node must occur before any outgoing edges. - for (ptrdiff_t j = 0; j < dims[1]; j++) { - for (ptrdiff_t i = 0; i < dims[0]; i++) { - ptrdiff_t src = source[j * dims[0] + i]; - ptrdiff_t src_i = src % dims[0]; - ptrdiff_t src_j = src / dims[0]; + int32_t flag = 0; - uint8_t flowdir = direction[src]; + Grid outdegree = GridCreate(GridU8, NULL, fd.cellsize, fd.dims); + for (ptrdiff_t e = 0; e < fd.count; e++) { + ptrdiff_t u = fd.source[e]; + ptrdiff_t v = fd.target[e]; - if (flowdir == 0) { - // src is a sink - continue; - } + ((uint8_t *)outdegree.data)[u] += 1; - uint8_t v = flowdir; - uint8_t r = 0; - while (v >>= 1) { - r++; - } - ptrdiff_t neighbor_i = src_i + i_offset[r]; - ptrdiff_t neighbor_j = src_j + j_offset[r]; - - ptrdiff_t neighbor_src = neighbor_j * dims[0] + neighbor_i; - // Check that we haven't seen neighbor in the topological sort yet - assert(marks[neighbor_src] != 0xff); - - marks[src] = 0xff; + if (((uint8_t *)outdegree.data)[v] > 0) { + flag = 1; } } - return 0; + assert(flag == 0); + GridFree(&outdegree); + return flag; } -int32_t test_flow_accumulation_max(float *acc, ptrdiff_t dims[2]) { - for (ptrdiff_t j = 0; j < dims[1]; j++) { - for (ptrdiff_t i = 0; i < dims[0]; i++) { - assert(acc[j * dims[0] + i] <= (dims[0] * dims[1])); - } - } - return 0; -} +int32_t test_flow_accumulation_max(Grid acc) { + Grid max = GridFill(acc.dims[0] * acc.dims[1], acc.cellsize, acc.dims); + Grid sub = GridSubtract(max, acc); + assert(GridAllNonnegative(sub)); -int32_t test_flow_accumulation_multimethod(float *acc1, float *acc2, - ptrdiff_t dims[2]) { - for (ptrdiff_t j = 0; j < dims[1]; j++) { - for (ptrdiff_t i = 0; i < dims[0]; i++) { - assert(acc1[j * dims[0] + i] == acc2[j * dims[0] + i]); - } - } + GridFree(&sub); + GridFree(&max); return 0; } @@ -529,7 +498,7 @@ int32_t test_stream_distance(float *node_distance, ptrdiff_t *stream_grid, A source pixel and its downstream neighbor should be part of the same basin, which should have an index greater than 0. */ -int32_t test_drainagebasins(ptrdiff_t *basins, ptrdiff_t *source, +int32_t test_drainagebasins(uint32_t *basins, ptrdiff_t *source, ptrdiff_t *target, ptrdiff_t edge_count, ptrdiff_t dims[2]) { for (ptrdiff_t e = 0; e < edge_count; e++) { @@ -579,9 +548,8 @@ int32_t test_lowerenv_convex(float *g, float *z, float *d, uint8_t *knickpoints, // drainagebasins should return the same drainage basins computed // using an appropriately initialized traverse_up_u32_or_and. -int32_t test_db_traverse(ptrdiff_t *basins, ptrdiff_t *source, - ptrdiff_t *target, ptrdiff_t edge_count, - ptrdiff_t dims[2]) { +int32_t test_db_traverse(uint32_t *basins, ptrdiff_t *source, ptrdiff_t *target, + ptrdiff_t edge_count, ptrdiff_t dims[2]) { std::vector ones(edge_count, 0xffffffff); // Edge weights std::vector traverse_basins(dims[0] * dims[1], 0); // Basin labels std::vector indegree(dims[0] * dims[1], 0); @@ -616,45 +584,38 @@ struct FlowRoutingData { float cellsize; // input data - std::vector dem; - std::vector bc; - std::vector fraction; + Grid dem; + Grid bc; + Grid fraction; // fillsinks - std::vector filled_dem; - // fillsinks_hybrid - std::vector queue; + Grid filled_dem; // identifyflats - std::vector flats; + Grid flats; // gwdt_compute_costs - std::vector conncomps; - std::vector costs; + Grid costs; + Grid conncomps; // gwdt - std::vector dist; - std::vector heap; - std::vector back; - std::vector prev; + Grid dist; + Grid heap; + Grid back; + Grid prev; // gradient8 - std::vector gradient; - std::vector gradient_mp; + Grid gradient; + Grid gradient_mp; // flow_routing_d8_carve - std::vector nodes; - std::vector direction; - std::vector marks; + Grid direction; // flow_routing_targets - std::vector source; - std::vector target; - ptrdiff_t edge_count; + Flow fd; // flow_accumulation - std::vector accum; - std::vector accum2; + Grid accum; // stream network std::vector stream_source; @@ -662,112 +623,114 @@ struct FlowRoutingData { std::vector stream_weight; std::vector stream_grid; ptrdiff_t stream_node_count; - std::vector basins; + + Grid basins; FlowRoutingData(ptrdiff_t input_dims[2], float cs, uint32_t seed) : dims({input_dims[0], input_dims[1]}), cellsize(cs), - dem(dims[0] * dims[1]), - bc(dims[0] * dims[1]), - fraction(dims[0] * dims[1]), - filled_dem(dims[0] * dims[1]), - queue(dims[0] * dims[1]), - flats(dims[0] * dims[1]), - conncomps(dims[0] * dims[1]), - costs(dims[0] * dims[1]), - dist(dims[0] * dims[1]), - heap(dims[0] * dims[1]), - back(dims[0] * dims[1]), - prev(dims[0] * dims[1]), - gradient(dims[0] * dims[1]), - gradient_mp(dims[0] * dims[1]), - nodes(dims[0] * dims[1]), - direction(dims[0] * dims[1]), - marks(dims[0] * dims[1]), - source(dims[0] * dims[1]), - target(dims[0] * dims[1]), - accum(dims[0] * dims[1]), - accum2(dims[0] * dims[1]), - stream_grid(dims[0] * dims[1]), - basins(dims[0] * dims[1]) { + stream_grid(dims[0] * dims[1]) { // Initialize DEM, boundary conditions for fillsinks and fraction - for (uint32_t col = 0; col < dims[1]; col++) { - for (uint32_t row = 0; row < dims[0]; row++) { - dem[col * dims[0] + row] = 100.0f * pcg4d(row, col, seed, 1); + dem = GridRandom(seed, cellsize, dims.data()); + bc = GridCreate(GridU8, NULL, cellsize, dims.data()); + fraction = GridFill(1.0f, cellsize, dims.data()); + for (ptrdiff_t col = 0; col < dims[1]; col++) { + for (ptrdiff_t row = 0; row < dims[0]; row++) { if (row == 0 || row == dims[0] - 1 || col == 0 || col == dims[1] - 1) { - bc[col * dims[0] + row] = 1; + ((uint8_t *)bc.data)[col * dims[0] + row] = 1; } else { - bc[col * dims[0] + row] = 0; + ((uint8_t *)bc.data)[col * dims[0] + row] = 0; } - - fraction[col * dims[0] + row] = 1.0f; } } } - void route_flow() { + ~FlowRoutingData() { + GridFree(&dem); + GridFree(&bc); + GridFree(&fraction); + GridFree(&filled_dem); + GridFree(&flats); + GridFree(&costs); + GridFree(&conncomps); + GridFree(&dist); + GridFree(&prev); + GridFree(&back); + GridFree(&heap); + GridFree(&gradient); + GridFree(&gradient_mp); + GridFree(&direction); + GridFree(&accum); + GridFree(&basins); + FlowFree(&fd); + } + + void fillsinks_hybrid() { ProfileFunction(prof); + filled_dem = GridFillsinks(dem); + } - tt::fillsinks(filled_dem.data(), dem.data(), bc.data(), dims.data()); - - tt::identifyflats(flats.data(), filled_dem.data(), dims.data()); - - tt::gwdt_computecosts(costs.data(), conncomps.data(), flats.data(), - dem.data(), filled_dem.data(), dims.data()); - - tt::gwdt(dist.data(), prev.data(), costs.data(), flats.data(), heap.data(), - back.data(), dims.data()); + void fillsinks() { + ProfileFunction(prof); + filled_dem = GridCreate(GridF32, NULL, dem.cellsize, dem.dims); - tt::flow_routing_d8_carve(nodes.data(), direction.data(), filled_dem.data(), - dist.data(), flats.data(), dims.data(), 0); + tt::fillsinks((float *)filled_dem.data, (float *)dem.data, + (uint8_t *)bc.data, dims.data()); + } - edge_count = - tt::flow_routing_d8_edgelist(source.data(), target.data(), nodes.data(), - direction.data(), dims.data(), 0); + void route_flow(bool hybrid) { + ProfileFunction(prof); - tt::flow_accumulation_edgelist(accum2.data(), source.data(), target.data(), - fraction.data(), NULL, edge_count, - dims.data()); + if (hybrid) { + fillsinks_hybrid(); + } else { + fillsinks(); + } - // Close timer - } + flats = GridIdentifyFlats(filled_dem); - void route_flow_hybrid() { - ProfileFunction(prof); + costs = GridZerosLike(dem); + conncomps = GridCreate(GridIdx, NULL, dem.cellsize, dem.dims); - tt::fillsinks_hybrid(filled_dem.data(), queue.data(), dem.data(), bc.data(), - dims.data()); + tt::gwdt_computecosts((float *)costs.data, (ptrdiff_t *)conncomps.data, + (int32_t *)flats.data, (float *)dem.data, + (float *)filled_dem.data, dims.data()); - tt::identifyflats(flats.data(), filled_dem.data(), dims.data()); + dist = GridZerosLike(costs); + prev = GridCreate(GridIdx, NULL, dist.cellsize, dist.dims); + heap = GridCreate(GridIdx, NULL, dist.cellsize, dist.dims); + back = GridCreate(GridIdx, NULL, dist.cellsize, dist.dims); + tt::gwdt((float *)dist.data, (ptrdiff_t *)prev.data, (float *)costs.data, + (int32_t *)flats.data, (ptrdiff_t *)heap.data, + (ptrdiff_t *)back.data, dims.data()); - tt::gwdt_computecosts(costs.data(), conncomps.data(), flats.data(), - dem.data(), filled_dem.data(), dims.data()); + fd = FlowInit(dem); - tt::gwdt(dist.data(), prev.data(), costs.data(), flats.data(), heap.data(), - back.data(), dims.data()); + direction = GridCreate(GridU8, NULL, dem.cellsize, dem.dims); + tt::flow_routing_d8_carve((ptrdiff_t *)fd.stream, (uint8_t *)direction.data, + (float *)filled_dem.data, (float *)dist.data, + (int32_t *)flats.data, dims.data(), 0); - tt::flow_routing_d8_carve(nodes.data(), direction.data(), filled_dem.data(), - dist.data(), flats.data(), dims.data(), 0); + fd.count = tt::flow_routing_d8_edgelist( + (ptrdiff_t *)fd.source, (ptrdiff_t *)fd.target, (ptrdiff_t *)fd.stream, + (uint8_t *)direction.data, dims.data(), 0); - edge_count = - tt::flow_routing_d8_edgelist(source.data(), target.data(), nodes.data(), - direction.data(), dims.data(), 0); + for (ptrdiff_t e = 0; e < fd.count; e++) { + fd.fraction[e] = 1.0; + } - tt::flow_accumulation_edgelist(accum2.data(), source.data(), target.data(), - fraction.data(), NULL, edge_count, - dims.data()); + accum = FlowAccumulation(fd); } void gradient8() { ProfileFunction(prof); - tt::gradient8(gradient.data(), dem.data(), cellsize, 0, dims.data()); + gradient = GridGradient8(dem, 0); } void gradient8_mp() { ProfileFunction(prof); - - tt::gradient8(gradient_mp.data(), dem.data(), cellsize, 1, dims.data()); + gradient_mp = GridGradient8(dem, 1); } void streamnetwork(float threshold) { @@ -779,8 +742,8 @@ struct FlowRoutingData { // the node attribute list. for (ptrdiff_t j = 0; j < dims[1]; j++) { for (ptrdiff_t i = 0; i < dims[0]; i++) { - ptrdiff_t u = nodes[j * dims[0] + i]; - if (accum2[u] >= threshold) { + ptrdiff_t u = ((ptrdiff_t *)fd.stream)[j * dims[0] + i]; + if (((float *)accum.data)[u] >= threshold) { // Pixel u is in the stream network stream_grid[u] = node_index++; } else { @@ -790,9 +753,9 @@ struct FlowRoutingData { } // Find edges that are part of the stream network - for (ptrdiff_t e = 0; e < edge_count; e++) { - ptrdiff_t u = source[e]; - ptrdiff_t v = target[e]; + for (ptrdiff_t e = 0; e < fd.count; e++) { + ptrdiff_t u = fd.source[e]; + ptrdiff_t v = fd.target[e]; // Compute distance between source and target pixels. If the // difference between u and v is 1 or dims[0], then they are @@ -801,7 +764,7 @@ struct FlowRoutingData { ? 1.0f * cellsize : SQRT2f * cellsize; - if (accum2[u] >= threshold) { + if (((float *)accum.data)[u] >= threshold) { // Pixel u is in the stream network // Add this edge to the stream network, but map u and v to @@ -817,58 +780,46 @@ struct FlowRoutingData { void drainagebasins() { ProfileFunction(prof); - tt::drainagebasins(basins.data(), source.data(), target.data(), edge_count, - dims.data()); + basins = FlowDrainageBasins(fd); } void runtests(bool hybrid) { - if (hybrid) { - route_flow_hybrid(); - } else { - route_flow(); - } + route_flow(hybrid); + + test_fillsinks_ge(dem, filled_dem); + test_fillsinks_filled(filled_dem); - test_fillsinks_ge(dem.data(), filled_dem.data(), dims.data()); - test_fillsinks_filled(filled_dem.data(), dims.data()); + test_identifyflats_flats(flats, filled_dem); + test_identifyflats_sills(flats, filled_dem); + test_identifyflats_presills(flats, filled_dem); - test_identifyflats_flats(flats.data(), filled_dem.data(), dims.data()); - test_identifyflats_sills(flats.data(), filled_dem.data(), dims.data()); - test_identifyflats_presills(flats.data(), filled_dem.data(), dims.data()); + test_gwdt_costs(costs, flats); + test_gwdt_conncomps(conncomps, flats); - test_gwdt_costs(costs.data(), flats.data(), dims.data()); - test_gwdt_conncomps(conncomps.data(), flats.data(), dims.data()); + test_gwdt(dist, prev, costs, flats); // route_flow and route_flow_hybrid do not compute gradients gradient8(); gradient8_mp(); - test_gradient8(gradient.data(), dem.data(), cellsize, dims.data()); - test_gradient8_mp(gradient.data(), gradient_mp.data(), dims.data()); + test_gradient8(gradient, dem); + test_gradient8_mp(gradient, gradient_mp); - test_routeflowd8_direction(direction.data(), filled_dem.data(), dist.data(), - flats.data(), dims.data()); - test_routeflowd8_tsort(marks.data(), nodes.data(), direction.data(), - dims.data()); + test_routeflowd8_direction(direction, filled_dem); - test_flow_routing_targets(target.data(), source.data(), direction.data(), - edge_count, dims.data()); + test_tsort(fd); + + test_flow_routing_targets((ptrdiff_t *)fd.target, (ptrdiff_t *)fd.source, + (uint8_t *)direction.data, fd.count, dims.data()); // Compute and test drainage basins drainagebasins(); - test_drainagebasins(basins.data(), source.data(), target.data(), edge_count, - dims.data()); - test_db_traverse(basins.data(), source.data(), target.data(), edge_count, - dims.data()); - - // route_flow and route_flow_hybrid only run the edgelist variant - // of flow_accumulation, so we must run it explicitly. - tt::flow_accumulation(accum.data(), nodes.data(), direction.data(), NULL, - dims.data()); - - test_flow_accumulation_max(accum.data(), dims.data()); + test_drainagebasins((uint32_t *)basins.data, (ptrdiff_t *)fd.source, + (ptrdiff_t *)fd.target, fd.count, dims.data()); + test_db_traverse((uint32_t *)basins.data, (ptrdiff_t *)fd.source, + (ptrdiff_t *)fd.target, fd.count, dims.data()); - test_flow_accumulation_multimethod(accum.data(), accum2.data(), - dims.data()); + test_flow_accumulation_max(accum); // Generate stream network streamnetwork(dims[0] * dims[1] / 20.0f); @@ -880,8 +831,8 @@ struct FlowRoutingData { stream_source.data(), stream_target.data(), stream_weight.data(), stream_source.size()); test_stream_distance(integral.data(), stream_grid.data(), distance.data(), - source.data(), target.data(), cellsize, source.size(), - dims.data()); + (ptrdiff_t *)fd.source, (ptrdiff_t *)fd.target, + cellsize, fd.count, dims.data()); std::vector pvF32(stream_node_count, 0.0f); std::vector pvF64(stream_node_count, 0.0); @@ -935,7 +886,7 @@ struct FlowRoutingData { for (ptrdiff_t i = 0; i < dims[0]; i++) { ptrdiff_t node = j * dims[0] + i; if (stream_grid[node] >= 0) { - z[stream_grid[node]] = dem[node]; + z[stream_grid[node]] = ((float *)dem.data)[node]; d[stream_grid[node]] = distance[node]; } } From 59a39be780153dd3598bcabbc58ed5d91e3954f6 Mon Sep 17 00:00:00 2001 From: William Kearney Date: Thu, 2 Apr 2026 12:21:44 +0200 Subject: [PATCH 2/3] Only run the jq profiler check on Ubuntu The GDAL requirement for random_dem means that it currently is not run on macos or Windows. --- .github/workflows/ci.yaml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index 1abc5747..b0b08b0c 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -129,6 +129,7 @@ jobs: - name: Test static library run: ctest --test-dir ${{github.workspace}}/build/static -C Release --output-on-failure - name: Validate JSON output of profiler + if: matrix.os == 'ubuntu-latest' run: ${{github.workspace}}/build/static/test/random_dem | jq -e . - name: Configure shared library run: > From c78993f15f0f8baf36cc9f059f2256585e9b0224 Mon Sep 17 00:00:00 2001 From: William Kearney Date: Thu, 2 Apr 2026 12:25:02 +0200 Subject: [PATCH 3/3] Only attempt to format random_dem if it has been created --- test/CMakeLists.txt | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 16963d75..c32392c9 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -166,7 +166,6 @@ if (TT_FORMAT_CHECK AND CLANG_FORMAT_EXE) set(FORMAT_TARGETS topotoolbox versioninfo - random_dem excesstopography filters swaths @@ -176,6 +175,10 @@ if (TT_FORMAT_CHECK AND CLANG_FORMAT_EXE) list(APPEND FORMAT_TARGETS snapshot) endif() + if (TARGET random_dem) + list(APPEND FORMAT_TARGETS random_dem) + endif() + set(ALL_SOURCES) foreach(FORMAT_TARGET ${FORMAT_TARGETS})