← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to vector/v.kernel/main.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
vector/v.kernel/main.c
19
19
#include <math.h>
20
20
#include <stdio.h>
21
21
#include <stdlib.h>
 
22
#include <float.h>
22
23
#include <string.h>
23
 
#include <float.h>
24
24
#include <grass/gis.h>
 
25
#include <grass/raster.h>
25
26
#include <grass/glocale.h>
26
27
#include <grass/gmath.h>
27
 
#include <grass/Vect.h>
 
28
#include <grass/vector.h>
28
29
#include "global.h"
29
30
 
30
31
 
74
75
int main(int argc, char **argv)
75
76
{
76
77
    struct Option *in_opt, *net_opt, *out_opt;
77
 
    struct Option *stddev_opt, *dsize_opt, *segmax_opt, *netmax_opt,
 
78
    struct Option *radius_opt, *dsize_opt, *segmax_opt, *netmax_opt,
78
79
        *multip_opt, *node_opt, *kernel_opt;
79
 
    struct Flag *flag_o, *flag_q, *flag_normalize, *flag_multiply, *flag_v;
 
80
    struct Flag *flag_o, *flag_q, *flag_normalize, *flag_multiply;
 
81
    char *desc;
80
82
 
81
 
    char *mapset;
82
83
    struct Map_info In, Net, Out;
 
84
    int overwrite;
83
85
    int fdout = -1, maskfd = -1;
84
86
    int node_method, kernel_function;
85
87
    int row, col;
89
91
    CELL *mask = NULL;
90
92
    DCELL *output_cell = NULL;
91
93
    double sigma, dmax, segmax, netmax, multip;
 
94
    char *tmpstr1, *tmpstr2;
92
95
 
93
96
    double **coordinate;
94
97
    double sigmaOptimal;
95
98
    struct GModule *module;
96
99
    double dsize;
97
 
    double term;
 
100
    double term = 0;
98
101
 
99
102
    double gausmax = 0;
100
103
    int notreachable = 0;
103
106
    G_gisinit(argv[0]);
104
107
 
105
108
    module = G_define_module();
106
 
    module->keywords = _("vector, kernel density");
107
 
    module->description =
108
 
        _("Generates a raster density map from vector point data using a moving kernel or "
109
 
         "optionally generates a vector density map on a vector network.");
 
109
    G_add_keyword(_("vector"));
 
110
    G_add_keyword(_("kernel density"));
 
111
    module->label =
 
112
        _("Generates a raster density map from vector points map.");
 
113
    module->description = _("Density is computed using a moving kernel. "
 
114
                            "Optionally generates a vector density map on a vector network.");
110
115
 
111
116
    in_opt = G_define_standard_option(G_OPT_V_INPUT);
112
 
    in_opt->description = _("Input vector with training points");
113
 
 
 
117
    in_opt->label = _("Name of input vector map with training points");
 
118
    in_opt->description = NULL;
 
119
    
114
120
    net_opt = G_define_standard_option(G_OPT_V_INPUT);
115
121
    net_opt->key = "net";
116
 
    net_opt->description = _("Input network vector map");
 
122
    net_opt->label = _("Name of input network vector map");
 
123
    net_opt->description = NULL;
117
124
    net_opt->required = NO;
 
125
    net_opt->guisection = _("Network");
118
126
 
119
127
    out_opt = G_define_option();
120
128
    out_opt->key = "output";
121
129
    out_opt->type = TYPE_STRING;
122
130
    out_opt->key_desc = "name";
123
131
    out_opt->required = YES;
124
 
    out_opt->description = _("Output raster/vector map");
 
132
    out_opt->label = _("Name for output raster/vector map");
 
133
    out_opt->description = _("Outputs vector map if network map is given, otherwise raster map");
125
134
 
126
 
    stddev_opt = G_define_option();
127
 
    stddev_opt->key = "stddeviation";
128
 
    stddev_opt->type = TYPE_DOUBLE;
129
 
    stddev_opt->required = YES;
130
 
    stddev_opt->description = _("Standard deviation in map units");
 
135
    radius_opt = G_define_option();
 
136
    radius_opt->key = "radius";
 
137
    radius_opt->type = TYPE_DOUBLE;
 
138
    radius_opt->required = YES;
 
139
    radius_opt->description = _("Kernel radius in map units");
131
140
 
132
141
    dsize_opt = G_define_option();
133
142
    dsize_opt->key = "dsize";
142
151
    segmax_opt->required = NO;
143
152
    segmax_opt->description = _("Maximum length of segment on network");
144
153
    segmax_opt->answer = "100.";
 
154
    segmax_opt->guisection = _("Network");
145
155
 
146
156
    netmax_opt = G_define_option();
147
157
    netmax_opt->key = "distmax";
149
159
    netmax_opt->required = NO;
150
160
    netmax_opt->description = _("Maximum distance from point to network");
151
161
    netmax_opt->answer = "100.";
 
162
    netmax_opt->guisection = _("Network");
152
163
 
153
164
    multip_opt = G_define_option();
154
 
    multip_opt->key = "mult";
 
165
    multip_opt->key = "multiplier";
155
166
    multip_opt->type = TYPE_DOUBLE;
156
167
    multip_opt->required = NO;
157
168
    multip_opt->description = _("Multiply the density result by this number");
164
175
    node_opt->description = _("Node method");
165
176
    node_opt->options = "none,split";
166
177
    node_opt->answer = "none";
167
 
    node_opt->descriptions =
168
 
        _("none;No method applied at nodes with more than 2 arcs;"
169
 
          "split;Equal split (Okabe 2009) applied at nodes;");
 
178
    desc = NULL;
 
179
    G_asprintf(&desc,
 
180
               "none;%s;split;%s",
 
181
               _("No method applied at nodes with more than 2 arcs"),
 
182
               _("Equal split (Okabe 2009) applied at nodes"));
 
183
    node_opt->descriptions = desc;
 
184
    node_opt->guisection = _("Network");
170
185
 
171
186
    kernel_opt = G_define_option();
172
187
    kernel_opt->key = "kernel";
191
206
    flag_normalize->key = 'n';
192
207
    flag_normalize->description =
193
208
        _("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");
 
209
    flag_normalize->guisection = _("Network");
194
210
 
195
211
    flag_multiply = G_define_flag();
196
212
    flag_multiply->key = 'm';
197
213
    flag_multiply->description =
198
 
        _("In network mode, multiply the result by number of input points.");
199
 
 
200
 
    flag_v = G_define_flag();
201
 
    flag_v->key = 'v';
202
 
    flag_v->description =
203
 
        _("Verbose module output (retained for backwards compatibility)");
204
 
 
 
214
        _("In network mode, multiply the result by number of input points");
 
215
    flag_multiply->guisection = _("Network");
 
216
    
 
217
    overwrite = G_check_overwrite(argc, argv);
205
218
    if (G_parser(argc, argv))
206
219
        exit(EXIT_FAILURE);
207
220
 
208
 
 
209
 
    /* remove for grass7 */
210
 
    if (flag_v->answer)
211
 
        G_set_verbose(G_verbose_max());
 
221
    if (net_opt->answer) {
 
222
        if (G_find_vector2(out_opt->answer, G_mapset()))
 
223
            if (overwrite)
 
224
                G_warning(_("Vector map <%s> already exists and will be overwritten"),
 
225
                          out_opt->answer);
 
226
            else
 
227
                G_fatal_error(_("Vector map <%s> already exists"),
 
228
                              out_opt->answer);
 
229
    }
 
230
    else {
 
231
        if (G_find_raster(out_opt->answer, G_mapset()))
 
232
            if (overwrite)
 
233
                G_warning(_("Raster map <%s> already exists and will be overwritten"),
 
234
                          out_opt->answer);
 
235
            else
 
236
                G_fatal_error(_("Raster map <%s> already exists"),
 
237
                              out_opt->answer);
 
238
    }
212
239
 
213
240
    /*read options */
214
 
    sigma = atof(stddev_opt->answer);
 
241
    dmax = atof(radius_opt->answer);
 
242
    sigma = dmax;
215
243
    dsize = atof(dsize_opt->answer);
216
244
    segmax = atof(segmax_opt->answer);
217
245
    netmax = atof(netmax_opt->answer);
260
288
 
261
289
    if (net_opt->answer) {
262
290
        Vect_check_input_output_name(in_opt->answer, out_opt->answer,
263
 
                                     GV_FATAL_EXIT);
 
291
                                     G_FATAL_EXIT);
264
292
        Vect_check_input_output_name(net_opt->answer, out_opt->answer,
265
 
                                     GV_FATAL_EXIT);
 
293
                                     G_FATAL_EXIT);
