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