1 /*!
2 * \file lib/vector/Vlib/net_analyze.c
3 *
4 * \brief Vector library - related fns for vector network analyses
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 static int From_node; /* from node set in SP and used by clipper for first arc */
22
clipper(dglGraph_s * pgraph,dglSPClipInput_s * pargIn,dglSPClipOutput_s * pargOut,void * pvarg)23 static int clipper(dglGraph_s * pgraph,
24 dglSPClipInput_s * pargIn,
25 dglSPClipOutput_s * pargOut, void *pvarg)
26 { /* caller's pointer */
27 dglInt32_t cost;
28 dglInt32_t from;
29
30 G_debug(3, "Net: clipper()");
31
32 from = dglNodeGet_Id(pgraph, pargIn->pnNodeFrom);
33
34 G_debug(3, " Edge = %d NodeFrom = %d NodeTo = %d edge cost = %d",
35 (int)dglEdgeGet_Id(pgraph, pargIn->pnEdge),
36 (int)from, (int)dglNodeGet_Id(pgraph, pargIn->pnNodeTo),
37 (int)pargOut->nEdgeCost);
38
39 if (from != From_node) { /* do not clip first */
40 if (dglGet_NodeAttrSize(pgraph) > 0) {
41 memcpy(&cost, dglNodeGet_Attr(pgraph, pargIn->pnNodeFrom),
42 sizeof(cost));
43 if (cost == -1) { /* closed, cannot go from this node except it is 'from' node */
44 G_debug(3, " closed node");
45 return 1;
46 }
47 else {
48 G_debug(3, " EdgeCost += %d (node)", (int)cost);
49 pargOut->nEdgeCost += cost;
50 }
51 }
52 }
53 else {
54 G_debug(3, " don't clip first node");
55 }
56
57 return 0;
58 }
59
60 /*!
61 \brief Converts shortest path result, which is calculated by DGLib on newtwork without turntable, into output format.
62 */
convert_dgl_shortest_path_result(struct Map_info * Map,dglSPReport_s * pSPReport,struct ilist * List)63 static int convert_dgl_shortest_path_result(struct Map_info *Map,
64 dglSPReport_s * pSPReport,
65 struct ilist *List)
66 {
67 int i, line;
68
69 Vect_reset_list(List);
70
71 for (i = 0; i < pSPReport->cArc; i++) {
72 line =
73 dglEdgeGet_Id(&(Map->dgraph.graph_s), pSPReport->pArc[i].pnEdge);
74 G_debug(2, "From %ld to %ld - cost %ld user %d distance %ld", pSPReport->pArc[i].nFrom, pSPReport->pArc[i].nTo, dglEdgeGet_Cost(&(Map->dgraph.graph_s), pSPReport->pArc[i].pnEdge) / Map->dgraph.cost_multip, /* this is the cost from clip() */
75 line, pSPReport->pArc[i].nDistance);
76 Vect_list_append(List, line);
77 }
78
79 return 0;
80 }
81
82 /*!
83 \brief Converts shortest path result, which is calculated by DGLib on newtwork with turntable, into output format.
84 */
ttb_convert_dgl_shortest_path_result(struct Map_info * Map,dglSPReport_s * pSPReport,int tucfield,struct ilist * List)85 static int ttb_convert_dgl_shortest_path_result(struct Map_info *Map,
86 dglSPReport_s * pSPReport,
87 int tucfield,
88 struct ilist *List)
89 {
90 int i, line_id, type, tucfield_idx;
91 int line_ucat;
92
93 Vect_reset_list(List);
94
95 tucfield_idx = Vect_cidx_get_field_index(Map, tucfield);
96
97 for (i = 0; i < pSPReport->cArc; i++) {
98 dglEdgeGet_Id(&(Map->dgraph.graph_s), pSPReport->pArc[i].pnEdge);
99
100 line_ucat =
101 dglNodeGet_Id(&(Map->dgraph.graph_s),
102 dglEdgeGet_Head(&(Map->dgraph.graph_s),
103 pSPReport->pArc[i].pnEdge));
104
105 /* get standard ucat numbers (DGLib does not like negative node numbers) */
106 if (line_ucat % 2 == 1)
107 line_ucat = ((line_ucat - 1) / -2);
108 else
109 line_ucat = (line_ucat) / 2;
110
111 /* skip virtual nodes */
112 if (Vect_cidx_find_next
113 (Map, tucfield_idx, abs(line_ucat), GV_LINE, 0, &type,
114 &line_id) == -1)
115 continue;
116
117 if (line_ucat < 0)
118 line_id *= -1;
119
120 G_debug(2, "From %ld to %ld - cost %ld user %d distance %ld", pSPReport->pArc[i].nFrom, pSPReport->pArc[i].nTo, dglEdgeGet_Cost(&(Map->dgraph.graph_s), pSPReport->pArc[i].pnEdge) / Map->dgraph.cost_multip, /* this is the cost from clip() */
121 line_ucat, pSPReport->pArc[i].nDistance);
122
123 Vect_list_append(List, line_id);
124 }
125
126 return 0;
127 }
128
129 /*!
130 \brief Finds shortest path on network using DGLib
131
132
133 \param Map vector map with build DGLib graph (see Vect_net_ttb_build_graph and Vect_net_build_graph)
134 \param from from node id in build the network
135 \param to to node in build the network
136 \param UseTtb the graph is build with/without turntable
137 \param tucfield layer with unique cats for turntable (relevant only when UseTtb = 1)
138 */
find_shortest_path(struct Map_info * Map,int from,int to,struct ilist * List,double * cost,int UseTtb,int tucfield)139 static int find_shortest_path(struct Map_info *Map, int from, int to,
140 struct ilist *List, double *cost, int UseTtb,
141 int tucfield)
142 {
143 int *pclip, cArc, nRet;
144 dglSPReport_s *pSPReport;
145 dglInt32_t nDistance;
146 int use_cache = 1; /* set to 0 to disable dglib cache */
147
148 G_debug(3, "find_shortest_path(): from = %d, to = %d", from, to);
149
150 /* Note : if from == to dgl goes to nearest node and returns back (dgl feature) =>
151 * check here for from == to */
152
153 if (List != NULL)
154 Vect_reset_list(List);
155
156 /* Check if from and to are identical, otherwise dglib returns path to neares node and back! */
157 if (from == to) {
158 if (cost != NULL)
159 *cost = 0;
160 return 0;
161 }
162
163 From_node = from;
164 pclip = NULL;
165 if (List != NULL) {
166 if (use_cache) {
167 nRet =
168 dglShortestPath(&(Map->dgraph.graph_s), &pSPReport,
169 (dglInt32_t) from, (dglInt32_t) to, clipper,
170 pclip, &(Map->dgraph.spCache));
171 }
172 else {
173 nRet =
174 dglShortestPath(&(Map->dgraph.graph_s), &pSPReport,
175 (dglInt32_t) from, (dglInt32_t) to, clipper,
176 pclip, NULL);
177 }
178 }
179 else {
180 if (use_cache) {
181 nRet =
182 dglShortestDistance(&(Map->dgraph.graph_s), &nDistance,
183 (dglInt32_t) from, (dglInt32_t) to,
184 clipper, pclip, &(Map->dgraph.spCache));
185 }
186 else {
187 nRet =
188 dglShortestDistance(&(Map->dgraph.graph_s), &nDistance,
189 (dglInt32_t) from, (dglInt32_t) to,
190 clipper, pclip, NULL);
191 }
192 }
193
194 if (nRet == 0) {
195 /* G_warning("Destination node %d is unreachable from node %d\n" , to , from); */
196 if (cost != NULL)
197 *cost = PORT_DOUBLE_MAX;
198 return -1;
199 }
200 else if (nRet < 0) {
201 G_warning(_("dglShortestPath error: %s"),
202 dglStrerror(&(Map->dgraph.graph_s)));
203 return -1;
204 }
205
206 if (List != NULL) {
207 if (UseTtb)
208 ttb_convert_dgl_shortest_path_result(Map, pSPReport, tucfield,
209 List);
210 else
211 convert_dgl_shortest_path_result(Map, pSPReport, List);
212 }
213
214 if (cost != NULL) {
215 if (List != NULL)
216 *cost = (double)pSPReport->nDistance / Map->dgraph.cost_multip;
217 else
218 *cost = (double)nDistance / Map->dgraph.cost_multip;
219 }
220
221 if (List != NULL) {
222 cArc = pSPReport->cArc;
223 dglFreeSPReport(&(Map->dgraph.graph_s), pSPReport);
224 }
225 else
226 cArc = 0;
227
228 return cArc;
229 }
230
231 /*!
232 \brief Find shortest path on network.
233
234 Costs for 'from' and 'to' nodes are not considered (SP found even if
235 'from' or 'to' are 'closed' (costs = -1) and costs of these
236 nodes are not added to SP costs result.
237
238 \param Map vector map with build graph (see Vect_net_ttb_build_graph and Vect_net_build_graph)
239 \param from start of the path
240 \param from_type if 0 - node id (intersection), if 1 - line unique cat
241 \param to end of the path
242 \param to_type if 0 - node id (intersection), if 1 - line unique cat
243 \param tucfield field with unique categories used in the turntable
244 \param[out] List list of line ids (path)
245 \param[out] cost costs value
246
247 \return number of segments
248 \return 0 is correct for from = to, or List == NULL ? sum of costs is better return value,
249 \return -1 : destination unreachable
250
251 */
252
253 int
Vect_net_ttb_shortest_path(struct Map_info * Map,int from,int from_type,int to,int to_type,int tucfield,struct ilist * List,double * cost)254 Vect_net_ttb_shortest_path(struct Map_info *Map, int from, int from_type,
255 int to, int to_type,
256 int tucfield, struct ilist *List, double *cost)
257 {
258 double x, y, z;
259 struct bound_box box;
260 struct boxlist *box_List;
261 struct line_cats *Cats;
262 int f, t;
263 int i_line, line, type, cfound;
264
265 box_List = Vect_new_boxlist(0);
266 Cats = Vect_new_cats_struct();
267
268 if (from_type == 0) { /* TODO duplicite code with to_type, move into function */
269 /* select points at node */
270 Vect_get_node_coor(Map, from, &x, &y, &z);
271 box.E = box.W = x;
272 box.N = box.S = y;
273 box.T = box.B = z;
274 Vect_select_lines_by_box(Map, &box, GV_POINT, box_List);
275
276 cfound = 0;
277
278 for (i_line = 0; i_line < box_List->n_values; i_line++) {
279 line = box_List->id[i_line];
280
281 type = Vect_read_line(Map, NULL, Cats, line);
282 if (!(type & GV_POINT))
283 continue;
284 if (Vect_cat_get(Cats, tucfield, &f)) {
285 ++cfound;
286 break;
287 }
288 }
289 if (!cfound)
290 G_fatal_error(_
291 ("Unable to find point with defined unique category for node <%d>."),
292 from);
293 else if (cfound > 1)
294 G_warning(_
295 ("There exists more than one point on node <%d> with unique category in field <%d>.\n"
296 "The unique category layer may not be valid."),
297 tucfield, from);
298
299 G_debug(2, "from node = %d, unique cat = %d ", from, f);
300 f = f * 2;
301 }
302 else {
303 if (from < 0)
304 f = from * -2 + 1;
305 else
306 f = from * 2;
307 G_debug(2, "from edge unique cat = %d", from);
308 }
309
310 if (to_type == 0) {
311 /* select points at node */
312 Vect_get_node_coor(Map, to, &x, &y, &z);
313 box.E = box.W = x;
314 box.N = box.S = y;
315 box.T = box.B = z;
316 Vect_select_lines_by_box(Map, &box, GV_POINT, box_List);
317
318 cfound = 0;
319
320 for (i_line = 0; i_line < box_List->n_values; i_line++) {
321 line = box_List->id[i_line];
322 type = Vect_read_line(Map, NULL, Cats, line);
323 if (!(type & GV_POINT))
324 continue;
325 if (Vect_cat_get(Cats, tucfield, &t)) {
326 cfound = 1;
327 break;
328 }
329 }
330 if (!cfound)
331 G_fatal_error(_
332 ("Unable to find point with defined unique category for node <%d>."),
333 to);
334 else if (cfound > 1)
335 G_warning(_
336 ("There exists more than one point on node <%d> with unique category in field <%d>.\n"
337 "The unique category layer may not be valid."),
338 tucfield, to);
339
340 G_debug(2, "to node = %d, unique cat = %d ", to, t);
341 t = t * 2 + 1;
342 }
343 else {
344 if (to < 0)
345 t = to * -2 + 1;
346 else
347 t = to * 2;
348 G_debug(2, "to edge unique cat = %d", to);
349 }
350
351 Vect_destroy_boxlist(box_List);
352 Vect_destroy_cats_struct(Cats);
353
354 return find_shortest_path(Map, f, t, List, cost, 1, tucfield);
355 }
356
357 /*!
358 \brief Find shortest path.
359
360 Costs for 'from' and 'to' nodes are not considered (SP found even if
361 'from' or 'to' are 'closed' (costs = -1) and costs of these
362 nodes are not added to SP costs result.
363
364 \param Map vector map with build graph (see Vect_net_ttb_build_graph and Vect_net_build_graph)
365 \param from from node
366 \param to to node
367 \param[out] List list of line ids (path)
368 \param[out] cost costs value
369
370 \return number of segments
371 \return 0 is correct for from = to, or List == NULL ? sum of costs is better return value,
372 \return -1 : destination unreachable
373
374 */
375 int
Vect_net_shortest_path(struct Map_info * Map,int from,int to,struct ilist * List,double * cost)376 Vect_net_shortest_path(struct Map_info *Map, int from, int to,
377 struct ilist *List, double *cost)
378 {
379 return find_shortest_path(Map, from, to, List, cost, 0, -1);
380 }
381
382 /*!
383 \brief Get graph structure
384
385 Graph is built by Vect_net_build_graph().
386
387 Returns NULL when graph is not built.
388
389 \param Map pointer to Map_info struct
390
391 \return pointer to dglGraph_s struct or NULL
392 */
Vect_net_get_graph(struct Map_info * Map)393 dglGraph_s *Vect_net_get_graph(struct Map_info * Map)
394 {
395 return &(Map->dgraph.graph_s);
396 }
397
398 /*!
399 \brief Returns in cost for given direction in *cost.
400
401 cost is set to -1 if closed.
402
403 \param Map vector map with build graph (see Vect_net_ttb_build_graph and Vect_net_build_graph)
404 \param line line id
405 \param direction direction (GV_FORWARD, GV_BACKWARD)
406 \param[out] cost
407
408 \return 1 OK
409 \return 0 does not exist (was not inserted)
410 */
411 int
Vect_net_get_line_cost(const struct Map_info * Map,int line,int direction,double * cost)412 Vect_net_get_line_cost(const struct Map_info *Map, int line, int direction,
413 double *cost)
414 {
415 /* dglInt32_t *pEdge; */
416
417 G_debug(5, "Vect_net_get_line_cost(): line = %d, dir = %d", line,
418 direction);
419
420 if (direction == GV_FORWARD) {
421 /* V1 has no index by line-id -> array used */
422 /*
423 pEdge = dglGetEdge(&(Map->dgraph.graph_s), line);
424 if (pEdge == NULL)
425 return 0;
426 *cost = (double) dglEdgeGet_Cost(&(Map->dgraph.graph_s), pEdge);
427 */
428 if (Map->dgraph.edge_fcosts[line] == -1) {
429 *cost = -1;
430 return 0;
431 }
432 else
433 *cost = Map->dgraph.edge_fcosts[line];
434 }
435 else if (direction == GV_BACKWARD) {
436 /*
437 pEdge = dglGetEdge(&(Map->dgraph.graph_s), -line);
438 if (pEdge == NULL)
439 return 0;
440 *cost = (double) dglEdgeGet_Cost(&(Map->dgraph.graph_s), pEdge);
441 */
442 if (Map->dgraph.edge_bcosts[line] == -1) {
443 *cost = -1;
444 return 0;
445 }
446 else
447 *cost = Map->dgraph.edge_bcosts[line];
448 G_debug(5, "Vect_net_get_line_cost(): edge_bcosts = %f",
449 Map->dgraph.edge_bcosts[line]);
450 }
451 else {
452 G_fatal_error(_("Wrong line direction in Vect_net_get_line_cost()"));
453 }
454
455 return 1;
456 }
457
458 /*!
459 \brief Get cost of node
460
461 \param Map vector map with build graph (see Vect_net_ttb_build_graph and Vect_net_build_graph)
462 \param node node id
463 \param[out] cost costs value
464
465 \return 1
466 */
Vect_net_get_node_cost(const struct Map_info * Map,int node,double * cost)467 int Vect_net_get_node_cost(const struct Map_info *Map, int node, double *cost)
468 {
469 G_debug(3, "Vect_net_get_node_cost(): node = %d", node);
470
471 *cost = Map->dgraph.node_costs[node];
472
473 G_debug(3, " -> cost = %f", *cost);
474
475 return 1;
476 }
477
478 /*!
479 \brief Find nearest node(s) on network.
480
481 \param Map vector map with build graph (see Vect_net_ttb_build_graph and Vect_net_build_graph)
482 \param x,y,z point coordinates (z coordinate NOT USED !)
483 \param direction (GV_FORWARD - from point to net, GV_BACKWARD - from net to point)
484 \param maxdist maximum distance to the network
485 \param[out] node1 pointer where to store the node number (or NULL)
486 \param[out] node2 pointer where to store the node number (or NULL)
487 \param[out] ln pointer where to store the nearest line number (or NULL)
488 \param[out] costs1 pointer where to store costs on nearest line to node1 (not costs from x,y,z to the line) (or NULL)
489 \param[out] costs2 pointer where to store costs on nearest line to node2 (not costs from x,y,z to the line) (or NULL)
490 \param[out] Points1 pointer to structure where to store vertices on nearest line to node1 (or NULL)
491 \param[out] Points2 pointer to structure where to store vertices on nearest line to node2 (or NULL)
492 \param[out] pointer where to distance to the line (or NULL)
493 \param[out] distance
494
495 \return number of nodes found (0,1,2)
496 */
Vect_net_nearest_nodes(struct Map_info * Map,double x,double y,double z,int direction,double maxdist,int * node1,int * node2,int * ln,double * costs1,double * costs2,struct line_pnts * Points1,struct line_pnts * Points2,double * distance)497 int Vect_net_nearest_nodes(struct Map_info *Map,
498 double x, double y, double z,
499 int direction, double maxdist,
500 int *node1, int *node2, int *ln, double *costs1,
501 double *costs2, struct line_pnts *Points1,
502 struct line_pnts *Points2, double *distance)
503 {
504 int line, n1, n2, nnodes;
505 int npoints;
506 int segment; /* nearest line segment (first is 1) */
507 static struct line_pnts *Points = NULL;
508 double cx, cy, cz, c1, c2;
509 double along; /* distance along the line to nearest point */
510 double length;
511
512 G_debug(3, "Vect_net_nearest_nodes() x = %f y = %f", x, y);
513
514 /* Reset */
515 if (node1)
516 *node1 = 0;
517 if (node2)
518 *node2 = 0;
519 if (ln)
520 *ln = 0;
521 if (costs1)
522 *costs1 = PORT_DOUBLE_MAX;
523 if (costs2)
524 *costs2 = PORT_DOUBLE_MAX;
525 if (Points1)
526 Vect_reset_line(Points1);
527 if (Points2)
528 Vect_reset_line(Points2);
529 if (distance)
530 *distance = PORT_DOUBLE_MAX;
531
532 if (!Points)
533 Points = Vect_new_line_struct();
534
535 /* Find nearest line */
536 line = Vect_find_line(Map, x, y, z, Map->dgraph.line_type, maxdist, 0, 0);
537
538 if (line < 1)
539 return 0;
540
541 Vect_read_line(Map, Points, NULL, line);
542 npoints = Points->n_points;
543 Vect_get_line_nodes(Map, line, &n1, &n2);
544
545 segment =
546 Vect_line_distance(Points, x, y, z, 0, &cx, &cy, &cz, distance, NULL,
547 &along);
548
549 G_debug(4, "line = %d n1 = %d n2 = %d segment = %d", line, n1, n2,
550 segment);
551
552 /* Check first or last point and return one node in that case */
553 G_debug(4, "cx = %f cy = %f first = %f %f last = %f %f", cx, cy,
554 Points->x[0], Points->y[0], Points->x[npoints - 1],
555 Points->y[npoints - 1]);
556
557 if (Points->x[0] == cx && Points->y[0] == cy) {
558 if (node1)
559 *node1 = n1;
560 if (ln)
561 *ln = line;
562 if (costs1)
563 *costs1 = 0;
564 if (Points1) {
565 Vect_append_point(Points1, x, y, z);
566 Vect_append_point(Points1, cx, cy, cz);
567 }
568 G_debug(3, "first node nearest");
569 return 1;
570 }
571 if (Points->x[npoints - 1] == cx && Points->y[npoints - 1] == cy) {
572 if (node1)
573 *node1 = n2;
574 if (ln)
575 *ln = line;
576 if (costs1)
577 *costs1 = 0;
578 if (Points1) {
579 Vect_append_point(Points1, x, y, z);
580 Vect_append_point(Points1, cx, cy, cz);
581 }
582 G_debug(3, "last node nearest");
583 return 1;
584 }
585
586 nnodes = 2;
587
588 /* c1 - costs to get from/to the first vertex */
589 /* c2 - costs to get from/to the last vertex */
590 if (direction == GV_FORWARD) { /* from point to net */
591 Vect_net_get_line_cost(Map, line, GV_BACKWARD, &c1);
592 Vect_net_get_line_cost(Map, line, GV_FORWARD, &c2);
593 }
594 else {
595 Vect_net_get_line_cost(Map, line, GV_FORWARD, &c1);
596 Vect_net_get_line_cost(Map, line, GV_BACKWARD, &c2);
597 }
598
599 if (c1 < 0)
600 nnodes--;
601 if (c2 < 0)
602 nnodes--;
603 if (nnodes == 0)
604 return 0; /* both directions closed */
605
606 length = Vect_line_length(Points);
607
608 if (ln)
609 *ln = line;
610
611 if (nnodes == 1 && c1 < 0) { /* first direction is closed, return node2 as node1 */
612 if (node1)
613 *node1 = n2;
614
615 if (costs1) { /* to node 2, i.e. forward */
616 *costs1 = c2 * (length - along) / length;
617 }
618
619 if (Points1) { /* to node 2, i.e. forward */
620 int i;
621
622 if (direction == GV_FORWARD) { /* from point to net */
623 Vect_append_point(Points1, x, y, z);
624 Vect_append_point(Points1, cx, cy, cz);
625 for (i = segment; i < npoints; i++)
626 Vect_append_point(Points1, Points->x[i], Points->y[i],
627 Points->z[i]);
628 }
629 else {
630 for (i = npoints - 1; i >= segment; i--)
631 Vect_append_point(Points1, Points->x[i], Points->y[i],
632 Points->z[i]);
633
634 Vect_append_point(Points1, cx, cy, cz);
635 Vect_append_point(Points1, x, y, z);
636 }
637 }
638 }
639 else {
640 if (node1)
641 *node1 = n1;
642 if (node2)
643 *node2 = n2;
644
645 if (costs1) { /* to node 1, i.e. backward */
646 *costs1 = c1 * along / length;
647 }
648
649 if (costs2) { /* to node 2, i.e. forward */
650 *costs2 = c2 * (length - along) / length;
651 }
652
653 if (Points1) { /* to node 1, i.e. backward */
654 int i;
655
656 if (direction == GV_FORWARD) { /* from point to net */
657 Vect_append_point(Points1, x, y, z);
658 Vect_append_point(Points1, cx, cy, cz);
659 for (i = segment - 1; i >= 0; i--)
660 Vect_append_point(Points1, Points->x[i], Points->y[i],
661 Points->z[i]);
662 }
663 else {
664 for (i = 0; i < segment; i++)
665 Vect_append_point(Points1, Points->x[i], Points->y[i],
666 Points->z[i]);
667
668 Vect_append_point(Points1, cx, cy, cz);
669 Vect_append_point(Points1, x, y, z);
670 }
671 }
672
673 if (Points2) { /* to node 2, i.e. forward */
674 int i;
675
676 if (direction == GV_FORWARD) { /* from point to net */
677 Vect_append_point(Points2, x, y, z);
678 Vect_append_point(Points2, cx, cy, cz);
679 for (i = segment; i < npoints; i++)
680 Vect_append_point(Points2, Points->x[i], Points->y[i],
681 Points->z[i]);
682 }
683 else {
684 for (i = npoints - 1; i >= segment; i--)
685 Vect_append_point(Points2, Points->x[i], Points->y[i],
686 Points->z[i]);
687
688 Vect_append_point(Points2, cx, cy, cz);
689 Vect_append_point(Points2, x, y, z);
690 }
691 }
692 }
693
694 return nnodes;
695 }
696
697 /*!
698 \brief Find shortest path on network between 2 points given by coordinates.
699
700 \param Map vector map with build graph (see Vect_net_ttb_build_graph and Vect_net_build_graph)
701 \param fx,fy,fz from point x coordinate (z ignored)
702 \param tx,ty,tz to point x coordinate (z ignored)
703 \param fmax maximum distance to the network from 'from'
704 \param tmax maximum distance to the network from 'to'
705 \param UseTtb the graph is build with/without turntable
706 \param tucfield field with unique categories used in the turntable
707 \param costs pointer where to store costs on the network (or NULL)
708 \param Points pointer to the structure where to store vertices of shortest path (or NULL)
709 \param List pointer to the structure where list of lines on the network is stored (or NULL)
710 \param NodesList pointer to the structure where list of nodes on the network is stored (or NULL)
711 \param FPoints pointer to the structure where to store line from 'from' to first network node (or NULL)
712 \param TPoints pointer to the structure where to store line from last network node to 'to' (or NULL)
713 \param fdist distance from 'from' to the net (or NULL)
714 \param tdist distance from 'to' to the net (or NULL)
715
716 \return 1 OK, 0 not reachable
717 */
718 static int
find_shortest_path_coor(struct Map_info * Map,double fx,double fy,double fz,double tx,double ty,double tz,double fmax,double tmax,int UseTtb,int tucfield,double * costs,struct line_pnts * Points,struct ilist * List,struct ilist * NodesList,struct line_pnts * FPoints,struct line_pnts * TPoints,double * fdist,double * tdist)719 find_shortest_path_coor(struct Map_info *Map,
720 double fx, double fy, double fz, double tx,
721 double ty, double tz, double fmax, double tmax,
722 int UseTtb, int tucfield,
723 double *costs, struct line_pnts *Points,
724 struct ilist *List, struct ilist *NodesList,
725 struct line_pnts *FPoints,
726 struct line_pnts *TPoints, double *fdist,
727 double *tdist)
728 {
729 int fnode[2], tnode[2]; /* nearest nodes, *node[1] is 0 if only one was found */
730 double fcosts[2], tcosts[2], cur_cst; /* costs to nearest nodes on the network */
731 int nfnodes, ntnodes, fline, tline;
732 static struct line_pnts *APoints, *SPoints, *fPoints[2], *tPoints[2];
733 static struct ilist *LList;
734 static int first = 1;
735 int reachable, shortcut;
736 int i, j, fn = 0, tn = 0;
737
738 /* from/to_point_node is set if from/to point projected to line
739 *falls exactly on node (shortcut -> fline == tline) */
740 int from_point_node = 0;
741 int to_point_node = 0;
742
743 G_debug(3, "Vect_net_shortest_path_coor()");
744
745 if (first) {
746 APoints = Vect_new_line_struct();
747 SPoints = Vect_new_line_struct();
748 fPoints[0] = Vect_new_line_struct();
749 fPoints[1] = Vect_new_line_struct();
750 tPoints[0] = Vect_new_line_struct();
751 tPoints[1] = Vect_new_line_struct();
752 LList = Vect_new_list();
753 first = 0;
754 }
755
756 /* Reset */
757 if (costs)
758 *costs = PORT_DOUBLE_MAX;
759 if (Points)
760 Vect_reset_line(Points);
761 if (fdist)
762 *fdist = 0;
763 if (tdist)
764 *tdist = 0;
765 if (List)
766 List->n_values = 0;
767 if (FPoints)
768 Vect_reset_line(FPoints);
769 if (TPoints)
770 Vect_reset_line(TPoints);
771 if (NodesList != NULL)
772 Vect_reset_list(NodesList);
773
774 /* Find nearest nodes */
775 fnode[0] = fnode[1] = tnode[0] = tnode[1] = 0;
776
777 nfnodes =
778 Vect_net_nearest_nodes(Map, fx, fy, fz, GV_FORWARD, fmax, &(fnode[0]),
779 &(fnode[1]), &fline, &(fcosts[0]),
780 &(fcosts[1]), fPoints[0], fPoints[1], fdist);
781 if (nfnodes == 0)
782 return 0;
783
784 if (nfnodes == 1 && fPoints[0]->n_points < 3) {
785 from_point_node = fnode[0];
786 }
787
788 ntnodes =
789 Vect_net_nearest_nodes(Map, tx, ty, tz, GV_BACKWARD, tmax,
790 &(tnode[0]), &(tnode[1]), &tline, &(tcosts[0]),
791 &(tcosts[1]), tPoints[0], tPoints[1], tdist);
792 if (ntnodes == 0)
793 return 0;
794
795 if (ntnodes == 1 && tPoints[0]->n_points < 3) {
796 to_point_node = tnode[0];
797 }
798
799
800 G_debug(3, "fline = %d tline = %d", fline, tline);
801
802 reachable = shortcut = 0;
803 cur_cst = PORT_DOUBLE_MAX;
804
805 /* It may happen, that 2 points are at the same line. */
806 /* TODO?: it could also happen that fline != tline but both points are on the same
807 * line if they fall on node but a different line was found. This case is correctly
808 * handled as normal non shortcut, but it could be added here. In that case
809 * NodesList collection must be changed */
810 if (fline == tline && (nfnodes > 1 || ntnodes > 1)) {
811 double len, flen, tlen, c, fseg, tseg;
812 double fcx, fcy, fcz, tcx, tcy, tcz;
813
814 Vect_read_line(Map, APoints, NULL, fline);
815 len = Vect_line_length(APoints);
816
817 /* distance along the line */
818 fseg =
819 Vect_line_distance(APoints, fx, fy, fz, 0, &fcx, &fcy, &fcz, NULL,
820 NULL, &flen);
821 tseg =
822 Vect_line_distance(APoints, tx, ty, tz, 0, &tcx, &tcy, &tcz, NULL,
823 NULL, &tlen);
824
825 Vect_reset_line(SPoints);
826 if (flen == tlen) {
827 cur_cst = 0;
828
829 Vect_append_point(SPoints, fx, fy, fz);
830 Vect_append_point(SPoints, fcx, fcy, fcz);
831 Vect_append_point(SPoints, tx, ty, tz);
832
833 reachable = shortcut = 1;
834 }
835 else if (flen < tlen) {
836 Vect_net_get_line_cost(Map, fline, GV_FORWARD, &c);
837 if (c >= 0) {
838 cur_cst = c * (tlen - flen) / len;
839
840 Vect_append_point(SPoints, fx, fy, fz);
841 Vect_append_point(SPoints, fcx, fcy, fcz);
842 for (i = fseg; i < tseg; i++)
843 Vect_append_point(SPoints, APoints->x[i], APoints->y[i],
844 APoints->z[i]);
845
846 Vect_append_point(SPoints, tcx, tcy, tcz);
847 Vect_append_point(SPoints, tx, ty, tz);
848
849 reachable = shortcut = 1;
850 }
851 }
852 else { /* flen > tlen */
853 Vect_net_get_line_cost(Map, fline, GV_BACKWARD, &c);
854 if (c >= 0) {
855 cur_cst = c * (flen - tlen) / len;
856
857 Vect_append_point(SPoints, fx, fy, fz);
858 Vect_append_point(SPoints, fcx, fcy, fcz);
859 for (i = fseg - 1; i >= tseg; i--)
860 Vect_append_point(SPoints, APoints->x[i], APoints->y[i],
861 APoints->z[i]);
862
863 Vect_append_point(SPoints, tcx, tcy, tcz);
864 Vect_append_point(SPoints, tx, ty, tz);
865
866 reachable = shortcut = 1;
867 }
868 }
869 }
870
871 /* Find the shortest variant from maximum 4 */
872 for (i = 0; i < nfnodes; i++) {
873 for (j = 0; j < ntnodes; j++) {
874 double ncst, cst;
875 int ret;
876
877 G_debug(3, "i = %d fnode = %d j = %d tnode = %d", i, fnode[i], j,
878 tnode[j]);
879
880 if (UseTtb)
881 ret =
882 Vect_net_ttb_shortest_path(Map, fnode[i], 0, tnode[j], 0,
883 tucfield, NULL, &ncst);
884 else
885 ret =
886 Vect_net_shortest_path(Map, fnode[i], tnode[j], NULL,
887 &ncst);
888 if (ret == -1)
889 continue; /* not reachable */
890
891 cst = fcosts[i] + ncst + tcosts[j];
892 if (reachable == 0 || cst < cur_cst) {
893 cur_cst = cst;
894 fn = i;
895 tn = j;
896 shortcut = 0;
897 }
898 reachable = 1;
899 }
900 }
901
902 G_debug(3, "reachable = %d shortcut = %d cur_cst = %f", reachable,
903 shortcut, cur_cst);
904 if (reachable) {
905 if (shortcut) {
906 if (Points)
907 Vect_append_points(Points, SPoints, GV_FORWARD);
908 if (NodesList) {
909 /* Check if from/to point projected to line falls on node and
910 *add it to the list */
911 if (from_point_node > 0)
912 Vect_list_append(NodesList, from_point_node);
913
914 if (to_point_node > 0)
915 Vect_list_append(NodesList, to_point_node);
916 }
917 }
918 else {
919 if (NodesList) {
920 /* it can happen that starting point falls on node but SP starts
921 * form the other node, add it in that case,
922 * similarly for to point below */
923 if (from_point_node > 0 && from_point_node != fnode[fn]) {
924 Vect_list_append(NodesList, from_point_node);
925 }
926
927 /* add starting net SP search node */
928 Vect_list_append(NodesList, fnode[fn]);
929 }
930
931 if (UseTtb)
932 Vect_net_ttb_shortest_path(Map, fnode[fn], 0, tnode[tn], 0,
933 tucfield, LList, NULL);
934 else
935 Vect_net_shortest_path(Map, fnode[fn], tnode[tn], LList,
936 NULL);
937
938 G_debug(3, "Number of lines %d", LList->n_values);
939
940 if (Points)
941 Vect_append_points(Points, fPoints[fn], GV_FORWARD);
942
943 if (FPoints)
944 Vect_append_points(FPoints, fPoints[fn], GV_FORWARD);
945
946 for (i = 0; i < LList->n_values; i++) {
947 int line;
948
949 line = LList->value[i];
950 G_debug(3, "i = %d line = %d", i, line);
951
952 if (Points) {
953 Vect_read_line(Map, APoints, NULL, abs(line));
954
955 if (line > 0)
956 Vect_append_points(Points, APoints, GV_FORWARD);
957 else
958 Vect_append_points(Points, APoints, GV_BACKWARD);
959 Points->n_points--;
960 }
961 if (NodesList) {
962 int node, node1, node2;
963
964 Vect_get_line_nodes(Map, abs(line), &node1, &node2);
965 /* add the second node, the first of first segmet was alread added */
966 if (line > 0)
967 node = node2;
968 else
969 node = node1;
970
971 Vect_list_append(NodesList, node);
972 }
973
974 if (List)
975 Vect_list_append(List, line);
976 }
977
978 if (Points) {
979 if (LList->n_values)
980 Points->n_points++;
981 Vect_append_points(Points, tPoints[tn], GV_FORWARD);
982 }
983
984 if (TPoints)
985 Vect_append_points(TPoints, tPoints[tn], GV_FORWARD);
986
987 if (NodesList) {
988 if (to_point_node > 0 && to_point_node != tnode[tn]) {
989 Vect_list_append(NodesList, to_point_node);
990 }
991 }
992 }
993
994 if (costs)
995 *costs = cur_cst;
996 if (Points)
997 Vect_line_prune(Points);
998 }
999
1000 return reachable;
1001 }
1002
1003 /*!
1004 \brief Find shortest path on network between 2 points given by coordinates.
1005
1006 \param Map vector map with build graph (see Vect_net_ttb_build_graph and Vect_net_build_graph)
1007 \param fx,fy,fz from point x coordinate (z ignored)
1008 \param tx,ty,tz to point x coordinate (z ignored)
1009 \param fmax maximum distance to the network from 'from'
1010 \param tmax maximum distance to the network from 'to'
1011 \param costs pointer where to store costs on the network (or NULL)
1012 \param Points pointer to the structure where to store vertices of shortest path (or NULL)
1013 \param List pointer to the structure where list of lines on the network is stored (or NULL)
1014 \param NodesList pointer to the structure where list of nodes on the network is stored (or NULL)
1015 \param FPoints pointer to the structure where to store line from 'from' to first network node (or NULL)
1016 \param TPoints pointer to the structure where to store line from last network node to 'to' (or NULL)
1017 \param fdist distance from 'from' to the net (or NULL)
1018 \param tdist distance from 'to' to the net (or NULL)
1019
1020 \return 1 OK, 0 not reachable
1021 */
1022 int
Vect_net_shortest_path_coor(struct Map_info * Map,double fx,double fy,double fz,double tx,double ty,double tz,double fmax,double tmax,double * costs,struct line_pnts * Points,struct ilist * List,struct ilist * NodesList,struct line_pnts * FPoints,struct line_pnts * TPoints,double * fdist,double * tdist)1023 Vect_net_shortest_path_coor(struct Map_info *Map,
1024 double fx, double fy, double fz, double tx,
1025 double ty, double tz, double fmax, double tmax,
1026 double *costs, struct line_pnts *Points,
1027 struct ilist *List, struct ilist *NodesList,
1028 struct line_pnts *FPoints,
1029 struct line_pnts *TPoints, double *fdist,
1030 double *tdist)
1031 {
1032 return find_shortest_path_coor(Map, fx, fy, fz, tx, ty, tz, fmax, tmax, 0,
1033 0, costs, Points, List, NodesList,
1034 FPoints, TPoints, fdist, tdist);
1035 }
1036
1037 /*!
1038 \brief Find shortest path on network with turntable between 2 points given by coordinates.
1039
1040 \param Map vector map with build graph (see Vect_net_ttb_build_graph and Vect_net_build_graph)
1041 \param fx,fy,fz from point x coordinate (z ignored)
1042 \param tx,ty,tz to point x coordinate (z ignored)
1043 \param fmax maximum distance to the network from 'from'
1044 \param tmax maximum distance to the network from 'to'
1045 \param tucfield field with unique categories used in the turntable
1046 \param costs pointer where to store costs on the network (or NULL)
1047 \param Points pointer to the structure where to store vertices of shortest path (or NULL)
1048 \param List pointer to the structure where list of lines on the network is stored (or NULL)
1049 \param NodesList pointer to the structure where list of nodes on the network is stored (or NULL)
1050 \param FPoints pointer to the structure where to store line from 'from' to first network node (or NULL)
1051 \param TPoints pointer to the structure where to store line from last network node to 'to' (or NULL)
1052 \param fdist distance from 'from' to the net (or NULL)
1053 \param tdist distance from 'to' to the net (or NULL)
1054
1055 \return 1 OK, 0 not reachable
1056 */
1057 int
Vect_net_ttb_shortest_path_coor(struct Map_info * Map,double fx,double fy,double fz,double tx,double ty,double tz,double fmax,double tmax,int tucfield,double * costs,struct line_pnts * Points,struct ilist * List,struct ilist * NodesList,struct line_pnts * FPoints,struct line_pnts * TPoints,double * fdist,double * tdist)1058 Vect_net_ttb_shortest_path_coor(struct Map_info *Map,
1059 double fx, double fy, double fz, double tx,
1060 double ty, double tz, double fmax,
1061 double tmax, int tucfield,
1062 double *costs, struct line_pnts *Points,
1063 struct ilist *List, struct ilist *NodesList,
1064 struct line_pnts *FPoints,
1065 struct line_pnts *TPoints, double *fdist,
1066 double *tdist)
1067 {
1068 return find_shortest_path_coor(Map, fx, fy, fz, tx, ty, tz, fmax, tmax,
1069 1, tucfield, costs, Points, List,
1070 NodesList, FPoints, TPoints, fdist, tdist);
1071 }
1072