r.neighbors: add option for exponential weighting (#597)

* add exponent weight

* add weights function option

* add weights function option

* avoid function duplication

* make choice none all lowercase

* Dont check for non-zero result from strcmp in if statement

* add weigth functions

* use function pointer to weigth functions

* apply grass_indent.sh

* update tests with the new options

* adjust documentation

* r.neighbors: Make weighting_function parameter descriptions translatable

Co-authored-by: Māris Nartišs <maris.gis@gmail.com>

raster/r.neighbors/local_proto.h

@@ -22,4 +22,6 @@ extern int null_cats(const char *);
 
 /* read_weights.c */
 extern void read_weights(const char *);
-extern void gaussian_weights(double);
+extern double gaussian(double, double);
+extern double exponential(double, double);
+extern void compute_weights(const char *, double);

raster/r.neighbors/main.c

@@ -42,6 +42,12 @@ struct menu
     char *text;			/* menu display - full description */
 };
 
+struct weight_functions
+{
+    char *name;			/* name  of the weight type */
+    char *text;			/* weight types display - full description */
+};
+
 enum out_type {
     T_FLOAT	= 1,
     T_INT	= 2,
@@ -148,7 +154,8 @@ int main(int argc, char *argv[])
 	struct Option *method, *size;
 	struct Option *title;
 	struct Option *weight;
-	struct Option *gauss;
+	struct Option *weighting_function;
+	struct Option *weighting_factor;
 	struct Option *quantile;
     } parm;
     struct
@@ -187,6 +194,14 @@ int main(int argc, char *argv[])
     parm.output = G_define_standard_option(G_OPT_R_OUTPUT);
     parm.output->multiple = YES;
 
+    parm.size = G_define_option();
+    parm.size->key = "size";
+    parm.size->type = TYPE_INTEGER;
+    parm.size->required = NO;
+    parm.size->description = _("Neighborhood size");
+    parm.size->answer = "3";
+    parm.size->guisection = _("Neighborhood");
+
     parm.method = G_define_option();
     parm.method->key = "method";
     parm.method->type = TYPE_STRING;
@@ -205,32 +220,36 @@ int main(int argc, char *argv[])
     parm.method->multiple = YES;
     parm.method->guisection = _("Neighborhood");
 
-    parm.size = G_define_option();
-    parm.size->key = "size";
-    parm.size->type = TYPE_INTEGER;
-    parm.size->required = NO;
-    parm.size->description = _("Neighborhood size");
-    parm.size->answer = "3";
-    parm.size->guisection = _("Neighborhood");
-
-    parm.title = G_define_option();
-    parm.title->key = "title";
-    parm.title->key_desc = "phrase";
-    parm.title->type = TYPE_STRING;
-    parm.title->required = NO;
-    parm.title->description = _("Title for output raster map");
+    parm.weighting_function = G_define_option();
+    parm.weighting_function->key = "weighting_function";
+    parm.weighting_function->type = TYPE_STRING;
+    parm.weighting_function->required = NO;
+    parm.weighting_function->answer = "none";
+    parm.weighting_function->options = "none,gaussian,exponential,file";
+    G_asprintf((char **)&(parm.weighting_function->descriptions),
+               "none;%s;"
+               "gaussian;%s;"
+               "exponential;%s;"
+               "file;%s;",
+               _("No weighting"),
+               _("Gaussian weighting function"),
+               _("Exponential weighting function"),
+               _("File with a custom weighting matrix"));
+    parm.weighting_function->description = _("Weighting function");
+    parm.weighting_function->multiple = NO;
+
+    parm.weighting_factor = G_define_option();
+    parm.weighting_factor->key = "weighting_factor";
+    parm.weighting_factor->type = TYPE_DOUBLE;
+    parm.weighting_factor->required = NO;
+    parm.weighting_factor->multiple = NO;
+    parm.weighting_factor->description = _("Factor used in the selected weighting function (ignored for none and file)");
 
     parm.weight = G_define_standard_option(G_OPT_F_INPUT);
     parm.weight->key = "weight";
     parm.weight->required = NO;
     parm.weight->description = _("Text file containing weights");
 
