GRASS Programmer's Manual  6.4.2(2012)
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
net.c
Go to the documentation of this file.
1 
20 #include <stdlib.h>
21 #include <string.h>
22 #include <fcntl.h>
23 #include <grass/gis.h>
24 #include <grass/dbmi.h>
25 #include <grass/Vect.h>
26 #include <grass/glocale.h>
27 
28 static int From_node; /* from node set in SP and used by clipper for first arc */
29 
30 static int clipper(dglGraph_s * pgraph,
31  dglSPClipInput_s * pargIn,
32  dglSPClipOutput_s * pargOut, void *pvarg)
33 { /* caller's pointer */
34  dglInt32_t cost;
35  dglInt32_t from;
36 
37  G_debug(3, "Net: clipper()");
38 
39  from = dglNodeGet_Id(pgraph, pargIn->pnNodeFrom);
40 
41  G_debug(3, " Edge = %d NodeFrom = %d NodeTo = %d edge cost = %d",
42  (int)dglEdgeGet_Id(pgraph, pargIn->pnEdge),
43  (int)from, (int)dglNodeGet_Id(pgraph, pargIn->pnNodeTo),
44  (int)pargOut->nEdgeCost);
45 
46  if (from != From_node) { /* do not clip first */
47  if (dglGet_NodeAttrSize(pgraph) > 0) {
48  memcpy(&cost, dglNodeGet_Attr(pgraph, pargIn->pnNodeFrom),
49  sizeof(cost));
50  if (cost == -1) { /* closed, cannot go from this node except it is 'from' node */
51  G_debug(3, " closed node");
52  return 1;
53  }
54  else {
55  G_debug(3, " EdgeCost += %d (node)", (int)cost);
56  pargOut->nEdgeCost += cost;
57  }
58  }
59  }
60  else {
61  G_debug(3, " don't clip first node");
62  }
63 
64  return 0;
65 }
66 
90 int
91 Vect_net_build_graph(struct Map_info *Map,
92  int ltype,
93  int afield,
94  int nfield,
95  const char *afcol,
96  const char *abcol,
97  const char *ncol, int geo, int algorithm)
98 {
99  int i, j, from, to, line, nlines, nnodes, ret, type, cat, skipped, cfound;
100  int dofw, dobw;
101  struct line_pnts *Points;
102  struct line_cats *Cats;
103  double dcost, bdcost, ll;
104  int cost, bcost;
105  dglGraph_s *gr;
106  dglInt32_t opaqueset[16] =
107  { 360000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
108  struct field_info *Fi;
109  dbDriver *driver = NULL;
110  dbHandle handle;
111  dbString stmt;
112  dbColumn *Column;
113  dbCatValArray fvarr, bvarr;
114  int fctype = 0, bctype = 0, nrec;
115 
116  /* TODO int costs -> double (waiting for dglib) */
117  G_debug(1, "Vect_build_graph(): ltype = %d, afield = %d, nfield = %d",
118  ltype, afield, nfield);
119  G_debug(1, " afcol = %s, abcol = %s, ncol = %s", afcol, abcol, ncol);
120 
121  G_message(_("Building graph..."));
122 
123  Map->graph_line_type = ltype;
124 
125  Points = Vect_new_line_struct();
126  Cats = Vect_new_cats_struct();
127 
128  ll = 0;
129  if (G_projection() == 3)
130  ll = 1; /* LL */
131 
132  if (afcol == NULL && ll && !geo)
133  Map->cost_multip = 1000000;
134  else
135  Map->cost_multip = 1000;
136 
137  nlines = Vect_get_num_lines(Map);
138  nnodes = Vect_get_num_nodes(Map);
139 
140  gr = &(Map->graph);
141 
142  /* Allocate space for costs, later replace by functions reading costs from graph */
143  Map->edge_fcosts = (double *)G_malloc((nlines + 1) * sizeof(double));
144  Map->edge_bcosts = (double *)G_malloc((nlines + 1) * sizeof(double));
145  Map->node_costs = (double *)G_malloc((nnodes + 1) * sizeof(double));
146  /* Set to -1 initially */
147  for (i = 1; i <= nlines; i++) {
148  Map->edge_fcosts[i] = -1; /* forward */
149  Map->edge_bcosts[i] = -1; /* backward */
150  }
151  for (i = 1; i <= nnodes; i++) {
152  Map->node_costs[i] = 0;
153  }
154 
155  if (ncol != NULL)
156  dglInitialize(gr, (dglByte_t) 1, sizeof(dglInt32_t), (dglInt32_t) 0,
157  opaqueset);
158  else
159  dglInitialize(gr, (dglByte_t) 1, (dglInt32_t) 0, (dglInt32_t) 0,
160  opaqueset);
161 
162  if (gr == NULL)
163  G_fatal_error(_("Unable to build network graph"));
164 
165  db_init_handle(&handle);
166  db_init_string(&stmt);
167 
168  if (abcol != NULL && afcol == NULL)
169  G_fatal_error(_("Forward costs column not specified"));
170 
171  /* --- Add arcs --- */
172  /* Open db connection */
173  if (afcol != NULL) {
174  /* Get field info */
175  if (afield < 1)
176  G_fatal_error(_("Arc field < 1"));
177  Fi = Vect_get_field(Map, afield);
178  if (Fi == NULL)
179  G_fatal_error(_("Database connection not defined for layer %d"),
180  afield);
181 
182  /* Open database */
183  driver = db_start_driver_open_database(Fi->driver, Fi->database);
184  if (driver == NULL)
185  G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
186  Fi->database, Fi->driver);
187 
188  /* Load costs to array */
189  if (db_get_column(driver, Fi->table, afcol, &Column) != DB_OK)
190  G_fatal_error(_("Column <%s> not found in table <%s>"),
191  afcol, Fi->table);
192 
193  fctype = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
194 
195  if (fctype != DB_C_TYPE_INT && fctype != DB_C_TYPE_DOUBLE)
196  G_fatal_error(_("Data type of column <%s> not supported (must be numeric)"),
197  afcol);
198 
199  db_CatValArray_init(&fvarr);
200  nrec =
201  db_select_CatValArray(driver, Fi->table, Fi->key, afcol, NULL,
202  &fvarr);
203  G_debug(1, "forward costs: nrec = %d", nrec);
204 
205  if (abcol != NULL) {
206  if (db_get_column(driver, Fi->table, abcol, &Column) != DB_OK)
207  G_fatal_error(_("Column <%s> not found in table <%s>"),
208  abcol, Fi->table);
209 
210  bctype = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
211 
212  if (bctype != DB_C_TYPE_INT && bctype != DB_C_TYPE_DOUBLE)
213  G_fatal_error(_("Data type of column <%s> not supported (must be numeric)"),
214  abcol);
215 
216  db_CatValArray_init(&bvarr);
217  nrec =
218  db_select_CatValArray(driver, Fi->table, Fi->key, abcol, NULL,
219  &bvarr);
220  G_debug(1, "backward costs: nrec = %d", nrec);
221  }
222  }
223 
224  skipped = 0;
225 
226  G_message(_("Registering arcs..."));
227 
228  for (i = 1; i <= nlines; i++) {
229  G_percent(i, nlines, 1); /* must be before any continue */
230  dofw = dobw = 1;
231  Vect_get_line_nodes(Map, i, &from, &to);
232  type = Vect_read_line(Map, Points, Cats, i);
233  if (!(type & ltype & (GV_LINE | GV_BOUNDARY)))
234  continue;
235 
236  if (afcol != NULL) {
237  if (!(Vect_cat_get(Cats, afield, &cat))) {
238  G_debug(2,
239  "Category of field %d not attached to the line %d -> line skipped",
240  afield, i);
241  skipped += 2; /* Both directions */
242  continue;
243  }
244  else {
245  if (fctype == DB_C_TYPE_INT) {
246  ret = db_CatValArray_get_value_int(&fvarr, cat, &cost);
247  dcost = cost;
248  }
249  else { /* DB_C_TYPE_DOUBLE */
250  ret =
251  db_CatValArray_get_value_double(&fvarr, cat, &dcost);
252  }
253  if (ret != DB_OK) {
254  G_warning(_("Database record for line %d (cat = %d, "
255  "forward/both direction(s)) not found "
256  "(forward/both direction(s) of line skipped)"),
257  i, cat);
258  dofw = 0;
259  }
260 
261  if (abcol != NULL) {
262  if (bctype == DB_C_TYPE_INT) {
263  ret =
264  db_CatValArray_get_value_int(&bvarr, cat, &bcost);
265  bdcost = bcost;
266  }
267  else { /* DB_C_TYPE_DOUBLE */
268  ret =
270  &bdcost);
271  }
272  if (ret != DB_OK) {
273  G_warning(_("Database record for line %d (cat = %d, "
274  "backword direction) not found"
275  "(direction of line skipped)"), i, cat);
276  dobw = 0;
277  }
278  }
279  else {
280  if (dofw)
281  bdcost = dcost;
282  else
283  dobw = 0;
284  }
285  }
286  }
287  else {
288  if (ll) {
289  if (geo)
290  dcost = Vect_line_geodesic_length(Points);
291  else
292  dcost = Vect_line_length(Points);
293  }
294  else
295  dcost = Vect_line_length(Points);
296 
297  bdcost = dcost;
298  }
299  if (dofw && dcost != -1) {
300  cost = (dglInt32_t) Map->cost_multip * dcost;
301  G_debug(5, "Add arc %d from %d to %d cost = %d", i, from, to,
302  cost);
303  ret =
304  dglAddEdge(gr, (dglInt32_t) from, (dglInt32_t) to,
305  (dglInt32_t) cost, (dglInt32_t) i);
306  Map->edge_fcosts[i] = dcost;
307  if (ret < 0)
308  G_fatal_error("Cannot add network arc");
309  }
310 
311  G_debug(5, "bdcost = %f edge_bcosts = %f", bdcost,
312  Map->edge_bcosts[i]);
313  if (dobw && bdcost != -1) {
314  bcost = (dglInt32_t) Map->cost_multip * bdcost;
315  G_debug(5, "Add arc %d from %d to %d bcost = %d", -i, to, from,
316  bcost);
317  ret =
318  dglAddEdge(gr, (dglInt32_t) to, (dglInt32_t) from,
319  (dglInt32_t) bcost, (dglInt32_t) - i);
320  Map->edge_bcosts[i] = bdcost;
321  if (ret < 0)
322  G_fatal_error(_("Cannot add network arc"));
323  }
324  }
325 
326  if (afcol != NULL && skipped > 0)
327  G_debug(2, "%d lines missing category of field %d skipped", skipped,
328  afield);
329 
330  if (afcol != NULL) {
332  db_CatValArray_free(&fvarr);
333 
334  if (abcol != NULL) {
335  db_CatValArray_free(&bvarr);
336  }
337  }
338 
339  /* Set node attributes */
340  G_debug(2, "Register nodes");
341  if (ncol != NULL) {
342  G_debug(2, "Set nodes' costs");
343  if (nfield < 1)
344  G_fatal_error("Node field < 1");
345 
346  G_message(_("Setting node costs..."));
347 
348  Fi = Vect_get_field(Map, nfield);
349  if (Fi == NULL)
350  G_fatal_error(_("Database connection not defined for layer %d"),
351  nfield);
352 
353  driver = db_start_driver_open_database(Fi->driver, Fi->database);
354  if (driver == NULL)
355  G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
356  Fi->database, Fi->driver);
357 
358  /* Load costs to array */
359  if (db_get_column(driver, Fi->table, ncol, &Column) != DB_OK)
360  G_fatal_error(_("Column <%s> not found in table <%s>"),
361  ncol, Fi->table);
362 
363  fctype = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
364 
365  if (fctype != DB_C_TYPE_INT && fctype != DB_C_TYPE_DOUBLE)
366  G_fatal_error(_("Data type of column <%s> not supported (must be numeric)"),
367  ncol);
368 
369  db_CatValArray_init(&fvarr);
370  nrec =
371  db_select_CatValArray(driver, Fi->table, Fi->key, ncol, NULL,
372  &fvarr);
373  G_debug(1, "node costs: nrec = %d", nrec);
374 
375  for (i = 1; i <= nnodes; i++) {
376  /* TODO: what happens if we set attributes of not existing node (skipped lines,
377  * nodes without lines) */
378 
379  nlines = Vect_get_node_n_lines(Map, i);
380  G_debug(2, " node = %d nlines = %d", i, nlines);
381  cfound = 0;
382  dcost = 0;
383  for (j = 0; j < nlines; j++) {
384  line = Vect_get_node_line(Map, i, j);
385  G_debug(2, " line (%d) = %d", j, line);
386  type = Vect_read_line(Map, NULL, Cats, line);
387  if (!(type & GV_POINT))
388  continue;
389  if (Vect_cat_get(Cats, nfield, &cat)) { /* point with category of field found */
390  /* Set costs */
391  if (fctype == DB_C_TYPE_INT) {
392  ret =
393  db_CatValArray_get_value_int(&fvarr, cat, &cost);
394  dcost = cost;
395  }
396  else { /* DB_C_TYPE_DOUBLE */
397  ret =
399  &dcost);
400  }
401  if (ret != DB_OK) {
402  G_warning(_("Database record for node %d (cat = %d) not found "
403  "(cost set to 0)"), i, cat);
404  }
405  cfound = 1;
406  break;
407  }
408  }
409  if (!cfound) {
410  G_debug(2,
411  "Category of field %d not attached to any points in node %d"
412  "(costs set to 0)", nfield, i);
413  }
414  if (dcost == -1) { /* closed */
415  cost = -1;
416  }
417  else {
418  cost = (dglInt32_t) Map->cost_multip * dcost;
419  }
420  G_debug(3, "Set node's cost to %d", cost);
422  (dglInt32_t *) (dglInt32_t) & cost);
423  Map->node_costs[i] = dcost;
424  }
426  db_CatValArray_free(&fvarr);
427  }
428 
429  G_message(_("Flattening the graph..."));
430  ret = dglFlatten(gr);
431  if (ret < 0)
432  G_fatal_error(_("GngFlatten error"));
433 
434  /* init SP cache */
435  /* disable to debug dglib cache */
436  dglInitializeSPCache(gr, &(Map->spCache));
437 
438  G_message(_("Graph was built"));
439 
440  return 0;
441 }
442 
443 
462 int
463 Vect_net_shortest_path(struct Map_info *Map, int from, int to,
464  struct ilist *List, double *cost)
465 {
466  int i, line, *pclip, cArc, nRet;
467  dglSPReport_s *pSPReport;
468  dglInt32_t nDistance;
469  int use_cache = 1; /* set to 0 to disable dglib cache */
470 
471  G_debug(3, "Vect_net_shortest_path(): from = %d, to = %d", from, to);
472 
473  /* Note : if from == to dgl goes to nearest node and returns back (dgl feature) =>
474  * check here for from == to */
475 
476  if (List != NULL)
477  Vect_reset_list(List);
478 
479  /* Check if from and to are identical, otherwise dglib returns path to neares node and back! */
480  if (from == to) {
481  if (cost != NULL)
482  *cost = 0;
483  return 0;
484  }
485 
486  From_node = from;
487 
488  pclip = NULL;
489  if (List != NULL) {
490  if (use_cache) {
491  nRet =
492  dglShortestPath(&(Map->graph), &pSPReport, (dglInt32_t) from,
493  (dglInt32_t) to, clipper, pclip, &(Map->spCache));
494  }
495  else {
496  nRet =
497  dglShortestPath(&(Map->graph), &pSPReport, (dglInt32_t) from,
498  (dglInt32_t) to, clipper, pclip, NULL);
499  }
500  }
501  else {
502  if (use_cache) {
503  nRet =
504  dglShortestDistance(&(Map->graph), &nDistance, (dglInt32_t) from,
505  (dglInt32_t) to, clipper, pclip, &(Map->spCache));
506  }
507  else {
508  nRet =
509  dglShortestDistance(&(Map->graph), &nDistance, (dglInt32_t) from,
510  (dglInt32_t) to, clipper, pclip, NULL);
511  }
512  }
513 
514  if (nRet == 0) {
515  /* G_warning( "Destination node %d is unreachable from node %d\n" , to , from ); */
516  if (cost != NULL)
517  *cost = PORT_DOUBLE_MAX;
518  return -1;
519  }
520  else if (nRet < 0) {
521  G_warning(_("dglShortestPath error: %s"), dglStrerror(&(Map->graph)));
522  return -1;
523  }
524 
525  if (List != NULL) {
526  for (i = 0; i < pSPReport->cArc; i++) {
527  line = dglEdgeGet_Id(&(Map->graph), pSPReport->pArc[i].pnEdge);
528  G_debug(2, "From %ld to %ld - cost %ld user %d distance %ld", pSPReport->pArc[i].nFrom, pSPReport->pArc[i].nTo, dglEdgeGet_Cost(&(Map->graph), pSPReport->pArc[i].pnEdge) / Map->cost_multip, /* this is the cost from clip() */
529  line, pSPReport->pArc[i].nDistance);
530  Vect_list_append(List, line);
531  }
532  }
533 
534  if (cost != NULL) {
535  if (List != NULL)
536  *cost = (double)pSPReport->nDistance / Map->cost_multip;
537  else
538  *cost = (double)nDistance / Map->cost_multip;
539  }
540 
541  if (List != NULL) {
542  cArc = pSPReport->cArc;
543  dglFreeSPReport(&(Map->graph), pSPReport);
544  }
545  else
546  cArc = 0;
547 
548  return (cArc);
549 }
550 
563 int
564 Vect_net_get_line_cost(struct Map_info *Map, int line, int direction,
565  double *cost)
566 {
567  /* dglInt32_t *pEdge; */
568 
569  G_debug(5, "Vect_net_get_line_cost(): line = %d, dir = %d", line,
570  direction);
571 
572  if (direction == GV_FORWARD) {
573  /* V1 has no index by line-id -> array used */
574  /*
575  pEdge = dglGetEdge ( &(Map->graph), line );
576  if ( pEdge == NULL ) return 0;
577  *cost = (double) dglEdgeGet_Cost ( &(Map->graph), pEdge );
578  */
579  if (Map->edge_fcosts[line] == -1) {
580  *cost = -1;
581  return 0;
582  }
583  else
584  *cost = Map->edge_fcosts[line];
585  }
586  else if (direction == GV_BACKWARD) {
587  /*
588  pEdge = dglGetEdge ( &(Map->graph), -line );
589  if ( pEdge == NULL ) return 0;
590  *cost = (double) dglEdgeGet_Cost ( &(Map->graph), pEdge );
591  */
592  if (Map->edge_bcosts[line] == -1) {
593  *cost = -1;
594  return 0;
595  }
596  else
597  *cost = Map->edge_bcosts[line];
598  G_debug(5, "Vect_net_get_line_cost(): edge_bcosts = %f",
599  Map->edge_bcosts[line]);
600  }
601  else {
602  G_fatal_error(_("Wrong line direction in Vect_net_get_line_cost()"));
603  }
604 
605  return 1;
606 }
607 
617 int Vect_net_get_node_cost(struct Map_info *Map, int node, double *cost)
618 {
619  G_debug(3, "Vect_net_get_node_cost(): node = %d", node);
620 
621  *cost = Map->node_costs[node];
622 
623  G_debug(3, " -> cost = %f", *cost);
624 
625  return 1;
626 }
627 
646 int Vect_net_nearest_nodes(struct Map_info *Map,
647  double x, double y, double z,
648  int direction, double maxdist,
649  int *node1, int *node2, int *ln, double *costs1,
650  double *costs2, struct line_pnts *Points1,
651  struct line_pnts *Points2, double *distance)
652 {
653  int line, n1, n2, nnodes;
654  int npoints;
655  int segment; /* nearest line segment (first is 1) */
656  static struct line_pnts *Points = NULL;
657  double cx, cy, cz, c1, c2;
658  double along; /* distance along the line to nearest point */
659  double length;
660 
661  G_debug(3, "Vect_net_nearest_nodes() x = %f y = %f", x, y);
662 
663  /* Reset */
664  if (node1)
665  *node1 = 0;
666  if (node2)
667  *node2 = 0;
668  if (ln)
669  *ln = 0;
670  if (costs1)
671  *costs1 = PORT_DOUBLE_MAX;
672  if (costs2)
673  *costs2 = PORT_DOUBLE_MAX;
674  if (Points1)
675  Vect_reset_line(Points1);
676  if (Points2)
677  Vect_reset_line(Points2);
678  if (distance)
679  *distance = PORT_DOUBLE_MAX;
680 
681  if (!Points)
682  Points = Vect_new_line_struct();
683 
684  /* Find nearest line */
685  line = Vect_find_line(Map, x, y, z, Map->graph_line_type, maxdist, 0, 0);
686 
687  if (line < 1)
688  return 0;
689 
690  Vect_read_line(Map, Points, NULL, line);
691  npoints = Points->n_points;
692  Vect_get_line_nodes(Map, line, &n1, &n2);
693 
694  segment =
695  Vect_line_distance(Points, x, y, z, 0, &cx, &cy, &cz, distance, NULL,
696  &along);
697 
698  G_debug(4, "line = %d n1 = %d n2 = %d segment = %d", line, n1, n2,
699  segment);
700 
701  /* Check first or last point and return one node in that case */
702  G_debug(4, "cx = %f cy = %f first = %f %f last = %f %f", cx, cy,
703  Points->x[0], Points->y[0], Points->x[npoints - 1],
704  Points->y[npoints - 1]);
705 
706  if (Points->x[0] == cx && Points->y[0] == cy) {
707  if (node1)
708  *node1 = n1;
709  if (ln)
710  *ln = line;
711  if (costs1)
712  *costs1 = 0;
713  if (Points1) {
714  Vect_append_point(Points1, x, y, z);
715  Vect_append_point(Points1, cx, cy, cz);
716  }
717  G_debug(3, "first node nearest");
718  return 1;
719  }
720  if (Points->x[npoints - 1] == cx && Points->y[npoints - 1] == cy) {
721  if (node1)
722  *node1 = n2;
723  if (ln)
724  *ln = line;
725  if (costs1)
726  *costs1 = 0;
727  if (Points1) {
728  Vect_append_point(Points1, x, y, z);
729  Vect_append_point(Points1, cx, cy, cz);
730  }
731  G_debug(3, "last node nearest");
732  return 1;
733  }
734 
735  nnodes = 2;
736 
737  /* c1 - costs to get from/to the first vertex */
738  /* c2 - costs to get from/to the last vertex */
739  if (direction == GV_FORWARD) { /* from point to net */
740  Vect_net_get_line_cost(Map, line, GV_BACKWARD, &c1);
741  Vect_net_get_line_cost(Map, line, GV_FORWARD, &c2);
742  }
743  else {
744  Vect_net_get_line_cost(Map, line, GV_FORWARD, &c1);
745  Vect_net_get_line_cost(Map, line, GV_BACKWARD, &c2);
746  }
747 
748  if (c1 < 0)
749  nnodes--;
750  if (c2 < 0)
751  nnodes--;
752  if (nnodes == 0)
753  return 0; /* both directions closed */
754 
755  length = Vect_line_length(Points);
756 
757  if (ln)
758  *ln = line;
759 
760  if (nnodes == 1 && c1 < 0) { /* first direction is closed, return node2 as node1 */
761  if (node1)
762  *node1 = n2;
763 
764  if (costs1) { /* to node 2, i.e. forward */
765  *costs1 = c2 * (length - along) / length;
766  }
767 
768  if (Points1) { /* to node 2, i.e. forward */
769  int i;
770 
771  if (direction == GV_FORWARD) { /* from point to net */
772  Vect_append_point(Points1, x, y, z);
773  Vect_append_point(Points1, cx, cy, cz);
774  for (i = segment; i < npoints; i++)
775  Vect_append_point(Points1, Points->x[i], Points->y[i],
776  Points->z[i]);
777  }
778  else {
779  for (i = npoints - 1; i >= segment; i--)
780  Vect_append_point(Points1, Points->x[i], Points->y[i],
781  Points->z[i]);
782 
783  Vect_append_point(Points1, cx, cy, cz);
784  Vect_append_point(Points1, x, y, z);
785  }
786  }
787  }
788  else {
789  if (node1)
790  *node1 = n1;
791  if (node2)
792  *node2 = n2;
793 
794  if (costs1) { /* to node 1, i.e. backward */
795  *costs1 = c1 * along / length;
796  }
797 
798  if (costs2) { /* to node 2, i.e. forward */
799  *costs2 = c2 * (length - along) / length;
800  }
801 
802  if (Points1) { /* to node 1, i.e. backward */
803  int i;
804 
805  if (direction == GV_FORWARD) { /* from point to net */
806  Vect_append_point(Points1, x, y, z);
807  Vect_append_point(Points1, cx, cy, cz);
808  for (i = segment - 1; i >= 0; i--)
809  Vect_append_point(Points1, Points->x[i], Points->y[i],
810  Points->z[i]);
811  }
812  else {
813  for (i = 0; i < segment; i++)
814  Vect_append_point(Points1, Points->x[i], Points->y[i],
815  Points->z[i]);
816 
817  Vect_append_point(Points1, cx, cy, cz);
818  Vect_append_point(Points1, x, y, z);
819  }
820  }
821 
822  if (Points2) { /* to node 2, i.e. forward */
823  int i;
824 
825  if (direction == GV_FORWARD) { /* from point to net */
826  Vect_append_point(Points2, x, y, z);
827  Vect_append_point(Points2, cx, cy, cz);
828  for (i = segment; i < npoints; i++)
829  Vect_append_point(Points2, Points->x[i], Points->y[i],
830  Points->z[i]);
831  }
832  else {
833  for (i = npoints - 1; i >= segment; i--)
834  Vect_append_point(Points2, Points->x[i], Points->y[i],
835  Points->z[i]);
836 
837  Vect_append_point(Points2, cx, cy, cz);
838  Vect_append_point(Points2, x, y, z);
839  }
840  }
841  }
842 
843  return nnodes;
844 }
845 
864 int
865 Vect_net_shortest_path_coor(struct Map_info *Map,
866  double fx, double fy, double fz, double tx,
867  double ty, double tz, double fmax, double tmax,
868  double *costs, struct line_pnts *Points,
869  struct ilist *List, struct line_pnts *FPoints,
870  struct line_pnts *TPoints, double *fdist,
871  double *tdist)
872 {
873  return Vect_net_shortest_path_coor2 ( Map, fx, fy, fz, tx, ty, tz, fmax, tmax,
874  costs, Points, List, NULL, FPoints, TPoints, fdist, tdist );
875 }
876 
896 int
897 Vect_net_shortest_path_coor2(struct Map_info *Map,
898  double fx, double fy, double fz, double tx,
899  double ty, double tz, double fmax, double tmax,
900  double *costs, struct line_pnts *Points,
901  struct ilist *List, struct ilist *NodesList,
902  struct line_pnts *FPoints,
903  struct line_pnts *TPoints, double *fdist,
904  double *tdist)
905 {
906  int fnode[2], tnode[2]; /* nearest nodes, *node[1] is 0 if only one was found */
907  double fcosts[2], tcosts[2], cur_cst; /* costs to nearest nodes on the network */
908  int nfnodes, ntnodes, fline, tline;
909  static struct line_pnts *APoints, *SPoints, *fPoints[2], *tPoints[2];
910  static struct ilist *LList;
911  static int first = 1;
912  int reachable, shortcut;
913  int i, j, fn = 0, tn = 0;
914 
915  /* from/to_point_node is set if from/to point projected to line
916  *falls exactly on node (shortcut -> fline == tline) */
917  int from_point_node = 0;
918  int to_point_node = 0;
919 
920  G_debug(3, "Vect_net_shortest_path_coor()");
921 
922  if (first) {
923  APoints = Vect_new_line_struct();
924  SPoints = Vect_new_line_struct();
925  fPoints[0] = Vect_new_line_struct();
926  fPoints[1] = Vect_new_line_struct();
927  tPoints[0] = Vect_new_line_struct();
928  tPoints[1] = Vect_new_line_struct();
929  LList = Vect_new_list();
930  first = 0;
931  }
932 
933  /* Reset */
934  if (costs)
935  *costs = PORT_DOUBLE_MAX;
936  if (Points)
937  Vect_reset_line(Points);
938  if (fdist)
939  *fdist = 0;
940  if (tdist)
941  *tdist = 0;
942  if (List)
943  List->n_values = 0;
944  if (FPoints)
945  Vect_reset_line(FPoints);
946  if (TPoints)
947  Vect_reset_line(TPoints);
948  if (NodesList != NULL)
949  Vect_reset_list(NodesList);
950 
951  /* Find nearest nodes */
952  fnode[0] = fnode[1] = tnode[0] = tnode[1] = 0;
953 
954  nfnodes =
955  Vect_net_nearest_nodes(Map, fx, fy, fz, GV_FORWARD, fmax, &(fnode[0]),
956  &(fnode[1]), &fline, &(fcosts[0]),
957  &(fcosts[1]), fPoints[0], fPoints[1], fdist);
958  if (nfnodes == 0)
959  return 0;
960 
961  if ( nfnodes == 1 && fPoints[0]->n_points < 3 ) {
962  from_point_node = fnode[0];
963  }
964 
965  ntnodes =
966  Vect_net_nearest_nodes(Map, tx, ty, tz, GV_BACKWARD, tmax,
967  &(tnode[0]), &(tnode[1]), &tline, &(tcosts[0]),
968  &(tcosts[1]), tPoints[0], tPoints[1], tdist);
969  if (ntnodes == 0)
970  return 0;
971 
972  if ( ntnodes == 1 && tPoints[0]->n_points < 3 ) {
973  to_point_node = tnode[0];
974  }
975 
976 
977  G_debug(3, "fline = %d tline = %d", fline, tline);
978 
979  reachable = shortcut = 0;
980  cur_cst = PORT_DOUBLE_MAX;
981 
982  /* It may happen, that 2 points are at the same line. */
983  /* TODO?: it could also happen that fline != tline but both points are on the same
984  * line if they fall on node but a different line was found. This case is correctly
985  * handled as normal non shortcut, but it could be added here. In that case
986  * NodesList collection must be changed */
987  if (fline == tline && (nfnodes > 1 || ntnodes > 1)) {
988  double len, flen, tlen, c, fseg, tseg;
989  double fcx, fcy, fcz, tcx, tcy, tcz;
990 
991  Vect_read_line(Map, APoints, NULL, fline);
992  len = Vect_line_length(APoints);
993 
994  /* distance along the line */
995  fseg =
996  Vect_line_distance(APoints, fx, fy, fz, 0, &fcx, &fcy, &fcz, NULL,
997  NULL, &flen);
998  tseg =
999  Vect_line_distance(APoints, tx, ty, tz, 0, &tcx, &tcy, &tcz, NULL,
1000  NULL, &tlen);
1001 
1002  Vect_reset_line(SPoints);
1003  if (flen == tlen) {
1004  cur_cst = 0;
1005  reachable = shortcut = 1;
1006  }
1007  else if (flen < tlen) {
1008  Vect_net_get_line_cost(Map, fline, GV_FORWARD, &c);
1009  if (c >= 0) {
1010  cur_cst = c * (tlen - flen) / len;
1011 
1012  Vect_append_point(SPoints, fx, fy, fz);
1013  Vect_append_point(SPoints, fcx, fcy, fcz);
1014  for (i = fseg; i < tseg; i++)
1015  Vect_append_point(SPoints, APoints->x[i], APoints->y[i],
1016  APoints->z[i]);
1017 
1018  Vect_append_point(SPoints, tcx, tcy, tcz);
1019  Vect_append_point(SPoints, tx, ty, tz);
1020 
1021  reachable = shortcut = 1;
1022  }
1023  }
1024  else { /* flen > tlen */
1025  Vect_net_get_line_cost(Map, fline, GV_BACKWARD, &c);
1026  if (c >= 0) {
1027  cur_cst = c * (flen - tlen) / len;
1028 
1029  Vect_append_point(SPoints, fx, fy, fz);
1030  Vect_append_point(SPoints, fcx, fcy, fcz);
1031  for (i = fseg - 1; i >= tseg; i--)
1032  Vect_append_point(SPoints, APoints->x[i], APoints->y[i],
1033  APoints->z[i]);
1034 
1035  Vect_append_point(SPoints, tcx, tcy, tcz);
1036  Vect_append_point(SPoints, tx, ty, tz);
1037 
1038  reachable = shortcut = 1;
1039  }
1040  }
1041  }
1042 
1043  /* Find the shortest variant from maximum 4 */
1044  for (i = 0; i < nfnodes; i++) {
1045  for (j = 0; j < ntnodes; j++) {
1046  double ncst, cst;
1047  int ret;
1048 
1049  G_debug(3, "i = %d fnode = %d j = %d tnode = %d", i, fnode[i], j,
1050  tnode[j]);
1051 
1052  ret =
1053  Vect_net_shortest_path(Map, fnode[i], tnode[j], NULL, &ncst);
1054  if (ret == -1)
1055  continue; /* not reachable */
1056 
1057  cst = fcosts[i] + ncst + tcosts[j];
1058  if (reachable == 0 || cst < cur_cst) {
1059  cur_cst = cst;
1060  fn = i;
1061  tn = j;
1062  shortcut = 0;
1063  }
1064  reachable = 1;
1065  }
1066  }
1067 
1068  G_debug(3, "reachable = %d shortcut = %d cur_cst = %f", reachable,
1069  shortcut, cur_cst);
1070  if (reachable) {
1071  int ret;
1072 
1073  if (shortcut) {
1074  if (Points)
1075  Vect_append_points(Points, SPoints, GV_FORWARD);
1076  if (NodesList) {
1077  /* Check if from/to point projected to line falls on node and
1078  *add it to the list */
1079  if (from_point_node > 0)
1080  Vect_list_append(NodesList, from_point_node);
1081 
1082  if (to_point_node > 0)
1083  Vect_list_append(NodesList, to_point_node);
1084  }
1085  }
1086  else {
1087  if (NodesList) {
1088  /* it can happen that starting point falls on node but SP starts
1089  * form the other node, add it in that case,
1090  * similarly for to point below */
1091  if (from_point_node > 0 && from_point_node != fnode[fn]) {
1092  Vect_list_append(NodesList, from_point_node);
1093  }
1094 
1095  /* add starting net SP search node */
1096  Vect_list_append(NodesList, fnode[fn]);
1097  }
1098  ret =
1099  Vect_net_shortest_path(Map, fnode[fn], tnode[tn], LList,
1100  NULL);
1101  G_debug(3, "Number of lines %d", LList->n_values);
1102 
1103  if (Points)
1104  Vect_append_points(Points, fPoints[fn], GV_FORWARD);
1105 
1106  if (FPoints)
1107  Vect_append_points(FPoints, fPoints[fn], GV_FORWARD);
1108 
1109  for (i = 0; i < LList->n_values; i++) {
1110  int line;
1111 
1112  line = LList->value[i];
1113  G_debug(3, "i = %d line = %d", i, line);
1114 
1115  if (Points) {
1116  Vect_read_line(Map, APoints, NULL, abs(line));
1117 
1118  if (line > 0)
1119  Vect_append_points(Points, APoints, GV_FORWARD);
1120  else
1121  Vect_append_points(Points, APoints, GV_BACKWARD);
1122  }
1123  if (NodesList) {
1124  int node, node1, node2;
1125  Vect_get_line_nodes(Map, abs(line), &node1, &node2);
1126  /* add the second node, the first of first segmet was alread added */
1127  if (line > 0)
1128  node = node2;
1129  else
1130  node = node1;
1131 
1132  Vect_list_append(NodesList, node);
1133  }
1134 
1135  if (List)
1136  Vect_list_append(List, line);
1137  }
1138 
1139  if (Points)
1140  Vect_append_points(Points, tPoints[tn], GV_FORWARD);
1141 
1142  if (TPoints)
1143  Vect_append_points(TPoints, tPoints[tn], GV_FORWARD);
1144 
1145  if (NodesList) {
1146  if (to_point_node > 0 && to_point_node != tnode[tn]) {
1147  Vect_list_append(NodesList, to_point_node);
1148  }
1149  }
1150  }
1151 
1152  if (costs)
1153  *costs = cur_cst;
1154  }
1155 
1156  return reachable;
1157 }