From 3f8a4b36aef813f17be6f662f97b5e173985c4b6 Mon Sep 17 00:00:00 2001 From: William Kearney Date: Thu, 16 Apr 2026 16:44:31 +0200 Subject: [PATCH] Implement libtopotoolbox-based flow routing Resolves #22 Resolves #15 This PR finishes implementing a pure libtopotoolbox workflow for flow routing for D8 with the 'carve' preprocessing option. New bindings are implemented for libtopotoolbox's `gwdt`, `flow_routing_d8_carve` and `flow_routing_d8_edgelist`. The latter two functions are combined in a single binding, `tt_flow_routing_d8_carve` because there is no real reason to call them separately. The flow routing interface in libtopotoolbox is expected to change in the near future, so this should minimize the amount of changes needed in the MATLAB code. createAuxiliaryTopo and the FLOWobj constructor now call only libtopotoolbox functions if uselibtt=true. This required a few changes to pass around the correct flats array. libtopotoolbox's `gwdt` expects a 32-bit integer with an encoding of flats, sills and presills, so we need to use libtopotoolbox's identifyflats explicitly to generate this. `flow_routing_d8_carve` also needs the flats array, but it is okay if they are just labels for the flats, so that array is also passed out of createAuxiliaryTopo for use in the FLOWobj constructor. The resulting FLOWobj is not exactly identical to that without libtopotoolbox because of TopoToolbox/libtopotoolbox#122. It is, however, identical to the pytopotoolbox output. The libtopotoolbox implementation is slightly faster than the non-libtopotoolbox one. The snapshot tests are updated to compare the auxiliary topography with a tolerance because libtopotoolbox does not give exactly identical answers as MATLAB for the gray-weighted distance transform. Signed-off-by: William Kearney --- bindings/CMakeLists.txt | 19 +++++++- bindings/tt_flow_routing_d8_carve.c | 47 ++++++++++++++++++++ bindings/tt_gwdt.c | 60 ++++++++++++++++++++++++++ tests/testSnapshot.m | 2 +- toolbox/@FLOWobj/FLOWobj.m | 20 ++++++--- toolbox/internal/createAuxiliaryTopo.m | 40 +++++++++++------ 6 files changed, 168 insertions(+), 20 deletions(-) create mode 100644 bindings/tt_flow_routing_d8_carve.c create mode 100644 bindings/tt_gwdt.c diff --git a/bindings/CMakeLists.txt b/bindings/CMakeLists.txt index 39b785b..96a7dd1 100644 --- a/bindings/CMakeLists.txt +++ b/bindings/CMakeLists.txt @@ -48,6 +48,22 @@ matlab_add_mex( LINK_TO topotoolbox ) +matlab_add_mex( + NAME tt_gwdt + MODULE + SRC tt_gwdt.c + R2018a + LINK_TO topotoolbox +) + +matlab_add_mex( + NAME tt_flow_routing_d8_carve + MODULE + SRC tt_flow_routing_d8_carve.c + R2018a + LINK_TO topotoolbox +) + matlab_add_mex( NAME tt_lowerenv MODULE @@ -113,7 +129,8 @@ matlab_add_mex( ) install( - TARGETS tt_has_topotoolbox tt_fillsinks tt_identifyflats tt_gwdt_computecosts + TARGETS tt_has_topotoolbox tt_fillsinks tt_identifyflats + tt_gwdt_computecosts tt_gwdt tt_flow_routing_d8_carve tt_lowerenv tt_hillshade tt_hillshade_fused tt_graphflood tt_gradient8 tt_excesstopography_fsm2d diff --git a/bindings/tt_flow_routing_d8_carve.c b/bindings/tt_flow_routing_d8_carve.c new file mode 100644 index 0000000..01b52fe --- /dev/null +++ b/bindings/tt_flow_routing_d8_carve.c @@ -0,0 +1,47 @@ +/* + + [SOURCE, TARGET, COUNT] = tt_flow_routing_d8_carve(DEMF, DIST, FLATS); + */ + +#include "matrix.h" +#include "mex.h" +#include "topotoolbox.h" + +#include +#include +#include + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { + // Validate input and output arguments + if (nrhs != 3) { + mexErrMsgIdAndTxt("tt3:tt_flow_routing_d8_carve:nrhs", "Three inputs required"); + } + if (nlhs != 3) { + mexErrMsgIdAndTxt("tt3:tt_flow_routing_d8_carve:nlhs", "Three outputs required"); + } + + float *demf = mxGetSingles(prhs[0]); + float *costs = mxGetSingles(prhs[1]); + int32_t *flats = mxGetInt32s(prhs[2]); + + ptrdiff_t dims[2] = {mxGetM(prhs[0]), mxGetN(prhs[0])}; + + plhs[0] = mxCreateNumericMatrix(dims[0] * dims[1], 1, mxINT64_CLASS, mxREAL); + int64_t *source = mxGetInt64s(plhs[0]); + + plhs[1] = mxCreateNumericMatrix(dims[0] * dims[1], 1, mxINT64_CLASS, mxREAL); + int64_t *target = mxGetInt64s(plhs[1]); + + ptrdiff_t *node = mxMalloc(sizeof(ptrdiff_t) * dims[0] * dims[1]); + uint8_t *direction = mxMalloc(sizeof(uint8_t) * dims[0] * dims[1]); + + flow_routing_d8_carve(node, direction, demf, costs, flats, dims, 0); + + ptrdiff_t edge_count = + flow_routing_d8_edgelist(source, target, node, direction, dims, 0); + + plhs[2] = mxCreateDoubleScalar(edge_count); + + mxFree(node); + mxFree(direction); +} diff --git a/bindings/tt_gwdt.c b/bindings/tt_gwdt.c new file mode 100644 index 0000000..0f18455 --- /dev/null +++ b/bindings/tt_gwdt.c @@ -0,0 +1,60 @@ +/* + +tt_gwdt.c + +[AUXTOPO] = tt_gwdt(DEM, DEMF, FLATS, COSTS) + +FLATS has to be the 32 bit encoded array returned by tt_identifyflats. + +*/ + +#include "matrix.h" +#include "mex.h" +#include "topotoolbox.h" + +#include +#include +#include + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { + // Validate input and output arguments + if (nrhs != 4) { + mexErrMsgIdAndTxt("tt3:tt_gwdt:nrhs", "Four inputs required"); + } + if (nlhs != 1) { + mexErrMsgIdAndTxt("tt3:tt_gwdt:nlhs", "One outputs required"); + } + // Extract input and output array data and dimensions + float *dem = mxGetSingles(prhs[0]); + float *demf = mxGetSingles(prhs[1]); + int32_t *flats = mxGetInt32s(prhs[2]); + float *costs = mxGetSingles(prhs[3]); + + // Create the dimensions. + ptrdiff_t dims[2] = {mxGetM(prhs[0]), mxGetN(prhs[0])}; + + // Create output and intermediate arrays + plhs[0] = mxCreateNumericMatrix(dims[0], dims[1], mxSINGLE_CLASS, mxREAL); + float *auxtopo = mxGetSingles(plhs[0]); + + // Allocate necessary intermediate arrays + if (dims[0] > 0 && dims[1] > PTRDIFF_MAX / dims[0]) { + mexErrMsgIdAndTxt("tt3:tt_gwdt:mxMalloc", + "Element count overflows"); + } + ptrdiff_t count = dims[0] * dims[1]; + + if (count > PTRDIFF_MAX / sizeof(ptrdiff_t)) { + mexErrMsgIdAndTxt("tt3:tt_gwdt:mxMalloc", + "Intermediate array size overflows"); + } + + ptrdiff_t *heap = mxMalloc(sizeof(ptrdiff_t) * count); + ptrdiff_t *back = mxMalloc(sizeof(ptrdiff_t) * count); + + gwdt(auxtopo, NULL, costs, flats, heap, back, dims); + + // Destroy intermediate arrays + mxFree(heap); + mxFree(back); +} diff --git a/tests/testSnapshot.m b/tests/testSnapshot.m index f44ef4d..188d03b 100644 --- a/tests/testSnapshot.m +++ b/tests/testSnapshot.m @@ -190,7 +190,7 @@ function auxtopo(testCase, dataset, uselibtt) D.GRIDobj2geotiff(result_file); else D_result = GRIDobj(result_file); - testCase.verifyEqual(D_result.Z, D.Z); + testCase.verifyEqual(D_result.Z, D.Z, RelTol=single(1e-5)); end end diff --git a/toolbox/@FLOWobj/FLOWobj.m b/toolbox/@FLOWobj/FLOWobj.m index a550597..bab0dd5 100644 --- a/toolbox/@FLOWobj/FLOWobj.m +++ b/toolbox/@FLOWobj/FLOWobj.m @@ -157,7 +157,7 @@ % If there are input arguments, the first input must be a GRIDobj assert(isa(DEM,'GRIDobj'),'TopoToolbox:WrongInput', ... 'First input argument must be a GRIDobj') - type = validatestring(type,{'single','multi','dinf'},'FLOWobj','type',2); + type = validatestring(type,{'single','libtt', 'multi','dinf'},'FLOWobj','type',2); % Transfer DEM (GRIDobj) properties to FLOWobj @@ -169,7 +169,7 @@ nrc = numel(DEM.Z); % Compute auxiliary topography - [D,DEM,SILLS] = createAuxiliaryTopo(DEM,... + [D, DEM, FLATS, SILLS] = createAuxiliaryTopo(DEM,... 'internaldrainage',options.internaldrainage,... 'preprocess',options.preprocess,... 'verbose',options.verbose,... @@ -182,9 +182,18 @@ % SILLS = logical matrix with locations of sills switch type - case 'single' - - if options.mex + case 'single' + if options.uselibtt + D(~FLATS) = 0.0; + [ix, ixc, count] = tt_flow_routing_d8_carve( ... + single(DEM.Z), single(D), int32(FLATS)); + + % libtopotoolbox ix and ixc arrays are 0-based. Add 1 + % to make them 1-based for MATLAB. + FD.ix = uint32(ix(1:count)) + 1; + FD.ixc = uint32(ixc(1:count)) + 1; + return + elseif options.mex % mexed version: % identifies steepest downward neighbors in the DEM and % the distance grid obtained from graydist and performs @@ -292,7 +301,6 @@ disp([char(datetime("now")) ' -- Ordered topology established']) end - case 'multi' %% Multiple flow direction D(SILLS) = 0; diff --git a/toolbox/internal/createAuxiliaryTopo.m b/toolbox/internal/createAuxiliaryTopo.m index e225260..aed5904 100644 --- a/toolbox/internal/createAuxiliaryTopo.m +++ b/toolbox/internal/createAuxiliaryTopo.m @@ -1,4 +1,4 @@ -function [D,DEMF,SILLS] = createAuxiliaryTopo(DEM,options) +function [D,DEMF, I, SILLS] = createAuxiliaryTopo(DEM,options) %CREATEAUXILIARYTOPO Calculate auxiliary topography from DEM % @@ -78,14 +78,25 @@ % internally drained basin (e.g. a flat lake surface). if options.internaldrainage [Iobj,SILLSobj,IntBasin] = identifyflats(DEMF,'uselibtt',options.uselibtt); + I = Iobj.Z; + SILLS = SILLSobj.Z; + + clear Iobj SILLSobj else - [Iobj,SILLSobj] = identifyflats(DEMF,'uselibtt',options.uselibtt); + if options.uselibtt + FLATS = tt_identifyflats(single(DEMF.Z)); + I = bitget(FLATS, 1) == 1; + SILLS = bitget(FLATS, 2) == 1; + else + [Iobj,SILLSobj] = identifyflats(DEMF,'uselibtt',options.uselibtt); + I = Iobj.Z; + SILLS = SILLSobj.Z; + + clear Iobj SILLSobj + end end -I = Iobj.Z; -SILLS = SILLSobj.Z; -clear Iobj SILLSobj % calculate sills for internal lake basins. These should be % located within the lake to force convergent flows @@ -140,7 +151,9 @@ % of pixels. It is correct to supply an int32 version of % the logical I array, since that will have 1 for all flats and % 0 for all other pixels. - D = tt_gwdt_computecosts(single(DEM.Z), single(DEMF.Z), int32(I)); + D = tt_gwdt_computecosts(single(DEM.Z), single(DEMF.Z), int32(FLATS)); + D = tt_gwdt(single(DEM.Z), single(DEMF.Z), int32(FLATS), single(D)); + D(~I) = -inf; else D = DEMF.Z-DEM.Z; @@ -191,12 +204,15 @@ disp([char(datetime("now")) ' -- Weights for graydist calculated']) end -% Here we calculate the auxiliary topography. That is, the -% cost surface seeded at socalled PreSillPixels, i.e. the -% pixel immediately upstream to sill pixels. -D(I) = inf; -D = graydist(double(D),double(PreSillPixel),'q') + 1; -D(I) = -inf; +if ~options.uselibtt + % Here we calculate the auxiliary topography. That is, the + % cost surface seeded at socalled PreSillPixels, i.e. the + % pixel immediately upstream to sill pixels. + D(I) = inf; + D = graydist(double(D),double(PreSillPixel),'q') + 1; + D(I) = -inf; +end +I = ~I; if options.verbose disp([char(datetime("now")) ' -- Auxiliary topography in flats calculated'])