r.geomorphon: parallelize with OpenMP using memory-bounded banding

raster/CMakeLists.txt

@@ -244,6 +244,7 @@ build_program_in_subdir(
 build_program_in_subdir(
     r.geomorphon
     DEPENDS ${LIBM} grass_gis grass_gmath grass_parson grass_raster
+    OPTIONAL_DEPENDS OpenMP::OpenMP_C
 )
 
 build_program_in_subdir(r.grow.distance DEPENDS grass_gis grass_raster ${LIBM})

raster/r.geomorphon/Makefile

@@ -3,7 +3,10 @@ MODULE_TOPDIR = ../..
 PGM = r.geomorphon
 
 LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB) $(PARSONLIB)
+EXTRA_LIBS = $(OPENMP_LIBPATH) $(OPENMP_LIB)
 DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+EXTRA_CFLAGS = $(OPENMP_CFLAGS)
+EXTRA_INC = $(OPENMP_INCPATH)
 
 include $(MODULE_TOPDIR)/include/Make/Module.make
 

raster/r.geomorphon/benchmark/benchmark_r_geomorphon.py

@@ -0,0 +1,96 @@
+"""Benchmarking of r.geomorphon
+raster (2D)
+"""
+
+from grass.exceptions import CalledModuleError, GrassError
+from grass.pygrass.modules import Module
+import grass.benchmark as bm
+
+# Baselines held fixed while one dimension is swept.
+BASE_MAPSIZE = 50e6  # cells
+BASE_SEARCH = 25  # outer search radius in cells
+BASE_MEMORY = 300  # MB
+MAPSIZES = [10e6, 50e6, 100e6]
+SEARCHES = [10, 25, 50]
+METRICS = ["time", "speedup", "efficiency"]
+
+
+def main():
+    # Sweep raster size at the baseline search radius.
+    results = []
+    for mapsize in MAPSIZES:
+        benchmark(
+            size=int(mapsize**0.5),
+            search=BASE_SEARCH,
+            memory=BASE_MEMORY,
+            label=f"r.geomorphon_{int(mapsize / 1e6)}M",
+            results=results,
+        )
+    plot(results, "rastersize")
+
+    # Sweep search radius at the baseline raster size.
+    results = []
+    for search in SEARCHES:
+        benchmark(
+            size=int(BASE_MAPSIZE**0.5),
+            search=search,
+            memory=BASE_MEMORY,
+            label=f"r.geomorphon_search_{search}",
+            results=results,
+        )
+    plot(results, "search")
+
+
+def benchmark(size, search, memory, label, results):
+    reference = "benchmark_r_geomorphon_reference"
+    output = "benchmark_r_geomorphon"
+    generate_map(rows=size, cols=size, fname=reference)
+    module = Module(
+        "r.geomorphon",
+        elevation=reference,
+        forms=output,
+        search=search,
+        memory=memory,
+        run_=False,
+        overwrite=True,
+    )
+    results.append(
+        bm.benchmark_nprocs(
+            module,
+            label=label,
+            max_nprocs=8,
+            repeat=3,
+        )
+    )
+    Module(
+        "g.remove",
+        quiet=True,
+        flags="f",
+        type="raster",
+        pattern="benchmark_r_geomorphon*",
+    )
+
+
+def plot(results, sweep):
+    for metric in METRICS:
+        bm.nprocs_plot(
+            results,
+            filename=f"r_geomorphon_{sweep}_{metric}.svg",
+            title=f"r.geomorphon {sweep} {metric}",
+            metric=metric,
+        )
+
+
+def generate_map(rows, cols, fname):
+    Module("g.region", flags="p", n=rows, e=cols, res=1, w=0, s=0)
+    # Generate using r.random.surface if r.surf.fractal fails
+    try:
+        print("Generating reference map using r.surf.fractal...")
+        Module("r.surf.fractal", output=fname, overwrite=True)
+    except (CalledModuleError, GrassError):
+        print("r.surf.fractal fails, using r.random.surface instead...")
+        Module("r.random.surface", output=fname, overwrite=True)
+
+
+if __name__ == "__main__":
+    main()

raster/r.geomorphon/local_proto.h

@@ -32,7 +32,6 @@
 typedef struct {
     char elevname[150];
     RASTER_MAP_TYPE raster_type;
-    FCELL **elev;
     int fd; /* file descriptor */
 } MAPS;
 
@@ -64,11 +63,24 @@ typedef enum {
     CNT   /* counter */
 } FORMS;
 
+/* Invariant compute_forms inputs shared by every cell in a run. */
+struct geomorphon_config {
+    double search_dist, skip_dist, flat_dist, max_resolution;
+    int extended, oneoff;
+};
+
+/* Per-cell compute_forms outputs. */
+struct geomorphon_result {
+    PATTERN *pattern;
+    int pattern_size;
+    FORMS cur_form, orig_form;
+    double eff_search, eff_skip, eff_flat;
+};
+
 /* main */
 GLOBAL MAPS elevation;
 GLOBAL int ncols, row_radius_size, row_buffer_size;
 GLOBAL int skip_cells;
-GLOBAL double search_distance, flat_distance;
 GLOBAL double flat_threshold, flat_threshold_height;
 GLOBAL struct Cell_head window;
 
@@ -77,15 +89,16 @@ GLOBAL enum { ANGLEV1, ANGLEV2, ANGLEV2_DISTANCE } compmode;
 
 /* memory */
 int open_map(MAPS *rast);
-int shift_buffers(int row);
-int free_map(FCELL **map, int n);
+int load_strip(int fd, RASTER_MAP_TYPE rtype, void *tmp_buf, FCELL **rows,
+               int abs_first, int count);
 int write_form_cat_colors(char *raster);
 int write_contrast_colors(char *);
 const char *form_short_name(const FORMS);
 const char *form_long_name(const FORMS);
 
 /* pattern */
-int calc_pattern(PATTERN *pattern, int row, int cur_row, int col, const int);
+int calc_pattern(PATTERN *pattern, int row, int cur_row, int col, const int,
+                 double search_distance, double flat_distance, FCELL **rows);
 extern const char *direction_name[];
 
 /* geom */