266
294
    }
267
295
 
268
296
    G_get_window(&window);
269
297
 
270
 
    G_message("STDDEV: %f\nRES: %f\tROWS: %d\tCOLS: %d",
271
 
              sigma, window.ew_res, window.rows, window.cols);
272
 
 
 
298
    G_verbose_message(_("Standard deviation: %f"), sigma);
 
299
    G_asprintf(&tmpstr1, _n("%d row", "%d rows", window.rows), window.rows);
 
300
    G_asprintf(&tmpstr2, _n("%d column", "%d columns", window.cols), window.cols);
 
301
    /* GTC First argument is resolution, second - number of rows as a text, third - number of columns as a text. */
 
302
    G_verbose_message(_("Output raster map: resolution: %f\t%s\t%s"),
 
303
                      window.ew_res, tmpstr1, tmpstr2);
 
304
    G_free(tmpstr1);
 
305
    G_free(tmpstr2); 
 
306
    
273
307
    /* Open input vector */
274
 
    if ((mapset = G_find_vector2(in_opt->answer, "")) == NULL)
275
 
        G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
276
 
 
277
308
    Vect_set_open_level(2);
278
 
    Vect_open_old(&In, in_opt->answer, mapset);
279
 
 
 
309
    if (Vect_open_old(&In, in_opt->answer, "") < 0)
 
