/****************************************************************************
*
* 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 kernel function or
*               optionally generates a vector density map on vector network 
*               with a 1D kernel
* 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
*   	    	for details.
*
*****************************************************************************/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <string.h>
#include <grass/gis.h>
#include <grass/raster.h>
#include <grass/glocale.h>
#include <grass/gmath.h>
#include <grass/vector.h>
#include "global.h"


static int ndists;		/* number of distances in dists */
static double *dists;		/* array of all distances < dmax */
static int npoints;
int net = 0;
static double dimension = 2.;


/* define score function L(window size) */
double L(double smooth)
{
    int ii;
    double resL, n, term;

    n = npoints;
    resL = 0.;
    term = 1. / pow((2. * M_PI), dimension / 2.);

    for (ii = 0; ii < ndists; ii++) {
	/*    resL+= gaussianFunction(dists[ii]/smooth,2.,dimension) - 2. * gaussianKernel(dists[ii]/smooth,term); */
	resL +=
	    gaussianFunction(dists[ii] / smooth, 2.,
			     dimension) -
	    2. * gaussianFunction(dists[ii] / smooth, 1., dimension);
    }

    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. * gaussianKernel(0.,term)) ) + (2./(n*pow(smooth,dimension)))*gaussianKernel(0.,term);   */
    G_debug(3, "smooth = %e resL = %e", smooth, resL);
    G_message(_("\tScore Value=%f\tsmoothing parameter (standard deviation)=%f"),
	      resL, smooth);

    return (resL);
}


