r.flowaccumulation: Support weighted flow accumulation (#1048)

* r.flowaccumulation: Support weighted flow accumulation

* Promote accum_map type if it's less than weight_map type

* Allow type for weighted accumulation

* Support weight map whose type is different from accumulation

src/raster/r.flowaccumulation/accumulate.c

@@ -1,8 +1,9 @@
 #include <grass/raster.h>
 #include "global.h"
 
-void accumulate(struct raster_map *dir_map, struct raster_map *accum_map,
-                int check_overflow, int use_less_memory, int use_zero)
+void accumulate(struct raster_map *dir_map, struct raster_map *weight_map,
+                struct raster_map *accum_map, int check_overflow,
+                int use_less_memory, int use_zero)
 {
     switch (accum_map->type) {
     case CELL_TYPE:
@@ -24,12 +25,16 @@ void accumulate(struct raster_map *dir_map, struct raster_map *accum_map,
             if (use_less_memory) {
                 if (use_zero)
                     accumulate_cmz(dir_map, accum_map);
+                else if (weight_map)
+                    accumulate_cmw(dir_map, weight_map, accum_map);
                 else
                     accumulate_cm(dir_map, accum_map);
             }
             else {
                 if (use_zero)
                     accumulate_cz(dir_map, accum_map);
+                else if (weight_map)
+                    accumulate_cw(dir_map, weight_map, accum_map);
                 else
                     accumulate_c(dir_map, accum_map);
             }
@@ -54,12 +59,16 @@ void accumulate(struct raster_map *dir_map, struct raster_map *accum_map,
             if (use_less_memory) {
                 if (use_zero)
                     accumulate_fmz(dir_map, accum_map);
+                else if (weight_map)
+                    accumulate_fmw(dir_map, weight_map, accum_map);
                 else
                     accumulate_fm(dir_map, accum_map);
             }
             else {
                 if (use_zero)
                     accumulate_fz(dir_map, accum_map);
+                else if (weight_map)
+                    accumulate_fw(dir_map, weight_map, accum_map);
                 else
                     accumulate_f(dir_map, accum_map);
             }
@@ -84,12 +93,16 @@ void accumulate(struct raster_map *dir_map, struct raster_map *accum_map,
             if (use_less_memory) {
                 if (use_zero)
                     accumulate_dmz(dir_map, accum_map);
+                else if (weight_map)
+                    accumulate_dmw(dir_map, weight_map, accum_map);
                 else
                     accumulate_dm(dir_map, accum_map);
             }
             else {
                 if (use_zero)
                     accumulate_dz(dir_map, accum_map);
+                else if (weight_map)
+                    accumulate_dw(dir_map, weight_map, accum_map);
                 else
                     accumulate_d(dir_map, accum_map);
             }

src/raster/r.flowaccumulation/accumulate_cmw.c

@@ -0,0 +1,4 @@
+#define ACCUM_RAST_TYPE CELL_TYPE
+#define USE_LESS_MEMORY
+#define USE_WEIGHT
+#include "accumulate_funcs.h"

src/raster/r.flowaccumulation/accumulate_cw.c

@@ -0,0 +1,3 @@
+#define ACCUM_RAST_TYPE CELL_TYPE
+#define USE_WEIGHT
+#include "accumulate_funcs.h"

src/raster/r.flowaccumulation/accumulate_dmw.c

@@ -0,0 +1,4 @@
+#define ACCUM_RAST_TYPE DCELL_TYPE
+#define USE_LESS_MEMORY
+#define USE_WEIGHT
+#include "accumulate_funcs.h"

src/raster/r.flowaccumulation/accumulate_dw.c

@@ -0,0 +1,3 @@
+#define ACCUM_RAST_TYPE DCELL_TYPE
+#define USE_WEIGHT
+#include "accumulate_funcs.h"

src/raster/r.flowaccumulation/accumulate_fmw.c

@@ -0,0 +1,4 @@
+#define ACCUM_RAST_TYPE FCELL_TYPE
+#define USE_LESS_MEMORY
+#define USE_WEIGHT
+#include "accumulate_funcs.h"

src/raster/r.flowaccumulation/accumulate_funcs.h

@@ -9,29 +9,48 @@
 #endif
 
 #if ACCUM_RAST_TYPE == CELL_TYPE
-#define ACCUMULATE(o, m, z)  accumulate_c##o##m##z
-#define NULLIFY_ZERO         nullify_zero_c
-#define ACCUM_MAP_CELLS      accum_map->cells.c
-#define ACCUM_TYPE           int
-#define SET_ACCUM_NULL(cell) Rast_set_c_null_value(cell, 1)
-#define IS_ACCUM_NULL(cell)  Rast_is_c_null_value(cell)
+#define ACCUMULATE(o, m, z, w) accumulate_c##o##m##z##w
+#define NULLIFY_ZERO           nullify_zero_c
+#define ACCUM_MAP_CELLS        accum_map->cells.c
+#define ACCUM_TYPE             int
+#define SET_ACCUM_NULL(cell)   Rast_set_c_null_value(cell, 1)
+#define IS_ACCUM_NULL(cell)    Rast_is_c_null_value(cell)
 #elif ACCUM_RAST_TYPE == FCELL_TYPE
-#define ACCUMULATE(o, m, z)  accumulate_f##o##m##z
-#define NULLIFY_ZERO         nullify_zero_f
-#define ACCUM_MAP_CELLS      accum_map->cells.f
-#define ACCUM_TYPE           float
-#define SET_ACCUM_NULL(cell) Rast_set_f_null_value(cell, 1)
-#define IS_ACCUM_NULL(cell)  Rast_is_f_null_value(cell)
+#define ACCUMULATE(o, m, z, w) accumulate_f##o##m##z##w
+#define NULLIFY_ZERO           nullify_zero_f
+#define ACCUM_MAP_CELLS        accum_map->cells.f
+#define ACCUM_TYPE             float
+#define SET_ACCUM_NULL(cell)   Rast_set_f_null_value(cell, 1)
+#define IS_ACCUM_NULL(cell)    Rast_is_f_null_value(cell)
 #else
-#define ACCUMULATE(o, m, z)  accumulate_d##o##m##z
-#define NULLIFY_ZERO         nullify_zero_d
-#define ACCUM_MAP_CELLS      accum_map->cells.d
-#define ACCUM_TYPE           double
-#define SET_ACCUM_NULL(cell) Rast_set_d_null_value(cell, 1)
-#define IS_ACCUM_NULL(cell)  Rast_is_d_null_value(cell)
+#define ACCUMULATE(o, m, z, w) accumulate_d##o##m##z##w
+#define NULLIFY_ZERO           nullify_zero_d
+#define ACCUM_MAP_CELLS        accum_map->cells.d
+#define ACCUM_TYPE             double
+#define SET_ACCUM_NULL(cell)   Rast_set_d_null_value(cell, 1)
+#define IS_ACCUM_NULL(cell)    Rast_is_d_null_value(cell)
+#endif
+
+#define ACCUM(row, col) ACCUM_MAP_CELLS[INDEX(row, col)]
+
+#ifdef USE_WEIGHT
+#define WEIGHT(row, col)                               \
+    (weight_map->type == CELL_TYPE                     \
+         ? weight_map->cells.c[INDEX(row, col)]        \
+         : (weight_map->type == FCELL_TYPE             \
+                ? weight_map->cells.f[INDEX(row, col)] \
+                : weight_map->cells.d[INDEX(row, col)]))
+#define IS_WEIGHT_NULL(row, col)                                              \
+    (weight_map->type == CELL_TYPE                                            \
+         ? Rast_is_c_null_value(&weight_map->cells.c[INDEX(row, col)])        \
+         : (weight_map->type == FCELL_TYPE                                    \
+                ? Rast_is_f_null_value(&weight_map->cells.f[INDEX(row, col)]) \
+                : Rast_is_d_null_value(                                       \
+                      &weight_map->cells.d[INDEX(row, col)])))
+#else
+#define WEIGHT(row, col) 1
 #endif
 
-#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) |                        \
@@ -48,14 +67,18 @@
 #ifdef USE_LESS_MEMORY
 #define UP(row, col) FIND_UP(row, col)
 #else
-#define UP(row, col) up_cells[(size_t)(row)*ncols + (col)]
+#define UP(row, col) up_cells[INDEX(row, col)]
 static unsigned char *up_cells;
 #endif
 
 static int nrows, ncols;
+static ACCUM_TYPE accum_not_ready;
 
-static void trace_down(struct raster_map *, struct raster_map *, int, int,
-                       ACCUM_TYPE);
+static void trace_down(struct raster_map *,
+#ifdef USE_WEIGHT
+                       struct raster_map *,
+#endif
+                       struct raster_map *, int, int, ACCUM_TYPE);
 static ACCUM_TYPE sum_up(struct raster_map *, int, int, int);
 
 void ACCUMULATE(
@@ -70,13 +93,27 @@ void ACCUMULATE(
 #ifdef USE_ZERO
     z
 #endif
-    )(struct raster_map *dir_map, struct raster_map *accum_map)
+    ,
+#ifdef USE_WEIGHT
+    w
+#endif
+    )(struct raster_map *dir_map,
+#ifdef USE_WEIGHT
+      struct raster_map *weight_map,
+#endif
+      struct raster_map *accum_map)
 {
     int row, col;
 
     nrows = dir_map->nrows;
     ncols = dir_map->ncols;
 
+#ifdef USE_WEIGHT
+    SET_ACCUM_NULL(&accum_not_ready);
+#else
+    accum_not_ready = 0;
+#endif
+
 #ifndef USE_ZERO
 #pragma omp parallel for schedule(dynamic)
     for (row = 0; row < nrows; row++)
@@ -101,16 +138,25 @@ void ACCUMULATE(
         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) && !UP(row, col))
-                trace_down(dir_map, accum_map, row, col, 1);
+            if (DIR(row, col) && !UP(row, col)
+#ifdef USE_WEIGHT
+                && !IS_WEIGHT_NULL(row, col)
+#endif
+            )
+                trace_down(dir_map,
+#ifdef USE_WEIGHT
+                           weight_map,
+#endif
+                           accum_map, row, col, WEIGHT(row, col));
     }
 
 #ifndef USE_LESS_MEMORY
     free(up_cells);
 #endif
 }
 
-#if !defined CHECK_OVERFLOW && !defined USE_LESS_MEMORY && !defined USE_ZERO
+#if !defined CHECK_OVERFLOW && !defined USE_LESS_MEMORY && \
+    !defined USE_ZERO && !defined USE_WEIGHT
 void NULLIFY_ZERO(struct raster_map *accum_map)
 {
     int row, col;
@@ -126,8 +172,12 @@ void NULLIFY_ZERO(struct raster_map *accum_map)
 }
 #endif
 
-static void trace_down(struct raster_map *dir_map, struct raster_map *accum_map,
-                       int row, int col, ACCUM_TYPE accum)
+static void trace_down(struct raster_map *dir_map,
+#ifdef USE_WEIGHT
+                       struct raster_map *weight_map,
+#endif
+                       struct raster_map *accum_map, int row, int col,
+                       ACCUM_TYPE accum)
 {
     int up;
     ACCUM_TYPE accum_up = 0;
@@ -189,12 +239,27 @@ static void trace_down(struct raster_map *dir_map, struct raster_map *accum_map,
     /* 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) ||
-        !(up = UP(row, col)) || !(accum_up = sum_up(accum_map, row, col, up)))
+        !(up = UP(row, col))
+#ifndef USE_WEIGHT
+        || (accum_up = sum_up(accum_map, row, col, up)) == accum_not_ready
+#endif
+    )
         return;
+#ifdef USE_WEIGHT
+    else {
+        accum_up = sum_up(accum_map, row, col, up);
+        if (IS_ACCUM_NULL(&accum_up))
+            return;
+    }
+#endif
 
     /* use gcc -O2 or -O3 flags for tail-call optimization
      * (-foptimize-sibling-calls) */
-    trace_down(dir_map, accum_map, row, col, accum_up + 1);
+    trace_down(dir_map,
+#ifdef USE_WEIGHT
+               weight_map,
+#endif
+               accum_map, row, col, accum_up + WEIGHT(row, col));
 }
 
 /* if any upstream cells have never been visited, 0 is returned; otherwise, the
@@ -211,7 +276,7 @@ static ACCUM_TYPE sum_up(struct raster_map *accum_map, int row, int col, int up)
         accum = ACCUM(row - 1, col - 1);
         if (IS_ACCUM_NULL(&accum))
 #endif
-            return 0;
+            return accum_not_ready;
         sum += accum;
     }
     if (up & N) {
@@ -221,7 +286,7 @@ static ACCUM_TYPE sum_up(struct raster_map *accum_map, int row, int col, int up)
         accum = ACCUM(row - 1, col);
         if (IS_ACCUM_NULL(&accum))
 #endif
-            return 0;
+            return accum_not_ready;
         sum += accum;
     }
     if (up & NE) {
@@ -231,7 +296,7 @@ static ACCUM_TYPE sum_up(struct raster_map *accum_map, int row, int col, int up)
         accum = ACCUM(row - 1, col + 1);
         if (IS_ACCUM_NULL(&accum))
 #endif
-            return 0;
+            return accum_not_ready;
         sum += accum;
     }
     if (up & W) {
@@ -241,7 +306,7 @@ static ACCUM_TYPE sum_up(struct raster_map *accum_map, int row, int col, int up)
         accum = ACCUM(row, col - 1);
         if (IS_ACCUM_NULL(&accum))
 #endif
-            return 0;
+            return accum_not_ready;
         sum += accum;
     }
     if (up & E) {
@@ -251,7 +316,7 @@ static ACCUM_TYPE sum_up(struct raster_map *accum_map, int row, int col, int up)
         accum = ACCUM(row, col + 1);
         if (IS_ACCUM_NULL(&accum))
 #endif
-            return 0;
+            return accum_not_ready;
         sum += accum;
     }
     if (up & SW) {
@@ -261,7 +326,7 @@ static ACCUM_TYPE sum_up(struct raster_map *accum_map, int row, int col, int up)
         accum = ACCUM(row + 1, col - 1);
         if (IS_ACCUM_NULL(&accum))
 #endif
-            return 0;
+            return accum_not_ready;
         sum += accum;
     }
     if (up & S) {
@@ -271,7 +336,7 @@ static ACCUM_TYPE sum_up(struct raster_map *accum_map, int row, int col, int up)
         accum = ACCUM(row + 1, col);
         if (IS_ACCUM_NULL(&accum))
 #endif
-            return 0;
+            return accum_not_ready;
         sum += accum;
     }
     if (up & SE) {
@@ -281,7 +346,7 @@ static ACCUM_TYPE sum_up(struct raster_map *accum_map, int row, int col, int up)
         accum = ACCUM(row + 1, col + 1);
         if (IS_ACCUM_NULL(&accum))
 #endif
-            return 0;
+            return accum_not_ready;
         sum += accum;
     }
 

src/raster/r.flowaccumulation/accumulate_fw.c

@@ -0,0 +1,3 @@
+#define ACCUM_RAST_TYPE FCELL_TYPE
+#define USE_WEIGHT
+#include "accumulate_funcs.h"