310
        G_fatal_error(_("Unable to open vector map <%s>"), in_opt->answer);
280
311
 
281
312
    if (net_opt->answer) {
282
313
        int nlines, line;
286
317
        net = 1;
287
318
        dimension = 1.;
288
319
        /* Open input network */
289
 
        if ((mapset = G_find_vector2(net_opt->answer, "")) == NULL)
290
 
            G_fatal_error(_("Network input map <%s> not found"),
291
 
                          net_opt->answer);
292
 
 
293
320
        Vect_set_open_level(2);
294
 
        Vect_open_old(&Net, net_opt->answer, mapset);
 
321
 
 
322
        if (Vect_open_old(&Net, net_opt->answer, "") < 0)
 
323
            G_fatal_error(_("Unable to open vector map <%s>"), net_opt->answer);
 
324
 
295
325
        Vect_net_build_graph(&Net, GV_LINES, 0, 0, NULL, NULL, NULL, 0, 0);
296
326
 
297
327
        if (!flag_q->answer) {
298
 
            Vect_open_new(&Out, out_opt->answer, 0);
 
328
            if (Vect_open_new(&Out, out_opt->answer, 0) < 0)
 
329
                G_fatal_error(_("Unable to create vector map <%s>"),
 
330
                                out_opt->answer);
299
331
            Vect_hist_command(&Out);
300
332
        }
301
333
 
314
346
        }
315
347
 
316
348
        if (notreachable > 0)
317
 
            G_warning(_("%d points outside threshold"), notreachable);
 
349
            G_warning(_n("%d point outside threshold",
 
350
                         "%d points outside threshold",
 
351
                         notreachable), notreachable);
318
352
    }
