[GRASS-SVN] r53956 - grass/trunk/vector/v.kernel
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Nov 21 12:57:47 PST 2012
Author: mmetz
Date: 2012-11-21 12:57:47 -0800 (Wed, 21 Nov 2012)
New Revision: 53956
Modified:
grass/trunk/vector/v.kernel/function.c
grass/trunk/vector/v.kernel/global.h
grass/trunk/vector/v.kernel/main.c
Log:
v.kernel optimization
Modified: grass/trunk/vector/v.kernel/function.c
===================================================================
--- grass/trunk/vector/v.kernel/function.c 2012-11-21 20:56:52 UTC (rev 53955)
+++ grass/trunk/vector/v.kernel/function.c 2012-11-21 20:57:47 UTC (rev 53956)
@@ -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 @@
/* 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 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 @@
/* 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");
Modified: grass/trunk/vector/v.kernel/global.h
===================================================================
--- grass/trunk/vector/v.kernel/global.h 2012-11-21 20:56:52 UTC (rev 53955)
+++ grass/trunk/vector/v.kernel/global.h 2012-11-21 20:57:47 UTC (rev 53956)
@@ -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_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);
Modified: grass/trunk/vector/v.kernel/main.c
===================================================================
--- grass/trunk/vector/v.kernel/main.c 2012-11-21 20:56:52 UTC (rev 53955)
+++ grass/trunk/vector/v.kernel/main.c 2012-11-21 20:57:47 UTC (rev 53956)
@@ -95,7 +95,7 @@
double sigmaOptimal;
struct GModule *module;
double dsize;
- double term;
+ double term = 0;
double gausmax = 0;
int notreachable = 0;
@@ -248,8 +248,6 @@
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 @@
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,12 +374,17 @@
}
}
- 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;
struct line_pnts *Points, *SPoints;
@@ -426,7 +428,7 @@
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 @@
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 @@
/* 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 @@
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_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 @@
if (dist <= dmax)
/* *gaussian += gaussianKernel(dist / sigma, term); */
- *gaussian += kernelFunction(2, sigma, dist);
+ *gaussian += kernelFunction(term, sigma, dist);
}
}
More information about the grass-commit
mailing list