v.kernel optimization

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

vector/v.kernel/function.c

@@ -4,7 +4,7 @@
 #include <grass/glocale.h>
 #include "global.h"
 
-static double (*kernelfn)(int dimension, double bandwidth, double x);
+static double (*kernelfn)(double term, double bandwidth, double x);
 
 /*********************** Gaussian ****************************/
 /* probability for gaussian distribution */
@@ -27,9 +27,9 @@ double gaussianFunction(double x, double sigma, double dimension)
 
 
 /* probability for gaussian distribution */
-double gaussianKernel(double x, double term)
+double gaussianKernel(double x, double termx)
 {
-    return (term * exp(-(x * x) / 2.));
+    return (termx * exp(-(x * x) / 2.));
 }
 
 
@@ -101,11 +101,16 @@ double euclidean_distance(double *x, double *y, int n)
 }
 
 
-/********************************************************************/
-double gaussianKernel2(int dimension, double bandwidth, double x)
+/*****************kernel density functions******************************/
+
+double gaussianKernel2(double term, double bandwidth, double x)
 {
+    /* term is set by setKernelFunction */
+
+    /*
     double term =
 	1. / (pow(bandwidth, dimension) * pow((2. * M_PI), dimension / 2.));
+    */
 
     x /= bandwidth;
 
@@ -114,136 +119,188 @@ double gaussianKernel2(int dimension, double bandwidth, double x)
 
 /* Note: these functions support currently only 1D and 2D, consider this for example 
  * before using them for 3d grid */
-double uniformKernel(int dimension, double bandwidth, double x)
+double uniformKernel(double term, double bandwidth, double x)
 {
-    double k;
+    /* term is set by setKernelFunction */
 
     if (x > bandwidth)
 	return 0;
 
+    /*
     if (dimension == 2)
-	k = 2. / (M_PI * pow(bandwidth, 2));
+	term = 2. / (M_PI * pow(bandwidth, 2));
     else
-	k = 1. / bandwidth;
+	term = 1. / bandwidth;
+    term *= (1. / 2);
 
     x /= bandwidth;
+    */
 
-    return k * (1. / 2);
+    return term;
 }
 
-double triangularKernel(int dimension, double bandwidth, double x)
+double triangularKernel(double term, double bandwidth, double x)
 {
-    double k;
+    /* term is set by setKernelFunction */
 
     if (x > bandwidth)
 	return 0;
 
+    /*
     if (dimension == 2)
-	k = 3. / (M_PI * pow(bandwidth, 2));
+	term = 3. / (M_PI * pow(bandwidth, 2));
     else
-	k = 1. / bandwidth;
+	term = 1. / bandwidth;
+    */
 
     x /= bandwidth;
 
-    return k * (1 - x);
+    return term * (1 - x);
 }
 
-double epanechnikovKernel(int dimension, double bandwidth, double x)
+double epanechnikovKernel(double term, double bandwidth, double x)
 {
-    double k;
+    /* term is set by setKernelFunction */
 
     if (x > bandwidth)
 	return 0;
 
+    /*
     if (dimension == 2)
-	k = 8. / (M_PI * 3. * pow(bandwidth, 2));
+	term = 8. / (M_PI * 3. * pow(bandwidth, 2));
     else
-	k = 1. / bandwidth;
+	term = 1. / bandwidth;
+    term *= (3. / 4.);
+    */
 
     x /= bandwidth;
 
-    return k * (3. / 4. * (1 - x * x));
+    /* return term * (3. / 4. * (1 - x * x)); */
+    return term * (1 - x * x);
 }
 
-double quarticKernel(int dimension, double bandwidth, double x)
+double quarticKernel(double term, double bandwidth, double x)
 {
-    double k;
+    /* term is set by setKernelFunction */
 
     if (x > bandwidth)
 	return 0;
 
+    /*
     if (dimension == 2)
-	k = 16. / (M_PI * 5. * pow(bandwidth, 2));
+	term = 16. / (M_PI * 5. * pow(bandwidth, 2));
     else
-	k = 1. / bandwidth;
+	term = 1. / bandwidth;
+    term *= (15. / 16.);
+    */
 
     x /= bandwidth;
 
-    return k * (15. / 16. * pow(1 - x * x, 2));
+    /* return term * (15. / 16. * pow(1 - x * x, 2)); */
+    return term * pow(1 - x * x, 2);
 }
 
-double triweightKernel(int dimension, double bandwidth, double x)
+double triweightKernel(double term, double bandwidth, double x)
 {
-    double k;
+    /* term is set by setKernelFunction */
 
     if (x > bandwidth)
 	return 0;
 
+    /*
     if (dimension == 2)
-	k = 128. / (M_PI * 35. * pow(bandwidth, 2));
+	term = 128. / (M_PI * 35. * pow(bandwidth, 2));
     else
-	k = 1. / bandwidth;
+	term = 1. / bandwidth;
+    term *= (35. / 32);
+    */
 
     x /= bandwidth;
 
-    return k * (35. / 32 * pow(1 - x * x, 3));
+    /* return term * (35. / 32 * pow(1 - x * x, 3)); */
+    return term * pow(1 - x * x, 3);
 }
 
-double cosineKernel(int dimension, double bandwidth, double x)
+double cosineKernel(double term, double bandwidth, double x)
 {
-    double k;
+    /* term is set by setKernelFunction */
 
     if (x > bandwidth)
 	return 0;
 
+    /*
     if (dimension == 2)
-	k = 1. / (2 * (M_PI / 2 - 1) * pow(bandwidth, 2));
+	term = 1. / (2 * (M_PI / 2 - 1) * pow(bandwidth, 2));
     else
-	k = 1. / bandwidth;
+	term = 1. / bandwidth;
+    term *= (M_PI / 4.);
+    */
 
     x /= bandwidth;
 
-    return k * (M_PI / 4. * cos(M_PI / 2. * x));
+    /* return term * (M_PI / 4. * cos(M_PI / 2. * x)); */
+    return term * cos(M_PI / 2. * x);
 }
 
-double kernelFunction(int dimension, double bandwidth, double x)
+double kernelFunction(double term, double bandwidth, double x)
 {
-    return kernelfn(dimension, bandwidth, x);
+    return kernelfn(term, bandwidth, x);
 }
 
-void setKernelFunction(int function)
+void setKernelFunction(int function, int dimension, double bandwidth, double *term)
 {
     switch (function) {
     case KERNEL_UNIFORM:
 	kernelfn = uniformKernel;
+	if (dimension == 2)
+	    *term = 2. / (M_PI * pow(bandwidth, 2));
+	else
+	    *term = 1. / bandwidth;
+	*term *= (1. / 2);
 	break;
     case KERNEL_TRIANGULAR:
 	kernelfn = triangularKernel;
+	if (dimension == 2)
+	    *term = 3. / (M_PI * pow(bandwidth, 2));
+	else
+	    *term = 1. / bandwidth;
 	break;
     case KERNEL_EPANECHNIKOV:
 	kernelfn = epanechnikovKernel;
+	if (dimension == 2)
+	    *term = 8. / (M_PI * 3. * pow(bandwidth, 2));
+	else
+	    *term = 1. / bandwidth;
+	*term *= (3. / 4.);
 	break;
     case KERNEL_QUARTIC:
 	kernelfn = quarticKernel;
+	if (dimension == 2)
+	    *term = 16. / (M_PI * 5. * pow(bandwidth, 2));
+	else
+	    *term = 1. / bandwidth;
+	*term *= (15. / 16.);
 	break;
     case KERNEL_TRIWEIGHT:
 	kernelfn = triweightKernel;
+	if (dimension == 2)
+	    *term = 128. / (M_PI * 35. * pow(bandwidth, 2));
+	else
+	    *term = 1. / bandwidth;
+	*term *= (35. / 32);
 	break;
     case KERNEL_GAUSSIAN:
 	kernelfn = gaussianKernel2;
+	*term =
+	    1. / (pow(bandwidth, dimension) * pow((2. * M_PI), dimension / 2.));
 	break;
     case KERNEL_COSINE:
 	kernelfn = cosineKernel;
+	if (dimension == 2)
+	    *term = 1. / (2 * (M_PI / 2 - 1) * pow(bandwidth, 2));
+	else
+	    *term = 1. / bandwidth;
+	*term *= (M_PI / 4.);
 	break;
     default:
 	G_fatal_error("Unknown kernel function");

vector/v.kernel/global.h

@@ -15,8 +15,8 @@
 #define KERNEL_COSINE       6
 
 
-void setKernelFunction(int function);
-double kernelFunction(int dimension, double bandwidth, double x);
+void setKernelFunction(int function, int dimension, double bandwidth, double *term);
+double kernelFunction(double term, double bandwidth, double x);
 
 double euclidean_distance(double *x, double *y, int n);
 double gaussian2dBySigma(double d, double sigma);
@@ -37,7 +37,7 @@ double compute_all_distances(double **coordinate, double **dists, int n,
 double compute_all_net_distances(struct Map_info *In, struct Map_info *Net,
 				 double netmax, double **dists, double dmax);
 void compute_distance(double N, double E, struct Map_info *In, double sigma,
-		      double *gaussian, double dmax);
+		      double term, double *gaussian, double dmax);
 void compute_net_distance(double x, double y, struct Map_info *In,
 			  struct Map_info *Net, double netmax, double sigma,
-			  double *gaussian, double dmax, int node_method);
+			  double term, double *gaussian, double dmax, int node_method);

vector/v.kernel/main.c

@@ -95,7 +95,7 @@ int main(int argc, char **argv)
     double sigmaOptimal;
     struct GModule *module;
     double dsize;
-    double term;
+    double term = 0;
 
     double gausmax = 0;
     int notreachable = 0;
@@ -248,8 +248,6 @@ int main(int argc, char **argv)
 	    G_fatal_error(_("Optimal standard deviation calculation is supported only for kernel function 'gaussian'."));
 	}
     }
