Replace flow_accumulation_edgelist with traverse_down_f32_add_mul

R/flow_accumulation.R

@@ -42,21 +42,19 @@ flow_accumulation <- function(fd,
   edge_count <- length(fd$source)
 
   # Compute flow routing using libtopotoolbox
-  output <- single(prod(dims))
+  output <- as.single(ezgetnal(fd, weights))
   result <- .C(
-    "wrap_flow_accumulation_edgelist",
-    accR = as.single(output), # float
+    "wrap_traverse_down_f32_add_mul",
+    accR = output, # float
+    fractionR = as.single(rep(1, edge_count)), # float
     sourceR = as.integer(fd$source), # ptrdiff_t
     targetR = as.integer(fd$target), # ptrdiff_t
-    fractionR = as.single(rep(1, edge_count)), # float
-    weightsR = as.single(ezgetnal(fd, weights)), # float
     edge_countR = as.integer(edge_count), # ptrdiff_t
-    dimsR = as.integer(dims), # ptrdiff_t
     NAOK = TRUE
   )$accR
 
   # Write result into SpatRaster (for now)
   fa <- fd$raster
   terra::values(fa) <- result
   GRIDobj(fa)
-}
+}

src/wrap_streamquad.c

@@ -0,0 +1,32 @@
+#include <stddef.h>
+#include <stdint.h>
+#include <R.h>
+
+#include "topotoolbox.h"
+#include "topotoolboxr.h"
+
+void wrap_traverse_down_f32_add_mul(float *accR,      // output
+                                    float *fractionR, // input
+                                    int *sourceR,     // ptrdiff_t
+                                    int *targetR,     // ptrdiff_t
+                                    int *edge_countR  // ptrdiff_t
+                                    ) {
+
+  // Transformation of integers and array allocation
+  ptrdiff_t edge_count = edge_countR[0];
+  ptrdiff_t *source = R_Calloc(edge_count, ptrdiff_t);
+  ptrdiff_t *target = R_Calloc(edge_count, ptrdiff_t);
+
+  // Convert sourceR and targetR to ptrdiff_t
+  for (ptrdiff_t idx = 0; idx < edge_count; idx++) {
+    source[idx] = (ptrdiff_t)sourceR[idx];
+    target[idx] = (ptrdiff_t)targetR[idx];
+  }
+
+  // Flow accumulation computation using libtopotoolbox
+  traverse_down_f32_add_mul(accR, fractionR, source, target, edge_count);
+
+  // Free memory
+  R_Free(source);
+  R_Free(target);
+}