source: grass/trunk/vector/v.kernel/main.c

Last change on this file was 66269, checked in by neteler, 9 months ago

v.kernel: always show target raster resolution; fix dangling curly braces -Wdangling-else

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id
  • Property svn:mime-type set to text/x-csrc
File size: 25.7 KB
Line 
1
2/****************************************************************************
3*
4* MODULE:       v.kernel
5*
6* AUTHOR(S):    Stefano Menegon, ITC-irst, Trento, Italy
7*               Radim Blazek (additional kernel functions, network part)
8* PURPOSE:      Generates a raster density map from vector points data using
9*               a moving kernel function or
10*               optionally generates a vector density map on vector network
11*               with a 1D kernel
12* COPYRIGHT:    (C) 2004-2011 by the GRASS Development Team
13*
14*               This program is free software under the GNU General Public
15*               License (>=v2). Read the file COPYING that comes with GRASS
16*               for details.
17*
18*****************************************************************************/
19#include <math.h>
20#include <stdio.h>
21#include <stdlib.h>
22#include <float.h>
23#include <string.h>
24#include <grass/gis.h>
25#include <grass/raster.h>
26#include <grass/glocale.h>
27#include <grass/gmath.h>
28#include <grass/vector.h>
29#include "global.h"
30
31
32static int ndists;              /* number of distances in dists */
33static double *dists;           /* array of all distances < dmax */
34static int npoints;
35int net = 0;
36static double dimension = 2.;
37
38
39/* define score function L(window size) */
40double L(double smooth)
41{
42    int ii;
43    double resL, n, term;
44
45    n = npoints;
46    resL = 0.;
47    term = 1. / pow((2. * M_PI), dimension / 2.);
48
49    for (ii = 0; ii < ndists; ii++) {
50        /*    resL+= gaussianFunction(dists[ii]/smooth,2.,dimension) - 2. * gaussianKernel(dists[ii]/smooth,term); */
51        resL +=
52            gaussianFunction(dists[ii] / smooth, 2.,
53                             dimension) -
54            2. * gaussianFunction(dists[ii] / smooth, 1., dimension);
55    }
56
57    if (!net)
58        resL *= 2.;
59
60    resL = (1. / (pow(n, 2.) * pow(smooth, dimension))) *
61           (resL + n * (gaussianFunction(0., 2., dimension) -
62           2. * gaussianFunction(0., 1., dimension))) + 
63           (2. / (n * pow(smooth, dimension))) *
64           gaussianFunction(0., 1., dimension);
65
66    /* 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);   */
67    G_debug(3, "smooth = %e resL = %e", smooth, resL);
68    G_message(_("\tScore Value=%f\tsmoothing parameter (standard deviation)=%f"),
69              resL, smooth);
70
71    return (resL);
72}
73
74
75int main(int argc, char **argv)
76{
77    struct Option *in_opt, *net_opt, *out_opt;
78    struct Option *radius_opt, *dsize_opt, *segmax_opt, *netmax_opt,
79        *multip_opt, *node_opt, *kernel_opt;
80    struct Flag *flag_o, *flag_q, *flag_normalize, *flag_multiply;
81    char *desc;
82
83    struct Map_info In, Net, Out;
84    int overwrite;
85    int fdout = -1, maskfd = -1;
86    int node_method, kernel_function;
87    int row, col;
88    struct Cell_head window;
89    double gaussian;
90    double N, E;
91    CELL *mask = NULL;
92    DCELL *output_cell = NULL;
93    double sigma, dmax, segmax, netmax, multip;
94    char *tmpstr1, *tmpstr2;
95
96    double **coordinate;
97    double sigmaOptimal;
98    struct GModule *module;
99    double dsize;
100    double term = 0;
101
102    double gausmax = 0;
103    int notreachable = 0;
104
105    /* Initialize the GIS calls */
106    G_gisinit(argv[0]);
107
108    module = G_define_module();
109    G_add_keyword(_("vector"));
110    G_add_keyword(_("kernel density"));
111    G_add_keyword(_("point density"));
112    G_add_keyword(_("heatmap"));
113    G_add_keyword(_("hotspot"));
114    module->label =
115        _("Generates a raster density map from vector points map.");
116    module->description = _("Density is computed using a moving kernel. "
117                            "Optionally generates a vector density map on a vector network.");
118
119    in_opt = G_define_standard_option(G_OPT_V_INPUT);
120    in_opt->label = _("Name of input vector map with training points");
121    in_opt->description = NULL;
122   
123    net_opt = G_define_standard_option(G_OPT_V_INPUT);
124    net_opt->key = "net";
125    net_opt->label = _("Name of input network vector map");
126    net_opt->description = NULL;
127    net_opt->required = NO;
128    net_opt->guisection = _("Network");
129
130    out_opt = G_define_option();
131    out_opt->key = "output";
132    out_opt->type = TYPE_STRING;
133    out_opt->key_desc = "name";
134    out_opt->required = YES;
135    out_opt->label = _("Name for output raster/vector map");
136    out_opt->description = _("Outputs vector map if network map is given, otherwise raster map");
137
138    radius_opt = G_define_option();
139    radius_opt->key = "radius";
140    radius_opt->type = TYPE_DOUBLE;
141    radius_opt->required = YES;
142    radius_opt->description = _("Kernel radius in map units");
143
144    dsize_opt = G_define_option();
145    dsize_opt->key = "dsize";
146    dsize_opt->type = TYPE_DOUBLE;
147    dsize_opt->required = NO;
148    dsize_opt->description = _("Discretization error in map units");
149    dsize_opt->answer = "0.";
150
151    segmax_opt = G_define_option();
152    segmax_opt->key = "segmax";
153    segmax_opt->type = TYPE_DOUBLE;
154    segmax_opt->required = NO;
155    segmax_opt->description = _("Maximum length of segment on network");
156    segmax_opt->answer = "100.";
157    segmax_opt->guisection = _("Network");
158
159    netmax_opt = G_define_option();
160    netmax_opt->key = "distmax";
161    netmax_opt->type = TYPE_DOUBLE;
162    netmax_opt->required = NO;
163    netmax_opt->description = _("Maximum distance from point to network");
164    netmax_opt->answer = "100.";
165    netmax_opt->guisection = _("Network");
166
167    multip_opt = G_define_option();
168    multip_opt->key = "multiplier";
169    multip_opt->type = TYPE_DOUBLE;
170    multip_opt->required = NO;
171    multip_opt->description = _("Multiply the density result by this number");
172    multip_opt->answer = "1.";
173
174    node_opt = G_define_option();
175    node_opt->key = "node";
176    node_opt->type = TYPE_STRING;
177    node_opt->required = NO;
178    node_opt->description = _("Node method");
179    node_opt->options = "none,split";
180    node_opt->answer = "none";
181    desc = NULL;
182    G_asprintf(&desc,
183               "none;%s;split;%s",
184               _("No method applied at nodes with more than 2 arcs"),
185               _("Equal split (Okabe 2009) applied at nodes"));
186    node_opt->descriptions = desc;
187    node_opt->guisection = _("Network");
188
189    kernel_opt = G_define_option();
190    kernel_opt->key = "kernel";
191    kernel_opt->type = TYPE_STRING;
192    kernel_opt->required = NO;
193    kernel_opt->description = _("Kernel function");
194    kernel_opt->options =
195        "uniform,triangular,epanechnikov,quartic,triweight,gaussian,cosine";
196    kernel_opt->answer = "gaussian";
197
198    flag_o = G_define_flag();
199    flag_o->key = 'o';
200    flag_o->description =
201        _("Try to calculate an optimal radius with given 'radius' taken as maximum (experimental)");
202
203    flag_q = G_define_flag();
204    flag_q->key = 'q';
205    flag_q->description =
206        _("Only calculate optimal radius and exit (no map is written)");
207
208    flag_normalize = G_define_flag();
209    flag_normalize->key = 'n';
210    flag_normalize->description =
211        _("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");
212    flag_normalize->guisection = _("Network");
213
214    flag_multiply = G_define_flag();
215    flag_multiply->key = 'm';
216    flag_multiply->description =
217        _("In network mode, multiply the result by number of input points");
218    flag_multiply->guisection = _("Network");
219   
220    overwrite = G_check_overwrite(argc, argv);
221    if (G_parser(argc, argv))
222        exit(EXIT_FAILURE);
223
224    if (net_opt->answer) {
225        if (G_find_vector2(out_opt->answer, G_mapset())) {
226            if (overwrite)
227                G_warning(_("Vector map <%s> already exists and will be overwritten"),
228                          out_opt->answer);
229            else
230                G_fatal_error(_("Vector map <%s> already exists"),
231                              out_opt->answer);
232        }
233    } else {
234        if (G_find_raster(out_opt->answer, G_mapset())) {
235            if (overwrite)
236                G_warning(_("Raster map <%s> already exists and will be overwritten"),
237                          out_opt->answer);
238            else
239                G_fatal_error(_("Raster map <%s> already exists"),
240                              out_opt->answer);
241        }
242    }
243
244    /*read options */
245    dmax = atof(radius_opt->answer);
246    sigma = dmax;
247    dsize = atof(dsize_opt->answer);
248    segmax = atof(segmax_opt->answer);
249    netmax = atof(netmax_opt->answer);
250    multip = atof(multip_opt->answer);
251
252    if (strcmp(node_opt->answer, "none") == 0)
253        node_method = NODE_NONE;
254    else if (strcmp(node_opt->answer, "split") == 0)
255        node_method = NODE_EQUAL_SPLIT;
256    else
257        G_fatal_error(_("Unknown node method"));
258
259    kernel_function = KERNEL_GAUSSIAN;
260    if (strcmp(kernel_opt->answer, "uniform") == 0)
261        kernel_function = KERNEL_UNIFORM;
262    else if (strcmp(kernel_opt->answer, "triangular") == 0)
263        kernel_function = KERNEL_TRIANGULAR;
264    else if (strcmp(kernel_opt->answer, "epanechnikov") == 0)
265        kernel_function = KERNEL_EPANECHNIKOV;
266    else if (strcmp(kernel_opt->answer, "quartic") == 0)
267        kernel_function = KERNEL_QUARTIC;
268    else if (strcmp(kernel_opt->answer, "triweight") == 0)
269        kernel_function = KERNEL_TRIWEIGHT;
270    else if (strcmp(kernel_opt->answer, "gaussian") == 0)
271        kernel_function = KERNEL_GAUSSIAN;
272    else if (strcmp(kernel_opt->answer, "cosine") == 0)
273        kernel_function = KERNEL_COSINE;
274    else
275        G_fatal_error(_("Unknown kernel function"));
276
277    if (flag_o->answer) {
278        if (net_opt->answer) {
279            if (node_method != NODE_NONE ||
280                kernel_function != KERNEL_GAUSSIAN) {
281                G_fatal_error(_("Optimal standard deviation calculation is supported only for node method 'none' and kernel function 'gaussian'."));
282            }
283        }
284        else if (kernel_function != KERNEL_GAUSSIAN) {
285            G_fatal_error(_("Optimal standard deviation calculation is supported only for kernel function 'gaussian'."));
286        }
287    }
288
289    if (flag_q->answer) {
290        flag_o->answer = 1;
291    }
292
293    if (net_opt->answer) {
294        Vect_check_input_output_name(in_opt->answer, out_opt->answer,
295                                     G_FATAL_EXIT);
296        Vect_check_input_output_name(net_opt->answer, out_opt->answer,
297                                     G_FATAL_EXIT);
298    }
299
300    G_get_window(&window);
301
302    G_verbose_message(_("Standard deviation: %f"), sigma);
303    G_asprintf(&tmpstr1, n_("%d row", "%d rows", window.rows), window.rows);
304    G_asprintf(&tmpstr2, n_("%d column", "%d columns", window.cols), window.cols);
305    /* GTC First argument is resolution, second - number of rows as a text, third - number of columns as a text. */
306    if (G_verbose() > G_verbose_std()) {
307       G_verbose_message(_("Output raster map: resolution: %f\t%s\t%s"), /* mhh, this assumes square pixels */
308                      window.ew_res, tmpstr1, tmpstr2);
309    } else {
310       G_message(_("Output raster map resolution: %f"), window.ew_res); /* mhh, this assumes square pixels */
311    }
312
313    G_free(tmpstr1);
314    G_free(tmpstr2); 
315   
316    /* Open input vector */
317    Vect_set_open_level(2);
318    if (Vect_open_old(&In, in_opt->answer, "") < 0)
319        G_fatal_error(_("Unable to open vector map <%s>"), in_opt->answer);
320
321    if (net_opt->answer) {
322        int nlines, line;
323        struct line_pnts *Points;
324
325        Points = Vect_new_line_struct();
326        net = 1;
327        dimension = 1.;
328        /* Open input network */
329        Vect_set_open_level(2);
330
331        if (Vect_open_old(&Net, net_opt->answer, "") < 0)
332            G_fatal_error(_("Unable to open vector map <%s>"), net_opt->answer);
333
334        Vect_net_build_graph(&Net, GV_LINES, 0, 0, NULL, NULL, NULL, 0, 0);
335
336        if (!flag_q->answer) {
337            if (Vect_open_new(&Out, out_opt->answer, 0) < 0)
338                G_fatal_error(_("Unable to create vector map <%s>"),
339                                out_opt->answer);
340            Vect_hist_command(&Out);
341        }
342
343        /* verify not reachable points */
344        nlines = Vect_get_num_lines(&In);
345        for (line = 1; line <= nlines; line++) {
346            int ltype;
347
348            ltype = Vect_read_line(&In, Points, NULL, line);
349            if (!(ltype & GV_POINTS))
350                continue;
351            if (Vect_find_line
352                (&Net, Points->x[0], Points->y[0], 0.0, GV_LINES, netmax, 0,
353                 0) == 0)
354                notreachable++;
355        }
356
357        if (notreachable > 0)
358            G_warning(n_("%d point outside threshold",
359                         "%d points outside threshold",
360                         notreachable), notreachable);
361    }
362    else {
363        /* check and open the name of output map */
364        if (!flag_q->answer) {
365            fdout = Rast_open_new(out_opt->answer, DCELL_TYPE);
366
367            /* open mask file */
368            if ((maskfd = Rast_maskfd()) >= 0)
369                mask = Rast_allocate_c_buf();
370            else
371                mask = NULL;
372
373            /* allocate output raster */
374            output_cell = Rast_allocate_buf(DCELL_TYPE);
375        }
376    }
377
378    /* valutazione distanza ottimale */
379    if (flag_o->answer) {
380        /* Note: sigmaOptimal calculates using ALL points (also those outside the region) */
381        G_message(_("Automatic choice of smoothing parameter (radius), maximum possible "
382                   "value of radius is set to %f"), sigma);
383
384        /* maximum distance 4*sigma (3.9*sigma ~ 1.0000), keep it small, otherwise it takes
385         * too much points and calculation on network becomes slow */
386        dmax = 4 * sigma;       /* used as maximum value */
387
388        G_message(_("Using maximum distance between points: %f"), dmax);
389
390        if (net_opt->answer) {
391            npoints = Vect_get_num_primitives(&In, GV_POINTS);
392            /* Warning: each distance is registered twice (both directions) */
393            ndists =
394                compute_all_net_distances(&In, &Net, netmax, &dists, dmax);
395        }
396        else {
397            /* Read points */
398            npoints = read_points(&In, &coordinate, dsize);
399            ndists = compute_all_distances(coordinate, &dists, npoints, dmax);
400        }
401
402        G_message(_("Number of input points: %d."), npoints);
403        G_message(n_("%d distance read from the map.",
404                     "%d distances read from the map.",
405                     ndists), ndists);
406
407        if (ndists == 0)
408            G_fatal_error(_("Distances between all points are beyond %e (4 * "
409                            "standard deviation), unable to calculate optimal value."),
410                          dmax);
411
412        /*  double iii;
413           for ( iii = 1.; iii <= 10000; iii++){
414           fprintf(stderr,"i=%f v=%.16f \n",iii,R(iii));
415           } */
416
417        /* sigma is used in brent as maximum possible value for sigmaOptimal */
418        sigmaOptimal = brent_iterate(L, 0.0, sigma, 1000);
419        G_message(_("Optimal smoothing parameter (standard deviation): %f."),
420                  sigmaOptimal);
421
422        /* Reset sigma to calculated optimal value */
423        sigma = sigmaOptimal;
424
425        if (flag_q->answer) {
426            Vect_close(&In);
427            if (net_opt->answer)
428                Vect_close(&Net);
429
430            exit(EXIT_SUCCESS);
431        }
432    }
433
434    if (kernel_function == KERNEL_GAUSSIAN)
435        sigma /= 4.;
436
437    if (net_opt->answer) {
438        setKernelFunction(kernel_function, 1, sigma, &term);
439    }
440    else {
441        setKernelFunction(kernel_function, 2, sigma, &term);
442    }
443
444    if (net) {
445        int line, nlines;
446        struct line_pnts *Points, *SPoints;
447        struct line_cats *SCats;
448        double total = 0.0;
449
450        G_verbose_message(_("Writing output vector map using smooth parameter %f"),
451                          sigma);
452        G_verbose_message(_("Normalising factor %f"),
453                          1. / gaussianFunction(sigma / 4., sigma, dimension));
454
455        /* Divide lines to segments and calculate gaussian for center of each segment */
456        Points = Vect_new_line_struct();
457        SPoints = Vect_new_line_struct();
458        SCats = Vect_new_cats_struct();
459
460        nlines = Vect_get_num_lines(&Net);
461        G_debug(3, "net nlines = %d", nlines);
462
463        for (line = 1; line <= nlines; line++) {
464            int seg, nseg, ltype;
465            double llength, length, x, y;
466
467            G_percent(line, nlines, 5);
468            ltype = Vect_read_line(&Net, Points, NULL, line);
469            if (!(ltype & GV_LINES))
470                continue;
471
472            llength = Vect_line_length(Points);
473            nseg = (int)(1 + llength / segmax);
474            length = llength / nseg;
475
476            G_debug(3, "net line = %d, nseg = %d, seg length = %f", line,
477                    nseg, length);
478
479            for (seg = 0; seg < nseg; seg++) {
480                double offset1, offset2;
481
482                offset1 = (seg + 0.5) * length;
483                Vect_point_on_line(Points, offset1, &x, &y, NULL, NULL, NULL);
484
485                G_debug(3, "  segment = %d, offset = %f, xy = %f %f", seg,
486                        offset1, x, y);
487
488                compute_net_distance(x, y, &In, &Net, netmax, sigma, term,
489                                     &gaussian, dmax, node_method);
490                gaussian *= multip;
491                if (gaussian > gausmax)
492                    gausmax = gaussian;
493
494                G_debug(3, "  gaussian = %f", gaussian);
495
496                /* Write segment */
497                if (gaussian > 0) {
498                    offset1 = seg * length;
499                    offset2 = (seg + 1) * length;
500                    if (offset2 > llength)
501                        offset2 = llength;
502                    Vect_line_segment(Points, offset1, offset2, SPoints);
503
504                    /* TODO!!! remove later
505                       if ( SPoints->n_points > 0 )
506                       Vect_append_point( SPoints, SPoints->x[SPoints->n_points-1],
507                       SPoints->y[SPoints->n_points-1], 0 );
508                     */
509                    Vect_reset_cats(SCats);
510                    Vect_cat_set(SCats, 1, (int)gaussian);
511
512                    Vect_write_line(&Out, GV_LINE, SPoints, SCats);
513
514                    total += length * gaussian;
515                }
516            }
517        }
518
519        if (flag_normalize->answer || flag_multiply->answer) {
520            double m = multip;
521
522            if (flag_normalize->answer) {
523                m /= total;
524            }
525            if (flag_multiply->answer) {
526                m *= (Vect_get_num_primitives(&In, GV_POINT) - notreachable);
527            }
528
529            Vect_build(&Out);
530
531            gausmax = 0.0;
532            nlines = Vect_get_num_lines(&Out);
533            for (line = 1; line <= nlines; line++) {
534                int cat;
535                double gaussian;
536
537                Vect_read_line(&Out, SPoints, SCats, line);
538
539                Vect_cat_get(SCats, 1, &cat);
540                gaussian = m * cat;
541                Vect_reset_cats(SCats);
542                Vect_cat_set(SCats, 1, (int)gaussian);
543                Vect_rewrite_line(&Out, line, GV_LINE, SPoints, SCats);
544                if (gaussian > gausmax)
545                    gausmax = gaussian;
546            }
547            Vect_build_partial(&Out, GV_BUILD_NONE);    /* to force rebuild */
548        }
549
550        Vect_close(&Net);
551
552        Vect_build(&Out);
553        Vect_close(&Out);
554    }
555    else {
556        /* spatial index handling, borrowed from lib/vector/Vlib/find.c */
557        struct bound_box box;
558        struct boxlist *NList = Vect_new_boxlist(1);
559
560        G_verbose_message(_("Writing output raster map using smooth parameter %f"),
561                          sigma);
562        G_verbose_message(_("Normalising factor %f"),
563                          1. / gaussianFunction(sigma / 4., sigma, dimension));
564
565        for (row = 0; row < window.rows; row++) {
566            G_percent(row, window.rows, 2);
567            if (mask)
568                Rast_get_c_row(maskfd, mask, row);
569
570            for (col = 0; col < window.cols; col++) {
571                /* don't interpolate outside of the mask */
572                if (mask && mask[col] == 0) {
573                    Rast_set_d_null_value(&output_cell[col], 1);
574                    continue;
575                }
576
577                N = Rast_row_to_northing(row + 0.5, &window);
578                E = Rast_col_to_easting(col + 0.5, &window);
579
580                if ((col & 31) == 0) {
581
582                    /* create bounding box 32x2*dmax size from the current cell center */
583                    box.N = N + dmax;
584                    box.S = N - dmax;
585                    box.E = E + dmax + 32 * window.ew_res;
586                    box.W = E - dmax;
587                    box.T = HUGE_VAL;
588                    box.B = -HUGE_VAL;
589
590                    Vect_select_lines_by_box(&In, &box, GV_POINT, NList);                   
591                }
592                box.N = N + dmax;
593                box.S = N - dmax;
594                box.E = E + dmax;
595                box.W = E - dmax;
596                box.T = HUGE_VAL;
597                box.B = -HUGE_VAL;
598                /* compute_distance(N, E, &In, sigma, term, &gaussian, dmax); */
599                compute_distance(N, E, sigma, term, &gaussian, dmax, &box, NList);
600
601                output_cell[col] = multip * gaussian;
602                if (gaussian > gausmax)
603                    gausmax = gaussian;
604            }
605            Rast_put_row(fdout, output_cell, DCELL_TYPE);
606        }
607        G_percent(1, 1, 1);
608       
609        Rast_close(fdout);
610    }
611
612    G_done_msg(_("Maximum value in output: %e."), multip * gausmax);
613
614    Vect_close(&In);
615
616    exit(EXIT_SUCCESS);
617}
618
619
620/* Read points to array return number of points */
621int read_points(struct Map_info *In, double ***coordinate, double dsize)
622{
623    int line, nlines, npoints, ltype, i = 0;
624    double **xySites;
625    static struct line_pnts *Points = NULL;
626
627    if (!Points)
628        Points = Vect_new_line_struct();
629
630    /* Allocate array of pointers */
631    npoints = Vect_get_num_primitives(In, GV_POINT);
632    xySites = (double **)G_calloc(npoints, sizeof(double *));
633
634    nlines = Vect_get_num_lines(In);
635
636    for (line = 1; line <= nlines; line++) {
637        ltype = Vect_read_line(In, Points, NULL, line);
638        if (!(ltype & GV_POINT))
639            continue;
640
641        xySites[i] = (double *)G_calloc((size_t) 2, sizeof(double));
642
643        xySites[i][0] = Points->x[0];
644        xySites[i][1] = Points->y[0];
645        i++;
646    }
647
648    *coordinate = xySites;
649
650    return (npoints);
651}
652
653
654/* Calculate distances < dmax between all sites in coordinate
655 * Return: number of distances in dists */
656double compute_all_distances(double **coordinate, double **dists, int n,
657                             double dmax)
658{
659    int ii, jj, kk;
660    size_t nn;
661
662    nn = n * (n - 1) / 2;
663    *dists = (double *)G_calloc(nn, sizeof(double));
664    kk = 0;
665
666    for (ii = 0; ii < n - 1; ii++) {
667        for (jj = ii + 1; jj < n; jj++) {
668            double dist;
669
670            dist = euclidean_distance(coordinate[ii], coordinate[jj], 2);
671            G_debug(3, "dist = %f", dist);
672
673            if (dist <= dmax) {
674                (*dists)[kk] = dist;
675                kk++;
676            }
677        }
678    }
679
680    return (kk);
681}
682
683
684/* Calculate distances < dmax between all sites in coordinate
685 * Return: number of distances in dists */
686double compute_all_net_distances(struct Map_info *In, struct Map_info *Net,
687                                 double netmax, double **dists, double dmax)
688{
689    int nn, kk, nalines, aline;
690    double dist;
691    struct line_pnts *APoints, *BPoints;
692    struct bound_box box;
693    struct boxlist *List;
694
695    APoints = Vect_new_line_struct();
696    BPoints = Vect_new_line_struct();
697    List = Vect_new_boxlist(0);
698
699    nn = Vect_get_num_primitives(In, GV_POINTS);
700    nn = nn * (nn - 1);
701    *dists = (double *)G_calloc(nn, sizeof(double));
702    kk = 0;
703
704    nalines = Vect_get_num_lines(In);
705    for (aline = 1; aline <= nalines; aline++) {
706        int i, altype;
707
708        G_debug(3, "  aline = %d", aline);
709
710        altype = Vect_read_line(In, APoints, NULL, aline);
711        if (!(altype & GV_POINTS))
712            continue;
713
714        box.E = APoints->x[0] + dmax;
715        box.W = APoints->x[0] - dmax;
716        box.N = APoints->y[0] + dmax;
717        box.S = APoints->y[0] - dmax;
718        box.T = PORT_DOUBLE_MAX;
719        box.B = -PORT_DOUBLE_MAX;
720
721        Vect_select_lines_by_box(In, &box, GV_POINT, List);
722        G_debug(3, "  %d points in box", List->n_values);
723
724        for (i = 0; i < List->n_values; i++) {
725            int bline, ret;
726
727            bline = List->id[i];
728
729            if (bline == aline)
730                continue;
731
732            G_debug(3, "    bline = %d", bline);
733            Vect_read_line(In, BPoints, NULL, bline);
734
735            ret =
736                Vect_net_shortest_path_coor(Net, APoints->x[0], APoints->y[0],
737                                            0.0, BPoints->x[0], BPoints->y[0],
738                                            0.0, netmax, netmax, &dist, NULL,
739                                            NULL, NULL, NULL, NULL, NULL, NULL);
740
741            G_debug(3, "  SP: %f %f -> %f %f", APoints->x[0], APoints->y[0],
742                    BPoints->x[0], BPoints->y[0]);
743
744            if (ret == 0) {
745                G_debug(3, "not reachable");
746                continue;       /* Not reachable */
747            }
748
749            G_debug(3, "  dist = %f", dist);
750
751            if (dist <= dmax) {
752                (*dists)[kk] = dist;
753                kk++;
754            }
755            G_debug(3, "  kk = %d", kk);
756        }
757    }
758
759    return (kk);
760}
761
762/* get number of arcs for a node */
763int count_node_arcs(struct Map_info *Map, int node)
764{
765    int i, n, line, type;
766    int count = 0;
767
768    n = Vect_get_node_n_lines(Map, node);
769    for (i = 0; i < n; i++) {
770        line = Vect_get_node_line(Map, node, i);
771        type = Vect_get_line_type(Map, abs(line));
772        if (type & GV_LINES)
773            count++;
774    }
775    return count;
776}
777
778/* Compute gausian for x, y along Net, using all points in In */
779void compute_net_distance(double x, double y, struct Map_info *In,
780                          struct Map_info *Net, double netmax, double sigma,
781                          double term, double *gaussian, double dmax, int node_method)
782{
783    int i;
784    double dist, kernel;
785    static struct line_pnts *FPoints = NULL;
786    struct bound_box box;
787    static struct boxlist *PointsList = NULL;
788    static struct ilist *NodesList = NULL;
789
790    if (!PointsList)
791        PointsList = Vect_new_boxlist(1);
792
793    if (node_method == NODE_EQUAL_SPLIT) {
794        if (!NodesList)
795            NodesList = Vect_new_list();
796
797        if (!FPoints)
798            FPoints = Vect_new_line_struct();
799    }
800
801    *gaussian = .0;
802
803    /* The network is usually much bigger than dmax and to calculate shortest path is slow
804     * -> use spatial index to select points
805     * enlarge the box by netmax (max permitted distance between a point and net) */
806    box.E = x + dmax + netmax;
807    box.W = x - dmax - netmax;
808    box.N = y + dmax + netmax;
809    box.S = y - dmax - netmax;
810    box.T = PORT_DOUBLE_MAX;
811    box.B = -PORT_DOUBLE_MAX;
812
813    Vect_select_lines_by_box(In, &box, GV_POINT, PointsList);
814    G_debug(3, "  %d points in box", PointsList->n_values);
815
816    for (i = 0; i < PointsList->n_values; i++) {
817        int line, ret;
818
819        line = PointsList->id[i];
820
821        G_debug(3, "  SP: %f %f -> %f %f", x, y, PointsList->box[i].E, PointsList->box[i].N);
822        /*ret = Vect_net_shortest_path_coor(Net, x, y, 0.0, Points->x[0], */
823        /*Points->y[0], 0.0, netmax, netmax, */
824        /*&dist, NULL, NULL, NULL, NULL, NULL, NULL, */
825        /*NULL); */
826        ret = Vect_net_shortest_path_coor(Net,
827                                           PointsList->box[i].E,
828                                           PointsList->box[i].N, 0.0,
829                                           x, y, 0.0, netmax, 1.0,
830                                           &dist, NULL,
831                                           NULL, NodesList, FPoints, NULL,
832                                           NULL, NULL);
833
834        if (ret == 0) {
835            G_debug(3, "not reachable");
836            continue;           /* Not reachable */
837        }
838
839        /* if (dist <= dmax)
840         *gaussian += gaussianKernel(dist / sigma, term); */
841
842        if (dist > dmax)
843            continue;
844
845        /* kernel = gaussianKernel(dist / sigma, term); */
846        kernel = kernelFunction(term, sigma, dist);
847
848        if (node_method == NODE_EQUAL_SPLIT) {
849            int j, node;
850            double ndiv = 1.;
851            int start = 0;
852
853            /* Count the nodes and arcs on path (n1-1)*(n2-1)* ... (ns-1) */
854
855            for (j = start; j < NodesList->n_values; j++) {
856                node = NodesList->value[j];
857
858                /* Divide into 2/n if point falls on a node */
859                if (j == 0 && FPoints->n_points < 3) {
860                    ndiv *= count_node_arcs(Net, node) / 2.;
861                }
862                else {
863                    ndiv *= count_node_arcs(Net, node) - 1;
864                }
865            }
866            kernel /= ndiv;
867        }
868        *gaussian += kernel;
869        G_debug(3, "  dist = %f gaussian = %f", dist, *gaussian);
870    }
871}
872
873void compute_distance(double N, double E, double sigma, double term,
874                      double *gaussian, double dmax, struct bound_box *box,
875                      struct boxlist *NList)
876{
877    int line, nlines;
878    double a[2], b[2];
879    double dist;
880
881    a[0] = E;
882    a[1] = N;
883
884    /* number of lines within dmax box  */
885    nlines = NList->n_values;
886
887    *gaussian = .0;
888
889    for (line = 0; line < nlines; line++) {
890
891        b[0] = NList->box[line].E;
892        b[1] = NList->box[line].N;
893
894        if (b[0] <= box->E && b[0] >= box->W &&
895            b[1] <= box->N && b[1] >= box->S)  {
896            dist = euclidean_distance(a, b, 2);
897
898            if (dist <= dmax)
899                /* *gaussian += gaussianKernel(dist / sigma, term); */
900                *gaussian += kernelFunction(term, sigma, dist);
901        }
902    }
903}
Note: See TracBrowser for help on using the repository browser.