GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-44ea597914
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) \param x,y,z point coordinates (z coordinate NOT USED
484  !) \param direction (GV_FORWARD - from point to net, GV_BACKWARD - from net
485  to point) \param maxdist maximum distance to the network \param[out] node1
486  pointer where to store the node number (or NULL) \param[out] node2 pointer
487  where to store the node number (or NULL) \param[out] ln pointer where to
488  store the nearest line number (or NULL) \param[out] costs1 pointer where to
489  store costs on nearest line to node1 (not costs from x,y,z to the line) (or
490  NULL) \param[out] costs2 pointer where to store costs on nearest line to
491  node2 (not costs from x,y,z to the line) (or NULL) \param[out] Points1
492  pointer to structure where to store vertices on nearest line to node1 (or
493  NULL) \param[out] Points2 pointer to structure where to store vertices on
494  nearest line to node2 (or NULL) \param[out] pointer where to distance to the
495  line (or NULL) \param[out] distance
496 
497  \return number of nodes found (0,1,2)
498  */
499 int Vect_net_nearest_nodes(struct Map_info *Map, double x, double y, double z,
500  int direction, double maxdist, int *node1,
501  int *node2, int *ln, double *costs1, double *costs2,
502  struct line_pnts *Points1, struct line_pnts *Points2,
503  double *distance)
504 {
505  int line, n1, n2, nnodes;
506  int npoints;
507  int segment; /* nearest line segment (first is 1) */
508  static struct line_pnts *Points = NULL;
509  double cx, cy, cz, c1, c2;
510  double along; /* distance along the line to nearest point */
511  double length;
512 
513  G_debug(3, "Vect_net_nearest_nodes() x = %f y = %f", x, y);
514 
515  /* Reset */
516  if (node1)
517  *node1 = 0;
518  if (node2)
519  *node2 = 0;
520  if (ln)
521  *ln = 0;
522  if (costs1)
523  *costs1 = PORT_DOUBLE_MAX;
524  if (costs2)
525  *costs2 = PORT_DOUBLE_MAX;
526  if (Points1)
527  Vect_reset_line(Points1);
528  if (Points2)
529  Vect_reset_line(Points2);
530  if (distance)
531  *distance = PORT_DOUBLE_MAX;
532 
533  if (!Points)
534  Points = Vect_new_line_struct();
535 
536  /* Find nearest line */
537  line = Vect_find_line(Map, x, y, z, Map->dgraph.line_type, maxdist, 0, 0);
538 
539  if (line < 1)
540  return 0;
541 
542  Vect_read_line(Map, Points, NULL, line);
543  npoints = Points->n_points;
544  Vect_get_line_nodes(Map, line, &n1, &n2);
545 
546  segment = Vect_line_distance(Points, x, y, z, 0, &cx, &cy, &cz, distance,
547  NULL, &along);
548 
549  G_debug(4, "line = %d n1 = %d n2 = %d segment = %d", line, n1, n2, segment);
550 
551  /* Check first or last point and return one node in that case */
552  G_debug(4, "cx = %f cy = %f first = %f %f last = %f %f", cx, cy,
553  Points->x[0], Points->y[0], Points->x[npoints - 1],
554  Points->y[npoints - 1]);
555 
556  if (Points->x[0] == cx && Points->y[0] == cy) {
557  if (node1)
558  *node1 = n1;
559  if (ln)
560  *ln = line;
561  if (costs1)
562  *costs1 = 0;
563  if (Points1) {
564  Vect_append_point(Points1, x, y, z);
565  Vect_append_point(Points1, cx, cy, cz);
566  }
567  G_debug(3, "first node nearest");
568  return 1;
569  }
570  if (Points->x[npoints - 1] == cx && Points->y[npoints - 1] == cy) {
571  if (node1)
572  *node1 = n2;
573  if (ln)
574  *ln = line;
575  if (costs1)
576  *costs1 = 0;
577  if (Points1) {
578  Vect_append_point(Points1, x, y, z);
579  Vect_append_point(Points1, cx, cy, cz);
580  }
581  G_debug(3, "last node nearest");
582  return 1;
583  }
584 
585  nnodes = 2;
586 
587  /* c1 - costs to get from/to the first vertex */
588  /* c2 - costs to get from/to the last vertex */
589  if (direction == GV_FORWARD) { /* from point to net */
590  Vect_net_get_line_cost(Map, line, GV_BACKWARD, &c1);
591  Vect_net_get_line_cost(Map, line, GV_FORWARD, &c2);
592  }
593  else {
594  Vect_net_get_line_cost(Map, line, GV_FORWARD, &c1);
595  Vect_net_get_line_cost(Map, line, GV_BACKWARD, &c2);
596  }
597 
598  if (c1 < 0)
599  nnodes--;
600  if (c2 < 0)
601  nnodes--;
602  if (nnodes == 0)
603  return 0; /* both directions closed */
604 
605  length = Vect_line_length(Points);
606 
607  if (ln)
608  *ln = line;
609 
610  if (nnodes == 1 &&
611  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
701  Vect_net_build_graph) \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
709  path (or NULL) \param List pointer to the structure where list of lines on
710  the network is stored (or NULL) \param NodesList pointer to the structure
711  where list of nodes on the network is stored (or NULL) \param FPoints pointer
712  to the structure where to store line from 'from' to first network node (or
713  NULL) \param TPoints pointer to the structure where to store line from last
714  network node to 'to' (or NULL) \param fdist distance from 'from' to the net
715  (or NULL) \param tdist distance from 'to' to the net (or NULL)
716 
717  \return 1 OK, 0 not reachable
718  */
719 static int
720 find_shortest_path_coor(struct Map_info *Map, double fx, double fy, double fz,
721  double tx, double ty, double tz, double fmax,
722  double tmax, int UseTtb, int tucfield, double *costs,
723  struct line_pnts *Points, struct ilist *List,
724  struct ilist *NodesList, struct line_pnts *FPoints,
725  struct line_pnts *TPoints, double *fdist, double *tdist)
726 {
727  int fnode[2],
728  tnode[2]; /* nearest nodes, *node[1] is 0 if only one was found */
729  double fcosts[2], tcosts[2],
730  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 = Vect_net_nearest_nodes(
778  Map, fx, fy, fz, GV_FORWARD, fmax, &(fnode[0]), &(fnode[1]), &fline,
779  &(fcosts[0]), &(fcosts[1]), fPoints[0], fPoints[1], fdist);
780  if (nfnodes == 0)
781  return 0;
782 
783  if (nfnodes == 1 && fPoints[0]->n_points < 3) {
784  from_point_node = fnode[0];
785  }
786 
787  ntnodes = Vect_net_nearest_nodes(
788  Map, tx, ty, tz, GV_BACKWARD, tmax, &(tnode[0]), &(tnode[1]), &tline,
789  &(tcosts[0]), &(tcosts[1]), tPoints[0], tPoints[1], tdist);
790  if (ntnodes == 0)
791  return 0;
792 
793  if (ntnodes == 1 && tPoints[0]->n_points < 3) {
794  to_point_node = tnode[0];
795  }
796 
797  G_debug(3, "fline = %d tline = %d", fline, tline);
798 
799  reachable = shortcut = 0;
800  cur_cst = PORT_DOUBLE_MAX;
801 
802  /* It may happen, that 2 points are at the same line. */
803  /* TODO?: it could also happen that fline != tline but both points are on
804  * the same line if they fall on node but a different line was found. This
805  * case is correctly handled as normal non shortcut, but it could be added
806  * here. In that case NodesList collection must be changed */
807  if (fline == tline && (nfnodes > 1 || ntnodes > 1)) {
808  double len, flen, tlen, c, fseg, tseg;
809  double fcx, fcy, fcz, tcx, tcy, tcz;
810 
811  Vect_read_line(Map, APoints, NULL, fline);
812  len = Vect_line_length(APoints);
813 
814  /* distance along the line */
815  fseg = Vect_line_distance(APoints, fx, fy, fz, 0, &fcx, &fcy, &fcz,
816  NULL, NULL, &flen);
817  tseg = Vect_line_distance(APoints, tx, ty, tz, 0, &tcx, &tcy, &tcz,
818  NULL, NULL, &tlen);
819 
820  Vect_reset_line(SPoints);
821  if (flen == tlen) {
822  cur_cst = 0;
823 
824  Vect_append_point(SPoints, fx, fy, fz);
825  Vect_append_point(SPoints, fcx, fcy, fcz);
826  Vect_append_point(SPoints, tx, ty, tz);
827 
828  reachable = shortcut = 1;
829  }
830  else if (flen < tlen) {
831  Vect_net_get_line_cost(Map, fline, GV_FORWARD, &c);
832  if (c >= 0) {
833  cur_cst = c * (tlen - flen) / len;
834 
835  Vect_append_point(SPoints, fx, fy, fz);
836  Vect_append_point(SPoints, fcx, fcy, fcz);
837  for (i = fseg; i < tseg; i++)
838  Vect_append_point(SPoints, APoints->x[i], APoints->y[i],
839  APoints->z[i]);
840 
841  Vect_append_point(SPoints, tcx, tcy, tcz);
842  Vect_append_point(SPoints, tx, ty, tz);
843 
844  reachable = shortcut = 1;
845  }
846  }
847  else { /* flen > tlen */
848  Vect_net_get_line_cost(Map, fline, GV_BACKWARD, &c);
849  if (c >= 0) {
850  cur_cst = c * (flen - tlen) / len;
851 
852  Vect_append_point(SPoints, fx, fy, fz);
853  Vect_append_point(SPoints, fcx, fcy, fcz);
854  for (i = fseg - 1; i >= tseg; i--)
855  Vect_append_point(SPoints, APoints->x[i], APoints->y[i],
856  APoints->z[i]);
857 
858  Vect_append_point(SPoints, tcx, tcy, tcz);
859  Vect_append_point(SPoints, tx, ty, tz);
860 
861  reachable = shortcut = 1;
862  }
863  }
864  }
865 
866  /* Find the shortest variant from maximum 4 */
867  for (i = 0; i < nfnodes; i++) {
868  for (j = 0; j < ntnodes; j++) {
869  double ncst, cst;
870  int ret;
871 
872  G_debug(3, "i = %d fnode = %d j = %d tnode = %d", i, fnode[i], j,
873  tnode[j]);
874 
875  if (UseTtb)
876  ret = Vect_net_ttb_shortest_path(Map, fnode[i], 0, tnode[j], 0,
877  tucfield, NULL, &ncst);
878  else
879  ret = Vect_net_shortest_path(Map, fnode[i], tnode[j], NULL,
880  &ncst);
881  if (ret == -1)
882  continue; /* not reachable */
883 
884  cst = fcosts[i] + ncst + tcosts[j];
885  if (reachable == 0 || cst < cur_cst) {
886  cur_cst = cst;
887  fn = i;
888  tn = j;
889  shortcut = 0;
890  }
891  reachable = 1;
892  }
893  }
894 
895  G_debug(3, "reachable = %d shortcut = %d cur_cst = %f", reachable, shortcut,
896  cur_cst);
897  if (reachable) {
898  if (shortcut) {
899  if (Points)
900  Vect_append_points(Points, SPoints, GV_FORWARD);
901  if (NodesList) {
902  /* Check if from/to point projected to line falls on node and
903  *add it to the list */
904  if (from_point_node > 0)
905  Vect_list_append(NodesList, from_point_node);
906 
907  if (to_point_node > 0)
908  Vect_list_append(NodesList, to_point_node);
909  }
910  }
911  else {
912  if (NodesList) {
913  /* it can happen that starting point falls on node but SP starts
914  * form the other node, add it in that case,
915  * similarly for to point below */
916  if (from_point_node > 0 && from_point_node != fnode[fn]) {
917  Vect_list_append(NodesList, from_point_node);
918  }
919 
920  /* add starting net SP search node */
921  Vect_list_append(NodesList, fnode[fn]);
922  }
923 
924  if (UseTtb)
925  Vect_net_ttb_shortest_path(Map, fnode[fn], 0, tnode[tn], 0,
926  tucfield, LList, NULL);
927  else
928  Vect_net_shortest_path(Map, fnode[fn], tnode[tn], LList, NULL);
929 
930  G_debug(3, "Number of lines %d", LList->n_values);
931 
932  if (Points)
933  Vect_append_points(Points, fPoints[fn], GV_FORWARD);
934 
935  if (FPoints)
936  Vect_append_points(FPoints, fPoints[fn], GV_FORWARD);
937 
938  for (i = 0; i < LList->n_values; i++) {
939  int line;
940 
941  line = LList->value[i];
942  G_debug(3, "i = %d line = %d", i, line);
943 
944  if (Points) {
945  Vect_read_line(Map, APoints, NULL, abs(line));
946 
947  if (line > 0)
948  Vect_append_points(Points, APoints, GV_FORWARD);
949  else
950  Vect_append_points(Points, APoints, GV_BACKWARD);
951  Points->n_points--;
952  }
953  if (NodesList) {
954  int node, node1, node2;
955 
956  Vect_get_line_nodes(Map, abs(line), &node1, &node2);
957  /* add the second node, the first of first segmet was
958  * already added */
959  if (line > 0)
960  node = node2;
961  else
962  node = node1;
963 
964  Vect_list_append(NodesList, node);
965  }
966 
967  if (List)
968  Vect_list_append(List, line);
969  }
970 
971  if (Points) {
972  if (LList->n_values)
973  Points->n_points++;
974  Vect_append_points(Points, tPoints[tn], GV_FORWARD);
975  }
976 
977  if (TPoints)
978  Vect_append_points(TPoints, tPoints[tn], GV_FORWARD);
979 
980  if (NodesList) {
981  if (to_point_node > 0 && to_point_node != tnode[tn]) {
982  Vect_list_append(NodesList, to_point_node);
983  }
984  }
985  }
986 
987  if (costs)
988  *costs = cur_cst;
989  if (Points)
990  Vect_line_prune(Points);
991  }
992 
993  return reachable;
994 }
995 
996 /*!
997  \brief Find shortest path on network between 2 points given by coordinates.
998 
999  \param Map vector map with build graph (see Vect_net_ttb_build_graph and
1000  Vect_net_build_graph) \param fx,fy,fz from point x coordinate (z ignored)
1001  \param tx,ty,tz to point x coordinate (z ignored)
1002  \param fmax maximum distance to the network from 'from'
1003  \param tmax maximum distance to the network from 'to'
1004  \param costs pointer where to store costs on the network (or NULL)
1005  \param Points pointer to the structure where to store vertices of shortest
1006  path (or NULL) \param List pointer to the structure where list of lines on
1007  the network is stored (or NULL) \param NodesList pointer to the structure
1008  where list of nodes on the network is stored (or NULL) \param FPoints pointer
1009  to the structure where to store line from 'from' to first network node (or
1010  NULL) \param TPoints pointer to the structure where to store line from last
1011  network node to 'to' (or NULL) \param fdist distance from 'from' to the net
1012  (or NULL) \param tdist distance from 'to' to the net (or NULL)
1013 
1014  \return 1 OK, 0 not reachable
1015  */
1016 int Vect_net_shortest_path_coor(struct Map_info *Map, double fx, double fy,
1017  double fz, double tx, double ty, double tz,
1018  double fmax, double tmax, double *costs,
1019  struct line_pnts *Points, struct ilist *List,
1020  struct ilist *NodesList,
1021  struct line_pnts *FPoints,
1022  struct line_pnts *TPoints, double *fdist,
1023  double *tdist)
1024 {
1025  return find_shortest_path_coor(Map, fx, fy, fz, tx, ty, tz, fmax, tmax, 0,
1026  0, costs, Points, List, NodesList, FPoints,
1027  TPoints, fdist, tdist);
1028 }
1029 
1030 /*!
1031  \brief Find shortest path on network with turntable between 2 points given by
1032  coordinates.
1033 
1034  \param Map vector map with build graph (see Vect_net_ttb_build_graph and
1035  Vect_net_build_graph) \param fx,fy,fz from point x coordinate (z ignored)
1036  \param tx,ty,tz to point x coordinate (z ignored)
1037  \param fmax maximum distance to the network from 'from'
1038  \param tmax maximum distance to the network from 'to'
1039  \param tucfield field with unique categories used in the turntable
1040  \param costs pointer where to store costs on the network (or NULL)
1041  \param Points pointer to the structure where to store vertices of shortest
1042  path (or NULL) \param List pointer to the structure where list of lines on
1043  the network is stored (or NULL) \param NodesList pointer to the structure
1044  where list of nodes on the network is stored (or NULL) \param FPoints pointer
1045  to the structure where to store line from 'from' to first network node (or
1046  NULL) \param TPoints pointer to the structure where to store line from last
1047  network node to 'to' (or NULL) \param fdist distance from 'from' to the net
1048  (or NULL) \param tdist distance from 'to' to the net (or NULL)
1049 
1050  \return 1 OK, 0 not reachable
1051  */
1052 int Vect_net_ttb_shortest_path_coor(struct Map_info *Map, double fx, double fy,
1053  double fz, double tx, double ty, double tz,
1054  double fmax, double tmax, int tucfield,
1055  double *costs, struct line_pnts *Points,
1056  struct ilist *List, struct ilist *NodesList,
1057  struct line_pnts *FPoints,
1058  struct line_pnts *TPoints, double *fdist,
1059  double *tdist)
1060 {
1061  return find_shortest_path_coor(Map, fx, fy, fz, tx, ty, tz, fmax, tmax, 1,
1062  tucfield, costs, Points, List, NodesList,
1063  FPoints, TPoints, fdist, tdist);
1064 }
#define NULL
Definition: ccmath.h:32
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:303
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:47
#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:1052
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:1016
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:499
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:706
int n_values
Number of values in the list.
Definition: gis.h:714
int * value
Array of values.
Definition: gis.h:710
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