r.flowaccumulation: Add a new addon that implements the Memory-Effici…

…ent Flow Accumulation parallel algorithm (Cho 2023) (#954)

* r.accum: Add a new addon that implements the Memory-Efficient Flow Accumulation parallel algorithm (Cho 2023)

* Fix format check

* Add ArcGIS format

* Rename to r.flowaccumulation

* Rename to r.flowaccumulation

* Fix link to r.accumulate manual

* Update src/raster/r.flowaccumulation/r.flowaccumulation.html

Co-authored-by: Markus Neteler <neteler@osgeo.org>

* Update src/raster/r.flowaccumulation/r.flowaccumulation.html

Co-authored-by: Markus Neteler <neteler@osgeo.org>

* Minor formatting

---------

Co-authored-by: Markus Neteler <neteler@osgeo.org>

src/raster/r.flowaccumulation/Makefile

@@ -0,0 +1,13 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.flowaccumulation
+
+LIBES = $(RASTERLIB) $(VECTORLIB) $(DBMILIB) $(GISLIB) $(MATHLIB)
+DEPENDENCIES = $(RASTERDEP) $(VECTORDEP) $(DBMIDEP) $(GISDEP)
+EXTRA_LIBS = $(OPENMP_LIBPATH) $(OPENMP_LIB)
+EXTRA_INC = $(VECT_INC) $(OPENMP_INCPATH)
+EXTRA_CFLAGS = $(VECT_CFLAGS) $(OPENMP_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd

src/raster/r.flowaccumulation/accumulate.c

@@ -0,0 +1,164 @@
+#include <stdlib.h>
+#include "global.h"
+
+#define DIR_NULL        0x80000000
+#define ACCUM(row, col) accum_map->cells[(size_t)(row)*ncols + (col)]
+#define FIND_UP(row, col)                                                     \
+    ((row > 0 ? (col > 0 && DIR(row - 1, col - 1) == SE ? NW : 0) |           \
+                    (DIR(row - 1, col) == S ? N : 0) |                        \
+                    (col < ncols - 1 && DIR(row - 1, col + 1) == SW ? NE : 0) \
+              : 0) |                                                          \
+     (col > 0 && DIR(row, col - 1) == E ? W : 0) |                            \
+     (col < ncols - 1 && DIR(row, col + 1) == W ? E : 0) |                    \
+     (row < nrows - 1                                                         \
+          ? (col > 0 && DIR(row + 1, col - 1) == NE ? SW : 0) |               \
+                (DIR(row + 1, col) == N ? S : 0) |                            \
+                (col < ncols - 1 && DIR(row + 1, col + 1) == NW ? SE : 0)     \
+          : 0))
+
+#ifdef USE_LESS_MEMORY
+#define UP(row, col) FIND_UP(row, col)
+#else
+#define UP(row, col) up_cells[(size_t)(row)*ncols + (col)]
+static unsigned char *up_cells;
+#endif
+
+static int nrows, ncols;
+
+static void trace_down(struct raster_map *, struct raster_map *, int, int, int);
+static int sum_up(struct raster_map *, int, int, int);
+
+void accumulate(struct raster_map *dir_map, struct raster_map *accum_map)
+{
+    int row, col;
+
+    nrows = dir_map->nrows;
+    ncols = dir_map->ncols;
+
+#ifndef USE_LESS_MEMORY
+    up_cells = calloc((size_t)nrows * ncols, sizeof *up_cells);
+
+#pragma omp parallel for schedule(dynamic) private(col)
+    for (row = 0; row < nrows; row++) {
+        for (col = 0; col < ncols; col++)
+            if (DIR(row, col) != DIR_NULL)
+                UP(row, col) = FIND_UP(row, col);
+    }
+#endif
+
+#pragma omp parallel for schedule(dynamic) private(col)
+    for (row = 0; row < nrows; row++) {
+        for (col = 0; col < ncols; col++)
+            /* if the current cell is not null and has no upstream cells, start
+             * tracing down */
+            if (DIR(row, col) != DIR_NULL && !UP(row, col))
+                trace_down(dir_map, accum_map, row, col, 1);
+    }
+
+#ifndef USE_LESS_MEMORY
+    free(up_cells);
+#endif
+}
+
+static void trace_down(struct raster_map *dir_map, struct raster_map *accum_map,
+                       int row, int col, int accum)
+{
+    int up, accum_up = 0;
+
+    /* accumulate the current cell itself */
+    ACCUM(row, col) = accum;
+
+    /* find the downstream cell */
+    switch (DIR(row, col)) {
+    case NW:
+        row--;
+        col--;
+        break;
+    case N:
+        row--;
+        break;
+    case NE:
+        row--;
+        col++;
+        break;
+    case W:
+        col--;
+        break;
+    case E:
+        col++;
+        break;
+    case SW:
+        row++;
+        col--;
+        break;
+    case S:
+        row++;
+        break;
+    case SE:
+        row++;
+        col++;
+        break;
+    }
+
+    /* if the downstream cell is null or any upstream cells of the downstream
+     * cell have never been visited, stop tracing down */
+    if (row < 0 || row >= nrows || col < 0 || col >= ncols ||
+        DIR(row, col) == DIR_NULL || !(up = UP(row, col)) ||
+        !(accum_up = sum_up(accum_map, row, col, up)))
+        return;
+
+    /* use gcc -O2 or -O3 flags for tail-call optimization
+     * (-foptimize-sibling-calls) */
+    trace_down(dir_map, accum_map, row, col, accum_up + 1);
+}
+
+/* if any upstream cells have never been visited, 0 is returned; otherwise, the
+ * sum of upstream accumulation is returned */
+static int sum_up(struct raster_map *accum_map, int row, int col, int up)
+{
+    int sum = 0, accum;
+
+#pragma omp flush(accum_map)
+    if (up & NW) {
+        if (!(accum = ACCUM(row - 1, col - 1)))
+            return 0;
+        sum += accum;
+    }
+    if (up & N) {
+        if (!(accum = ACCUM(row - 1, col)))
+            return 0;
+        sum += accum;
+    }
+    if (up & NE) {
+        if (!(accum = ACCUM(row - 1, col + 1)))
+            return 0;
+        sum += accum;
+    }
+    if (up & W) {
+        if (!(accum = ACCUM(row, col - 1)))
+            return 0;
+        sum += accum;
+    }
+    if (up & E) {
+        if (!(accum = ACCUM(row, col + 1)))
+            return 0;
+        sum += accum;
+    }
+    if (up & SW) {
+        if (!(accum = ACCUM(row + 1, col - 1)))
+            return 0;
+        sum += accum;
+    }
+    if (up & S) {
+        if (!(accum = ACCUM(row + 1, col)))
+            return 0;
+        sum += accum;
+    }
+    if (up & SE) {
+        if (!(accum = ACCUM(row + 1, col + 1)))
+            return 0;
+        sum += accum;
+    }
+
+    return sum;
+}

src/raster/r.flowaccumulation/gettimeofday.c

@@ -0,0 +1,28 @@
+#ifdef _MSC_VER
+#include <stdint.h>
+#include <winsock2.h>
+
+// https://stackoverflow.com/a/26085827/16079666
+int gettimeofday(struct timeval *tp, struct timezone *tzp)
+{
+    // Note: some broken versions only have 8 trailing zero's, the correct
+    // epoch has 9 trailing zero's This magic number is the number of 100
+    // nanosecond intervals since January 1, 1601 (UTC) until 00:00:00 January
+    // 1, 1970
+    static const uint64_t EPOCH = ((uint64_t)116444736000000000ULL);
+
+    SYSTEMTIME system_time;
+    FILETIME file_time;
+    uint64_t time;
+
+    GetSystemTime(&system_time);
+    SystemTimeToFileTime(&system_time, &file_time);
+    time = ((uint64_t)file_time.dwLowDateTime);
+    time += ((uint64_t)file_time.dwHighDateTime) << 32;
+
+    tp->tv_sec = (long)((time - EPOCH) / 10000000L);
+    tp->tv_usec = (long)(system_time.wMilliseconds * 1000);
+
+    return 0;
+}
+#endif

src/raster/r.flowaccumulation/global.h

@@ -0,0 +1,35 @@
+#ifndef _GLOBAL_H_
+#define _GLOBAL_H_
+
+#ifdef _MSC_VER
+#include <winsock2.h>
+/* gettimeofday.c */
+int gettimeofday(struct timeval *, struct timezone *);
+#else
+#include <sys/time.h>
+#endif
+#include <grass/raster.h>
+
+#define E             1
+#define SE            2
+#define S             4
+#define SW            8
+#define W             16
+#define NW            32
+#define N             64
+#define NE            128
+
+#define DIR(row, col) dir_map->cells[(size_t)(row)*ncols + (col)]
+
+struct raster_map {
+    int nrows, ncols;
+    CELL *cells;
+};
+
+/* timeval_diff.c */
+long long timeval_diff(struct timeval *, struct timeval *, struct timeval *);
+
+/* accumulate.c */
+void accumulate(struct raster_map *, struct raster_map *);
+
+#endif