GRASS GIS 7 Programmer's Manual  7.5.svn(2018)-r72846
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
net_analyze.c
Go to the documentation of this file.
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 
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  */
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  */
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),
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 */
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  */
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
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)
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)
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);
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
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  */
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
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  */
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  */
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
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
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
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 }
struct boxlist * Vect_new_boxlist(int have_boxes)
Creates and initializes a struct boxlist.
int Vect_select_lines_by_box(struct Map_info *Map, const struct bound_box *Box, int type, struct boxlist *list)
Select lines with bounding boxes by box.
Definition: sindex.c:34
Bounding box.
Definition: dig_structs.h:65
#define GV_FORWARD
Line direction indicator forward/backward.
Definition: dig_defines.h:178
int * id
Array of ids.
Definition: dig_structs.h:1755
int Vect_cidx_find_next(const struct Map_info *Map, int field_index, int cat, int type_mask, int start_index, int *type, int *id)
Find next line/area id for given category, start_index and type_mask.
Definition: Vlib/cindex.c:325
double W
West.
Definition: dig_structs.h:82
int Vect_net_shortest_path(struct Map_info *Map, int from, int to, struct ilist *List, double *cost)
Find shortest path.
Definition: net_analyze.c:376
int Vect_get_node_coor(const struct Map_info *Map, int num, double *x, double *y, double *z)
Get node coordinates.
Definition: level_two.c:278
dglInt32_t * pnEdge
Definition: graph.h:215
#define GV_BACKWARD
Definition: dig_defines.h:179
struct line_pnts * Vect_new_line_struct()
Creates and initializes a line_pnts structure.
Definition: line.c:45
struct ilist * Vect_new_list(void)
Creates and initializes a struct ilist.
dglInt32_t * pnEdge
Definition: graph.h:96
dglInt32_t dglNodeGet_Id(dglGraph_s *pGraph, dglInt32_t *pnNode)
int n_points
Number of points.
Definition: dig_structs.h:1692
dglInt32_t nDistance
Definition: graph.h:216
void Vect_destroy_boxlist(struct boxlist *list)
Frees all memory associated with a struct boxlist, including the struct itself.
dglInt32_t dglEdgeGet_Cost(dglGraph_s *pGraph, dglInt32_t *pnEdge)
int Vect_cidx_get_field_index(const struct Map_info *Map, int field)
Get layer index for given layer number.
Definition: Vlib/cindex.c:113
int n_values
Number of values in the list.
Definition: gis.h:676
double * edge_fcosts
Forward costs used for graph.
Definition: dig_structs.h:1238
dglInt32_t * dglEdgeGet_Head(dglGraph_s *pGraph, dglInt32_t *pnEdge)
double E
East.
Definition: dig_structs.h:78
void Vect_destroy_cats_struct(struct line_cats *p)
Frees all memory associated with line_cats structure, including the struct itself.
#define NULL
Definition: ccmath.h:32
#define x
double * edge_bcosts
backward costs used for graph
Definition: dig_structs.h:1242
int n_values
Number of items in the list.
Definition: dig_structs.h:1767
#define GV_POINT
Feature types used in memory on run time (may change)
Definition: dig_defines.h:182
double Vect_line_length(const struct line_pnts *Points)
Calculate line length, 3D-length in case of 3D vector line.
Definition: line.c:576
int Vect_append_point(struct line_pnts *Points, double x, double y, double z)
Appends one point to the end of a line.
Definition: line.c:149
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:160
Feature category info.
Definition: dig_structs.h:1702
dglInt32_t nTo
Definition: graph.h:214
int Vect_append_points(struct line_pnts *Points, const struct line_pnts *APoints, int direction)
Appends points to the end of a line.
Definition: line.c:337
#define GV_LINE
Definition: dig_defines.h:183
double N
North.
Definition: dig_structs.h:70
struct Graph_info dgraph
Graph info (built for network analysis)
Definition: dig_structs.h:1398
int type
Definition: gis.h:512
double * x
Array of X coordinates.
Definition: dig_structs.h:1680
Feature geometry info - coordinates.
Definition: dig_structs.h:1675
double t
Definition: r_raster.c:39
char * dglStrerror(dglGraph_s *pgraph)
dglGraph_s graph_s
Graph structure.
Definition: dig_structs.h:1228
int Vect_cat_get(const struct line_cats *Cats, int field, int *cat)
Get first found category of given field.
int Vect_line_distance(const struct line_pnts *points, double ux, double uy, double uz, int with_z, double *px, double *py, double *pz, double *dist, double *spdist, double *lpdist)
Calculate distance of point to line.
Definition: line.c:651
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
dglSPArc_s * pArc
Definition: graph.h:229
int Vect_net_get_line_cost(const struct Map_info *Map, int line, int direction, double *cost)
Returns in cost for given direction in *cost.
Definition: net_analyze.c:412
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)
Find shortest path on network between 2 points given by coordinates.
Definition: net_analyze.c:1023
double B
Bottom.
Definition: dig_structs.h:90
int Vect_reset_list(struct ilist *list)
Reset ilist structure.
long dglInt32_t
Definition: type.h:37
int Vect_list_append(struct ilist *list, int val)
Append new item to the end of list if not yet present.
double T
Top.
Definition: dig_structs.h:86
dglInt32_t nFrom
Definition: graph.h:213
int dglShortestDistance(dglGraph_s *pGraph, dglInt32_t *pnDistance, dglInt32_t nStart, dglInt32_t nDestination, dglSPClip_fn fnClip, void *pvClipArg, dglSPCache_s *pCache)
dglSPCache_s spCache
Shortest path cache.
Definition: dig_structs.h:1232
dglInt32_t nEdgeCost
Definition: graph.h:104
Vector map info.
Definition: dig_structs.h:1259
#define PORT_DOUBLE_MAX
Limits for portable types.
Definition: dig_defines.h:66
double * y
Array of Y coordinates.
Definition: dig_structs.h:1684
int Vect_get_line_nodes(const struct Map_info *Map, int line, int *n1, int *n2)
Get line nodes.
Definition: level_two.c:307
struct line_cats * Vect_new_cats_struct()
Creates and initializes line_cats structure.
int Vect_line_prune(struct line_pnts *Points)
Remove duplicate points, i.e. zero length segments.
Definition: line.c:281
dglInt32_t * pnNodeFrom
Definition: graph.h:95
dglInt32_t * pnNodeTo
Definition: graph.h:97
List of bounding boxes with id.
Definition: dig_structs.h:1750
int cost_multip
Edge and node costs multiplicator.
Definition: dig_structs.h:1250
double S
South.
Definition: dig_structs.h:74
int Vect_net_get_node_cost(const struct Map_info *Map, int node, double *cost)
Get cost of node.
Definition: net_analyze.c:467
dglInt32_t cArc
Definition: graph.h:228
dglInt32_t * dglNodeGet_Attr(dglGraph_s *pGraph, dglInt32_t *pnNode)
double * z
Array of Z coordinates.
Definition: dig_structs.h:1688
#define _(str)
Definition: glocale.h:13
List of integers.
Definition: gis.h:667
int Vect_read_line(const struct Map_info *Map, struct line_pnts *line_p, struct line_cats *line_c, int line)
Read vector feature (topological level required)
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)
Find shortest path on network with turntable between 2 points given by coordinates.
Definition: net_analyze.c:1058
void Vect_reset_line(struct line_pnts *Points)
Reset line.
Definition: line.c:130
int * value
Array of values.
Definition: gis.h:672
int 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)
Find nearest node(s) on network.
Definition: net_analyze.c:497
double * node_costs
Node costs used for graph.
Definition: dig_structs.h:1246
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)
Find shortest path on network.
Definition: net_analyze.c:254
dglGraph_s * Vect_net_get_graph(struct Map_info *Map)
Get graph structure.
Definition: net_analyze.c:393
int line_type
Line type used to build the graph.
Definition: dig_structs.h:1224
int Vect_find_line(struct Map_info *map, double ux, double uy, double uz, int type, double maxdist, int with_z, int exclude)
Find the nearest line.
dglInt32_t nDistance
Definition: graph.h:227
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204
dglInt32_t dglEdgeGet_Id(dglGraph_s *pGraph, dglInt32_t *pnEdge)
void dglFreeSPReport(dglGraph_s *pgraph, dglSPReport_s *pSPReport)
int dglGet_NodeAttrSize(dglGraph_s *pgraph)
int dglShortestPath(dglGraph_s *pGraph, dglSPReport_s **ppReport, dglInt32_t nStart, dglInt32_t nDestination, dglSPClip_fn fnClip, void *pvClipArg, dglSPCache_s *pCache)