Implement libtopotoolbox-based flow routing

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

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 <limits.h>
+#include <stddef.h>
+#include <stdint.h>
+
+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);
+}

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 <limits.h>
+#include <stddef.h>
+#include <stdint.h>
+
+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);
+}

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
 

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;

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