int main(int argc, char **argv)
{
    struct Option *in_opt, *net_opt, *out_opt, *net_out_opt;
    struct Option *radius_opt, *dsize_opt, *segmax_opt, *netmax_opt,
	*multip_opt, *node_opt, *kernel_opt;
    struct Flag *flag_o, *flag_q, *flag_normalize, *flag_multiply;
    char *desc;

    struct Map_info In, Net, Out;
    int fdout = -1, maskfd = -1;
    int node_method, kernel_function;
    int row, col;
    struct Cell_head window;
    double gaussian;
    double N, E;
    CELL *mask = NULL;
    DCELL *output_cell = NULL;
    double sigma, dmax, segmax, netmax, multip;
    char *tmpstr1, *tmpstr2;
    struct History history;

    double **coordinate;
    double sigmaOptimal;
    struct GModule *module;
    double dsize;
    double term = 0;

    double gausmax = 0;
    int notreachable = 0;

    /* Initialize the GIS calls */
    G_gisinit(argv[0]);

    module = G_define_module();
    G_add_keyword(_("vector"));
    G_add_keyword(_("kernel density"));
    G_add_keyword(_("point density"));
    G_add_keyword(_("heatmap"));
    G_add_keyword(_("hotspot"));
    module->label =
        _("Generates a raster density map from vector points map.");
    module->description = _("Density is computed using a moving kernel. "
                            "Optionally generates a vector density map on a vector network.");

    in_opt = G_define_standard_option(G_OPT_V_INPUT);
    in_opt->label = _("Name of input vector map with training points");
    in_opt->description = NULL;
    in_opt->guisection = _("Basic");
    
    net_opt = G_define_standard_option(G_OPT_V_INPUT);
    net_opt->key = "net";
    net_opt->label = _("Name of input network vector map");
    net_opt->description = NULL;
    net_opt->required = NO;
    net_opt->guisection = _("Network");

    out_opt = G_define_standard_option(G_OPT_R_OUTPUT);
    out_opt->key = "output";
    out_opt->required = NO;
    out_opt->label = _("Name for output raster map");
    out_opt->description = NULL;
    out_opt->guisection = _("Basic");

    net_out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
    net_out_opt->key = "net_output";
    net_out_opt->required = NO;
    net_out_opt->label = _("Name for output vector density map");
    net_out_opt->description = _("Outputs vector map if network map is given");
    net_out_opt->guisection = _("Network");

    radius_opt = G_define_option();
    radius_opt->key = "radius";
    radius_opt->type = TYPE_DOUBLE;
    radius_opt->required = YES;
    radius_opt->description = _("Kernel radius in map units");
    radius_opt->guisection = _("Basic");

    dsize_opt = G_define_option();
    dsize_opt->key = "dsize";
    dsize_opt->type = TYPE_DOUBLE;
    dsize_opt->required = NO;
    dsize_opt->description = _("Discretization error in map units");
    dsize_opt->answer = "0.";

    segmax_opt = G_define_option();
    segmax_opt->key = "segmax";
    segmax_opt->type = TYPE_DOUBLE;
    segmax_opt->required = NO;
    segmax_opt->description = _("Maximum length of segment on network");
    segmax_opt->answer = "100.";
    segmax_opt->guisection = _("Network");

    netmax_opt = G_define_option();
    netmax_opt->key = "distmax";
    netmax_opt->type = TYPE_DOUBLE;
    netmax_opt->required = NO;
    netmax_opt->description = _("Maximum distance from point to network");
    netmax_opt->answer = "100.";
    netmax_opt->guisection = _("Network");

    multip_opt = G_define_option();
    multip_opt->key = "multiplier";
    multip_opt->type = TYPE_DOUBLE;
    multip_opt->required = NO;
    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 = NO;
    node_opt->description = _("Node method");
    node_opt->options = "none,split";
    node_opt->answer = "none";
    desc = NULL;
    G_asprintf(&desc,
	       "none;%s;split;%s",
	       _("No method applied at nodes with more than 2 arcs"),
	       _("Equal split (Okabe 2009) applied at nodes"));
    node_opt->descriptions = desc;
    node_opt->guisection = _("Network");

    kernel_opt = G_define_option();
    kernel_opt->key = "kernel";
    kernel_opt->type = TYPE_STRING;
    kernel_opt->required = NO;
    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 =
	_("Try to calculate an optimal radius with given 'radius' taken as maximum (experimental)");

    flag_q = G_define_flag();
    flag_q->key = 'q';
    flag_q->description =
	_("Only calculate optimal radius and exit (no map is written)");

    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 * multiplier");
    flag_normalize->guisection = _("Network");

    flag_multiply = G_define_flag();
    flag_multiply->key = 'm';
    flag_multiply->description =
	_("In network mode, multiply the result by number of input points");
    flag_multiply->guisection = _("Network");

    G_option_required(out_opt, net_out_opt, NULL);
    G_option_exclusive(out_opt, net_out_opt, NULL);
    /* TODO: this should be activated for GRASS 8
    G_option_requires(net_opt, net_out_opt, NULL);
    */
    if (G_parser(argc, argv))
	exit(EXIT_FAILURE);


    if (net_opt->answer && out_opt->answer) {
	G_warning(_("Use option net_output if you compute network density. "
		"Name provided in option output will be used for net_output."));
	net_out_opt->answer = out_opt->answer;
	out_opt->answer = NULL;
    }

    /*read options */
    dmax = atof(radius_opt->answer);
    sigma = dmax;
    dsize = atof(dsize_opt->answer);
    segmax = atof(segmax_opt->answer);
    netmax = atof(netmax_opt->answer);
    multip = atof(multip_opt->answer);

    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;
    }

    if (net_opt->answer) {
	Vect_check_input_output_name(in_opt->answer, net_out_opt->answer,
				     G_FATAL_EXIT);
	Vect_check_input_output_name(net_opt->answer, net_out_opt->answer,
				     G_FATAL_EXIT);
    }

    G_get_window(&window);

    G_verbose_message(_("Standard deviation: %f"), sigma);
    G_asprintf(&tmpstr1, n_("%d row", "%d rows", window.rows), window.rows);
    G_asprintf(&tmpstr2, n_("%d column", "%d columns", window.cols), window.cols);
    /* GTC First argument is resolution, second - number of rows as a text, third - number of columns as a text. */
    if (G_verbose() > G_verbose_std()) {
       G_verbose_message(_("Output raster map: resolution: %f\t%s\t%s"), /* mhh, this assumes square pixels */
                      window.ew_res, tmpstr1, tmpstr2);
    } else {
       G_message(_("Output raster map resolution: %f"), window.ew_res); /* mhh, this assumes square pixels */
    }

    G_free(tmpstr1);
    G_free(tmpstr2); 
    
    /* Open input vector */
    Vect_set_open_level(2);
    if (Vect_open_old(&In, in_opt->answer, "") < 0)
	G_fatal_error(_("Unable to open vector map <%s>"), in_opt->answer);

    if (net_opt->answer) {
	int nlines, line;
	struct line_pnts *Points;

	Points = Vect_new_line_struct();
	net = 1;
	dimension = 1.;
	/* Open input network */
	Vect_set_open_level(2);

	if (Vect_open_old(&Net, net_opt->answer, "") < 0)
	    G_fatal_error(_("Unable to open vector map <%s>"), net_opt->answer);

	Vect_net_build_graph(&Net, GV_LINES, 0, 0, NULL, NULL, NULL, 0, 0);

	if (!flag_q->answer) {
	    if (Vect_open_new(&Out, net_out_opt->answer, 0) < 0)
		G_fatal_error(_("Unable to create vector map <%s>"),
				net_out_opt->answer);
	    Vect_hist_command(&Out);
	}

	/* verify not reachable points */
	nlines = Vect_get_num_lines(&In);
	for (line = 1; line <= nlines; line++) {
	    int ltype;

	    ltype = Vect_read_line(&In, Points, NULL, line);
	    if (!(ltype & GV_POINTS))
		continue;
	    if (Vect_find_line
		(&Net, Points->x[0], Points->y[0], 0.0, GV_LINES, netmax, 0,
		 0) == 0)
		notreachable++;
	}

	if (notreachable > 0)
	    G_warning(n_("%d point outside threshold",
                         "%d points outside threshold",
                         notreachable), notreachable);
    }
    else {
	/* check and open the name of output map */
	if (!flag_q->answer) {
	    fdout = Rast_open_new(out_opt->answer, DCELL_TYPE);

	    /* open mask file */
	    if ((maskfd = Rast_maskfd()) >= 0)
		mask = Rast_allocate_c_buf();
	    else
		mask = NULL;

	    /* allocate output raster */
	    output_cell = Rast_allocate_buf(DCELL_TYPE);
	}
    }

    /* valutazione distanza ottimale */
    if (flag_o->answer) {
	/* Note: sigmaOptimal calculates using ALL points (also those outside the region) */
	G_message(_("Automatic choice of smoothing parameter (radius), maximum possible "
		   "value of radius 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 */
	dmax = 4 * sigma;	/* used as maximum value */

	G_message(_("Using maximum distance between points: %f"), dmax);

	if (net_opt->answer) {
	    npoints = Vect_get_num_primitives(&In, GV_POINTS);
	    /* Warning: each distance is registered twice (both directions) */
	    ndists =
		compute_all_net_distances(&In, &Net, netmax, &dists, dmax);
	}
	else {
	    /* Read points */
	    npoints = read_points(&In, &coordinate, dsize);
	    ndists = compute_all_distances(coordinate, &dists, npoints, dmax);
	}

	G_message(_("Number of input points: %d."), npoints);
	G_message(n_("%d distance read from the map.",
                     "%d distances read from the map.",
                     ndists), ndists);

	if (ndists == 0)
	    G_fatal_error(_("Distances between all points are beyond %e (4 * "
			    "standard deviation), unable to calculate optimal value."),
			  dmax);

	/*  double iii;
	   for ( iii = 1.; iii <= 10000; iii++){
	   fprintf(stderr,"i=%f v=%.16f \n",iii,R(iii));
	   } */

	/* sigma is used in brent as maximum possible value for sigmaOptimal */
	sigmaOptimal = brent_iterate(L, 0.0, sigma, 1000);
	G_message(_("Optimal smoothing parameter (standard deviation): %f."),
		  sigmaOptimal);

	/* Reset sigma to calculated optimal value */
	sigma = sigmaOptimal;

	if (flag_q->answer) {
	    Vect_close(&In);
	    if (net_opt->answer)
		Vect_close(&Net);

	    exit(EXIT_SUCCESS);
	}
    }

    if (kernel_function == KERNEL_GAUSSIAN)
	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;
	struct line_cats *SCats;
	double total = 0.0;

	G_verbose_message(_("Writing output vector map using smooth parameter %f"),
                          sigma);
	G_verbose_message(_("Normalising factor %f"),
                          1. / gaussianFunction(sigma / 4., sigma, dimension));

	/* Divide lines to segments and calculate gaussian for center of each segment */
	Points = Vect_new_line_struct();
	SPoints = Vect_new_line_struct();
	SCats = Vect_new_cats_struct();

	nlines = Vect_get_num_lines(&Net);
	G_debug(3, "net nlines = %d", nlines);

	for (line = 1; line <= nlines; line++) {
	    int seg, nseg, ltype;
	    double llength, length, x, y;

            G_percent(line, nlines, 5);
	    ltype = Vect_read_line(&Net, Points, NULL, line);
	    if (!(ltype & GV_LINES))
		continue;

	    llength = Vect_line_length(Points);
	    nseg = (int)(1 + llength / segmax);
	    length = llength / nseg;

	    G_debug(3, "net line = %d, nseg = %d, seg length = %f", line,
		    nseg, length);

	    for (seg = 0; seg < nseg; seg++) {
		double offset1, offset2;

		offset1 = (seg + 0.5) * length;
		Vect_point_on_line(Points, offset1, &x, &y, NULL, NULL, NULL);

		G_debug(3, "  segment = %d, offset = %f, xy = %f %f", seg,
			offset1, x, y);

		compute_net_distance(x, y, &In, &Net, netmax, sigma, term,
				     &gaussian, dmax, node_method);
		gaussian *= multip;
		if (gaussian > gausmax)
		    gausmax = gaussian;

		G_debug(3, "  gaussian = %f", gaussian);

		/* Write segment */
		if (gaussian > 0) {
		    offset1 = seg * length;
		    offset2 = (seg + 1) * length;
		    if (offset2 > llength)
			offset2 = llength;
		    Vect_line_segment(Points, offset1, offset2, SPoints);

		    /* TODO!!! remove later
		       if ( SPoints->n_points > 0 ) 
		       Vect_append_point( SPoints, SPoints->x[SPoints->n_points-1], 
		       SPoints->y[SPoints->n_points-1], 0 );
		     */
		    Vect_reset_cats(SCats);
		    Vect_cat_set(SCats, 1, (int)gaussian);

		    Vect_write_line(&Out, GV_LINE, SPoints, SCats);

		    total += length * gaussian;
		}
	    }
	}

	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);
	Vect_close(&Out);
    }
    else {
	/* spatial index handling, borrowed from lib/vector/Vlib/find.c */
	struct bound_box box;
	struct boxlist *NList = Vect_new_boxlist(1);

	G_verbose_message(_("Writing output raster map using smooth parameter %f"),
                          sigma);
	G_verbose_message(_("Normalising factor %f"),
                          1. / gaussianFunction(sigma / 4., sigma, dimension));

	for (row = 0; row < window.rows; row++) {
	    G_percent(row, window.rows, 2);
	    if (mask)
		Rast_get_c_row(maskfd, mask, row);

	    for (col = 0; col < window.cols; col++) {
		/* don't interpolate outside of the mask */
		if (mask && Rast_is_c_null_value(&mask[col])) {
		    Rast_set_d_null_value(&output_cell[col], 1);
		    continue;
		}

		N = Rast_row_to_northing(row + 0.5, &window);
		E = Rast_col_to_easting(col + 0.5, &window);

		if ((col & 31) == 0) {

		    /* create bounding box 32x2*dmax size from the current cell center */
		    box.N = N + dmax;
		    box.S = N - dmax;
		    box.E = E + dmax + 32 * window.ew_res;
		    box.W = E - dmax;
		    box.T = HUGE_VAL;
		    box.B = -HUGE_VAL;

		    Vect_select_lines_by_box(&In, &box, GV_POINT, NList);		    
		}
		box.N = N + dmax;
		box.S = N - dmax;
		box.E = E + dmax;
		box.W = E - dmax;
		box.T = HUGE_VAL;
		box.B = -HUGE_VAL;
		/* compute_distance(N, E, &In, sigma, term, &gaussian, dmax); */
		compute_distance(N, E, sigma, term, &gaussian, dmax, &box, NList);

		output_cell[col] = multip * gaussian;
		if (gaussian > gausmax)
		    gausmax = gaussian;
	    }
	    Rast_put_row(fdout, output_cell, DCELL_TYPE);
	}
        G_percent(1, 1, 1);
        
	Rast_close(fdout);

    Rast_short_history(out_opt->answer, "raster", &history);
    Rast_command_history(&history);
    Rast_write_history(out_opt->answer, &history);
    }

    G_done_msg(_("Maximum value in output: %e."), multip * gausmax);

    Vect_close(&In);

    exit(EXIT_SUCCESS);
}


