1 /*!
2 * \file lib/vector/Vlib/net_build.c
3 *
4 * \brief Vector library - related fns for vector network building
5 *
6 * Higher level functions for reading/writing/manipulating vectors.
7 *
8 * (C) 2001-2009, 2014 by the GRASS Development Team
9 *
10 * This program is free software under the GNU General Public License
11 * (>=v2). Read the file COPYING that comes with GRASS for details.
12 *
13 * \author Radim Blazek
14 * \author Stepan Turek stepan.turek seznam.cz (turns support)
15 */
16
17 #include <grass/dbmi.h>
18 #include <grass/vector.h>
19 #include <grass/glocale.h>
20
21 /*!
22 \brief Build network graph with turntable.
23
24 Internal format for edge costs is integer, costs are multiplied
25 before conversion to int by 1000 and for lengths LL without geo flag by 1000000.
26 The same multiplication factor is used for nodes.
27 Costs in database column may be 'integer' or 'double precision' number >= 0
28 or -1 for infinity i.e. arc or node is closed and cannot be traversed
29 If record in table is not found for arcs, costs for arc are set to 0.
30 If record in table is not found for node, costs for node are set to 0.
31
32 \param Map vector map
33 \param ltype line type for arcs
34 \param afield arc costs field (if 0, use length)
35 \param nfield node costs field (if 0, do not use node costs)
36 \param tfield field where turntable is attached
37 \param tucfield field with unique categories used in the turntable
38 \param afcol column with forward costs for arc
39 \param abcol column with backward costs for arc (if NULL, back costs = forward costs),
40 \param ncol column with costs for nodes (if NULL, do not use node costs),
41 \param geo use geodesic calculation for length (LL),
42 \param algorithm not used (in future code for algorithm)
43
44 \return 0 on success, 1 on error
45 */
46
47 int
Vect_net_ttb_build_graph(struct Map_info * Map,int ltype,int afield,int nfield,int tfield,int tucfield,const char * afcol,const char * abcol,const char * ncol,int geo,int algorithm)48 Vect_net_ttb_build_graph(struct Map_info *Map,
49 int ltype,
50 int afield,
51 int nfield,
52 int tfield,
53 int tucfield,
54 const char *afcol,
55 const char *abcol,
56 const char *ncol, int geo, int algorithm)
57 {
58 /* TODO very long function, split into smaller ones */
59 int i, j, from, to, line, nlines, nnodes, ret, type, cat, skipped, cfound;
60 struct line_pnts *Points;
61 struct line_cats *Cats;
62 double dcost, bdcost, ll;
63 int cost, bcost;
64 dglGraph_s *gr;
65 dglInt32_t opaqueset[16] =
66 { 360000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
67 struct field_info *Fi;
68 dbDriver *driver = NULL;
69 dbDriver *ttbdriver = NULL;
70 dbHandle handle;
71 dbString stmt;
72 dbColumn *Column;
73 dbCatValArray fvarr, bvarr;
74 int fctype = 0, bctype = 0, nrec, nturns;
75
76 int ln_cat, nnode_lns, i_line, line_id, i_virt_edge;
77 struct line_cats *ln_Cats;
78 double x, y, z;
79 struct bound_box box;
80 struct boxlist *List;
81
82 dglInt32_t dgl_cost = cost;
83
84 /*TODO attributes of turntable should be stored in one place */
85 const char *tcols[] = { "cat", "ln_from", "ln_to", "cost", "isec", NULL
86 };
87 dbCatValArray tvarrs[5];
88 int tctype[5] = { 0 };
89 int tucfield_idx;
90
91 int t, f;
92 int node_pt_id, turn_cat, tucfound;
93 int isec;
94
95 /* TODO int costs -> double (waiting for dglib) */
96 G_debug(1, "Vect_net_ttb_build_graph(): "
97 "ltype = %d, afield = %d, nfield = %d, tfield = %d, tucfield = %d ",
98 ltype, afield, nfield, tfield, tucfield);
99 G_debug(1, " afcol = %s, abcol = %s, ncol = %s", afcol, abcol, ncol);
100
101 G_message(_("Building graph..."));
102
103 Map->dgraph.line_type = ltype;
104
105 Points = Vect_new_line_struct();
106 Cats = Vect_new_cats_struct();
107
108 ll = 0;
109 if (G_projection() == 3)
110 ll = 1; /* LL */
111
112 if (afcol == NULL && ll && !geo)
113 Map->dgraph.cost_multip = 1000000;
114 else
115 Map->dgraph.cost_multip = 1000;
116
117 nlines = Vect_get_num_lines(Map);
118 nnodes = Vect_get_num_nodes(Map);
119
120 gr = &(Map->dgraph.graph_s);
121
122 /* Allocate space for costs, later replace by functions reading costs from graph */
123 Map->dgraph.edge_fcosts =
124 (double *)G_malloc((nlines + 1) * sizeof(double));
125 Map->dgraph.edge_bcosts =
126 (double *)G_malloc((nlines + 1) * sizeof(double));
127 Map->dgraph.node_costs =
128 (double *)G_malloc((nnodes + 1) * sizeof(double));
129
130 /* Set to -1 initially */
131 for (i = 1; i <= nlines; i++) {
132 Map->dgraph.edge_fcosts[i] = -1; /* forward */
133 Map->dgraph.edge_bcosts[i] = -1; /* backward */
134 }
135 for (i = 1; i <= nnodes; i++) {
136 Map->dgraph.node_costs[i] = 0;
137 }
138
139 dglInitialize(gr, (dglByte_t) 1, sizeof(dglInt32_t), (dglInt32_t) 0,
140 opaqueset);
141
142 if (gr == NULL)
143 G_fatal_error(_("Unable to build network graph"));
144
145 db_init_handle(&handle);
146 db_init_string(&stmt);
147
148 if (abcol != NULL && afcol == NULL)
149 G_fatal_error(_("Forward costs column not specified"));
150
151 /* --- Add arcs --- */
152 /* Open db connection */
153
154 /* Get field info */
155 if (tfield < 1)
156 G_fatal_error(_("Turntable field < 1"));
157 Fi = Vect_get_field(Map, tfield);
158 if (Fi == NULL)
159 G_fatal_error(_("Database connection not defined for layer %d"),
160 tfield);
161
162 /* Open database */
163 ttbdriver = db_start_driver_open_database(Fi->driver, Fi->database);
164 if (ttbdriver == NULL)
165 G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
166 Fi->database, Fi->driver);
167
168 i = 0;
169 while (tcols[i]) {
170 /* Load costs to array */
171 if (db_get_column(ttbdriver, Fi->table, tcols[i], &Column) != DB_OK)
172 G_fatal_error(_("Turntable column <%s> not found in table <%s>"),
173 tcols[i], Fi->table);
174
175 tctype[i] = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
176
177 if ((tctype[i] == DB_C_TYPE_INT || tctype[i] == DB_C_TYPE_DOUBLE) &&
178 !strcmp(tcols[i], "cost")) ;
179 else if (tctype[i] == DB_C_TYPE_INT) ;
180 else
181 G_fatal_error(_
182 ("Data type of column <%s> not supported (must be numeric)"),
183 tcols[i]);
184
185 db_CatValArray_init(&tvarrs[i]);
186 nturns =
187 db_select_CatValArray(ttbdriver, Fi->table, Fi->key, tcols[i],
188 NULL, &tvarrs[i]);
189 ++i;
190 }
191
192 G_debug(1, "forward costs: nrec = %d", nturns);
193
194 /* Set node attributes */
195 G_message("Register nodes");
196 if (ncol != NULL) {
197
198 G_debug(2, "Set nodes' costs");
199 if (nfield < 1)
200 G_fatal_error("Node field < 1");
201
202 G_message(_("Setting node costs..."));
203
204 Fi = Vect_get_field(Map, nfield);
205 if (Fi == NULL)
206 G_fatal_error(_("Database connection not defined for layer %d"),
207 nfield);
208
209 driver = db_start_driver_open_database(Fi->driver, Fi->database);
210 if (driver == NULL)
211 G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
212 Fi->database, Fi->driver);
213
214 /* Load costs to array */
215 if (db_get_column(driver, Fi->table, ncol, &Column) != DB_OK)
216 G_fatal_error(_("Column <%s> not found in table <%s>"),
217 ncol, Fi->table);
218
219 fctype = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
220
221 if (fctype != DB_C_TYPE_INT && fctype != DB_C_TYPE_DOUBLE)
222 G_fatal_error(_
223 ("Data type of column <%s> not supported (must be numeric)"),
224 ncol);
225
226 db_CatValArray_init(&fvarr);
227
228 nrec =
229 db_select_CatValArray(driver, Fi->table, Fi->key, ncol, NULL,
230 &fvarr);
231 G_debug(1, "node costs: nrec = %d", nrec);
232
233 tucfield_idx = Vect_cidx_get_field_index(Map, tucfield);
234 }
235
236 List = Vect_new_boxlist(0);
237 ln_Cats = Vect_new_cats_struct();
238
239 G_message("Building turns graph...");
240
241 i_virt_edge = -1;
242 for (i = 1; i <= nnodes; i++) {
243 /* TODO: what happens if we set attributes of non existing node (skipped lines,
244 * nodes without lines) */
245
246 /* select points at node */
247 Vect_get_node_coor(Map, i, &x, &y, &z);
248 box.E = box.W = x;
249 box.N = box.S = y;
250 box.T = box.B = z;
251 Vect_select_lines_by_box(Map, &box, GV_POINT, List);
252
253 G_debug(2, " node = %d nlines = %d", i, List->n_values);
254 cfound = 0;
255 dcost = 0;
256 tucfound = 0;
257
258 for (j = 0; j < List->n_values; j++) {
259 line = List->id[j];
260 G_debug(2, " line (%d) = %d", j, line);
261 type = Vect_read_line(Map, NULL, Cats, line);
262 if (!(type & GV_POINT))
263 continue;
264 /* get node column costs */
265 if (ncol != NULL && !cfound && Vect_cat_get(Cats, nfield, &cat)) { /* point with category of field found */
266 /* Set costs */
267 if (fctype == DB_C_TYPE_INT) {
268 ret = db_CatValArray_get_value_int(&fvarr, cat, &cost);
269 dcost = cost;
270 }
271 else { /* DB_C_TYPE_DOUBLE */
272 ret =
273 db_CatValArray_get_value_double(&fvarr, cat, &dcost);
274 }
275 if (ret != DB_OK) {
276 G_warning(_
277 ("Database record for node %d (cat = %d) not found "
278 "(cost set to 0)"), i, cat);
279 }
280 cfound = 1;
281 Map->dgraph.node_costs[i] = dcost;
282 }
283
284 /* add virtual nodes and lines, which represents the intersections
285 there are added two nodes for every intersection, which are linked
286 with the nodes (edges in primal graph).
287 the positive node - when we are going from the intersection
288 the negative node - when we are going to the intersection
289
290 TODO There are more possible approaches in virtual nodes management.
291 We can also add and remove them dynamically as they are needed
292 for analysis when Vect_net_ttb_shortest_path is called
293 (problem of flattening graph).
294 Currently this static solution was chosen, because it cost
295 time only when graph is build. However it costs more memory space.
296 For Dijkstra algorithm this expansion should not be serious problem
297 because we can only get into positive node or go from the negative
298 node.
299
300 */
301
302 ret = Vect_cat_get(Cats, tucfield, &cat);
303 if (!tucfound && ret) { /* point with category of field found */
304 /* find lines which belongs to the intersection */
305 nnode_lns = Vect_get_node_n_lines(Map, i);
306
307 for (i_line = 0; i_line < nnode_lns; i_line++) {
308
309 line_id = Vect_get_node_line(Map, i, i_line);
310 Vect_read_line(Map, NULL, ln_Cats, abs(line_id));
311 Vect_cat_get(ln_Cats, tucfield, &ln_cat);
312
313 if (line_id < 0)
314 ln_cat *= -1;
315 f = cat * 2;
316
317 if (ln_cat < 0)
318 t = ln_cat * -2 + 1;
319 else
320 t = ln_cat * 2;
321
322 G_debug(5,
323 "Add arc %d for virtual node from %d to %d cost = %d",
324 i_virt_edge, f, t, 0);
325
326 /* positive, start virtual node */
327 ret = dglAddEdge(gr, (dglInt32_t) f, (dglInt32_t) t,
328 (dglInt32_t) 0,
329 (dglInt32_t) (i_virt_edge));
330 if (ret < 0)
331 G_fatal_error(_
332 ("Cannot add network arc for virtual node connection."));
333
334 t = cat * 2 + 1;
335 i_virt_edge--;
336
337 if (-ln_cat < 0)
338 f = ln_cat * 2 + 1;
339 else
340 f = ln_cat * -2;
341
342 G_debug(5,
343 "Add arc %d for virtual node from %d to %d cost = %d",
344 i_virt_edge, f, t, 0);
345
346 /* negative, destination virtual node */
347 ret = dglAddEdge(gr, (dglInt32_t) f, (dglInt32_t) t,
348 (dglInt32_t) 0,
349 (dglInt32_t) (i_virt_edge));
350 if (ret < 0)
351 G_fatal_error(_
352 ("Cannot add network arc for virtual node connection."));
353
354 i_virt_edge--;
355 }
356 tucfound++;
357 }
358 else if (ret)
359 tucfound++;
360 }
361
362 if (tucfound > 1)
363 G_warning(_
364 ("There exists more than one point of node <%d> with unique category field <%d>.\n"
365 "The unique categories layer is not valid therefore you will probably get incorrect results."),
366 tucfield, i);
367
368 if (ncol != NULL && !cfound)
369 G_debug(2,
370 "Category of field %d is not attached to any points in node %d"
371 "(costs set to 0)", nfield, i);
372 }
373
374 Vect_destroy_cats_struct(ln_Cats);
375
376 for (i = 1; i <= nturns; i++) {
377 /* select points at node */
378
379 /* TODO use cursors */
380 db_CatValArray_get_value_int(&tvarrs[0], i, &turn_cat);
381
382 db_CatValArray_get_value_int(&tvarrs[1], i, &from);
383 db_CatValArray_get_value_int(&tvarrs[2], i, &to);
384
385 db_CatValArray_get_value_int(&tvarrs[4], i, &isec);
386 dcost = 0.0;
387 if (ncol != NULL) {
388 /* TODO optimization do not do it for every turn in intersection again */
389 if (Vect_cidx_find_next
390 (Map, tucfield_idx, isec, GV_POINT, 0, &type,
391 &node_pt_id) == -1) {
392 G_warning(_
393 ("Unable to find point representing intersection <%d> in unique categories field <%d>.\n"
394 "Cost for the intersection was set to 0.\n"
395 "The unique categories layer is not valid therefore you will probably get incorrect results."),
396 isec, tucfield);
397 }
398 else {
399 Vect_read_line(Map, Points, Cats, node_pt_id);
400
401 node_pt_id =
402 Vect_find_node(Map, *Points->x, *Points->y, *Points->z,
403 0.0, WITHOUT_Z);
404
405 if (node_pt_id == 0) {
406 G_warning(_
407 ("Unable to find node for point representing intersection <%d> in unique categories field <%d>.\n"
408 "Cost for the intersection was set to 0.\n"
409 "The unique categories layer is not valid therefore you will probably get incorrect results."),
410 isec, tucfield);
411 }
412 else {
413 G_debug(2, " node = %d", node_pt_id);
414 dcost = Map->dgraph.node_costs[node_pt_id];
415 }
416 }
417 }
418
419 G_debug(2, "Set node's cost to %f", dcost);
420
421 if (dcost >= 0) {
422 /* Set costs from turntable */
423 if (tctype[3] == DB_C_TYPE_INT) {
424 ret = db_CatValArray_get_value_int(&tvarrs[3], i, &cost);
425 dcost = cost;
426 }
427 else /* DB_C_TYPE_DOUBLE */
428 ret = db_CatValArray_get_value_double(&tvarrs[3], i, &dcost);
429
430 if (ret != DB_OK) {
431 G_warning(_
432 ("Database record for turn with cat = %d in not found. "
433 "(The turn was skipped."), i);
434 continue;
435 }
436
437 if (dcost >= 0) {
438
439 if (ncol != NULL)
440 cost =
441 (Map->dgraph.node_costs[node_pt_id] +
442 dcost) * (dglInt32_t) Map->dgraph.cost_multip;
443 else
444 cost = dcost * (dglInt32_t) Map->dgraph.cost_multip;
445
446 /* dglib does not like negative id's of nodes */
447 if (from < 0)
448 f = from * -2 + 1;
449 else
450 f = from * 2;
451
452 if (to < 0)
453 t = to * -2 + 1;
454 else
455 t = to * 2;
456
457 G_debug(5, "Add arc/turn %d for turn from %d to %d cost = %d",
458 turn_cat, f, t, cost);
459
460 ret = dglAddEdge(gr, (dglInt32_t) f, (dglInt32_t) t,
461 (dglInt32_t) cost, (dglInt32_t) (turn_cat));
462
463 if (ret < 0)
464 G_fatal_error(_
465 ("Cannot add network arc representing turn."));
466 }
467 }
468 }
469
470 Vect_destroy_boxlist(List);
471
472 i = 0;
473 while (tcols[i]) {
474 db_CatValArray_free(&tvarrs[i]);
475 ++i;
476 }
477
478 if (ncol != NULL) {
479 db_close_database_shutdown_driver(driver);
480 db_CatValArray_free(&fvarr);
481 }
482
483 /* Open db connection */
484 if (afcol != NULL) {
485 /* Get field info */
486 if (afield < 1)
487 G_fatal_error(_("Arc field < 1"));
488 Fi = Vect_get_field(Map, afield);
489 if (Fi == NULL)
490 G_fatal_error(_("Database connection not defined for layer %d"),
491 afield);
492
493 /* Open database */
494 driver = db_start_driver_open_database(Fi->driver, Fi->database);
495 if (driver == NULL)
496 G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
497 Fi->database, Fi->driver);
498
499 /* Load costs to array */
500 if (db_get_column(driver, Fi->table, afcol, &Column) != DB_OK)
501 G_fatal_error(_("Column <%s> not found in table <%s>"),
502 afcol, Fi->table);
503
504 fctype = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
505
506 if (fctype != DB_C_TYPE_INT && fctype != DB_C_TYPE_DOUBLE)
507 G_fatal_error(_
508 ("Data type of column <%s> not supported (must be numeric)"),
509 afcol);
510
511 db_CatValArray_init(&fvarr);
512 nrec =
513 db_select_CatValArray(driver, Fi->table, Fi->key, afcol, NULL,
514 &fvarr);
515 G_debug(1, "forward costs: nrec = %d", nrec);
516
517 if (abcol != NULL) {
518 if (db_get_column(driver, Fi->table, abcol, &Column) != DB_OK)
519 G_fatal_error(_("Column <%s> not found in table <%s>"),
520 abcol, Fi->table);
521
522 bctype = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
523
524 if (bctype != DB_C_TYPE_INT && bctype != DB_C_TYPE_DOUBLE)
525 G_fatal_error(_
526 ("Data type of column <%s> not supported (must be numeric)"),
527 abcol);
528
529 db_CatValArray_init(&bvarr);
530 nrec =
531 db_select_CatValArray(driver, Fi->table, Fi->key, abcol, NULL,
532 &bvarr);
533 G_debug(1, "backward costs: nrec = %d", nrec);
534 }
535 }
536
537 skipped = 0;
538
539 G_message(_("Registering arcs..."));
540
541 for (i = 1; i <= nlines; i++) {
542 G_percent(i, nlines, 1); /* must be before any continue */
543
544 type = Vect_read_line(Map, Points, Cats, i);
545 if (!(type & ltype & (GV_LINE | GV_BOUNDARY)))
546 continue;
547
548 Vect_get_line_nodes(Map, i, &from, &to);
549
550 dcost = bdcost = 0;
551
552 cfound = Vect_cat_get(Cats, tucfield, &cat);
553 if (!cfound)
554 continue;
555
556 if (cfound > 1)
557 G_warning(_
558 ("Line with id <%d> has more unique categories defined in field <%d>.\n"
559 "The unique categories layer is not valid therefore you will probably get incorrect results."),
560 i, tucfield);
561
562 if (afcol != NULL) {
563 if (!(Vect_cat_get(Cats, afield, &cat))) {
564 G_debug(2,
565 "Category of field %d not attached to the line %d -> cost was set to 0",
566 afield, i);
567 skipped += 2; /* Both directions */
568 }
569 else {
570 if (fctype == DB_C_TYPE_INT) {
571 ret = db_CatValArray_get_value_int(&fvarr, cat, &cost);
572 dcost = cost;
573 }
574 else { /* DB_C_TYPE_DOUBLE */
575 ret =
576 db_CatValArray_get_value_double(&fvarr, cat, &dcost);
577 }
578 if (ret != DB_OK) {
579 G_warning(_("Database record for line %d (cat = %d, "
580 "forward/both direction(s)) not found "
581 "(cost was set to 0)"), i, cat);
582 }
583
584 if (abcol != NULL) {
585 if (bctype == DB_C_TYPE_INT) {
586 ret =
587 db_CatValArray_get_value_int(&bvarr, cat, &bcost);
588 bdcost = bcost;
589 }
590 else { /* DB_C_TYPE_DOUBLE */
591 ret =
592 db_CatValArray_get_value_double(&bvarr, cat,
593 &bdcost);
594 }
595 if (ret != DB_OK) {
596 G_warning(_("Database record for line %d (cat = %d, "
597 "backword direction) not found"
598 "(cost was set to 0)"), i, cat);
599 }
600 }
601 else
602 bdcost = dcost;
603
604 Vect_cat_get(Cats, tucfield, &cat);
605 }
606 }
607 else {
608 if (ll) {
609 if (geo)
610 dcost = Vect_line_geodesic_length(Points);
611 else
612 dcost = Vect_line_length(Points);
613 }
614 else
615 dcost = Vect_line_length(Points);
616
617 bdcost = dcost;
618 }
619
620 cost = (dglInt32_t) Map->dgraph.cost_multip * dcost;
621 dgl_cost = cost;
622
623 cat = cat * 2;
624
625 G_debug(5, "Setinng node %d cost: %d", cat, cost);
626 dglNodeSet_Attr(gr, dglGetNode(gr, (dglInt32_t) cat), &dgl_cost);
627
628 Map->dgraph.edge_fcosts[i] = dcost;
629
630 cost = (dglInt32_t) Map->dgraph.cost_multip * bdcost;
631 dgl_cost = cost;
632
633 ++cat;
634
635 G_debug(5, "Setinng node %d cost: %d", cat, cost);
636 dglNodeSet_Attr(gr, dglGetNode(gr, (dglInt32_t) cat), &dgl_cost);
637
638 Map->dgraph.edge_bcosts[i] = bdcost;
639 }
640
641 if (afcol != NULL && skipped > 0)
642 G_debug(2, "%d lines missing category of field %d skipped", skipped,
643 afield);
644
645 if (afcol != NULL) {
646 db_close_database_shutdown_driver(driver);
647 db_CatValArray_free(&fvarr);
648
649 if (abcol != NULL) {
650 db_CatValArray_free(&bvarr);
651 }
652 }
653 db_close_database_shutdown_driver(ttbdriver);
654
655 Vect_destroy_line_struct(Points);
656 Vect_destroy_cats_struct(Cats);
657
658 G_message(_("Flattening the graph..."));
659 ret = dglFlatten(gr);
660 if (ret < 0)
661 G_fatal_error(_("GngFlatten error"));
662
663 /* init SP cache */
664 /* disable to debug dglib cache */
665 dglInitializeSPCache(gr, &(Map->dgraph.spCache));
666
667 G_message(_("Graph was built"));
668 return 0;
669 }
670
671 /*!
672 \brief Build network graph.
673
674 Internal format for edge costs is integer, costs are multiplied
675 before conversion to int by 1000 and for lengths LL without geo flag by 1000000.
676 The same multiplication factor is used for nodes.
677 Costs in database column may be 'integer' or 'double precision' number >= 0
678 or -1 for infinity i.e. arc or node is closed and cannot be traversed
679 If record in table is not found for arcs, arc is skip.
680 If record in table is not found for node, costs for node are set to 0.
681
682 \param Map vector map
683 \param ltype line type for arcs
684 \param afield arc costs field (if 0, use length)
685 \param nfield node costs field (if 0, do not use node costs)
686 \param afcol column with forward costs for arc
687 \param abcol column with backward costs for arc (if NULL, back costs = forward costs),
688 \param ncol column with costs for nodes (if NULL, do not use node costs),
689 \param geo use geodesic calculation for length (LL),
690 \param version graph version to create (1, 2, 3)
691
692 \return 0 on success, 1 on error
693 */
694 int
Vect_net_build_graph(struct Map_info * Map,int ltype,int afield,int nfield,const char * afcol,const char * abcol,const char * ncol,int geo,int version)695 Vect_net_build_graph(struct Map_info *Map,
696 int ltype,
697 int afield,
698 int nfield,
699 const char *afcol,
700 const char *abcol,
701 const char *ncol, int geo, int version)
702 {
703 /* TODO long function, split into smaller ones */
704 int i, j, from, to, line, nlines, nnodes, ret, type, cat, skipped, cfound;
705 int dofw, dobw;
706 struct line_pnts *Points;
707 struct line_cats *Cats;
708 double dcost, bdcost, ll;
709 int cost, bcost;
710 dglGraph_s *gr;
711 dglInt32_t dgl_cost;
712 dglInt32_t opaqueset[16] =
713 { 360000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
714 struct field_info *Fi;
715 dbDriver *driver = NULL;
716 dbHandle handle;
717 dbString stmt;
718 dbColumn *Column;
719 dbCatValArray fvarr, bvarr;
720 int fctype = 0, bctype = 0, nrec;
721
722 /* TODO int costs -> double (waiting for dglib) */
723 G_debug(1, "Vect_net_build_graph(): ltype = %d, afield = %d, nfield = %d",
724 ltype, afield, nfield);
725 G_debug(1, " afcol = %s, abcol = %s, ncol = %s", afcol, abcol, ncol);
726
727 G_message(_("Building graph..."));
728
729 Map->dgraph.line_type = ltype;
730
731 Points = Vect_new_line_struct();
732 Cats = Vect_new_cats_struct();
733
734 ll = 0;
735 if (G_projection() == 3)
736 ll = 1; /* LL */
737
738 if (afcol == NULL && ll && !geo)
739 Map->dgraph.cost_multip = 1000000;
740 else
741 Map->dgraph.cost_multip = 1000;
742
743 nlines = Vect_get_num_lines(Map);
744 nnodes = Vect_get_num_nodes(Map);
745
746 gr = &(Map->dgraph.graph_s);
747
748 /* Allocate space for costs, later replace by functions reading costs from graph */
749 Map->dgraph.edge_fcosts =
750 (double *)G_malloc((nlines + 1) * sizeof(double));
751 Map->dgraph.edge_bcosts =
752 (double *)G_malloc((nlines + 1) * sizeof(double));
753 Map->dgraph.node_costs =
754 (double *)G_malloc((nnodes + 1) * sizeof(double));
755 /* Set to -1 initially */
756 for (i = 1; i <= nlines; i++) {
757 Map->dgraph.edge_fcosts[i] = -1; /* forward */
758 Map->dgraph.edge_bcosts[i] = -1; /* backward */
759 }
760 for (i = 1; i <= nnodes; i++) {
761 Map->dgraph.node_costs[i] = 0;
762 }
763
764 if (version < 1 || version > 3)
765 version = 1;
766
767 if (ncol != NULL)
768 dglInitialize(gr, (dglByte_t) version, sizeof(dglInt32_t), (dglInt32_t) 0,
769 opaqueset);
770 else
771 dglInitialize(gr, (dglByte_t) version, (dglInt32_t) 0, (dglInt32_t) 0,
772 opaqueset);
773
774 if (gr == NULL)
775 G_fatal_error(_("Unable to build network graph"));
776
777 db_init_handle(&handle);
778 db_init_string(&stmt);
779
780 if (abcol != NULL && afcol == NULL)
781 G_fatal_error(_("Forward costs column not specified"));
782
783 /* --- Add arcs --- */
784 /* Open db connection */
785 if (afcol != NULL) {
786 /* Get field info */
787 if (afield < 1)
788 G_fatal_error(_("Arc field < 1"));
789 Fi = Vect_get_field(Map, afield);
790 if (Fi == NULL)
791 G_fatal_error(_("Database connection not defined for layer %d"),
792 afield);
793
794 /* Open database */
795 driver = db_start_driver_open_database(Fi->driver, Fi->database);
796 if (driver == NULL)
797 G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
798 Fi->database, Fi->driver);
799
800 /* Load costs to array */
801 if (db_get_column(driver, Fi->table, afcol, &Column) != DB_OK)
802 G_fatal_error(_("Column <%s> not found in table <%s>"),
803 afcol, Fi->table);
804
805 fctype = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
806
807 if (fctype != DB_C_TYPE_INT && fctype != DB_C_TYPE_DOUBLE)
808 G_fatal_error(_
809 ("Data type of column <%s> not supported (must be numeric)"),
810 afcol);
811
812 db_CatValArray_init(&fvarr);
813 nrec =
814 db_select_CatValArray(driver, Fi->table, Fi->key, afcol, NULL,
815 &fvarr);
816 G_debug(1, "forward costs: nrec = %d", nrec);
817
818 if (abcol != NULL) {
819 if (db_get_column(driver, Fi->table, abcol, &Column) != DB_OK)
820 G_fatal_error(_("Column <%s> not found in table <%s>"),
821 abcol, Fi->table);
822
823 bctype = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
824
825 if (bctype != DB_C_TYPE_INT && bctype != DB_C_TYPE_DOUBLE)
826 G_fatal_error(_
827 ("Data type of column <%s> not supported (must be numeric)"),
828 abcol);
829
830 db_CatValArray_init(&bvarr);
831 nrec =
832 db_select_CatValArray(driver, Fi->table, Fi->key, abcol, NULL,
833 &bvarr);
834 G_debug(1, "backward costs: nrec = %d", nrec);
835 }
836 }
837
838 skipped = 0;
839
840 G_message(_("Registering arcs..."));
841
842 for (i = 1; i <= nlines; i++) {
843 G_percent(i, nlines, 1); /* must be before any continue */
844 dofw = dobw = 1;
845 type = Vect_read_line(Map, Points, Cats, i);
846 if (!(type & ltype & (GV_LINE | GV_BOUNDARY)))
847 continue;
848
849 Vect_get_line_nodes(Map, i, &from, &to);
850
851 if (afcol != NULL) {
852 if (!(Vect_cat_get(Cats, afield, &cat))) {
853 G_debug(2,
854 "Category of field %d not attached to the line %d -> line skipped",
855 afield, i);
856 skipped += 2; /* Both directions */
857 continue;
858 }
859 else {
860 if (fctype == DB_C_TYPE_INT) {
861 ret = db_CatValArray_get_value_int(&fvarr, cat, &cost);
862 dcost = cost;
863 }
864 else { /* DB_C_TYPE_DOUBLE */
865 ret =
866 db_CatValArray_get_value_double(&fvarr, cat, &dcost);
867 }
868 if (ret != DB_OK) {
869 G_warning(_("Database record for line %d (cat = %d, "
870 "forward/both direction(s)) not found "
871 "(forward/both direction(s) of line skipped)"),
872 i, cat);
873 dofw = 0;
874 }
875
876 if (abcol != NULL) {
877 if (bctype == DB_C_TYPE_INT) {
878 ret =
879 db_CatValArray_get_value_int(&bvarr, cat, &bcost);
880 bdcost = bcost;
881 }
882 else { /* DB_C_TYPE_DOUBLE */
883 ret =
884 db_CatValArray_get_value_double(&bvarr, cat,
885 &bdcost);
886 }
887 if (ret != DB_OK) {
888 G_warning(_("Database record for line %d (cat = %d, "
889 "backword direction) not found"
890 "(direction of line skipped)"), i, cat);
891 dobw = 0;
892 }
893 }
894 else {
895 if (dofw)
896 bdcost = dcost;
897 else
898 dobw = 0;
899 }
900 }
901 }
902 else {
903 if (ll) {
904 if (geo)
905 dcost = Vect_line_geodesic_length(Points);
906 else
907 dcost = Vect_line_length(Points);
908 }
909 else
910 dcost = Vect_line_length(Points);
911
912 bdcost = dcost;
913 }
914 if (dofw && dcost != -1) {
915 cost = (dglInt32_t) Map->dgraph.cost_multip * dcost;
916 G_debug(5, "Add arc %d from %d to %d cost = %d", i, from, to,
917 cost);
918 ret =
919 dglAddEdge(gr, (dglInt32_t) from, (dglInt32_t) to,
920 (dglInt32_t) cost, (dglInt32_t) i);
921 Map->dgraph.edge_fcosts[i] = dcost;
922 if (ret < 0)
923 G_fatal_error("Cannot add network arc");
924 }
925
926 G_debug(5, "bdcost = %f edge_bcosts = %f", bdcost,
927 Map->dgraph.edge_bcosts[i]);
928 if (dobw && bdcost != -1) {
929 bcost = (dglInt32_t) Map->dgraph.cost_multip * bdcost;
930 G_debug(5, "Add arc %d from %d to %d bcost = %d", -i, to, from,
931 bcost);
932 ret =
933 dglAddEdge(gr, (dglInt32_t) to, (dglInt32_t) from,
934 (dglInt32_t) bcost, (dglInt32_t) - i);
935 Map->dgraph.edge_bcosts[i] = bdcost;
936 if (ret < 0)
937 G_fatal_error(_("Cannot add network arc"));
938 }
939 }
940
941 if (afcol != NULL && skipped > 0)
942 G_debug(2, "%d lines missing category of field %d skipped", skipped,
943 afield);
944
945 if (afcol != NULL) {
946 db_close_database_shutdown_driver(driver);
947 db_CatValArray_free(&fvarr);
948
949 if (abcol != NULL) {
950 db_CatValArray_free(&bvarr);
951 }
952 }
953
954 /* Set node attributes */
955 G_debug(2, "Register nodes");
956 if (ncol != NULL) {
957 double x, y, z;
958 struct bound_box box;
959 struct boxlist *List;
960
961 List = Vect_new_boxlist(0);
962
963 G_debug(2, "Set nodes' costs");
964 if (nfield < 1)
965 G_fatal_error("Node field < 1");
966
967 G_message(_("Setting node costs..."));
968
969 Fi = Vect_get_field(Map, nfield);
970 if (Fi == NULL)
971 G_fatal_error(_("Database connection not defined for layer %d"),
972 nfield);
973
974 driver = db_start_driver_open_database(Fi->driver, Fi->database);
975 if (driver == NULL)
976 G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
977 Fi->database, Fi->driver);
978
979 /* Load costs to array */
980 if (db_get_column(driver, Fi->table, ncol, &Column) != DB_OK)
981 G_fatal_error(_("Column <%s> not found in table <%s>"),
982 ncol, Fi->table);
983
984 fctype = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
985
986 if (fctype != DB_C_TYPE_INT && fctype != DB_C_TYPE_DOUBLE)
987 G_fatal_error(_
988 ("Data type of column <%s> not supported (must be numeric)"),
989 ncol);
990
991 db_CatValArray_init(&fvarr);
992 nrec =
993 db_select_CatValArray(driver, Fi->table, Fi->key, ncol, NULL,
994 &fvarr);
995 G_debug(1, "node costs: nrec = %d", nrec);
996
997 for (i = 1; i <= nnodes; i++) {
998 /* TODO: what happens if we set attributes of not existing node (skipped lines,
999 * nodes without lines) */
1000
1001 /* select points at node */
1002 Vect_get_node_coor(Map, i, &x, &y, &z);
1003 box.E = box.W = x;
1004 box.N = box.S = y;
1005 box.T = box.B = z;
1006 Vect_select_lines_by_box(Map, &box, GV_POINT, List);
1007
1008 G_debug(2, " node = %d nlines = %d", i, List->n_values);
1009 cfound = 0;
1010 dcost = 0;
1011
1012 for (j = 0; j < List->n_values; j++) {
1013 line = List->id[j];
1014 G_debug(2, " line (%d) = %d", j, line);
1015 type = Vect_read_line(Map, NULL, Cats, line);
1016 if (!(type & GV_POINT))
1017 continue;
1018 if (Vect_cat_get(Cats, nfield, &cat)) { /* point with category of field found */
1019 /* Set costs */
1020 if (fctype == DB_C_TYPE_INT) {
1021 ret =
1022 db_CatValArray_get_value_int(&fvarr, cat, &cost);
1023 dcost = cost;
1024 }
1025 else { /* DB_C_TYPE_DOUBLE */
1026 ret =
1027 db_CatValArray_get_value_double(&fvarr, cat,
1028 &dcost);
1029 }
1030 if (ret != DB_OK) {
1031 G_warning(_
1032 ("Database record for node %d (cat = %d) not found "
1033 "(cost set to 0)"), i, cat);
1034 }
1035 cfound = 1;
1036 break;
1037 }
1038 }
1039 if (!cfound) {
1040 G_debug(2,
1041 "Category of field %d not attached to any points in node %d"
1042 "(costs set to 0)", nfield, i);
1043 }
1044 if (dcost == -1) { /* closed */
1045 cost = -1;
1046 }
1047 else {
1048 cost = (dglInt32_t) Map->dgraph.cost_multip * dcost;
1049 }
1050
1051 dgl_cost = cost;
1052 G_debug(3, "Set node's cost to %d", cost);
1053
1054 dglNodeSet_Attr(gr, dglGetNode(gr, (dglInt32_t) i), &dgl_cost);
1055
1056 Map->dgraph.node_costs[i] = dcost;
1057 }
1058 db_close_database_shutdown_driver(driver);
1059 db_CatValArray_free(&fvarr);
1060
1061 Vect_destroy_boxlist(List);
1062 }
1063
1064 G_message(_("Flattening the graph..."));
1065 ret = dglFlatten(gr);
1066 if (ret < 0)
1067 G_fatal_error(_("GngFlatten error"));
1068
1069 /* init SP cache */
1070 /* disable to debug dglib cache */
1071 dglInitializeSPCache(gr, &(Map->dgraph.spCache));
1072
1073 G_message(_("Graph was built"));
1074
1075 return 0;
1076 }
1077