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'])