-    parm.gauss = G_define_option();
-    parm.gauss->key = "gauss";
-    parm.gauss->type = TYPE_DOUBLE;
-    parm.gauss->required = NO;
-    parm.gauss->description = _("Sigma (in cells) for Gaussian filter");
-
     parm.quantile = G_define_option();
     parm.quantile->key = "quantile";
     parm.quantile->type = TYPE_DOUBLE;
@@ -238,6 +257,14 @@ int main(int argc, char *argv[])
     parm.quantile->multiple = YES;
     parm.quantile->description = _("Quantile to calculate for method=quantile");
     parm.quantile->options = "0.0-1.0";
+    parm.quantile->guisection = _("Neighborhood");
+
+    parm.title = G_define_option();
+    parm.title->key = "title";
+    parm.title->key_desc = "phrase";
+    parm.title->type = TYPE_STRING;
+    parm.title->required = NO;
+    parm.title->description = _("Title for output raster map");
 
     flag.align = G_define_flag();
     flag.align->key = 'a';
@@ -258,13 +285,19 @@ int main(int argc, char *argv[])
 	G_fatal_error(_("Neighborhood size must be odd"));
     ncb.dist = ncb.nsize / 2;
 
-    if (parm.weight->answer && flag.circle->answer)
+    if (strcmp(parm.weighting_function->answer, "none") && flag.circle->answer)
 	G_fatal_error(_("-%c and %s= are mutually exclusive"),
-			flag.circle->key, parm.weight->key);
+			flag.circle->key, parm.weighting_function->answer);
 
-    if (parm.weight->answer && parm.gauss->answer)
-	G_fatal_error(_("%s= and %s= are mutually exclusive"),
-			parm.weight->key, parm.gauss->key);
+    if (strcmp(parm.weighting_function->answer, "file") == 0 && !parm.weight->answer)
+	G_fatal_error(_("File with weighting matrix is missing."));
+
+    /* Check if weighting factor is given for all other weighting functions*/
+    if (strcmp(parm.weighting_function->answer, "none") &&
+        strcmp(parm.weighting_function->answer, "file") &&
+        !parm.weighting_factor->answer)
+	G_fatal_error(_("Weighting function '%s' requires a %s."),
+			parm.weighting_function->answer, parm.weighting_factor->key);
 
     ncb.oldcell = parm.input->answer;
 
@@ -299,12 +332,15 @@ int main(int argc, char *argv[])
     weights = 0;
     ncb.weights = NULL;
     ncb.mask = NULL;