/* Read points to array return number of points */
int read_points(struct Map_info *In, double ***coordinate, double dsize)
{
    int line, nlines, npoints, ltype, i = 0;
    double **xySites;
    static struct line_pnts *Points = NULL;

    if (!Points)
	Points = Vect_new_line_struct();

    /* Allocate array of pointers */
    npoints = Vect_get_num_primitives(In, GV_POINT);
    xySites = (double **)G_calloc(npoints, sizeof(double *));

    nlines = Vect_get_num_lines(In);

    for (line = 1; line <= nlines; line++) {
	ltype = Vect_read_line(In, Points, NULL, line);
	if (!(ltype & GV_POINT))
	    continue;

	xySites[i] = (double *)G_calloc((size_t) 2, sizeof(double));

	xySites[i][0] = Points->x[0];
	xySites[i][1] = Points->y[0];
	i++;
    }

    *coordinate = xySites;

    return (npoints);
}


/* Calculate distances < dmax between all sites in coordinate 
 * Return: number of distances in dists */
double compute_all_distances(double **coordinate, double **dists, int n,
			     double dmax)
{
    int ii, jj, kk;
    size_t nn;

    nn = n * (n - 1) / 2;
    *dists = (double *)G_calloc(nn, sizeof(double));
    kk = 0;

    for (ii = 0; ii < n - 1; ii++) {
	for (jj = ii + 1; jj < n; jj++) {
	    double dist;

	    dist = euclidean_distance(coordinate[ii], coordinate[jj], 2);
	    G_debug(3, "dist = %f", dist);

	    if (dist <= dmax) {
		(*dists)[kk] = dist;
		kk++;
	    }
	}
    }

    return (kk);
}


