r.horizon: Support parallel computing for the raster mode by OpenMP

raster/r.horizon/Makefile

@@ -2,9 +2,10 @@ MODULE_TOPDIR = ../..
 
 PGM = r.horizon
 
-LIBES = $(PARSONLIB) $(GPROJLIB) $(RASTERLIB) $(GISLIB) $(MATHLIB) $(PROJLIB)
+LIBES = $(PARSONLIB) $(GPROJLIB) $(RASTERLIB) $(GISLIB) $(MATHLIB) $(PROJLIB) $(OPENMP_LIBPATH) $(OPENMP_LIB)
 DEPENDENCIES = $(GPROJDEP) $(RASTERDEP) $(GISDEP)
-EXTRA_INC = $(PROJINC) $(GDALCFLAGS)
+EXTRA_CFLAGS = $(OPENMP_CFLAGS)
+EXTRA_INC = $(PROJINC) $(GDALCFLAGS) $(OPENMP_INCPATH)
 
 include $(MODULE_TOPDIR)/include/Make/Module.make
 

raster/r.horizon/benchmark/benchmark_rhorizon.py

@@ -0,0 +1,67 @@
+"""Benchmarking of r.horizon
+point mode, one direction
+
+@author Chung-Yuan Liang, 2024
+"""
+
+from grass.exceptions import CalledModuleError
+from grass.pygrass.modules import Module
+
+import grass.benchmark as bm
+
+
+def main():
+    results = []
+    metrics = ["time", "speedup", "efficiency"]
+    mapsizes = [1e6, 2e6, 4e6, 8e6]
+
+    # run benchmarks
+    for mapsize in mapsizes:
+        benchmark(
+            size=int(mapsize**0.5),
+            step=0,
+            label=f"r.horizon_{int(mapsize/1e6)}M",
+            results=results,
+        )
+
+    # plot results
+    for metric in metrics:
+        bm.nprocs_plot(
+            results,
+            title=f"r.horizon raster mode {metric}",
+            metric=metric,
+        )
+
+
+def benchmark(size, step, label, results):
+    reference = "benchmark_r_horizon_reference_map"
+    output = "benchmark_r_horizon"
+
+    generate_map(rows=size, cols=size, fname=reference)
+    module = Module(
+        "r.horizon",
+        elevation=reference,
+        output=output,
+        direction=0,
+        step=step,
+    )
+
+    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_horizon*"
+    )
+
+
+def generate_map(rows, cols, fname):
+    Module("g.region", flags="p", rows=rows, cols=cols, res=1)
+    # 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:
+        print("r.surf.fractal fails, using r.random.surface instead...")
+        Module("r.random.surface", output=fname, overwrite=True)
+
+
+if __name__ == "__main__":
+    main()

raster/r.horizon/main.c

@@ -31,6 +31,9 @@ Program was refactored by Anna Petrasova to remove most global variables.
  *   Free Software Foundation, Inc.,
  *   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  */
+#if defined(_OPENMP)
+#include <omp.h>
+#endif
 
 #include <stdio.h>
 #include <stdlib.h>
