[GRASS-SVN] r36894 - grass/branches/develbranch_6/raster/r.neighbors

[svn_grass at osgeo.org]

Author: neteler
Date: 2009-04-25 02:30:54 -0400 (Sat, 25 Apr 2009)
New Revision: 36894

Modified:
   grass/branches/develbranch_6/raster/r.neighbors/local_proto.h
   grass/branches/develbranch_6/raster/r.neighbors/main.c
   grass/branches/develbranch_6/raster/r.neighbors/readweights.c
Log:
glynn: Add Gaussian filter (merge from trunk, r36892)

Modified: grass/branches/develbranch_6/raster/r.neighbors/local_proto.h
===================================================================
--- grass/branches/develbranch_6/raster/r.neighbors/local_proto.h	2009-04-25 05:59:00 UTC (rev 36893)
+++ grass/branches/develbranch_6/raster/r.neighbors/local_proto.h	2009-04-25 06:30:54 UTC (rev 36894)
@@ -22,3 +22,4 @@
 
 /* read_weights.c */
 extern void read_weights(const char *);
+extern void gaussian_weights(double);

Modified: grass/branches/develbranch_6/raster/r.neighbors/main.c
===================================================================
--- grass/branches/develbranch_6/raster/r.neighbors/main.c	2009-04-25 05:59:00 UTC (rev 36893)
+++ grass/branches/develbranch_6/raster/r.neighbors/main.c	2009-04-25 06:30:54 UTC (rev 36894)
@@ -89,6 +89,7 @@
 	struct Option *method, *size;
 	struct Option *title;
 	struct Option *weight;
+	struct Option *gauss;
     } parm;
     struct
     {
@@ -153,6 +154,12 @@
     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");
@@ -189,6 +196,8 @@
 	G_fatal_error(_("<%s> is an illegal file name"), p);
     }
     ncb.newcell.mapset = G_mapset();
+    if (parm.weight->answer && parm.gauss->answer)
+	G_fatal_error(_("weight= and gauss= are mutually exclusive"));
 
     if (!flag.align->answer) {
 	if (G_get_cellhd(ncb.oldcell.name, ncb.oldcell.mapset, &cellhd) < 0)
@@ -240,6 +249,11 @@
 	if (!newvalue_w)
 	    weights_mask();
     }
+    else if (parm.gauss->answer) {
+	if (!newvalue_w)
+	    G_fatal_error(_("Method %s not compatible with Gaussian filter"), parm.method->answer);
+	gaussian_weights(atof(parm.gauss->answer));
+    }
     else
 	newvalue_w = NULL;
 

Modified: grass/branches/develbranch_6/raster/r.neighbors/readweights.c
===================================================================
--- grass/branches/develbranch_6/raster/r.neighbors/readweights.c	2009-04-25 05:59:00 UTC (rev 36893)
+++ grass/branches/develbranch_6/raster/r.neighbors/readweights.c	2009-04-25 06:30:54 UTC (rev 36894)
@@ -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 @@
 
     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);
+	}
+    }
+}