/* Calculate distances < dmax between all sites in coordinate 
 * Return: number of distances in dists */
double compute_all_net_distances(struct Map_info *In, struct Map_info *Net,
				 double netmax, double **dists, double dmax)
{
    int nn, kk, nalines, aline;
    double dist;
    struct line_pnts *APoints, *BPoints;
    struct bound_box box;
    struct boxlist *List;

    APoints = Vect_new_line_struct();
    BPoints = Vect_new_line_struct();
    List = Vect_new_boxlist(0);

    nn = Vect_get_num_primitives(In, GV_POINTS);
    nn = nn * (nn - 1);
    *dists = (double *)G_calloc(nn, sizeof(double));
    kk = 0;

    nalines = Vect_get_num_lines(In);
    for (aline = 1; aline <= nalines; aline++) {
	int i, altype;

	G_debug(3, "  aline = %d", aline);

	altype = Vect_read_line(In, APoints, NULL, aline);
	if (!(altype & GV_POINTS))
	    continue;

	box.E = APoints->x[0] + dmax;
	box.W = APoints->x[0] - dmax;
	box.N = APoints->y[0] + dmax;
	box.S = APoints->y[0] - dmax;
	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);

	for (i = 0; i < List->n_values; i++) {
	    int bline, ret;

	    bline = List->id[i];

	    if (bline == aline)
		continue;

	    G_debug(3, "    bline = %d", bline);
	    Vect_read_line(In, BPoints, NULL, bline);

	    ret =
		Vect_net_shortest_path_coor(Net, APoints->x[0], APoints->y[0],
					    0.0, BPoints->x[0], BPoints->y[0],
					    0.0, netmax, netmax, &dist, NULL,
					    NULL, NULL, NULL, NULL, NULL, NULL);

	    G_debug(3, "  SP: %f %f -> %f %f", APoints->x[0], APoints->y[0],
		    BPoints->x[0], BPoints->y[0]);

	    if (ret == 0) {
		G_debug(3, "not reachable");
		continue;	/* Not reachable */
	    }

	    G_debug(3, "  dist = %f", dist);

	    if (dist <= dmax) {
		(*dists)[kk] = dist;
		kk++;
	    }
	    G_debug(3, "  kk = %d", kk);
	}
    }

    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_get_line_type(Map, abs(line));
	if (type & GV_LINES)
	    count++;
    }
    return count;
}