@@ -163,7 +166,7 @@ int main(int argc, char *argv[])
     struct {
         struct Option *elevin, *dist, *coord, *direction, *horizon, *step,
             *start, *end, *bufferzone, *e_buff, *w_buff, *n_buff, *s_buff,
-            *maxdistance, *format, *output;
+            *maxdistance, *format, *output, *nprocs;
     } parm;
 
     struct {
@@ -175,6 +178,7 @@ int main(int argc, char *argv[])
     G_add_keyword(_("raster"));
     G_add_keyword(_("solar"));
     G_add_keyword(_("sun position"));
+    G_add_keyword(_("parallel"));
     module->label =
         _("Computes horizon angle height from a digital elevation model.");
     module->description =
@@ -309,6 +313,8 @@ int main(int argc, char *argv[])
         _("Name of file for output (use output=- for stdout)");
     parm.output->guisection = _("Point mode");
 
+    parm.nprocs = G_define_standard_option(G_OPT_M_NPROCS);
+
     flag.horizonDistance = G_define_flag();
     flag.horizonDistance->key = 'l';
     flag.horizonDistance->description =
@@ -328,6 +334,29 @@ int main(int argc, char *argv[])
     if (G_parser(argc, argv))
         exit(EXIT_FAILURE);
 
+    int threads = atoi(parm.nprocs->answer);
+#if defined(_OPENMP)
+    /* Set the number of threads */
+    if (threads < 0) {
+        threads += omp_get_num_procs() + 1;
+    }
+    else if (threads == 0) {
+        G_warning(_("At least one thread should be used, the number of threads "
+                    "will be set on <%d>"),
+                  1);
+        threads = 1;
+    }
+    G_message(_("Using %d threads for parallel computing."), threads);
+    omp_set_num_threads(threads);
+#else
+    /* The number of threads is always 1 when OpenMP is not available */
+    if (threads != 1) {
+        G_warning(_("GRASS GIS is not compiled with OpenMP support, parallel "
+                    "computation is disabled."));
+    }
+    threads = 1;
+#endif
+
     struct Cell_head cellhd;
     struct Cell_head new_cellhd;
     G_get_set_window(&cellhd);
@@ -857,6 +886,7 @@ void calculate_point_mode(const Settings *settings, const Geometry *geometry,
         horizons = json_value_get_array(horizons_value);
         break;
     }
+
     for (int i = 0; i < printCount; i++) {
         JSON_Value *value;
         JSON_Object *object;
@@ -1179,6 +1209,7 @@ void calculate_raster_mode(const Settings *settings, const Geometry *geometry,
             _("Calculating map %01d of %01d (angle %.2f, raster map <%s>)"),
             (k + 1), arrayNumInt, angle_deg, shad_filename);
 
+#pragma omp parallel for schedule(static, 1) default(shared)
         for (int j = hor_row_start; j < hor_row_end; j++) {
             G_percent(j - hor_row_start, hor_numrows - 1, 2);
             for (int i = hor_col_start; i < hor_col_end; i++) {
@@ -1219,12 +1250,11 @@ void calculate_raster_mode(const Settings *settings, const Geometry *geometry,
                     if (settings->degreeOutput) {
                         shadow_angle *= rad2deg;
                     }
-
                     horizon_raster[j - buffer_s][i - buffer_w] = shadow_angle;
 
                 } /* undefs */
-            }
-        }
+            }     /* end of loop over columns */
+        }         /* end of parallel section */
 
         G_debug(1, "OUTGR() starts...");
         OUTGR(settings, shad_filename, cellhd);
@@ -1265,4 +1295,9 @@ void calculate_raster_mode(const Settings *settings, const Geometry *geometry,
         Rast_write_history(shad_filename, &history);
         G_free(shad_filename);
     }
+
+    /* free memory */
+    for (int l = 0; l < hor_numrows; l++)
+        G_free(horizon_raster[l]);
+    G_free(horizon_raster);
 }

raster/r.horizon/r.horizon.html

@@ -145,6 +145,49 @@ <h2>METHOD</h2>
 <p> All horizon values are positive (or zero). While negative values are
 in theory possible, <b>r.horizon</b> currently does not support them.
 
+<h2>NOTES</h2>
+The parallelization is only implemented for the raster mode.
+Although the point mode has been parallelized,
+the run time is much shorter than the raster mode.
+
+<h3>Performance</h3>
+To enable parallel processing, the user can specify the number of threads
+to be used with the <b>nprocs</b> parameter (default 1).
+Figures below show benchmark results running on
+Intel® Core™ i5-13600K CPU @ 3.5GHz.
+See benchmark scripts
+in the source code for more details.
+
+<p>
+The iterations of the algorithm used in <b>r.horizon</b>
+depends on the topography. As a result, the benchmark results may vary
+depending on the topography of the study area.
+The benchmark results below are examples based on
+the North Carolina sample dataset.
+
+<div align="center" style="margin: 10px">
+  <img src="rhorizon_raster_time.png"
+  alt="time for r.horizon with different map sizes" border="0">
+  <br>
+  <i>Figure 1: Benchmark shows execution time for different
+    number of cells (1M, 2M, 4M, and 8M) on the raster mode.
+</div>
+<div align="center" style="margin: 10px">
+  <img src="rhorizon_raster_speedup.png"
+  alt="speedup for r.horizon with different map sizes" border="0">
+  <br>
+  <i>Figure 2: Benchmark shows speedup for different
+    numbers of cells (1M, 2M, 4M, and 8M) on the raster mode.
+</div>
+<div align="center" style="margin: 10px">
+  <img src="rhorizon_raster_efficiency.png"
+  alt="efficiency for r.horizon with different map sizes" border="0">
+  <br>
+  <i>Figure 3: Benchmark shows efficiency for different
+    numbers of cells (1M, 2M, 4M, and 8M) on the raster mode.
+</div>
+
+
 <h2>EXAMPLES</h2>
 
 The examples are intended for the North Carolina sample dataset.