r.texture: support parallel computing by OpenMP

raster/r.texture/Makefile

@@ -2,8 +2,10 @@ MODULE_TOPDIR = ../..
 
 PGM = r.texture
 
-LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB)
+LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB) $(OPENMP_LIBPATH) $(OPENMP_LIB)
 DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+EXTRA_CFLAGS = $(OPENMP_CFLAGS)
+EXTRA_INC = $(OPENMP_INCPATH)
 
 include $(MODULE_TOPDIR)/include/Make/Module.make
 

raster/r.texture/benchmark/benchmark_rtexture.py

@@ -0,0 +1,60 @@
+"""Benchmarking of r.texture
+
+@author Aaron Saw Min Sern
+        Chung-Yuan Liang
+"""
+
+from grass.exceptions import CalledModuleError
+from grass.pygrass.modules import Module
+from subprocess import DEVNULL
+
+import grass.benchmark as bm
+
+
+def main():
+    results = []
+    mapsizes = [1e6, 2e6, 4e6, 8e6]
+    metrics = ["time", "speedup", "efficiency"]
+
+    for mapsize in mapsizes:
+        benchmark(int(mapsize**0.5), f"r.texture_{int(mapsize/1e6)}M", results)
+
+    for metric in metrics:
+        bm.nprocs_plot(
+            results,
+            title=f"r.texture -a {metric}",
+            metric=metric,
+        )
+
+
+def benchmark(size, label, results):
+    reference = "r_texture_reference_map"
+    output = "benchmark_r_texture"
+    generate_map(rows=size, cols=size, fname=reference)
+    module = Module(
+        "r.texture",
+        input=reference,
+        output=output,
+        a=True,
+        run_=False,
+        stdout_=DEVNULL,
+        overwrite=True,
+    )
+    results.append(bm.benchmark_nprocs(module, label=label, max_nprocs=8, repeat=3))
+    Module("g.remove", quiet=True, flags="f", type="raster", name=reference)
+    Module("g.remove", quiet=True, flags="f", type="raster", pattern=f"{output}_*")
+
+
+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.texture/benchmark/benchmark_rtexture_window.py

@@ -0,0 +1,62 @@
+"""Benchmarking of r.texture
+
+@author Aaron Saw Min Sern
+        Chung-Yuan Liang
+"""
+
+from grass.exceptions import CalledModuleError
+from grass.pygrass.modules import Module
+from subprocess import DEVNULL
+
+import grass.benchmark as bm
+
+
+def main():
+    results = []
+    mapsize = 1e6
+    window_sizes = [3, 9, 15, 27]
+    metrics = ["time", "speedup", "efficiency"]
+
+    for window in window_sizes:
+        benchmark(int(mapsize**0.5), window, f"r.texture_{window}x{window}", results)
+
+    for metric in metrics:
+        bm.nprocs_plot(
+            results,
+            title=f"r.texture -a {metric}",
+            metric=metric,
+        )
+
+
+def benchmark(size, window, label, results):
+    reference = "r_texture_reference_map"
+    output = "benchmark_r_texture"
+    generate_map(rows=size, cols=size, fname=reference)
+    module = Module(
+        "r.texture",
+        input=reference,
+        output=output,
+        size=window,
+        a=True,
+        run_=False,
+        stdout_=DEVNULL,
+        overwrite=True,
+    )
+    results.append(bm.benchmark_nprocs(module, label=label, max_nprocs=8, repeat=3))
+    Module("g.remove", quiet=True, flags="f", type="raster", name=reference)
+    Module("g.remove", quiet=True, flags="f", type="raster", pattern=f"{output}_*")
+
+
+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.texture/execute.c

