r.mfilter: implement parallelization with OpenMP (#1708)

Co-authored-by: Anna Petrasova <kratochanna@gmail.com>

raster/r.mfilter/Makefile

@@ -2,8 +2,9 @@ MODULE_TOPDIR = ../..
 
 PGM = r.mfilter
 
-LIBES = $(ROWIOLIB) $(RASTERLIB) $(GISLIB)
+LIBES = $(ROWIOLIB) $(RASTERLIB) $(GISLIB) $(OMPLIB)
 DEPENDENCIES = $(ROWIODEP) $(RASTERDEP) $(GISDEP)
+EXTRA_CFLAGS = $(OMPCFLAGS)
 
 include $(MODULE_TOPDIR)/include/Make/Module.make
 

raster/r.mfilter/benchmark/benchmark_r_mfilter_nprocs.py

@@ -0,0 +1,73 @@
+"""Benchmarking of r.mfilter
+
+@author Aaron Saw Min Sern
+"""
+
+from grass.exceptions import CalledModuleError
+from grass.pygrass.modules import Module
+from grass.script import tempfile
+from subprocess import DEVNULL
+
+import grass.benchmark as bm
+
+
+def main():
+    results = []
+
+    # Users can add more or modify existing reference maps
+    benchmark(7071, "r.mfilter_50M", results)
+    benchmark(10000, "r.mfilter_100M", results)
+    benchmark(14142, "r.mfilter_200M", results)
+    benchmark(20000, "r.mfilter_400M", results)
+
+    bm.nprocs_plot(results)
+
+
+def benchmark(size, label, results):
+    reference = "r_mfilter_reference_map"
+    output = "benchmark_r_mfilter_nprocs"
+    filter = tempfile()
+    with open(filter, "w") as w:
+        w.write(
+            """MATRIX 9
+                   1 1 1 1 1 1 1 1 1
+                   1 2 1 2 1 2 1 2 1
+                   1 1 3 1 3 1 3 1 1
+                   1 2 1 4 1 4 1 2 1
+                   1 1 3 1 5 1 3 1 1
+                   1 2 1 4 1 4 1 2 1
+                   1 1 3 1 3 1 3 1 1
+                   1 2 1 2 1 2 1 2 1
+                   1 1 1 1 1 1 1 1 1
+                   DIVISOR 81
+                   TYPE    P"""
+        )
+
+    generate_map(rows=size, cols=size, fname=reference)
+    module = Module(
+        "r.mfilter",
+        input=reference,
+        output=output,
+        filter=filter,
+        run_=False,
+        stdout_=DEVNULL,
+        overwrite=True,
+    )
+    results.append(bm.benchmark_nprocs(module, label=label, max_nprocs=16, repeat=3))
+    Module("g.remove", quiet=True, flags="f", type="raster", name=reference)
+    Module("g.remove", quiet=True, flags="f", type="raster", name=output)
+
+
+def generate_map(rows, cols, fname):
+    Module("g.region", flags="p", s=0, n=rows, w=0, e=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)
+    except CalledModuleError:
+        print("r.surf.fractal fails, using r.random.surface instead...")
+        Module("r.random.surface", output=fname)
+
+
+if __name__ == "__main__":
+    main()

raster/r.mfilter/execute.c

@@ -1,12 +1,17 @@
+#if defined(_OPENMP)
+    #include <omp.h>
+#endif
+
 #include <unistd.h>
 #include <grass/rowio.h>
 #include <grass/raster.h>
 #include "glob.h"
 #include "filter.h"
 
-int execute_filter(ROWIO * r, int out, FILTER * filter, DCELL * cell)
+int execute_filter(ROWIO *r, int *out, FILTER *filter, DCELL **cell)
 {
     int i;
+    int t;
     int count;
     int size;
     int row, rcount;
@@ -14,12 +19,25 @@ int execute_filter(ROWIO * r, int out, FILTER * filter, DCELL * cell)
     int startx, starty;
     int dx, dy;
     int mid;
-    DCELL **bufs, **box, *cp;
+    int old_nprocs = 0;
+    DCELL ***bufs, ***box, *cp;
+
+    if (nprocs > 1 && filter->type == SEQUENTIAL) {
+        /* disable parallel temporarily */
+        old_nprocs = nprocs;
+        nprocs = 1;
+    }
+
 
     size = filter->size;
     mid = size / 2;
-    bufs = (DCELL **) G_malloc(size * sizeof(DCELL *));
-    box = (DCELL **) G_malloc(size * sizeof(DCELL *));
+    bufs = (DCELL ***) G_malloc(nprocs * sizeof(DCELL **));
+    box = (DCELL ***) G_malloc(nprocs * sizeof(DCELL **));
+
+    for (t = 0; t < nprocs; t++) {
+        bufs[t] = (DCELL **) G_malloc(size * sizeof(DCELL *));
+        box[t] = (DCELL **) G_malloc(size * sizeof(DCELL *));
+    }
 
     switch (filter->start) {
     case UR:
@@ -56,66 +74,94 @@ int execute_filter(ROWIO * r, int out, FILTER * filter, DCELL * cell)
     ccount = ncols - (size - 1);
 
     /* rewind output */
-    lseek(out, 0L, 0);
+    lseek(out[MASTER], 0L, SEEK_SET);
 
     /* copy border rows to output */
     row = starty;
     for (i = 0; i < mid; i++) {
-	cp = (DCELL *) Rowio_get(r, row);
-	write(out, cp, buflen);
-	row += dy;
+        cp = (DCELL *) Rowio_get(&r[MASTER], row);
+        if (write(out[MASTER], cp, buflen) < 0) 
+            G_fatal_error("Error writing temporary file");
+        row += dy;
     }
 
     /* for each row */
-    for (count = 0; count < rcount; count++) {
-	G_percent(count, rcount, 2);
+    int id = MASTER;
+    int start =  0;
+    int end = rcount;
+    int work = 0;
+    DCELL* cellp = cell[MASTER];
+
+    #pragma omp parallel firstprivate(starty, id, start, end, cellp) private(i, count, row, col, cp) if(nprocs > 1)
+    {
+    #if defined(_OPENMP)
+    if (nprocs > 1) {
+        id = omp_get_thread_num();
+        start = rcount * id/nprocs;
+        end = rcount * (id+1)/nprocs;
+        cellp = cell[id];
+        starty += start * dy;
+        lseek(out[id], (off_t) buflen * (mid + start), SEEK_SET);
+    }
+    #endif
+
+    for (count = start; count < end; count++) {
+	G_percent(work, rcount, 2);
 	row = starty;
 	starty += dy;
 	/* get "size" rows */
 	for (i = 0; i < size; i++) {
-	    bufs[i] = (DCELL *) Rowio_get(r, row);
-	    box[i] = bufs[i] + startx;
+	    bufs[id][i] = (DCELL *) Rowio_get(&r[id], row);
+	    box[id][i] = bufs[id][i] + startx;
 	    row += dy;
 	}
 	if (filter->type == SEQUENTIAL)
-	    cell = bufs[mid];
+	    cellp = bufs[id][mid];
 	/* copy border */
-	cp = cell;
+	cp = cellp;
 	for (i = 0; i < mid; i++)
-	    *cp++ = bufs[mid][i];
+	    *cp++ = bufs[id][mid][i];
 
 	/* filter row */
 	col = ccount;
 	while (col--) {
-	    if (null_only) {
-		if (Rast_is_d_null_value(&box[mid][mid]))
-		    *cp++ = apply_filter(filter, box);
-		else
-		    *cp++ = box[mid][mid];
-	    }
-	    else {
-		*cp++ = apply_filter(filter, box);
+	    if (null_only && !Rast_is_d_null_value(&box[id][mid][mid])) {
+		    *cp++ = box[id][mid][mid];
+	    } else {
+            *cp++ = apply_filter(filter, box[id]);
 	    }
 	    for (i = 0; i < size; i++)
-		box[i] += dx;
+		box[id][i] += dx;
 	}
 
 	/* copy border */
 	for (i = ncols - mid; i < ncols; i++)
-	    *cp++ = bufs[mid][i];
+	    *cp++ = bufs[id][mid][i];
 
 	/* write row */
-	write(out, cell, buflen);
+	write(out[id], cellp, buflen);
+    #pragma omp atomic update
+    work++;
     }
-    G_percent(count, rcount, 2);
+    }
+    G_percent(work, rcount, 2);
+    starty = rcount * dy;
+    lseek(out[MASTER], (off_t) buflen * (mid + rcount), SEEK_SET);
 
     /* copy border rows to output */
     row = starty + mid * dy;
     for (i = 0; i < mid; i++) {
-	cp = (DCELL *) Rowio_get(r, row);
-	write(out, cp, buflen);
+	cp = (DCELL *) Rowio_get(&r[MASTER], row);
+        if (write(out[MASTER], cp, buflen) < 0)
+            G_fatal_error("Error writing temporary file");
 	row += dy;
     }
 
+    if (old_nprocs != 0) {
+        /* restore parallel execution */
+        nprocs = old_nprocs;
+        old_nprocs = 0;
+    }
+
     return 0;
 }

raster/r.mfilter/filter.h

@@ -28,4 +28,4 @@ FILTER *get_filter(char *, int *, char *);
 int perform_filter(const char *, const char *, FILTER *, int, int);
 
 /* execute.c */
-int execute_filter(ROWIO *, int, FILTER *, DCELL *);
+int execute_filter(ROWIO *, int *, FILTER *, DCELL **);

raster/r.mfilter/glob.h

@@ -1,5 +1,6 @@
-
+extern const int MASTER;
 extern int nrows, ncols;
+extern int nprocs;
 extern int buflen;
 extern int direction;
 extern int null_only;

raster/r.mfilter/main.c

@@ -5,15 +5,19 @@
  * AUTHOR(S):    Michael Shapiro, CERL (original contributor)
  *               Roberto Flor <flor itc.it>, Markus Neteler <neteler itc.it>
  *               Glynn Clements <glynn gclements.plus.com>, Jachym Cepicky <jachym les-ejk.cz>,
- *               Jan-Oliver Wagner <jan intevation.de>
- * PURPOSE:      
- * COPYRIGHT:    (C) 1999-2006, 2010 by the GRASS Development Team
+ *               Jan-Oliver Wagner <jan intevation.de>,
+ *               Aaron Saw Min Sern
+ * PURPOSE:      Performs raster map matrix filter
+ * COPYRIGHT:    (C) 1999-2022 by the GRASS Development Team
  *
  *               This program is free software under the GNU General Public
  *               License (>=v2). Read the file COPYING that comes with GRASS
  *               for details.
  *
  *****************************************************************************/
+#if defined(_OPENMP)
+#include <omp.h>
+#endif
 
 #include <stdlib.h>
 #include <string.h>
@@ -25,7 +29,9 @@
 #include "filter.h"
 #include "glob.h"
 
+const int MASTER = 0;
 int nrows, ncols;
+int nprocs;
 int buflen;
 int direction;
 int null_only;
@@ -49,6 +55,7 @@ int main(int argc, char **argv)
     struct Option *opt3;
     struct Option *opt4;
     struct Option *opt5;
+    struct Option *opt6;
 
     G_gisinit(argv[0]);
 
@@ -57,6 +64,8 @@ int main(int argc, char **argv)
     G_add_keyword(_("algebra"));
     G_add_keyword(_("statistics"));
     G_add_keyword(_("filter"));
+    G_add_keyword(_("parallel"));
+
     module->description = _("Performs raster map matrix filter.");
 
     /* Define the different options */
@@ -85,6 +94,8 @@ int main(int argc, char **argv)
     opt5->required = NO;
     opt5->description = _("Output raster map title");
 
+    opt6 = G_define_standard_option(G_OPT_M_NPROCS);
+
     /* Define the different flags */
 
     /* this isn't implemented at all 
@@ -107,6 +118,20 @@ int main(int argc, char **argv)
     null_only = flag2->answer;
 
     sscanf(opt4->answer, "%d", &repeat);
+    sscanf(opt6->answer, "%d", &nprocs);
+    if (nprocs < 1)
+    {
+      G_fatal_error(_("<%d> is not valid number of threads."), nprocs);
+    }
+#if defined(_OPENMP)
+    omp_set_num_threads(nprocs);
+#else
+    if (nprocs != 1)
+        G_warning(_("GRASS is compiled without OpenMP support. Ignoring "
+                    "threads setting."));
+    nprocs = 1;
+#endif
+
     out_name = opt2->answer;
     filt_name = opt3->answer;
 

raster/r.mfilter/perform.c

@@ -10,18 +10,25 @@
 int perform_filter(const char *in_name, const char *out_name,
 		   FILTER * filter, int nfilters, int repeat)
 {
-    int in;
-    int out;
+    int *in;
+    int *out;
     int n;
     int pass;
-    ROWIO r;
+    ROWIO *r;
     char *tmp1, *tmp2;
     int count;
     int row;
-    DCELL *cell;
+    int t;
+    DCELL **cell;
 
+    cell = G_malloc(nprocs * sizeof(DCELL*));
+    for (t = 0; t < nprocs; t++) {
+        cell[t] = Rast_allocate_d_buf();
+    }
 
-    cell = Rast_allocate_d_buf();
+    in = G_malloc(nprocs * sizeof(int));
+    out = G_malloc(nprocs * sizeof(int));
+    r = G_malloc(nprocs * sizeof(ROWIO));
 
     count = 0;
     for (pass = 0; pass < repeat; pass++) {
@@ -30,58 +37,76 @@ int perform_filter(const char *in_name, const char *out_name,
 	    G_debug(1, "Filter %d", n + 1);
 
 	    if (count == 0) {
-		in = Rast_open_old(in_name, "");
-
-		G_debug(1, "Open raster map %s = %d", in_name, in);
-
-		close(creat(tmp1 = G_tempfile(), 0666));
-		out = open(tmp1, 2);
-		if (out < 0)
-		    G_fatal_error(_("Unable to create temporary file"));
+            for (t = 0; t < nprocs; t++) {
+                in[t] = Rast_open_old(in_name, "");
+
+                G_debug(1, "Open raster map %s = %d", in_name, in[t]);
+            }
+            close(creat(tmp1 = G_tempfile(), 0666));
+
+            for (t = 0; t < nprocs; t++) {
+                out[t] = open(tmp1, 2);
+                if (out[t] < 0) {
+                    G_fatal_error(_("Unable to create temporary file"));
+                }
+            }
 	    }
 	    else if (count == 1) {
 
-		G_debug(1, "Closing raster map");
-
-		Rast_close(in);
-		in = out;
-		close(creat(tmp2 = G_tempfile(), 0666));
-		out = open(tmp2, 2);
-		if (out < 0)
-		    G_fatal_error(_("Unable to create temporary file"));
+            G_debug(1, "Closing raster map");
+            for (t = 0; t < nprocs; t++) {
+                Rast_close(in[t]);
+                in[t] = out[t];
+            }
+            close(creat(tmp2 = G_tempfile(), 0666));
+
+            for (t = 0; t < nprocs; t++) {
+                out[t] = open(tmp2, 2);
+                if (out[t] < 0) {
+                    G_fatal_error(_("Unable to create temporary file"));
+                }
+            }
 	    }
 	    else {
-		int fd;
+            int fd;
 
-		G_debug(1, "Swap temp files");
+            G_debug(1, "Swap temp files");
 
-		fd = in;
-		in = out;
-		out = fd;
+            for (t = 0; t < nprocs; t++) {
+                fd = in[t];
+                in[t] = out[t];
+                out[t] = fd;
+            }
 	    }
 
-	    Rowio_setup(&r, in, filter[n].size, buflen,
-			count ? getrow : getmaprow, NULL);
+        for (t = 0; t < nprocs; t++) {
+            Rowio_setup(&r[t], in[t], filter[n].size, buflen,
+                count ? getrow : getmaprow, NULL);
+        }
 
-	    execute_filter(&r, out, &filter[n], cell);
+        execute_filter(r, out, &filter[n], cell);
 
-	    Rowio_release(&r);
+        for (t = 0; t < nprocs; t++) {
+            Rowio_release(&r[t]);
+        }
 	}
     }
 
     if (count == 1)
-	Rast_close(in);
+    for (t = 0; t < nprocs; t++)
+        Rast_close(in[t]);
     else if (count > 1)
-	close(in);
+    for (t = 0; t < nprocs; t++)
+        close(in[t]);
 
     /* copy final result to output raster map */
-    in = out;
-    out = Rast_open_fp_new(out_name);
+    in[MASTER] = out[MASTER];
+    out[MASTER] = Rast_open_fp_new(out_name);
 
     G_message(_("Writing raster map <%s>"), out_name);
     for (row = 0; row < nrows; row++) {
-	getrow(in, cell, row, buflen);
-	Rast_put_d_row(out, cell);
+        getrow(in[MASTER], cell[MASTER], row, buflen);
+        Rast_put_d_row(out[MASTER], cell[MASTER]);
     }
 
     /* remove the temporary files before closing so that the Rast_close()
@@ -91,7 +116,7 @@ int perform_filter(const char *in_name, const char *out_name,
 	unlink(tmp1);
     if (count > 1)
 	unlink(tmp2);
-    Rast_close(out);
+    Rast_close(out[MASTER]);
 
     return 0;
 }

raster/r.mfilter/r.mfilter.html

@@ -22,7 +22,7 @@ below, under FILTERS.
 The <b>repeat</b> parameter defines the number of times the <em>filter</em>
 is to be applied to the <em>input</em> data.
 
-<h2>FILTERS</h2>
+<h3>FILTERS</h3>
 
 The <em>filter</em> file is a normal UNIX ASCII file designed by the user.
 It has the following format:
@@ -83,7 +83,7 @@ For example, the following describes two filters:
 </dl>
 
 
-<h2>EXAMPLE FILTER FILE</h2>
+<h3>EXAMPLE FILTER FILE</h3>
 
 <div class="code"><pre>
       TITLE     3x3 average, non-null data only, followed by 5x5 average
@@ -104,7 +104,7 @@ For example, the following describes two filters:
      TYPE      P
 </pre></div>
 
-<h2>HOW THE FILTER WORKS</h2>
+<h3>HOW THE FILTER WORKS</h3>
 
 The filter process produces a new category value for each cell
 in the input raster map layer by multiplying the category values of the
@@ -126,6 +126,25 @@ then the next filter is applied to the intermediate result to
 produce another intermediate result;  and so on, until the
 final filter is applied.  Then the output cell is written.
 
+
+<h3>PERFORMANCE</h3>
+<p>By specifying the number of parallel processes with <b>nprocs</b> option,
+<em>r.mfilter</em> can run significantly faster, see benchmarks below.
+
+<div align="center" style="margin: 10px">
+     <img src="r_mfilter_benchmark_1.png" alt="benchmark for number of cells" border="0">
+     <img src="r_mfilter_benchmark_2.png" alt="benchmark for window size" border="0">
+     <br>
+     <i>Figure: Benchmark on the left shows execution time for different
+     number of cells for 9x9 matrix, benchmark on the right shows execution time
+     for 16 billion cells for different matrix sizes. (Intel Core i9-10940X CPU @ 3.30GHz x 28) </i>
+     </div>
+<p> Note that parallelization is implemented only for the parallel filter,
+not the sequential one.
+To take advantage of the parallelization, GRASS GIS
+needs to compiled with OpenMP enabled.
+
+
 <h2>NOTES</h2>
 
 If the resolution of the geographic region does not agree with the
@@ -139,14 +158,16 @@ is set properly.
 <a href="g.region.html">g.region</a>,
 <a href="r.clump.html">r.clump</a>,
 <a href="r.neighbors.html">r.neighbors</a>,
-<a href="r.resamp.filter.html">r.resamp.filter</a>
+<a href="r.resamp.filter.html">r.resamp.filter</a>,
+<a href="https://grasswiki.osgeo.org/wiki/Raster_Parallelization_with_OpenMP">Raster Parallelization with OpenMP</a>
 </em>
 
 <h2>AUTHOR</h2>
 
 Glynn Clements.
 Based upon r.mfilter, by Michael Shapiro,
-U.S.Army Construction Engineering Research Laboratory
+U.S.Army Construction Engineering Research Laboratory.<br>
+Aaron Saw Min Sern (OpenMP support).
 
 <!--
 <p>

raster/r.mfilter/testsuite/test_r_mfilter.py

@@ -0,0 +1,434 @@
+from tempfile import NamedTemporaryFile
+from grass.gunittest.case import TestCase
+from grass.gunittest.main import test
+
+
+class TestNeighbors(TestCase):
+    """
+
+    Used dataset: nc_spm_full_v2alphav2
+
+    Test cases:
+    test_sequential: Test output with sequential filter type
+    test_parallel: Test output with parallel filter type
+    test_sequential_null: Test output with sequential filter type
+        with null mode enabled
+    test_parallel_null: Test output with parallel filter type
+        with null mode enabled
+    """
+
+    test_results = {
+        "test_sequential": {
+            "n": 2025000,
+            "null_cells": 0,
+            "cells": 2025000,
+            "min": 55.5787925720215,
+            "max": 156.124557495117,
+            "range": 100.545764923096,
+            "mean": 110.374947784588,
+            "mean_of_abs": 110.374947784588,
+            "stddev": 20.2637480154117,
+            "variance": 410.6194836321,
+            "coeff_var": 18.3590102846157,
+            "sum": 223509269.26379,
+        },
+        "test_parallel": {
+            "n": 2025000,
+            "null_cells": 0,
+            "cells": 2025000,
+            "min": 55.5787925720215,
+            "max": 156.223892211914,
+            "range": 100.645099639893,
+            "mean": 110.375174079437,
+            "mean_of_abs": 110.375174079437,
+            "stddev": 20.2824544942723,
+            "variance": 411.377960312227,
+            "coeff_var": 18.3759207298509,
+            "sum": 223509727.51086,
+        },
+        "test_sequential_null": {
+            False: {
+                "n": 1789490,
+                "null_cells": 235510,
+                "cells": 2025000,
+                "min": 34301.08203125,
+                "max": 43599.45703125,
+                "range": 9298.375,
+                "mean": 39040.3073035648,
+                "mean_of_abs": 39040.3073035648,
+                "stddev": 338.861109540213,
+                "variance": 114826.851558824,
+                "coeff_var": 0.867977567147046,
+                "sum": 69862239516.6562,
+            },
+            True: {
+                "n": 1789490,
+                "null_cells": 235510,
+                "cells": 2025000,
+                "min": 34300,
+                "max": 43600,
+                "range": 9300,
+                "mean": 39041.4984470043,
+                "mean_of_abs": 39041.4984470043,
+                "stddev": 348.205753496913,
+                "variance": 121247.246768353,
+                "coeff_var": 0.891886242454486,
+                "sum": 69864371055.9297,
+            },
+        },
+        "test_parallel_null": {
+            False: {
+                "n": 46381,
+                "null_cells": 1978619,
+                "cells": 2025000,
+                "min": 34300,
+                "max": 43600,
+                "range": 9300,
+                "mean": 38988.3010371973,
+                "mean_of_abs": 38988.3010371973,
+                "stddev": 792.424493234645,
+                "variance": 627936.577478183,
+                "coeff_var": 2.03246736111589,
+                "sum": 1808316390.40625,
+            },
+            True: {
+                "n": 46381,
+                "null_cells": 1978619,
+                "cells": 2025000,
+                "min": 34300,
+                "max": 43600,
+                "range": 9300,
+                "mean": 38988.2992790859,
+                "mean_of_abs": 38988.2992790859,
+                "stddev": 798.805978312437,
+                "variance": 638090.990987689,
+                "coeff_var": 2.04883514562774,
+                "sum": 1808316308.86328,
+            },
+        },
+        "test_multiple_filters": {
+            "n": 2025000,
+            "null_cells": 0,
+            "cells": 2025000,
+            "min": 55.5787925720215,
+            "max": 155.957977294922,
+            "range": 100.3791847229,
+            "mean": 110.374591878668,
+            "mean_of_abs": 110.374591878668,
+            "stddev": 20.235686035401,
+            "variance": 409.482989323322,
+            "coeff_var": 18.3336451722925,
+            "sum": 223508548.554302,
+        },
+        "test_repeated_filters": {
+            "n": 2025000,
+            "null_cells": 0,
+            "cells": 2025000,
+            "min": 55.5787925720215,
+            "max": 155.465209960938,
+            "range": 99.886417388916,
+            "mean": 110.373036950725,
+            "mean_of_abs": 110.373036950725,
+            "stddev": 20.1413729988748,
+            "variance": 405.674906279802,
+            "coeff_var": 18.2484541110042,
+            "sum": 223505399.825218,
+        },
+    }
+
+    filter_options = {
+        "sequential": b"""MATRIX 5
+                          1 1 1 1 1
+                          1 1 1 1 1
+                          1 1 1 1 1
+                          1 1 1 1 1
+                          1 1 1 1 1
+                          DIVISOR 0
+                          TYPE    S""",
+        "parallel": b"""MATRIX 5
+                        1 1 1 1 1
+                        1 1 1 1 1
+                        1 1 1 1 1
+                        1 1 1 1 1
+                        1 1 1 1 1
+                        DIVISOR 0
+                        TYPE    P""",
+        "mix": b"""MATRIX 5
+                   1 1 1 1 1
+                   1 1 1 1 1
+                   1 1 1 1 1
+                   1 1 1 1 1
+                   1 1 1 1 1
+                   DIVISOR 0
+                   TYPE    P
+
+                   MATRIX 5
+                   1 1 1 1 1
+                   1 1 1 1 1
+                   1 1 1 1 1
+                   1 1 1 1 1
+                   1 1 1 1 1
+                   DIVISOR 0
+                   TYPE    S
+
+                   MATRIX 3
+                   1 1 1 
+                   1 1 1 
+                   1 1 1 
+                   DIVISOR 9
+                   TYPE    P""",
+    }
+
+    # TODO: replace by unified handing of maps
+    to_remove = []
+
+    def create_filter(self, options):
+        """Create a temporary filter file with the given name and options."""
+        f = NamedTemporaryFile()
+        f.write(options)
+        f.flush()
+        return f
+
+    @classmethod
+    def setUpClass(cls):
+        cls.use_temp_region()
+        cls.runModule("g.region", raster="elevation")
+
+    @classmethod
+    def tearDownClass(cls):
+        cls.del_temp_region()
+        if cls.to_remove:
+            cls.runModule(
+                "g.remove", flags="f", type="raster", name=",".join(cls.to_remove)
+            )
+
+    def test_sequential(self):
+        """Test output with sequential filter type."""
+        test_case = "test_sequential"
+        output = "{}_raster".format(test_case)
+        output_threaded = "{}_threaded_raster".format(test_case)
+        self.to_remove.extend([output, output_threaded])
+
+        filter = self.create_filter(self.filter_options["sequential"])
+        self.assertModule(
+            "r.mfilter",
+            input="elevation",
+            output=output,
+            filter=filter.name,
+        )
+        self.assertModule(
+            "r.mfilter",
+            input="elevation",
+            output=output_threaded,
+            filter=filter.name,
+            nprocs=4,
+        )
+        filter.close()
+        self.assertRasterFitsUnivar(
+            raster=output,
+            reference=self.test_results[test_case],
+            precision=1e-5,
+        )
+        self.assertRasterFitsUnivar(
+            raster=output_threaded,
+            reference=self.test_results[test_case],
+            precision=1e-5,
+        )
+
+    def test_parallel(self):
+        """Test output with parallel filter type."""
+        test_case = "test_parallel"
+        output = "{}_raster".format(test_case)
+        output_threaded = "{}_threaded_raster".format(test_case)
+        self.to_remove.extend([output, output_threaded])
+
+        filter = self.create_filter(self.filter_options["parallel"])
+        self.assertModule(
+            "r.mfilter",
+            input="elevation",
+            output=output,
+            filter=filter.name,
+        )
+        self.assertModule(
+            "r.mfilter",
+            input="elevation",
+            output=output_threaded,
+            filter=filter.name,
+            nprocs=4,
+        )
+        filter.close()
+        self.assertRasterFitsUnivar(
+            raster=output,
+            reference=self.test_results[test_case],
+            precision=1e-5,
+        )
+        self.assertRasterFitsUnivar(
+            raster=output_threaded,
+            reference=self.test_results[test_case],
+            precision=1e-5,
+        )
+
+    def test_sequential_null(self):
+        """Test output with sequential filter type with null mode enabled."""
+        test_case = "test_sequential_null"
+        output = "{}_raster".format(test_case)
+        output_z = "{}_z_raster".format(test_case)
+        self.to_remove.extend([output, output_z])
+
+        filter = self.create_filter(self.filter_options["sequential"])
+        self.assertModule(
+            "r.mfilter",
+            input="lakes",
+            output=output,
+            filter=filter.name,
+        )
+        self.assertModule(
+            "r.mfilter",
+            input="lakes",
+            output=output_z,
+            filter=filter.name,
+            flags="z",
+        )
+        filter.close()
+        self.assertRasterFitsUnivar(
+            raster=output,
+            reference=self.test_results[test_case][False],
+            precision=1e-5,
+        )
+        self.assertRasterFitsUnivar(
+            raster=output_z,
+            reference=self.test_results[test_case][True],
+            precision=1e-5,
+        )
+
+    def test_parallel_null(self):
+        """Test output with parallel filter type with null mode enabled."""
+        test_case = "test_parallel_null"
+        output = "{}_raster".format(test_case)
+        output_threaded = "{}_threaded_raster".format(test_case)
+        output_z = "{}_z_raster".format(test_case)
+        output_z_threaded = "{}_z_threaded_raster".format(test_case)
+        self.to_remove.extend([output, output_threaded, output_z, output_z_threaded])
+
+        filter = self.create_filter(self.filter_options["parallel"])
+        self.assertModule(
+            "r.mfilter",
+            input="lakes",
+            output=output,
+            filter=filter.name,
+        )
+        self.assertModule(
+            "r.mfilter",
+            input="lakes",
+            output=output_threaded,
+            filter=filter.name,
+            nprocs=4,
+        )
+        self.assertModule(
+            "r.mfilter",
+            input="lakes",
+            output=output_z,
+            filter=filter.name,
+            flags="z",
+        )
+        self.assertModule(
+            "r.mfilter",
+            input="lakes",
+            output=output_z_threaded,
+            filter=filter.name,
+            flags="z",
+            nprocs=4,
+        )
+        filter.close()
+        self.assertRasterFitsUnivar(
+            raster=output,
+            reference=self.test_results[test_case][False],
+            precision=1e-5,
+        )
+        self.assertRasterFitsUnivar(
+            raster=output_threaded,
+            reference=self.test_results[test_case][False],
+            precision=1e-5,
+        )
+        self.assertRasterFitsUnivar(
+            raster=output_z,
+            reference=self.test_results[test_case][True],
+            precision=1e-5,
+        )
+        self.assertRasterFitsUnivar(
+            raster=output_z_threaded,
+            reference=self.test_results[test_case][True],
+            precision=1e-5,
+        )
+
+    def test_multiple_filters(self):
+        """Test output with multiple filters."""
+        test_case = "test_multiple_filters"
+        output = "{}_raster".format(test_case)
+        output_threaded = "{}_threaded_raster".format(test_case)
+        self.to_remove.extend([output, output_threaded])
+
+        filter = self.create_filter(self.filter_options["mix"])
+        self.assertModule(
+            "r.mfilter",
+            input="elevation",
+            output=output,
+            filter=filter.name,
+        )
+        self.assertModule(
+            "r.mfilter",
+            input="elevation",
+            output=output_threaded,
+            filter=filter.name,
+            nprocs=4,
+        )
+        filter.close()
+        self.assertRasterFitsUnivar(
+            raster=output,
+            reference=self.test_results[test_case],
+            precision=1e-5,
+        )
+        self.assertRasterFitsUnivar(
+            raster=output_threaded,
+            reference=self.test_results[test_case],
+            precision=1e-5,
+        )
+
+    def test_repeated_filters(self):
+        """Test output with repeated filters."""
+        test_case = "test_repeated_filters"
+        output = "{}_raster".format(test_case)
+        output_threaded = "{}_threaded_raster".format(test_case)
+        self.to_remove.extend([output, output_threaded])
+
+        filter = self.create_filter(self.filter_options["mix"])
+        self.assertModule(
+            "r.mfilter",
+            input="elevation",
+            output=output,
+            filter=filter.name,
+            repeat=3,
+        )
+        self.assertModule(
+            "r.mfilter",
+            input="elevation",
+            output=output_threaded,
+            filter=filter.name,
+            repeat=3,
+            nprocs=4,
+        )
+        filter.close()
+        self.assertRasterFitsUnivar(
+            raster=output,
+            reference=self.test_results[test_case],
+            precision=1e-5,
+        )
+        self.assertRasterFitsUnivar(
+            raster=output_threaded,
+            reference=self.test_results[test_case],
+            precision=1e-5,
+        )
+
+
+if __name__ == "__main__":
+    test()