/* Compute Gaussian 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, int node_method)
{
    int i;
    double dist, kernel;
    static struct line_pnts *FPoints = NULL;
    struct bound_box box;
    static struct boxlist *PointsList = NULL;
    static struct ilist *NodesList = NULL;

    if (!PointsList)
	PointsList = Vect_new_boxlist(1);

    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
     * 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, PointsList);
    G_debug(3, "  %d points in box", PointsList->n_values);

    for (i = 0; i < PointsList->n_values; i++) {
	int line, ret;

	line = PointsList->id[i];

	G_debug(3, "  SP: %f %f -> %f %f", x, y, PointsList->box[i].E, PointsList->box[i].N);
	/*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, */
	/*NULL); */
	ret = Vect_net_shortest_path_coor(Net,
					   PointsList->box[i].E,
					   PointsList->box[i].N, 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)
	    continue;

	/* kernel = gaussianKernel(dist / sigma, term); */
	kernel = kernelFunction(term, 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, double sigma, double term,
                      double *gaussian, double dmax, struct bound_box *box,
		      struct boxlist *NList)
{
    int line, nlines;
    double a[2], b[2];
    double dist;

    a[0] = E;
    a[1] = N;

    /* number of lines within dmax box  */
    nlines = NList->n_values;

    *gaussian = .0;

    for (line = 0; line < nlines; line++) {

	b[0] = NList->box[line].E;
	b[1] = NList->box[line].N;

	if (b[0] <= box->E && b[0] >= box->W &&
	    b[1] <= box->N && b[1] >= box->S)  {
	    dist = euclidean_distance(a, b, 2);

	    if (dist <= dmax)
		/* *gaussian += gaussianKernel(dist / sigma, term); */
		*gaussian += kernelFunction(term, sigma, dist);
	}
    }
}