319
353
    else {
320
354
        /* check and open the name of output map */
321
355
        if (!flag_q->answer) {
322
 
            if (G_legal_filename(out_opt->answer) < 0)
323
 
                G_fatal_error(_("<%s> is an illegal file name"),
324
 
                              out_opt->answer);
325
 
 
326
 
            G_set_fp_type(DCELL_TYPE);
327
 
            if ((fdout = G_open_raster_new(out_opt->answer, DCELL_TYPE)) < 0)
328
 
                G_fatal_error(_("Unable to create raster map <%s>"),
329
 
                              out_opt->answer);
 
356
            fdout = Rast_open_new(out_opt->answer, DCELL_TYPE);
330
357
 
331
358
            /* open mask file */
332
 
            if ((maskfd = G_maskfd()) >= 0)
333
 
                mask = G_allocate_cell_buf();
 
359
            if ((maskfd = Rast_maskfd()) >= 0)
 
360
                mask = Rast_allocate_c_buf();
334
361
            else
335
362
                mask = NULL;
336
363
 
337
364
            /* allocate output raster */
338
 
            output_cell = G_allocate_raster_buf(DCELL_TYPE);
 
365
            output_cell = Rast_allocate_buf(DCELL_TYPE);
339
366
        }
340
367
    }
341
368
 
364
391
        }
365
392
 
366
393
        G_message(_("Number of input points: %d."), npoints);
367
 
        G_message(_("%d distances read from the map."), ndists);
 
394
        G_message(_n("%d distance read from the map.",
 
395
                     "%d distances read from the map.",
 
396
                     ndists), ndists);
368
397
 
369
398
        if (ndists == 0)
370
399
            G_fatal_error(_("Distances between all points are beyond %e (4 * "
393
422
        }
394
423
    }
395
424
 
396
 
    term = 1. / (pow(sigma, dimension) * pow((2. * M_PI), dimension / 2.));
397
 
 
398
 
    dmax = sigma;
399
425
    if (kernel_function == KERNEL_GAUSSIAN)
400
 
        dmax = sigma * 4.;
 
426
        sigma /= 4.;
 
427
 
 
428
    if (net_opt->answer) {
 
429
        setKernelFunction(kernel_function, 1, sigma, &term);
 
430
    }
 
431
    else {
 
432
        setKernelFunction(kernel_function, 2, sigma, &term);
 
433
    }
401
434
 
402
435
    if (net) {
403
436
        int line, nlines;
405
438
        struct line_cats *SCats;
406
439
        double total = 0.0;
407
440
 
408
 
        G_message(_("\nWriting output vector map using smooth parameter=%f."),
409
 
                  sigma);
410
 
        G_message(_("\nNormalising factor=%f."),
411
 
                  1. / gaussianFunction(sigma / 4., sigma, dimension));
 
441
        G_verbose_message(_("Writing output vector map using smooth parameter %f"),
 
442
                          sigma);
 
443
        G_verbose_message(_("Normalising factor %f"),
 
444
                          1. / gaussianFunction(sigma / 4., sigma, dimension));
412
445
 
413
446
        /* Divide lines to segments and calculate gaussian for center of each segment */
414
447
        Points = Vect_new_line_struct();
422
455
            int seg, nseg, ltype;
423
456
            double llength, length, x, y;
424
457
 
 
458
            G_percent(line, nlines, 5);
425
459
            ltype = Vect_read_line(&Net, Points, NULL, line);
426
460
            if (!(ltype & GV_LINES))
427
461
                continue;
443
477
                        offset1, x, y);
444
478
 
445
479
                compute_net_distance(x, y, &In, &Net, netmax, sigma, term,
446
 
                                     &gaussian, dmax, node_method,
447
 
                                     kernel_function);
 
480
                                     &gaussian, dmax, node_method);
448
481
                gaussian *= multip;
449
482
                if (gaussian > gausmax)
450
483
                    gausmax = gaussian;
472
505
                    total += length * gaussian;
473
506
                }
474
507
            }
475
 
            G_percent(line, nlines, 1);
476
508
        }
477
509
 
478
510
        if (flag_normalize->answer || flag_multiply->answer) {
512
544
        Vect_close(&Out);
513
545
    }