-    if (parm.weight->answer) {
+    if (strcmp(parm.weighting_function->answer, "file") == 0) {
 	read_weights(parm.weight->answer);
 	weights = 1;
     }
-    else if (parm.gauss->answer) {
-	gaussian_weights(atof(parm.gauss->answer));
+    else if (strcmp(parm.weighting_function->answer, "none")) {
+	G_verbose_message(_("Computing %s weights..."),
+			      parm.weighting_function->answer);
+	compute_weights(parm.weighting_function->answer,
+	                atof(parm.weighting_factor->answer));
 	weights = 1;
     }
     
@@ -325,19 +361,14 @@ int main(int argc, char *argv[])
 		out->method_fn_w = menu[method].method_w;
 	    }
 	    else {
-		if (parm.weight->answer) {
+		if (strcmp(parm.weighting_function->answer,"none")) {
 		    G_warning(_("Method %s not compatible with weighing window, using weight mask instead"),
 			      method_name);
 		    if (!have_weights_mask) {
 			weights_mask();
 			have_weights_mask = 1;
 		    }
-		}
-		else if (parm.gauss->answer) {
-		    G_warning(_("Method %s not compatible with Gaussian filter, using unweighed version instead"),
-			      method_name);
-		}
-		
+		}		
 		out->method_fn = menu[method].method;
 		out->method_fn_w = NULL;
 	    }

raster/r.neighbors/r.neighbors.html

@@ -1,11 +1,11 @@
 <h2>DESCRIPTION</h2>
 
 <em><b>r.neighbors</b></em> looks at each cell in a raster input
-file, and examines the values assigned to the
-cells in some user-defined "neighborhood" around it.  It
+map, and examines the values assigned to the
+cells in some user-defined "neighborhood" around it. It
 outputs a new raster map layer in which each cell is
 assigned a value that is some (user-specified)
-function of the values in that cell's neighborhood.  For
+function of the values in that cell's neighborhood. For
 example, each cell in the output layer might be assigned a
 value equal to the average of the values
 appearing in its 3 x 3 cell "neighborhood" in the input
@@ -58,15 +58,46 @@ Without using the selection map, the output map would look like this:
 1 1 1 1 1
 </pre></div>
 
-<p>Optionally, the user can also specify the <b>TITLE</b> to
-be assigned to the raster map layer <b>output</b>, elect
-to not align the resolution of the output with that of the
-input (the <b>-a</b> option), and run <em><b>r.neighbors</b></em>
-with a custom matrix weights with the <em>weight</em> option.
+<p>It is also possible to weigh cells within the neighborhood. This can
+be either done with a custom weights matrix or by specifying a weighting
+function.
+
+<p>In order to use a custom weights matrix, <i>file</i> needs to be
+specified as a <b>weighting_function</b> and a path to a text file
+containing the weights needs to be given in the <b>weight</b> option.
+
+<p>Alternatively, <i>gaussian</i> and <i>exponential</i> weighting
+functions can be selected as weighting function.
+
+<p>For the <i>gaussian</i> weighting function, the user specifies the
+sigma value (&sigma;) for the gauss filter in the <b>weighting_factor</b>
+option. The sigma value represents the standard deviation of the gaussian
+distribution, where the weighting formula for the gaussian filter is
+defined as follows:
+<p>exp(-(x&#42;x+y&#42;y)/(2&#42;&sigma;&#94;2))/(2&#42;&pi;&#42;&sigma;&#94;2)
+<p>Lower values for sigma result in a steeper curve, so that more weight
+is put on cells close to the center of the moving window with a steeper
+decrease in weights with distance from the center.
+
+<p>For the <i>exponential</i> weighting function, the user specifies a
+factor for an exponential kernel in the <b>weighting_factor</b>.
+Negative factors result in negative exponential decrease in weights from
+the center cell. The weighting formula for the exponential kernel is
+defined as follows:
+<p>exp(factor*sqrt(x&#42;x+y&#42;y))
+<p>Stronger negative values for the factor result in a steeper curve,
+so that more weight is put on cells close to the center of the moving
+window with a steeper decrease in weights with distance from the center.
+
+<p>Optionally, the user can also run <em><b>r.neighbors</b></em> specify
+the <b>TITLE</b> to be assigned to the raster map layer <b>output</b>,
+select to not align the resolution of the output with that of the
+input (the <b>-a</b> option).
 These options are described further below.
+
 <p>
 <em>Neighborhood Operation Methods:</em>
-The <b>neighborhood</b> operators determine what new 
+The <b>neighborhood</b> operators determine what new
 value a center cell in a neighborhood will have after examining
 values inside its neighboring cells.
 Each cell in a raster map layer becomes the center cell of a neighborhood 
@@ -188,13 +219,15 @@ For example, a size value of 3 will result in
 <p>
 <em>Matrix weights:</em>
 A custom matrix can be used if none of the neighborhood operation
-methods are desirable by using the <b>weight</b>.  This option must
+methods are desirable by using the <b>weight</b>. This option must
 be used in conjunction with the <b>size</b> option to specify the
-matrix size.  The weights desired are to be entered into a text file.
-For example, to calculate the focal mean with a matrix <b>size</b> of
-3,
+matrix size and <i>file</i> needs to be specified as
+<b>weighting_function</b>. The weights desired are to be entered into a
+text file. For example, to calculate the focal mean with a matrix
+<b>size</b> of 3,
 <div class="code"><pre>
-r.neigbors in=input.map out=output.map size=3 weight=weights.txt
+r.neigbors in=input.map out=output.map size=3 weighting_function=file \
+weight=weights.txt
 </pre></div>
 
 The contents of the weight.txt file:
@@ -215,7 +248,7 @@ This corresponds to the following 3x3 matrix:
 +-+-+-+
 </pre></div>
 
-To calculate an annulus shaped neighborhood the contents of weight.txt file 
+To calculate an annulus shaped neighborhood the contents of weight.txt file
 may be e.g. for size=5:
 <div class="code"><pre>
  0 1 1 1 0
@@ -401,7 +434,7 @@ r.neighbors input=random_cells output=counts method=count
 
 Optionally - exclude centre cell from the count (= only look around)
 <div class="code"><pre>
-r.mapcalc "cound_around = if( isnull(random_cells), counts, counts - 1)"
+r.mapcalc "count_around = if( isnull(random_cells), counts, counts - 1)"
 </pre></div>
 
 

raster/r.neighbors/readweights.c

@@ -1,3 +1,4 @@
+#include <string.h>
 #include <stdio.h>
 #include <math.h>
 #include <grass/gis.h>
@@ -12,33 +13,55 @@ void read_weights(const char *filename)
 
     ncb.weights = G_malloc(ncb.nsize * sizeof(DCELL *));
     for (i = 0; i < ncb.nsize; i++)
-	ncb.weights[i] = G_malloc(ncb.nsize * sizeof(DCELL));
+        ncb.weights[i] = G_malloc(ncb.nsize * sizeof(DCELL));
 
     if (!fp)
-	G_fatal_error(_("Unable to open weights file %s"), filename);
+        G_fatal_error(_("Unable to open weights file %s"), filename);
 
     for (i = 0; i < ncb.nsize; i++)
-	for (j = 0; j < ncb.nsize; j++)
-	    if (fscanf(fp, "%lf", &ncb.weights[i][j]) != 1)
-		G_fatal_error(_("Error reading weights file %s"), filename);
+        for (j = 0; j < ncb.nsize; j++)
+            if (fscanf(fp, "%lf", &ncb.weights[i][j]) != 1)
+                G_fatal_error(_("Error reading weights file %s"), filename);
 
     fclose(fp);
 }
 
-void gaussian_weights(double sigma)
+double gaussian(double factor, double squared_distance)
+{
+    double sigma2 = factor * factor;
+
+    return exp(-squared_distance / (2 * sigma2)) / (2 * M_PI * sigma2);
+}
+
+double exponential(double factor, double squared_distance)
+{
+    return exp(factor * sqrt(squared_distance));
+}
+
+void compute_weights(const char *function_type, double factor)
 {
-    double sigma2 = sigma * sigma;
     int i, j;
+    double (*weight) (double, double);
+
+    if (!strcmp(function_type, "gaussian")) {
+        weight = gaussian;
+    }
+    else if (!strcmp(function_type, "exponential")) {
+        weight = exponential;
+    }
+
 
     ncb.weights = G_malloc(ncb.nsize * sizeof(DCELL *));
     for (i = 0; i < ncb.nsize; i++)
-	ncb.weights[i] = G_malloc(ncb.nsize * sizeof(DCELL));
+        ncb.weights[i] = G_malloc(ncb.nsize * sizeof(DCELL));
 
     for (i = 0; i < ncb.nsize; i++) {
-	double y = i - ncb.dist;
-	for (j = 0; j < ncb.nsize; j++) {
-	    double x = j - ncb.dist;
-	    ncb.weights[i][j] = exp(-(x*x+y*y)/(2*sigma2))/(2*M_PI*sigma2);
-	}
+        double y = i - ncb.dist;
+
+        for (j = 0; j < ncb.nsize; j++) {
+            double x = j - ncb.dist;
+
+            ncb.weights[i][j] = weight(factor, x * x + y * y);
+        }
     }
 }

raster/r.neighbors/testsuite/test_r_neighbors.py

@@ -531,10 +531,10 @@ class TestNeighbors(TestCase):
         self.to_remove.extend(outputs)
         self.assertModule(
             "r.neighbors",
-            flags="c",
             input="elevation",
             size=7,
-            gauss=2.0,
+            weighting_function="gaussian",
+            weighting_factor=2.0,
             output=",".join(outputs),
             method=list(self.test_options[test_case].keys()),
         )
@@ -558,6 +558,7 @@ class TestNeighbors(TestCase):
             "r.neighbors",
             input="elevation",
             size=5,
+            weighting_function="file",
             output=",".join(outputs),
             method=list(self.test_options[test_case].keys()),
             weight=weights,
@@ -570,6 +571,7 @@ class TestNeighbors(TestCase):
             input="elevation",
             flags="c",
             size=5,
+            weighting_function="file",
             output="{}_fails".format(test_case),
             method="sum",
             weight=weights,