@@ -1,3 +1,7 @@
+#if defined(_OPENMP)
+#include <omp.h>
+#endif
+
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
@@ -17,8 +21,7 @@
  ***************************************************************************************************/
 int execute_texture(CELL **data, struct dimensions *dim,
                     struct menu *measure_menu, int *measure_idx,
-                    struct output_setting *out_set)
-
+                    struct output_setting *out_set, int threads)
 {
     int size = dim->size;
     int dist = dim->dist;
@@ -27,23 +30,34 @@ int execute_texture(CELL **data, struct dimensions *dim,
     int n_outputs = dim->n_outputs;
     int n_measures = dim->n_measures;
     int *outfd = out_set->outfd;
+
     RASTER_MAP_TYPE out_data_type = out_set->out_data_type;
     struct Flag *flag_null = out_set->flag_null;
     struct Flag *flag_ind = out_set->flag_ind;
 
     int offset = size / 2;
     int i, j, row, col, first_row, first_col, last_row, last_col;
+    int trow;
     int have_px, have_py, have_pxpys, have_pxpyd;
-    FCELL **fbuf;
-    FCELL measure; /* Containing measure done */
-    struct matvec *mv;
 
-    fbuf = G_malloc(n_outputs * sizeof(FCELL *));
-    for (i = 0; i < n_outputs; i++)
-        fbuf[i] = Rast_allocate_buf(out_data_type);
+    FCELL ***fbuf_threads; /* Buffer for each thread */
+    FCELL measure;         /* Containing measure done */
+    struct matvec **mvs;   /* matrices and vectors for each thread */
+
+    /* allocate memory*/
+    /* fbuf_threads[0] is used for writing out when program is not in loops*/
+    fbuf_threads = (FCELL ***)G_malloc(sizeof(FCELL **) * threads);
+    for (i = 0; i < threads; i++) {
+        fbuf_threads[i] = (FCELL **)G_malloc(n_outputs * sizeof(FCELL *));
+        for (j = 0; j < n_outputs; j++)
+            fbuf_threads[i][j] = Rast_allocate_buf(out_data_type);
+    }
 
-    mv = G_malloc(sizeof(struct matvec));
-    alloc_vars(size, mv);
+    mvs = (struct matvec **)G_malloc(sizeof(struct matvec *) * threads);
+    for (i = 0; i < threads; i++) {
+        mvs[i] = G_malloc(sizeof(struct matvec));
+        alloc_vars(size, mvs[i]);
+    }
 
     /* variables needed */
     if (measure_menu[2].useme || measure_menu[11].useme ||
@@ -76,11 +90,11 @@ int execute_texture(CELL **data, struct dimensions *dim,
         last_col = ncols;
     }
 
-    Rast_set_f_null_value(fbuf[0], ncols);
+    Rast_set_f_null_value(fbuf_threads[0][0], ncols);
 
     for (row = 0; row < first_row; row++) {
         for (i = 0; i < n_outputs; i++) {
-            Rast_put_row(outfd[i], fbuf[0], out_data_type);
+            Rast_put_row(outfd[i], fbuf_threads[0][0], out_data_type);
         }
     }
     if (n_measures > 1)
@@ -90,62 +104,80 @@ int execute_texture(CELL **data, struct dimensions *dim,
     else
         G_message(_("Calculating %s..."), measure_menu[measure_idx[0]].desc);
 
-    for (row = first_row; row < last_row; row++) {
-        G_percent(row, nrows, 2);
-        for (i = 0; i < n_outputs; i++)
-            Rast_set_f_null_value(fbuf[i], ncols);
-
-        /*process the data */
-        for (col = first_col; col < last_col; col++) {
-            if (!set_vars(mv, data, row, col, size, offset, dist,
-                          flag_null->answer)) {
-                for (i = 0; i < n_outputs; i++)
-                    Rast_set_f_null_value(&(fbuf[i][col]), 1);
-                continue;
-            }
-            /* for all angles (0, 45, 90, 135) */
-            for (i = 0; i < 4; i++) {
-                set_angle_vars(mv, i, have_px, have_py, have_pxpys, have_pxpyd);
-                /* for all requested textural measures */
-                for (j = 0; j < n_measures; j++) {
-
-                    measure =
-                        (FCELL)h_measure(measure_menu[measure_idx[j]].idx, mv);
-
-                    if (flag_ind->answer) {
-                        /* output for each angle separately */
-                        fbuf[j * 4 + i][col] = measure;
-                    }
-                    else {
-                        /* use average over all angles for each measure */
-                        if (i == 0)
-                            fbuf[j][col] = measure;
-                        else if (i < 3)
-                            fbuf[j][col] += measure;
-                        else
-                            fbuf[j][col] = (fbuf[j][col] + measure) / 4.0;
+#pragma omp parallel private(row, col, i, j, measure, trow) default(shared)
+    {
+#pragma omp for schedule(static, 1) ordered
+        for (row = first_row; row < last_row; row++) {
+            trow = row % threads; /* Obtain thread row id */
+            G_percent(row, nrows, 2);
+
+            /* initialize the output row */
+            for (i = 0; i < n_outputs; i++)
+                Rast_set_f_null_value(fbuf_threads[trow][i], ncols);
+
+            /*process the data */
+            for (col = first_col; col < last_col; col++) {
+                if (!set_vars(mvs[trow], data, row, col, size, offset, dist,
+                              flag_null->answer)) {
+                    for (i = 0; i < n_outputs; i++)
+                        Rast_set_f_null_value(&(fbuf_threads[trow][i][col]), 1);
+                    continue;
+                }
+                /* for all angles (0, 45, 90, 135) */
+                for (i = 0; i < 4; i++) {
+                    set_angle_vars(mvs[trow], i, have_px, have_py, have_pxpys,
+                                   have_pxpyd);
+                    /* for all requested textural measures */
+                    for (j = 0; j < n_measures; j++) {
+                        measure = (FCELL)h_measure(
+                            measure_menu[measure_idx[j]].idx, mvs[trow]);
+                        if (flag_ind->answer) {
+                            /* output for each angle separately */
+                            fbuf_threads[trow][j * 4 + i][col] = measure;
+                        }
+                        else {
+                            /* use average over all angles for each measure */
+                            if (i == 0)
+                                fbuf_threads[trow][j][col] = measure;
+                            else if (i < 3)
+                                fbuf_threads[trow][j][col] += measure;
+                            else
+                                fbuf_threads[trow][j][col] =
+                                    (fbuf_threads[trow][j][col] + measure) /
+                                    4.0;
+                        }
                     }
                 }
             }
+#pragma omp ordered
+            {
+                for (i = 0; i < n_outputs; i++)
+                    Rast_put_row(outfd[i], fbuf_threads[trow][i],
+                                 out_data_type);
+            }
         }
+    } /* end of parallel section */
 
-        for (i = 0; i < n_outputs; i++)
-            Rast_put_row(outfd[i], fbuf[i], out_data_type);
-    }
-
-    Rast_set_f_null_value(fbuf[0], ncols);
+    Rast_set_f_null_value(fbuf_threads[0][0], ncols);
     for (row = last_row; row < nrows; row++) {
         for (i = 0; i < n_outputs; i++) {
-            Rast_put_row(outfd[i], fbuf[0], out_data_type);
+            Rast_put_row(outfd[i], fbuf_threads[0][0], out_data_type);
         }
     }
     G_percent(nrows, nrows, 1);
 
-    for (i = 0; i < n_outputs; i++)
-        G_free(fbuf[i]);
-    G_free(fbuf);
-    dealloc_vars(mv);
-    G_free(mv);
+    for (i = 0; i < threads; i++) {
+        for (j = 0; j < n_outputs; j++)
+            G_free(fbuf_threads[i][j]);
+        G_free(fbuf_threads[i]);
+    }
+    G_free(fbuf_threads);
+
+    for (i = 0; i < threads; i++) {
+        dealloc_vars(mvs[i]);
+        G_free(mvs[i]);
+    }
+    G_free(mvs);
 
     return 0;
 }

raster/r.texture/execute.h

@@ -20,4 +20,4 @@ typedef struct output_setting {
 
 int execute_texture(CELL **data, struct dimensions *dim,
                     struct menu *measure_menu, int *measure_idx,
-                    struct output_setting *out_set);
+                    struct output_setting *out_set, int threads);