514
546
    else {
515
 
        G_message(_("\nWriting output raster map using smooth parameter=%f."),
516
 
                  sigma);
517
 
        G_message(_("\nNormalising factor=%f."),
518
 
                  1. / gaussianFunction(sigma / 4., sigma, dimension));
 
547
        /* spatial index handling, borrowed from lib/vector/Vlib/find.c */
 
548
        struct bound_box box;
 
549
        struct boxlist *NList = Vect_new_boxlist(1);
 
550
 
 
551
        G_verbose_message(_("Writing output raster map using smooth parameter %f"),
 
552
                          sigma);
 
553
        G_verbose_message(_("Normalising factor %f"),
 
554
                          1. / gaussianFunction(sigma / 4., sigma, dimension));
519
555
 
520
556
        for (row = 0; row < window.rows; row++) {
521
557
            G_percent(row, window.rows, 2);
522
 
            if (mask) {
523
 
                if (G_get_map_row(maskfd, mask, row) < 0)
524
 
                    G_fatal_error(_("Unable to read MASK"));
525
 
            }
 
558
            if (mask)
 
559
                Rast_get_c_row(maskfd, mask, row);
526
560
 
527
561
            for (col = 0; col < window.cols; col++) {
528
562
                /* don't interpolate outside of the mask */
529
563
                if (mask && mask[col] == 0) {
530
 
                    G_set_d_null_value(&output_cell[col], 1);
 
564
                    Rast_set_d_null_value(&output_cell[col], 1);
531
565
                    continue;
532
566
                }
533
567
 
534
 
                N = G_row_to_northing(row + 0.5, &window);
535
 
                E = G_col_to_easting(col + 0.5, &window);
536
 
 
 
568
                N = Rast_row_to_northing(row + 0.5, &window);
 
569
                E = Rast_col_to_easting(col + 0.5, &window);
 
570
 
 
571
                if ((col & 31) == 0) {
 
572
 
 
573
                    /* create bounding box 32x2*dmax size from the current cell center */
 
574
                    box.N = N + dmax;
 
575
                    box.S = N - dmax;
 
576
                    box.E = E + dmax + 32 * window.ew_res;
 
577
                    box.W = E - dmax;
 
578
                    box.T = HUGE_VAL;
 
579
                    box.B = -HUGE_VAL;
 
580
 
 
581
                    Vect_select_lines_by_box(&In, &box, GV_POINT, NList);                   
 
582
                }
 
583
                box.N = N + dmax;
 
584
                box.S = N - dmax;
 
585
                box.E = E + dmax;
 
586
                box.W = E - dmax;
 
587
                box.T = HUGE_VAL;
 
588
                box.B = -HUGE_VAL;
537
589
                /* compute_distance(N, E, &In, sigma, term, &gaussian, dmax); */
538
 
                compute_distance(N, E, &In, sigma, term, &gaussian, dmax,
539
 
                                 kernel_function);
 
590
                compute_distance(N, E, sigma, term, &gaussian, dmax, &box, NList);
540
591
 
541
592
                output_cell[col] = multip * gaussian;
542
593
                if (gaussian > gausmax)
543
594
                    gausmax = gaussian;
544
595
            }
545
 
            G_put_raster_row(fdout, output_cell, DCELL_TYPE);
 
596
            Rast_put_row(fdout, output_cell, DCELL_TYPE);
546
597
        }
547
 
 
548
 
        G_close_cell(fdout);
 
598
        G_percent(1, 1, 1);
 
599
        
 
600
        Rast_close(fdout);
549
601
    }
550
602
 
551
 
    G_message(_("Maximum value in output: %e."), multip * gausmax);
 
603
    G_done_msg(_("Maximum value in output: %e."), multip * gausmax);
552
604
 
553
605
    Vect_close(&In);
554
606
 