-    
-    setKernelFunction(kernel_function);
 
     if (flag_q->answer) {
 	flag_o->answer = 1;
@@ -309,7 +307,6 @@ int main(int argc, char **argv)
     else {
 	/* check and open the name of output map */
 	if (!flag_q->answer) {
-	    Rast_set_fp_type(DCELL_TYPE);
 	    fdout = Rast_open_new(out_opt->answer, DCELL_TYPE);
 
 	    /* open mask file */
@@ -377,11 +374,16 @@ int main(int argc, char **argv)
 	}
     }
 
-    term = 1. / (pow(sigma, dimension) * pow((2. * M_PI), dimension / 2.));
-
     dmax = sigma;
     if (kernel_function == KERNEL_GAUSSIAN)
-	dmax = sigma * 4.;
+	dmax = sigma * 4.;  /* should be sigma /= 4.; */
+
+    if (net_opt->answer) {
+	setKernelFunction(kernel_function, 1, sigma, &term);
+    }
+    else {
+	setKernelFunction(kernel_function, 2, sigma, &term);
+    }
 
     if (net) {
 	int line, nlines;
@@ -426,7 +428,7 @@ int main(int argc, char **argv)
 		G_debug(3, "  segment = %d, offset = %f, xy = %f %f", seg,
 			offset1, x, y);
 
-		compute_net_distance(x, y, &In, &Net, netmax, sigma,
+		compute_net_distance(x, y, &In, &Net, netmax, sigma, term,
 				     &gaussian, dmax, node_method);
 		gaussian *= multip;
 		if (gaussian > gausmax)
@@ -516,7 +518,7 @@ int main(int argc, char **argv)
 		E = Rast_col_to_easting(col + 0.5, &window);
 
 		/* compute_distance(N, E, &In, sigma, term, &gaussian, dmax); */
-		compute_distance(N, E, &In, sigma, &gaussian, dmax);
+		compute_distance(N, E, &In, sigma, term, &gaussian, dmax);
 
 		output_cell[col] = multip * gaussian;
 		if (gaussian > gausmax)
@@ -697,7 +699,7 @@ int count_node_arcs(struct Map_info *Map, int node)
 /* Compute gausian for x, y along Net, using all points in In */
 void compute_net_distance(double x, double y, struct Map_info *In,
 			  struct Map_info *Net, double netmax, double sigma,
-			  double *gaussian, double dmax, int node_method)
+			  double term, double *gaussian, double dmax, int node_method)
 {
     int i;
     double dist, kernel;
@@ -762,7 +764,7 @@ void compute_net_distance(double x, double y, struct Map_info *In,
 	    continue;
 
 	/* kernel = gaussianKernel(dist / sigma, term); */
-	kernel = kernelFunction(1, sigma, dist);
+	kernel = kernelFunction(term, sigma, dist);
 
 	if (node_method == NODE_EQUAL_SPLIT) {
 	    int j, node;
@@ -790,7 +792,7 @@ void compute_net_distance(double x, double y, struct Map_info *In,
 }
 
 void compute_distance(double N, double E, struct Map_info *In,
-		      double sigma, double *gaussian,
+		      double sigma, double term, double *gaussian,
 		      double dmax)
 {
     int line, nlines;
@@ -830,6 +832,6 @@ void compute_distance(double N, double E, struct Map_info *In,
 
 	if (dist <= dmax)
 	    /* *gaussian += gaussianKernel(dist / sigma, term); */
-	    *gaussian += kernelFunction(2, sigma, dist);
+	    *gaussian += kernelFunction(term, sigma, dist);
     }
 }