diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index 1abc574..b0b08b0 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: > diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index f30fb42..c32392c 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) @@ -163,7 +166,6 @@ if (TT_FORMAT_CHECK AND CLANG_FORMAT_EXE) set(FORMAT_TARGETS topotoolbox versioninfo - random_dem excesstopography filters swaths @@ -173,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}) diff --git a/test/flow.c b/test/flow.c new file mode 100644 index 0000000..7ea50c3 --- /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 0000000..2d7219e --- /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 0000000..6ffb780 --- /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 0000000..a2bcaf4 --- /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 568c502..af278e5 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]; } }