[GRASS-SVN] r45955 -
grass/branches/releasebranch_6_4/vector/v.kernel
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Apr 13 15:02:30 EDT 2011
Author: martinl
Date: 2011-04-13 12:02:30 -0700 (Wed, 13 Apr 2011)
New Revision: 45955
Modified:
grass/branches/releasebranch_6_4/vector/v.kernel/description.html
grass/branches/releasebranch_6_4/vector/v.kernel/function.c
grass/branches/releasebranch_6_4/vector/v.kernel/global.h
grass/branches/releasebranch_6_4/vector/v.kernel/main.c
Log:
Radim Blazek: additional kernel functions, equal split on networks (backport from trunk r45393, 45398, r45406)
(merge r45407 from devbr6)
Modified: grass/branches/releasebranch_6_4/vector/v.kernel/description.html
===================================================================
--- grass/branches/releasebranch_6_4/vector/v.kernel/description.html 2011-04-13 16:28:23 UTC (rev 45954)
+++ grass/branches/releasebranch_6_4/vector/v.kernel/description.html 2011-04-13 19:02:30 UTC (rev 45955)
@@ -1,9 +1,16 @@
<H2>DESCRIPTION</H2>
-Generates a raster density map from vector points data using
-a moving 2D isotropic Gaussian kernel or
-optionally generates a vector density map on vector network
-with a 1D kernel.
+<em>v.kernel</em> generates a raster density map from vector points data using
+a moving kernel. Available <a href="http://en.wikipedia.org/wiki/Kernel_(statistics)#Kernel_functions_in_common_use">kernel density functions</a> are <em>uniform,
+triangular, epanechnikov, quartic, triweight, gaussian, cosine</em>,
+default is <em>gaussian</em>.
+<p>
+The module can also generate a vector density map on a vector network.
+Conventional kernel functions produce biased estimates by overestimating
+the densities around network nodes, whereas the equal split method of
+Okabe et al. (2009) produces unbiased density estimates. The equal split
+method uses the kernel function selected with the <em>kernel</em> option
+and can be enabled with <em>node=split</em>.
<H2>NOTES</H2>
@@ -18,16 +25,23 @@
<H2>LIMITATIONS</H2>
-The modules only considers the presence of point, but not
+The module only considers the presence of points, but not
(yet) any attribute values.
<H2>SEE ALSO</H2>
<A HREF="v.surf.rst.html">v.surf.rst</A>
-<H2>AUTHOR</H2>
+<H2>REFERENCES</H2>
+Okabe, A., Satoh, T., Sugihara, K. (2009). <i>A kernel density estimation
+method for networks, its computational method and a GIS-based tool</i>.
+<b>International Journal of Geographical Information Science</b>, Vol 23(1),
+pp. 7-32.<br>
+DOI: <a href="http://dx.doi.org/10.1080/13658810802475491">10.1080/13658810802475491</a>
-Stefano Menegon, <a href=http://mpa.itc.it>ITC-irst</a>, Trento, Italy
+<H2>AUTHORS</H2>
+
+Stefano Menegon, <a href="http://mpa.itc.it/">ITC-irst</a>, Trento, Italy
<BR>
-Radim Blazek (network part)
+Radim Blazek (additional kernel density functions and network part)
<p><i>Last changed: $Date$</i>
Modified: grass/branches/releasebranch_6_4/vector/v.kernel/function.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/v.kernel/function.c 2011-04-13 16:28:23 UTC (rev 45954)
+++ grass/branches/releasebranch_6_4/vector/v.kernel/function.c 2011-04-13 19:02:30 UTC (rev 45955)
@@ -1,9 +1,10 @@
#include <math.h>
#include <grass/gis.h>
#include <grass/Vect.h>
+#include <grass/glocale.h>
#include "global.h"
-
+/*********************** Gaussian ****************************/
/* probability for gaussian distribution */
double gaussian2dBySigma(double d, double sigma)
{
@@ -63,10 +64,8 @@
res = 1. / (M_PI * (d * d + rs * rs));
}
else {
- res =
- segno(a) * (a / M_PI) * (pow(rs, 2. * a)) * (1 /
- pow(d * d + rs * rs,
- lambda));
+ res = segno(a) * (a / M_PI) * (pow(rs, 2. * a)) *
+ (1 / pow(d * d + rs * rs, lambda));
}
/* res=1./(M_PI*(d*d+rs*rs)); */
@@ -79,6 +78,7 @@
double d;
d = sqrt(-2 * sigma * sigma * log(prob * M_PI * 2 * sigma * sigma));
+
return (d);
}
@@ -97,3 +97,150 @@
return sqrt(out);
}
+
+
+/********************************************************************/
+double gaussianKernel2(int dimension, double bandwidth, double x)
+{
+ double term =
+ 1. / (pow(bandwidth, dimension) * pow((2. * M_PI), dimension / 2.));
+
+ x /= bandwidth;
+
+ return (term * exp(-(x * x) / 2.));
+}
+
+/* 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 k;
+
+ if (x > bandwidth)
+ return 0;
+
+ if (dimension == 2)
+ k = 2. / (M_PI * pow(bandwidth, 2));
+ else
+ k = 1. / bandwidth;
+
+ x /= bandwidth;
+
+ return k * (1. / 2);
+}
+
+double triangularKernel(int dimension, double bandwidth, double x)
+{
+ double k;
+
+ if (x > bandwidth)
+ return 0;
+
+ if (dimension == 2)
+ k = 3. / (M_PI * pow(bandwidth, 2));
+ else
+ k = 1. / bandwidth;
+
+ x /= bandwidth;
+
+ return k * (1 - x);
+}
+
+double epanechnikovKernel(int dimension, double bandwidth, double x)
+{
+ double k;
+
+ if (x > bandwidth)
+ return 0;
+
+ if (dimension == 2)
+ k = 8. / (M_PI * 3. * pow(bandwidth, 2));
+ else
+ k = 1. / bandwidth;
+
+ x /= bandwidth;
+
+ return k * (3. / 4. * (1 - x * x));
+}
+
+double quarticKernel(int dimension, double bandwidth, double x)
+{
+ double k;
+
+ if (x > bandwidth)
+ return 0;
+
+ if (dimension == 2)
+ k = 16. / (M_PI * 5. * pow(bandwidth, 2));
+ else
+ k = 1. / bandwidth;
+
+ x /= bandwidth;
+
+ return k * (15. / 16. * pow(1 - x * x, 2));
+}
+
+double triweightKernel(int dimension, double bandwidth, double x)
+{
+ double k;
+
+ if (x > bandwidth)
+ return 0;
+
+ if (dimension == 2)
+ k = 128. / (M_PI * 35. * pow(bandwidth, 2));
+ else
+ k = 1. / bandwidth;
+
+ x /= bandwidth;
+
+ return k * (35. / 32 * pow(1 - x * x, 3));
+}
+
+double cosineKernel(int dimension, double bandwidth, double x)
+{
+ double k;
+
+ if (x > bandwidth)
+ return 0;
+
+ if (dimension == 2)
+ k = 1. / (2 * (M_PI / 2 - 1) * pow(bandwidth, 2));
+ else
+ k = 1. / bandwidth;
+
+ x /= bandwidth;
+
+ return k * (M_PI / 4. * cos(M_PI / 2. * x));
+}
+
+double kernelFunction(int function, int dimension, double bandwidth, double x)
+{
+ if (dimension > 2 && function != KERNEL_GAUSSIAN) {
+ G_fatal_error(_("Dimension > 2 supported only by gaussian function"));
+ }
+ switch (function) {
+ case KERNEL_UNIFORM:
+ return uniformKernel(dimension, bandwidth, x);
+ break;
+ case KERNEL_TRIANGULAR:
+ return triangularKernel(dimension, bandwidth, x);
+ break;
+ case KERNEL_EPANECHNIKOV:
+ return epanechnikovKernel(dimension, bandwidth, x);
+ break;
+ case KERNEL_QUARTIC:
+ return quarticKernel(dimension, bandwidth, x);
+ break;
+ case KERNEL_TRIWEIGHT:
+ return triweightKernel(dimension, bandwidth, x);
+ break;
+ case KERNEL_GAUSSIAN:
+ return gaussianKernel2(dimension, bandwidth, x);
+ break;
+ case KERNEL_COSINE:
+ return cosineKernel(dimension, bandwidth, x);
+ break;
+ }
+ G_fatal_error("Unknown kernel function");
+}
Modified: grass/branches/releasebranch_6_4/vector/v.kernel/global.h
===================================================================
--- grass/branches/releasebranch_6_4/vector/v.kernel/global.h 2011-04-13 16:28:23 UTC (rev 45954)
+++ grass/branches/releasebranch_6_4/vector/v.kernel/global.h 2011-04-13 19:02:30 UTC (rev 45955)
@@ -1,5 +1,22 @@
#include <grass/gis.h>
+#define NODE_NONE 0
+/* According Okabe 2009 */
+#define NODE_SIMILAR 1
+#define NODE_EQUAL_SPLIT 2
+#define NODE_CONTINUOUS_EQUAL_SPLIT 3
+
+#define KERNEL_UNIFORM 0
+#define KERNEL_TRIANGULAR 1
+#define KERNEL_EPANECHNIKOV 2
+#define KERNEL_QUARTIC 3
+#define KERNEL_TRIWEIGHT 4
+#define KERNEL_GAUSSIAN 5
+#define KERNEL_COSINE 6
+
+
+double kernelFunction(int function, int dimension, double bandwidth, double x);
+
double euclidean_distance(double *x, double *y, int n);
double gaussian2dBySigma(double d, double sigma);
double gaussianFunction(double x, double sigma, double dimension);
@@ -18,7 +35,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 term, double *gaussian, double dmax);
+ double term, double *gaussian, double dmax, int kernel_function);
void compute_net_distance(double x, double y, struct Map_info *In,
struct Map_info *Net, double netmax, double sigma,
- double term, double *gaussian, double dmax);
+ double term, double *gaussian, double dmax, int node_method, int kernel_function);
Modified: grass/branches/releasebranch_6_4/vector/v.kernel/main.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/v.kernel/main.c 2011-04-13 16:28:23 UTC (rev 45954)
+++ grass/branches/releasebranch_6_4/vector/v.kernel/main.c 2011-04-13 19:02:30 UTC (rev 45955)
@@ -4,11 +4,12 @@
* MODULE: v.kernel
*
* AUTHOR(S): Stefano Menegon, ITC-irst, Trento, Italy
+* Radim Blazek (additional kernel functions, network part)
* PURPOSE: Generates a raster density map from vector points data using
-* a moving 2D isotropic Gaussian kernel or
+* a moving kernel function or
* optionally generates a vector density map on vector network
* with a 1D kernel
-* COPYRIGHT: (C) 2004 by the GRASS Development Team
+* COPYRIGHT: (C) 2004-2011 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
@@ -33,7 +34,7 @@
static double dimension = 2.;
- /* define score function L(window size) */
+/* define score function L(window size) */
double L(double smooth)
{
int ii;
@@ -54,18 +55,11 @@
if (!net)
resL *= 2.;
- resL =
- (1. / (pow(n, 2.) * pow(smooth, dimension))) * (resL +
- n *
- (gaussianFunction
- (0., 2.,
- dimension) -
- 2. *
- gaussianFunction(0.,
- 1.,
- dimension)))
- + (2. / (n * pow(smooth, dimension))) * gaussianFunction(0., 1.,
- dimension);
+ resL = (1. / (pow(n, 2.) * pow(smooth, dimension))) *
+ (resL + n * (gaussianFunction(0., 2., dimension) -
+ 2. * gaussianFunction(0., 1., dimension))) +
+ (2. / (n * pow(smooth, dimension))) *
+ gaussianFunction(0., 1., dimension);
/* resL = (1./(pow(n,2.)*pow(smooth,dimension))) * (resL + n*( gaussianFunction(0.,2.,dimension) - 2. * gaussianKernel(0.,term)) ) + (2./(n*pow(smooth,dimension)))*gaussianKernel(0.,term); */
G_debug(3, "smooth = %e resL = %e", smooth, resL);
@@ -80,12 +74,13 @@
{
struct Option *in_opt, *net_opt, *out_opt;
struct Option *stddev_opt, *dsize_opt, *segmax_opt, *netmax_opt,
- *multip_opt;
- struct Flag *flag_o, *flag_v, *flag_q;
+ *multip_opt, *node_opt, *kernel_opt;
+ struct Flag *flag_o, *flag_q, *flag_normalize, *flag_multiply;
char *mapset;
struct Map_info In, Net, Out;
int fdout = 0, maskfd = 0;
+ int node_method, kernel_function;
int row, col;
struct Cell_head window;
double gaussian;
@@ -101,6 +96,7 @@
double term;
double gausmax = 0;
+ int notreachable = 0;
/* Initialize the GIS calls */
G_gisinit(argv[0]);
@@ -108,8 +104,8 @@
module = G_define_module();
module->keywords = _("vector, kernel density");
module->description =
- _("Generates a raster density map from vector points data using a moving 2D isotropic Gaussian kernel or "
- "optionally generates a vector density map on vector network with a 1D kernel.");
+ _("Generates a raster density map from vector point data using a moving kernel or "
+ "optionally generates a vector density map on a vector network.");
in_opt = G_define_standard_option(G_OPT_V_INPUT);
in_opt->description = _("Input vector with training points");
@@ -160,6 +156,26 @@
multip_opt->description = _("Multiply the density result by this number");
multip_opt->answer = "1.";
+ node_opt = G_define_option();
+ node_opt->key = "node";
+ node_opt->type = TYPE_STRING;
+ node_opt->required = YES;
+ node_opt->description = _("Node method");
+ node_opt->options = "none,split";
+ node_opt->answer = "none";
+ node_opt->descriptions =
+ _("none;No method applied at nodes with more than 2 arcs;"
+ "split;Equal split (Okabe 2009) applied at nodes;");
+
+ kernel_opt = G_define_option();
+ kernel_opt->key = "kernel";
+ kernel_opt->type = TYPE_STRING;
+ kernel_opt->required = YES;
+ kernel_opt->description = _("Kernel function");
+ kernel_opt->options =
+ "uniform,triangular,epanechnikov,quartic,triweight,gaussian,cosine";
+ kernel_opt->answer = "gaussian";
+
flag_o = G_define_flag();
flag_o->key = 'o';
flag_o->description =
@@ -170,11 +186,16 @@
flag_q->description =
_("Only calculate optimal standard deviation and exit (no map is written)");
- /* please, remove before GRASS 7 released */
- flag_v = G_define_flag();
- flag_v->key = 'v';
- flag_v->description = _("Run verbosely");
+ flag_normalize = G_define_flag();
+ flag_normalize->key = 'n';
+ flag_normalize->description =
+ _("In network mode, normalize values by sum of density multiplied by length of each segment. Integral over the output map then gives 1.0 * mult");
+ flag_multiply = G_define_flag();
+ flag_multiply->key = 'm';
+ flag_multiply->description =
+ _("In network mode, multiply the result by number of input points.");
+
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
@@ -185,14 +206,43 @@
netmax = atof(netmax_opt->answer);
multip = atof(multip_opt->answer);
- /* please, remove before GRASS 7 released */
- if (flag_v->answer) {
- putenv("GRASS_VERBOSE=3");
- G_warning(_("The '-v' flag is superseded and will be removed "
- "in future. Please use '--verbose' instead."));
+ if (strcmp(node_opt->answer, "none") == 0)
+ node_method = NODE_NONE;
+ else if (strcmp(node_opt->answer, "split") == 0)
+ node_method = NODE_EQUAL_SPLIT;
+ else
+ G_fatal_error(_("Unknown node method"));
+
+ kernel_function = KERNEL_GAUSSIAN;
+ if (strcmp(kernel_opt->answer, "uniform") == 0)
+ kernel_function = KERNEL_UNIFORM;
+ else if (strcmp(kernel_opt->answer, "triangular") == 0)
+ kernel_function = KERNEL_TRIANGULAR;
+ else if (strcmp(kernel_opt->answer, "epanechnikov") == 0)
+ kernel_function = KERNEL_EPANECHNIKOV;
+ else if (strcmp(kernel_opt->answer, "quartic") == 0)
+ kernel_function = KERNEL_QUARTIC;
+ else if (strcmp(kernel_opt->answer, "triweight") == 0)
+ kernel_function = KERNEL_TRIWEIGHT;
+ else if (strcmp(kernel_opt->answer, "gaussian") == 0)
+ kernel_function = KERNEL_GAUSSIAN;
+ else if (strcmp(kernel_opt->answer, "cosine") == 0)
+ kernel_function = KERNEL_COSINE;
+ else
+ G_fatal_error(_("Unknown kernel function"));
+
+ if (flag_o->answer) {
+ if (net_opt->answer) {
+ if (node_method != NODE_NONE ||
+ kernel_function != KERNEL_GAUSSIAN) {
+ G_fatal_error(_("Optimal standard deviation calculation is supported only for node method 'none' and kernel function 'gaussian'."));
+ }
+ }
+ else if (kernel_function != KERNEL_GAUSSIAN) {
+ G_fatal_error(_("Optimal standard deviation calculation is supported only for kernel function 'gaussian'."));
+ }
}
-
if (flag_q->answer) {
flag_o->answer = 1;
}
@@ -219,7 +269,6 @@
if (net_opt->answer) {
int nlines, line;
- int notreachable = 0;
struct line_pnts *Points;
Points = Vect_new_line_struct();
@@ -282,8 +331,8 @@
/* valutazione distanza ottimale */
if (flag_o->answer) {
/* Note: sigmaOptimal calculates using ALL points (also those outside the region) */
- G_message(_("Automatic choose of smoothing parameter (standard deviation), maximum possible "
- "value of standard deviation is was set to %f"), sigma);
+ G_message(_("Automatic choice of smoothing parameter (standard deviation), maximum possible "
+ "value of standard deviation is set to %f"), sigma);
/* maximum distance 4*sigma (3.9*sigma ~ 1.0000), keep it small, otherwise it takes
* too much points and calculation on network becomes slow */
@@ -334,12 +383,16 @@
}
term = 1. / (pow(sigma, dimension) * pow((2. * M_PI), dimension / 2.));
- dmax = sigma * 4.;
+ dmax = sigma;
+ if (kernel_function == KERNEL_GAUSSIAN)
+ dmax = sigma * 4.;
+
if (net) {
int line, nlines;
struct line_pnts *Points, *SPoints;
struct line_cats *SCats;
+ double total = 0.0;
G_message(_("\nWriting output vector map using smooth parameter=%f."),
sigma);
@@ -379,7 +432,8 @@
offset1, x, y);
compute_net_distance(x, y, &In, &Net, netmax, sigma, term,
- &gaussian, dmax);
+ &gaussian, dmax, node_method,
+ kernel_function);
gaussian *= multip;
if (gaussian > gausmax)
gausmax = gaussian;
@@ -403,11 +457,44 @@
Vect_cat_set(SCats, 1, (int)gaussian);
Vect_write_line(&Out, GV_LINE, SPoints, SCats);
+
+ total += length * gaussian;
}
}
G_percent(line, nlines, 1);
}
+ if (flag_normalize->answer || flag_multiply->answer) {
+ double m = multip;
+
+ if (flag_normalize->answer) {
+ m /= total;
+ }
+ if (flag_multiply->answer) {
+ m *= (Vect_get_num_primitives(&In, GV_POINT) - notreachable);
+ }
+
+ Vect_build(&Out);
+
+ gausmax = 0.0;
+ nlines = Vect_get_num_lines(&Out);
+ for (line = 1; line <= nlines; line++) {
+ int cat;
+ double gaussian;
+
+ Vect_read_line(&Out, SPoints, SCats, line);
+
+ Vect_cat_get(SCats, 1, &cat);
+ gaussian = m * cat;
+ Vect_reset_cats(SCats);
+ Vect_cat_set(SCats, 1, (int)gaussian);
+ Vect_rewrite_line(&Out, line, GV_LINE, SPoints, SCats);
+ if (gaussian > gausmax)
+ gausmax = gaussian;
+ }
+ Vect_build_partial(&Out, GV_BUILD_NONE); /* to force rebuild */
+ }
+
Vect_close(&Net);
Vect_build(&Out);
@@ -436,7 +523,10 @@
N = G_row_to_northing(row + 0.5, &window);
E = G_col_to_easting(col + 0.5, &window);
- compute_distance(N, E, &In, sigma, term, &gaussian, dmax);
+ /* compute_distance(N, E, &In, sigma, term, &gaussian, dmax); */
+ compute_distance(N, E, &In, sigma, term, &gaussian, dmax,
+ kernel_function);
+
output_cell[col] = multip * gaussian;
if (gaussian > gausmax)
gausmax = gaussian;
@@ -447,7 +537,7 @@
G_close_cell(fdout);
}
- G_message(_("Maximum value in output: %e."), gausmax);
+ G_message(_("Maximum value in output: %e."), multip * gausmax);
Vect_close(&In);
@@ -597,66 +687,125 @@
return (kk);
}
+/* get number of arcs for a node */
+int count_node_arcs(struct Map_info *Map, int node)
+{
+ int i, n, line, type;
+ int count = 0;
+ n = Vect_get_node_n_lines(Map, node);
+ for (i = 0; i < n; i++) {
+ line = Vect_get_node_line(Map, node, i);
+ type = Vect_read_line(Map, NULL, NULL, abs(line));
+ if (type & GV_LINES)
+ count++;
+ }
+ return count;
+}
+
/* 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 term, double *gaussian, double dmax)
+ double term, double *gaussian, double dmax,
+ int node_method, int kernel_function)
{
int i;
- double dist;
+ double dist, kernel;
static struct line_pnts *Points = NULL;
+ static struct line_pnts *FPoints = NULL;
BOUND_BOX box;
- static struct ilist *List = NULL;
+ static struct ilist *PointsList = NULL;
+ static struct ilist *NodesList = NULL;
if (!Points)
Points = Vect_new_line_struct();
- if (!List)
- List = Vect_new_list();
+ if (!PointsList)
+ PointsList = Vect_new_list();
+ if (node_method == NODE_EQUAL_SPLIT) {
+ if (!NodesList)
+ NodesList = Vect_new_list();
+
+ if (!FPoints)
+ FPoints = Vect_new_line_struct();
+ }
+
*gaussian = .0;
/* The network is usually much bigger than dmax and to calculate shortest path is slow
- * -> use spatial index to select points */
- box.E = x + dmax;
- box.W = x - dmax;
- box.N = y + dmax;
- box.S = y - dmax;
+ * -> use spatial index to select points
+ * enlarge the box by netmax (max permitted distance between a point and net) */
+ box.E = x + dmax + netmax;
+ box.W = x - dmax - netmax;
+ box.N = y + dmax + netmax;
+ box.S = y - dmax - netmax;
box.T = PORT_DOUBLE_MAX;
box.B = -PORT_DOUBLE_MAX;
- Vect_select_lines_by_box(In, &box, GV_POINT, List);
- G_debug(3, " %d points in box", List->n_values);
+ Vect_select_lines_by_box(In, &box, GV_POINT, PointsList);
+ G_debug(3, " %d points in box", PointsList->n_values);
- for (i = 0; i < List->n_values; i++) {
+ for (i = 0; i < PointsList->n_values; i++) {
int line, ret;
- line = List->value[i];
+ line = PointsList->value[i];
Vect_read_line(In, Points, NULL, line);
G_debug(3, " SP: %f %f -> %f %f", x, y, Points->x[0], Points->y[0]);
- ret =
- Vect_net_shortest_path_coor(Net, x, y, 0.0, Points->x[0],
- Points->y[0], 0.0, netmax, netmax,
- &dist, NULL, NULL, NULL, NULL, NULL,
- NULL);
+ /*ret = Vect_net_shortest_path_coor(Net, x, y, 0.0, Points->x[0], */
+ /*Points->y[0], 0.0, netmax, netmax, */
+ /*&dist, NULL, NULL, NULL, NULL, NULL, */
+ /*NULL); */
+ ret = Vect_net_shortest_path_coor2(Net,
+ Points->x[0], Points->y[0], 0.0,
+ x, y, 0.0, netmax, 1.0,
+ &dist, NULL,
+ NULL, NodesList, FPoints, NULL,
+ NULL, NULL);
if (ret == 0) {
G_debug(3, "not reachable");
continue; /* Not reachable */
}
- if (dist <= dmax)
- *gaussian += gaussianKernel(dist / sigma, term);
+ /* if (dist <= dmax)
+ *gaussian += gaussianKernel(dist / sigma, term); */
+ if (dist > dmax)
+ continue;
+
+ /* kernel = gaussianKernel(dist / sigma, term); */
+ kernel = kernelFunction(kernel_function, 1, sigma, dist);
+
+ if (node_method == NODE_EQUAL_SPLIT) {
+ int j, node;
+ double ndiv = 1.;
+ int start = 0;
+
+ /* Count the nodes and arcs on path (n1-1)*(n2-1)* ... (ns-1) */
+
+ for (j = start; j < NodesList->n_values; j++) {
+ node = NodesList->value[j];
+
+ /* Divide into 2/n if point falls on a node */
+ if (j == 0 && FPoints->n_points < 3) {
+ ndiv *= count_node_arcs(Net, node) / 2.;
+ }
+ else {
+ ndiv *= count_node_arcs(Net, node) - 1;
+ }
+ }
+ kernel /= ndiv;
+ }
+ *gaussian += kernel;
G_debug(3, " dist = %f gaussian = %f", dist, *gaussian);
}
}
void compute_distance(double N, double E, struct Map_info *In,
double sigma, double term, double *gaussian,
- double dmax)
+ double dmax, int kernel_function)
{
int line, nlines;
double a[2], b[2];
@@ -698,7 +847,9 @@
dist = euclidean_distance(a, b, 2);
if (dist <= dmax)
- *gaussian += gaussianKernel(dist / sigma, term);
+ /* *gaussian += gaussianKernel(dist / sigma, term); */
+ *gaussian += kernelFunction(kernel_function, 2, sigma, dist);
+
}
}
More information about the grass-commit
mailing list