Add MFD implementation

include/topotoolbox.h

@@ -3391,4 +3391,108 @@ void resolve_flats_lcat(uint8_t *direction, uint8_t *resolved, float *aux,
 TOPOTOOLBOX_API
 void resolve_flats_lcat_weights(float *weight, ptrdiff_t count);
 
+/**
+   @brief Compute multiple flow directions
+
+   @details MFD distributes flow to all downstream neighbors with
+   weights proportional to the gradient.
+
+   @param[out] direction A bitmap edge set
+   @parblock
+
+   A uint8_t array of the same size as the DEM. Each of the 8 bits of
+   each pixel corresponds to an outgoing edge.
+
+   @endparblock
+
+   @param[out] totalgradient Sum of downslope gradients
+   @parblock
+
+   A float array of the same size as the DEM. This is the sum of the
+   gradients from each pixel to its downslope neighbors. This is not
+   particularly useful, but it is needed to normalize the weights and it
+   is easy to compute during the first pass.
+
+   @endparblock
+
+   @param[in] dem The input DEM
+   @parblock
+
+   A float array with size dims[0] x dims[1].
+
+   @endparblock
+
+   @param[in]  dims The dimensions of the DEM
+   @parblock
+   A pair of ptrdiff_t, fastest changing dimension first.
+   @endparblock
+
+   @param[in] order The memory order of the underlying arrays
+   @parblock
+   0 for column-major, 1 for row-major
+   @endparblock
+ */
+TOPOTOOLBOX_API
+void flow_routing_mfd_directions(uint8_t *direction, float *totalgradient,
+                                 float *dem, ptrdiff_t dims[2], int order);
+
+/**
+   @brief Compute and store the MFD edge weights
+
+   @details Once the flow directions and flow proportions have been
+   computed, the edge weights need to be stored in an array in an
+   order that will be used by functions like flow_routing_tsort.
+
+   @param[out] weight The weight array
+   @parblock
+
+   An array of floats that will be filled with the edge weights. This
+   should be preallocated with a size equal to the number of edges in
+   the direction bitmap. This size can be computed using
+   edgeset_count.
+
+   @endparblock
+
+   @param[in] direction The input direction bitmap
+   @parblock
+
+   An array of uint8_t that should come from the direction argument of
+   flow_routing_dinf_directions. This is a bitmap edge set encoded as
+   described elsewhere.
+
+   @endparblock
+
+   @param[out] totalgradient Sum of downslope gradients
+   @parblock
+
+   A float array of the same size as the DEM. This should come from
+   flow_routing_mfd_directions. It is the sum of the gradients from
+   each pixel to its downslope neighbors. It is temporary and can be
+   discarded after the weights are computed.
+
+   @endparblock
+
+
+   @param[in] dem The input DEM
+   @parblock
+
+   A float array with size dims[0] x dims[1].
+
+   @endparblock
+
+   @param[in]  dims The dimensions of the DEM
+   @parblock
+   A pair of ptrdiff_t, fastest changing dimension first.
+   @endparblock
+
+   @param[in] order The memory order of the underlying arrays
+   @parblock
+   0 for column-major, 1 for row-major
+   @endparblock
+ */
+TOPOTOOLBOX_API
+void flow_routing_mfd_weights(float *weight, float *totalgradient,
+                              uint8_t *direction, float *dem, ptrdiff_t dims[2],
+                              int order);
+
 #endif  // TOPOTOOLBOX_H

src/CMakeLists.txt

@@ -40,6 +40,7 @@ add_library(topotoolbox
   dinf.c
   d8.c
   lcat.c
+  mfd.c
 )
 
 # Define the include directory

src/mfd.c

@@ -0,0 +1,64 @@
+#define TOPOTOOLBOX_BUILD
+
+#include <assert.h>
+#include <math.h>
+#include <stddef.h>
+#include <stdint.h>
+
+#include "topotoolbox.h"
+
+#define SQRT2 1.41421356237309504880f
+
+TOPOTOOLBOX_API
+void flow_routing_mfd_directions(uint8_t *direction, float *totalgradient,
+                                 float *dem, ptrdiff_t dims[2], int order) {
+  ptrdiff_t e[2][8] = {{0, -1, -1, -1, 0, 1, 1, 1},
+                       {1, 1, 0, -1, -1, -1, 0, 1}};
+
+  float d[8] = {1.0, SQRT2, 1.0, SQRT2, 1.0, SQRT2, 1.0, SQRT2};
+
+  for (ptrdiff_t j = 0; j < dims[1]; j++) {
+    for (ptrdiff_t i = 0; i < dims[0]; i++) {
+      float z = dem[j * dims[0] + i];
+      for (int n = 0; n < 8; n++) {
+        ptrdiff_t i1 = i + e[order & 1][n];
+        ptrdiff_t j1 = j + e[(order ^ 1) & 1][n];
+
+        if (i1 < 0 || i1 >= dims[0] || j1 < 0 || j1 >= dims[1]) continue;
+
+        float z2 = dem[j1 * dims[0] + i1];
+        direction[j * dims[0] + i] |= (z > z2) << n;
+        totalgradient[j * dims[0] + i] += fmaxf(0.0, z - z2) / d[n];
+      }
+    }
+  }
+}
+
+TOPOTOOLBOX_API
+void flow_routing_mfd_weights(float *weight, float *totalgradient,
+                              uint8_t *direction, float *dem, ptrdiff_t dims[2],
+                              int order) {
+  ptrdiff_t e[2][8] = {{0, -1, -1, -1, 0, 1, 1, 1},
+                       {1, 1, 0, -1, -1, -1, 0, 1}};
+
+  float d[8] = {1.0, SQRT2, 1.0, SQRT2, 1.0, SQRT2, 1.0, SQRT2};
+
+  ptrdiff_t edge = 0;
+  for (ptrdiff_t j = 0; j < dims[1]; j++) {
+    for (ptrdiff_t i = 0; i < dims[0]; i++) {
+      ptrdiff_t idx = j * dims[0] + i;
+      float z = dem[idx];
+      float g = totalgradient[idx];
+
+      for (int n = 0; n < 8; n++) {
+        if (direction[idx] & (1 << n)) {
+          ptrdiff_t i1 = i + e[order & 1][n];
+          ptrdiff_t j1 = j + e[(order ^ 1) & 1][n];
+          float z2 = dem[j1 * dims[0] + i1];
+
+          weight[edge++] = ((z - z2) / d[n]) / g;
+        }
+      }
+    }
+  }
+}

test/CMakeLists.txt

@@ -176,6 +176,18 @@ if(GDAL_FOUND AND TT_SNAPSHOT_DIR)
   add_test(NAME d8 COMMAND d8 ${TT_SNAPSHOT_DIR})
   set_tests_properties(d8 PROPERTIES ENVIRONMENT_MODIFICATION
     "PATH=path_list_prepend:$<$<BOOL:${WIN32}>:$<TARGET_FILE_DIR:topotoolbox>>")
+
+  # TEST : mfd
+  #
+  add_executable(mfd mfd.c utils.c grid.c flow.c)
+  if(TT_SANITIZE AND NOT MSVC)
+    target_compile_options(mfd PRIVATE "$<$<CONFIG:DEBUG>:-fsanitize=address>")
+    target_link_options(mfd PRIVATE "$<$<CONFIG:DEBUG>:-fsanitize=address>")
+  endif()
+  target_link_libraries(mfd PRIVATE topotoolbox GDAL::GDAL)
+  add_test(NAME mfd COMMAND mfd ${TT_SNAPSHOT_DIR})
+  set_tests_properties(mfd PROPERTIES ENVIRONMENT_MODIFICATION
+    "PATH=path_list_prepend:$<$<BOOL:${WIN32}>:$<TARGET_FILE_DIR:topotoolbox>>")
 endif()