Add Gaussian filter

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@36892 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.neighbors/local_proto.h

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

raster/r.neighbors/main.c

@@ -89,6 +89,7 @@ int main(int argc, char *argv[])
 	struct Option *method, *size;
 	struct Option *title;
 	struct Option *weight;
+	struct Option *gauss;
     } parm;
     struct
     {
@@ -152,6 +153,12 @@ int main(int argc, char *argv[])
     parm.weight->gisprompt = "old_file,file,input";
     parm.weight->description = _("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");
+
     flag.align = G_define_flag();
     flag.align->key = 'a';
     flag.align->description = _("Do not align output with the input");
@@ -174,6 +181,9 @@ int main(int argc, char *argv[])
     if (parm.weight->answer && flag.circle->answer)
 	G_fatal_error(_("weight= and -c are mutually exclusive"));
 
+    if (parm.weight->answer && parm.gauss->answer)
+	G_fatal_error(_("weight= and gauss= are mutually exclusive"));
+
     ncb.oldcell = parm.input->answer;
     ncb.newcell = parm.output->answer;
 
@@ -227,6 +237,11 @@ int main(int argc, char *argv[])
 	if (!newvalue_w)
 	    weights_mask();
     }
+    else if (parm.gauss->answer) {
+	if (!newvalue_w)
+	    G_fatal_error(_("method %s not compatible with Gaussian filter"));
+	gaussian_weights(atof(parm.gauss->answer));
+    }
     else
 	newvalue_w = NULL;
 

raster/r.neighbors/readweights.c

@@ -1,4 +1,5 @@
 #include <stdio.h>
+#include <math.h>
 #include <grass/gis.h>
 #include <grass/glocale.h>
 #include "ncb.h"
@@ -23,3 +24,21 @@ void read_weights(const char *filename)
 
     fclose(fp);
 }
+
+void gaussian_weights(double sigma)
+{
+    double sigma2 = sigma * sigma;
+    int i, j;
+
+    ncb.weights = G_malloc(ncb.nsize * sizeof(DCELL *));
+    for (i = 0; i < ncb.nsize; i++)
+	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);
+	}
+    }
+}