627
679
{
628
680
    int nn, kk, nalines, aline;
629
681
    double dist;
630
 
    struct line_pnts *Points;
631
 
    BOUND_BOX box;
632
 
    struct ilist *List;
 
682
    struct line_pnts *APoints, *BPoints;
 
683
    struct bound_box box;
 
684
    struct boxlist *List;
633
685
 
634
 
    Points = Vect_new_line_struct();
635
 
    List = Vect_new_list();
 
686
    APoints = Vect_new_line_struct();
 
687
    BPoints = Vect_new_line_struct();
 
688
    List = Vect_new_boxlist(0);
636
689
 
637
690
    nn = Vect_get_num_primitives(In, GV_POINTS);
638
691
    nn = nn * (nn - 1);
645
698
 
646
699
        G_debug(3, "  aline = %d", aline);
647
700
 
648
 
        altype = Vect_read_line(In, Points, NULL, aline);
 
701
        altype = Vect_read_line(In, APoints, NULL, aline);
649
702
        if (!(altype & GV_POINTS))
650
703
            continue;
651
704
 
652
 
        box.E = Points->x[0] + dmax;
653
 
        box.W = Points->x[0] - dmax;
654
 
        box.N = Points->y[0] + dmax;
655
 
        box.S = Points->y[0] - dmax;
 
705
        box.E = APoints->x[0] + dmax;
 
706
        box.W = APoints->x[0] - dmax;
 
707
        box.N = APoints->y[0] + dmax;
 
708
        box.S = APoints->y[0] - dmax;
656
709
        box.T = PORT_DOUBLE_MAX;
657
710
        box.B = -PORT_DOUBLE_MAX;
658
711
 
662
715
        for (i = 0; i < List->n_values; i++) {
663
716
            int bline, ret;
664
717
 
665
 
            bline = List->value[i];
 
718
            bline = List->id[i];
666
719
 
667
720
            if (bline == aline)
668
721
                continue;
669
722
 
670
723
            G_debug(3, "    bline = %d", bline);
671
 
 
672
 
            Vect_get_line_box(In, bline, &box);
673
 
 
674
 
 
675
 
            G_debug(3, "  SP: %f %f -> %f %f", Points->x[0], Points->y[0],
676
 
                    box.E, box.N);
 
724
            Vect_read_line(In, BPoints, NULL, bline);
677
725
 
678
726
            ret =
679
 
                Vect_net_shortest_path_coor(Net, Points->x[0], Points->y[0],
680
 
                                            0.0, box.E, box.N,
 
727
                Vect_net_shortest_path_coor(Net, APoints->x[0], APoints->y[0],
 
728
                                            0.0, BPoints->x[0], BPoints->y[0],
681
729
                                            0.0, netmax, netmax, &dist, NULL,
682
730
                                            NULL, NULL, NULL, NULL, NULL);
 
731
 
 
732
            G_debug(3, "  SP: %f %f -> %f %f", APoints->x[0], APoints->y[0],
 
733
                    BPoints->x[0], BPoints->y[0]);
 
734
 
683
735
            if (ret == 0) {
684
736
                G_debug(3, "not reachable");
685
737
                continue;       /* Not reachable */
707
759
    n = Vect_get_node_n_lines(Map, node);
708
760
    for (i = 0; i < n; i++) {
709
761
        line = Vect_get_node_line(Map, node, i);
710
 
        type = Vect_read_line(Map, NULL, NULL, abs(line));
 
762
        type = Vect_get_line_type(Map, abs(line));
711
763
        if (type & GV_LINES)
712
764
            count++;
713
765
    }
717
769
/* Compute gausian for x, y along Net, using all points in In */
718
770
void compute_net_distance(double x, double y, struct Map_info *In,
719
771
                          struct Map_info *Net, double netmax, double sigma,
720
 
                          double term, double *gaussian, double dmax,
721
 
                          int node_method, int kernel_function)
 
772
                          double term, double *gaussian, double dmax, int node_method)
722
773
{
723
774
    int i;
724
775
    double dist, kernel;
725
776
    static struct line_pnts *FPoints = NULL;
726
 
    BOUND_BOX box;
727
 
    static struct ilist *PointsList = NULL;
 
777
    struct bound_box box;
 
778
    static struct boxlist *PointsList = NULL;
728
779
    static struct ilist *NodesList = NULL;
729
780
 
730
781
    if (!PointsList)
731
 
        PointsList = Vect_new_list();
 
782
        PointsList = Vect_new_boxlist(1);
732
783
 
733
784
    if (node_method == NODE_EQUAL_SPLIT) {
734
785
        if (!NodesList)
756
807
    for (i = 0; i < PointsList->n_values; i++) {
757
808
        int line, ret;
758
809
 
759
 
        line = PointsList->value[i];
760
 
 
761
 
        Vect_get_line_box(In, line, &box);
762
 
 
763
 
        G_debug(3, "  SP: %f %f -> %f %f", x, y, box.E, box.N);
764
 
 
 
810
        line = PointsList->id[i];
 
811
 
 
812
        G_debug(3, "  SP: %f %f -> %f %f", x, y, PointsList->box[i].E, PointsList->box[i].N);
 
813
        /*ret = Vect_net_shortest_path_coor(Net, x, y, 0.0, Points->x[0], */
 
814
        /*Points->y[0], 0.0, netmax, netmax, */
 
815
        /*&dist, NULL, NULL, NULL, NULL, NULL, */
 
816
        /*NULL); */
765
817
        ret = Vect_net_shortest_path_coor2(Net,
766
 
                                           box.E, box.N, 0.0,
 
818
                                           PointsList->box[i].E,
 
819
                                           PointsList->box[i].N, 0.0,
767
820
                                           x, y, 0.0, netmax, 1.0,
768
821
                                           &dist, NULL,
769
822
                                           NULL, NodesList, FPoints, NULL,
781
834
            continue;
782
835
 
783
836
        /* kernel = gaussianKernel(dist / sigma, term); */
784
 
        kernel = kernelFunction(kernel_function, 1, sigma, dist);
 
837
        kernel = kernelFunction(term, sigma, dist);
785
838
 
786
839
        if (node_method == NODE_EQUAL_SPLIT) {
787
840
            int j, node;
808
861
    }
809
862
}
810
863
 
811
 
void compute_distance(double N, double E, struct Map_info *In,
812
 
                      double sigma, double term, double *gaussian,
813
 
                      double dmax, int kernel_function)
 
864
void compute_distance(double N, double E, double sigma, double term,
 
865
                      double *gaussian, double dmax, struct bound_box *box,
 
866
                      struct boxlist *NList)
814
867
{
815
868
    int line, nlines;
816
869
    double a[2], b[2];
817
870
    double dist;
818
871
 
819
 
    /* spatial index handling, borrowed from lib/vector/Vlib/find.c */
820
 
    BOUND_BOX box;
821
 
    static struct ilist *NList = NULL;
822
 
 
823
872
    a[0] = E;
824
873
    a[1] = N;
825
874
 
826
 
    if (!NList)
827
 
        NList = Vect_new_list();
828
 
 
829
 
    /* create bounding box 2x2*dmax size from the current cell center */
830
 
    box.N = N + dmax;
831
 
    box.S = N - dmax;
832
 
    box.E = E + dmax;
833
 
    box.W = E - dmax;
834
 
    box.T = HUGE_VAL;
835
 
    box.B = -HUGE_VAL;
836
 
 
837
875
    /* number of lines within dmax box  */
838
 
    nlines = Vect_select_lines_by_box(In, &box, GV_POINT, NList);
 
876
    nlines = NList->n_values;
839
877
 
840
878
    *gaussian = .0;
841
879
 
842
880
    for (line = 0; line < nlines; line++) {
843
881
 
844
 
        Vect_get_line_box(In, NList->value[line], &box);
845
 
        b[0] = box.E;
846
 
        b[1] = box.N;
847
 
 
848
 
        dist = euclidean_distance(a, b, 2);
849
 
 
850
 
        if (dist <= dmax)
851
 
            /* *gaussian += gaussianKernel(dist / sigma, term); */
852
 
            *gaussian += kernelFunction(kernel_function, 2, sigma, dist);
 
882
        b[0] = NList->box[line].E;
 
883
        b[1] = NList->box[line].N;
 
884
 
 
885
        if (b[0] <= box->E && b[0] >= box->W &&
 
886
            b[1] <= box->N && b[1] >= box->S)  {
 
887
            dist = euclidean_distance(a, b, 2);
 
888
 
 
889
            if (dist <= dmax)
 
890
                /* *gaussian += gaussianKernel(dist / sigma, term); */
 
891
                *gaussian += kernelFunction(term, sigma, dist);
 
892
        }
853
893
    }
854
894
}