00001
00020 #include<stdlib.h>
00021 #include<string.h>
00022 #include<fcntl.h>
00023 #include<grass/gis.h>
00024 #include<grass/dbmi.h>
00025 #include<grass/Vect.h>
00026 #include<grass/glocale.h>
00027
00028 static int From_node;
00029
00030 static int clipper(dglGraph_s * pgraph,
00031 dglSPClipInput_s * pargIn,
00032 dglSPClipOutput_s * pargOut, void *pvarg)
00033 {
00034 dglInt32_t cost;
00035 dglInt32_t from;
00036
00037 G_debug(3, "Net: clipper()");
00038
00039 from = dglNodeGet_Id(pgraph, pargIn->pnNodeFrom);
00040
00041 G_debug(3, " Edge = %d NodeFrom = %d NodeTo = %d edge cost = %d",
00042 (int)dglEdgeGet_Id(pgraph, pargIn->pnEdge),
00043 (int)from, (int)dglNodeGet_Id(pgraph, pargIn->pnNodeTo),
00044 (int)pargOut->nEdgeCost);
00045
00046 if (from != From_node) {
00047 if (dglGet_NodeAttrSize(pgraph) > 0) {
00048 memcpy(&cost, dglNodeGet_Attr(pgraph, pargIn->pnNodeFrom),
00049 sizeof(cost));
00050 if (cost == -1) {
00051 G_debug(3, " closed node");
00052 return 1;
00053 }
00054 else {
00055 G_debug(3, " EdgeCost += %d (node)", (int)cost);
00056 pargOut->nEdgeCost += cost;
00057 }
00058 }
00059 }
00060 else {
00061 G_debug(3, " don't clip first node");
00062 }
00063
00064 return 0;
00065 }
00066
00090 int
00091 Vect_net_build_graph(struct Map_info *Map,
00092 int ltype,
00093 int afield,
00094 int nfield,
00095 const char *afcol,
00096 const char *abcol,
00097 const char *ncol, int geo, int algorithm)
00098 {
00099 int i, j, from, to, line, nlines, nnodes, ret, type, cat, skipped, cfound;
00100 int dofw, dobw;
00101 struct line_pnts *Points;
00102 struct line_cats *Cats;
00103 double dcost, bdcost, ll;
00104 int cost, bcost;
00105 dglGraph_s *gr;
00106 dglInt32_t opaqueset[16] =
00107 { 360000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
00108 struct field_info *Fi;
00109 dbDriver *driver = NULL;
00110 dbHandle handle;
00111 dbString stmt;
00112 dbColumn *Column;
00113 dbCatValArray fvarr, bvarr;
00114 int fctype = 0, bctype = 0, nrec;
00115
00116
00117 G_debug(1, "Vect_build_graph(): ltype = %d, afield = %d, nfield = %d",
00118 ltype, afield, nfield);
00119 G_debug(1, " afcol = %s, abcol = %s, ncol = %s", afcol, abcol, ncol);
00120
00121 G_message(_("Building graph..."));
00122
00123 Map->graph_line_type = ltype;
00124
00125 Points = Vect_new_line_struct();
00126 Cats = Vect_new_cats_struct();
00127
00128 ll = 0;
00129 if (G_projection() == 3)
00130 ll = 1;
00131
00132 if (afcol == NULL && ll && !geo)
00133 Map->cost_multip = 1000000;
00134 else
00135 Map->cost_multip = 1000;
00136
00137 nlines = Vect_get_num_lines(Map);
00138 nnodes = Vect_get_num_nodes(Map);
00139
00140 gr = &(Map->graph);
00141
00142
00143 Map->edge_fcosts = (double *)G_malloc((nlines + 1) * sizeof(double));
00144 Map->edge_bcosts = (double *)G_malloc((nlines + 1) * sizeof(double));
00145 Map->node_costs = (double *)G_malloc((nnodes + 1) * sizeof(double));
00146
00147 for (i = 1; i <= nlines; i++) {
00148 Map->edge_fcosts[i] = -1;
00149 Map->edge_bcosts[i] = -1;
00150 }
00151 for (i = 1; i <= nnodes; i++) {
00152 Map->node_costs[i] = 0;
00153 }
00154
00155 if (ncol != NULL)
00156 dglInitialize(gr, (dglByte_t) 1, sizeof(dglInt32_t), (dglInt32_t) 0,
00157 opaqueset);
00158 else
00159 dglInitialize(gr, (dglByte_t) 1, (dglInt32_t) 0, (dglInt32_t) 0,
00160 opaqueset);
00161
00162 if (gr == NULL)
00163 G_fatal_error(_("Unable to build network graph"));
00164
00165 db_init_handle(&handle);
00166 db_init_string(&stmt);
00167
00168 if (abcol != NULL && afcol == NULL)
00169 G_fatal_error(_("Forward costs column not specified"));
00170
00171
00172
00173 if (afcol != NULL) {
00174
00175 if (afield < 1)
00176 G_fatal_error(_("Arc field < 1"));
00177 Fi = Vect_get_field(Map, afield);
00178 if (Fi == NULL)
00179 G_fatal_error(_("Database connection not defined for layer %d"),
00180 afield);
00181
00182
00183 driver = db_start_driver_open_database(Fi->driver, Fi->database);
00184 if (driver == NULL)
00185 G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
00186 Fi->database, Fi->driver);
00187
00188
00189 if (db_get_column(driver, Fi->table, afcol, &Column) != DB_OK)
00190 G_fatal_error(_("Column <%s> not found in table <%s>"),
00191 afcol, Fi->table);
00192
00193 fctype = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
00194
00195 if (fctype != DB_C_TYPE_INT && fctype != DB_C_TYPE_DOUBLE)
00196 G_fatal_error(_("Data type of column <%s> not supported (must be numeric)"),
00197 afcol);
00198
00199 db_CatValArray_init(&fvarr);
00200 nrec =
00201 db_select_CatValArray(driver, Fi->table, Fi->key, afcol, NULL,
00202 &fvarr);
00203 G_debug(1, "forward costs: nrec = %d", nrec);
00204
00205 if (abcol != NULL) {
00206 if (db_get_column(driver, Fi->table, abcol, &Column) != DB_OK)
00207 G_fatal_error(_("Column <%s> not found in table <%s>"),
00208 abcol, Fi->table);
00209
00210 bctype = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
00211
00212 if (bctype != DB_C_TYPE_INT && bctype != DB_C_TYPE_DOUBLE)
00213 G_fatal_error(_("Data type of column <%s> not supported (must be numeric)"),
00214 abcol);
00215
00216 db_CatValArray_init(&bvarr);
00217 nrec =
00218 db_select_CatValArray(driver, Fi->table, Fi->key, abcol, NULL,
00219 &bvarr);
00220 G_debug(1, "backward costs: nrec = %d", nrec);
00221 }
00222 }
00223
00224 skipped = 0;
00225
00226 G_message(_("Registering arcs..."));
00227
00228 for (i = 1; i <= nlines; i++) {
00229 G_percent(i, nlines, 1);
00230 dofw = dobw = 1;
00231 Vect_get_line_nodes(Map, i, &from, &to);
00232 type = Vect_read_line(Map, Points, Cats, i);
00233 if (!(type & ltype & (GV_LINE | GV_BOUNDARY)))
00234 continue;
00235
00236 if (afcol != NULL) {
00237 if (!(Vect_cat_get(Cats, afield, &cat))) {
00238 G_debug(2,
00239 "Category of field %d not attached to the line %d -> line skipped",
00240 afield, i);
00241 skipped += 2;
00242 continue;
00243 }
00244 else {
00245 if (fctype == DB_C_TYPE_INT) {
00246 ret = db_CatValArray_get_value_int(&fvarr, cat, &cost);
00247 dcost = cost;
00248 }
00249 else {
00250 ret =
00251 db_CatValArray_get_value_double(&fvarr, cat, &dcost);
00252 }
00253 if (ret != DB_OK) {
00254 G_warning(_("Database record for line %d (cat = %d, "
00255 "forward/both direction(s)) not found "
00256 "(forward/both direction(s) of line skipped)"),
00257 i, cat);
00258 dofw = 0;
00259 }
00260
00261 if (abcol != NULL) {
00262 if (bctype == DB_C_TYPE_INT) {
00263 ret =
00264 db_CatValArray_get_value_int(&bvarr, cat, &bcost);
00265 bdcost = bcost;
00266 }
00267 else {
00268 ret =
00269 db_CatValArray_get_value_double(&bvarr, cat,
00270 &bdcost);
00271 }
00272 if (ret != DB_OK) {
00273 G_warning(_("Database record for line %d (cat = %d, "
00274 "backword direction) not found"
00275 "(direction of line skipped)"), i, cat);
00276 dobw = 0;
00277 }
00278 }
00279 else {
00280 if (dofw)
00281 bdcost = dcost;
00282 else
00283 dobw = 0;
00284 }
00285 }
00286 }
00287 else {
00288 if (ll) {
00289 if (geo)
00290 dcost = Vect_line_geodesic_length(Points);
00291 else
00292 dcost = Vect_line_length(Points);
00293 }
00294 else
00295 dcost = Vect_line_length(Points);
00296
00297 bdcost = dcost;
00298 }
00299 if (dofw && dcost != -1) {
00300 cost = (dglInt32_t) Map->cost_multip * dcost;
00301 G_debug(5, "Add arc %d from %d to %d cost = %d", i, from, to,
00302 cost);
00303 ret =
00304 dglAddEdge(gr, (dglInt32_t) from, (dglInt32_t) to,
00305 (dglInt32_t) cost, (dglInt32_t) i);
00306 Map->edge_fcosts[i] = dcost;
00307 if (ret < 0)
00308 G_fatal_error("Cannot add network arc");
00309 }
00310
00311 G_debug(5, "bdcost = %f edge_bcosts = %f", bdcost,
00312 Map->edge_bcosts[i]);
00313 if (dobw && bdcost != -1) {
00314 bcost = (dglInt32_t) Map->cost_multip * bdcost;
00315 G_debug(5, "Add arc %d from %d to %d bcost = %d", -i, to, from,
00316 bcost);
00317 ret =
00318 dglAddEdge(gr, (dglInt32_t) to, (dglInt32_t) from,
00319 (dglInt32_t) bcost, (dglInt32_t) - i);
00320 Map->edge_bcosts[i] = bdcost;
00321 if (ret < 0)
00322 G_fatal_error(_("Cannot add network arc"));
00323 }
00324 }
00325
00326 if (afcol != NULL && skipped > 0)
00327 G_debug(2, "%d lines missing category of field %d skipped", skipped,
00328 afield);
00329
00330 if (afcol != NULL) {
00331 db_close_database_shutdown_driver(driver);
00332 db_CatValArray_free(&fvarr);
00333
00334 if (abcol != NULL) {
00335 db_CatValArray_free(&bvarr);
00336 }
00337 }
00338
00339
00340 G_debug(2, "Register nodes");
00341 if (ncol != NULL) {
00342 G_debug(2, "Set nodes' costs");
00343 if (nfield < 1)
00344 G_fatal_error("Node field < 1");
00345
00346 G_message(_("Setting node costs..."));
00347
00348 Fi = Vect_get_field(Map, nfield);
00349 if (Fi == NULL)
00350 G_fatal_error(_("Database connection not defined for layer %d"),
00351 nfield);
00352
00353 driver = db_start_driver_open_database(Fi->driver, Fi->database);
00354 if (driver == NULL)
00355 G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
00356 Fi->database, Fi->driver);
00357
00358
00359 if (db_get_column(driver, Fi->table, ncol, &Column) != DB_OK)
00360 G_fatal_error(_("Column <%s> not found in table <%s>"),
00361 ncol, Fi->table);
00362
00363 fctype = db_sqltype_to_Ctype(db_get_column_sqltype(Column));
00364
00365 if (fctype != DB_C_TYPE_INT && fctype != DB_C_TYPE_DOUBLE)
00366 G_fatal_error(_("Data type of column <%s> not supported (must be numeric)"),
00367 ncol);
00368
00369 db_CatValArray_init(&fvarr);
00370 nrec =
00371 db_select_CatValArray(driver, Fi->table, Fi->key, ncol, NULL,
00372 &fvarr);
00373 G_debug(1, "node costs: nrec = %d", nrec);
00374
00375 for (i = 1; i <= nnodes; i++) {
00376
00377
00378
00379 nlines = Vect_get_node_n_lines(Map, i);
00380 G_debug(2, " node = %d nlines = %d", i, nlines);
00381 cfound = 0;
00382 dcost = 0;
00383 for (j = 0; j < nlines; j++) {
00384 line = Vect_get_node_line(Map, i, j);
00385 G_debug(2, " line (%d) = %d", j, line);
00386 type = Vect_read_line(Map, NULL, Cats, line);
00387 if (!(type & GV_POINT))
00388 continue;
00389 if (Vect_cat_get(Cats, nfield, &cat)) {
00390
00391 if (fctype == DB_C_TYPE_INT) {
00392 ret =
00393 db_CatValArray_get_value_int(&fvarr, cat, &cost);
00394 dcost = cost;
00395 }
00396 else {
00397 ret =
00398 db_CatValArray_get_value_double(&fvarr, cat,
00399 &dcost);
00400 }
00401 if (ret != DB_OK) {
00402 G_warning(_("Database record for node %d (cat = %d) not found "
00403 "(cost set to 0)"), i, cat);
00404 }
00405 cfound = 1;
00406 break;
00407 }
00408 }
00409 if (!cfound) {
00410 G_debug(2,
00411 "Category of field %d not attached to any points in node %d"
00412 "(costs set to 0)", nfield, i);
00413 }
00414 if (dcost == -1) {
00415 cost = -1;
00416 }
00417 else {
00418 cost = (dglInt32_t) Map->cost_multip * dcost;
00419 }
00420 G_debug(3, "Set node's cost to %d", cost);
00421 dglNodeSet_Attr(gr, dglGetNode(gr, (dglInt32_t) i),
00422 (dglInt32_t *) (dglInt32_t) & cost);
00423 Map->node_costs[i] = dcost;
00424 }
00425 db_close_database_shutdown_driver(driver);
00426 db_CatValArray_free(&fvarr);
00427 }
00428
00429 G_message(_("Flattening the graph..."));
00430 ret = dglFlatten(gr);
00431 if (ret < 0)
00432 G_fatal_error(_("GngFlatten error"));
00433
00434
00435
00436 dglInitializeSPCache(gr, &(Map->spCache));
00437
00438 G_message(_("Graph was built"));
00439
00440 return 0;
00441 }
00442
00443
00462 int
00463 Vect_net_shortest_path(struct Map_info *Map, int from, int to,
00464 struct ilist *List, double *cost)
00465 {
00466 int i, line, *pclip, cArc, nRet;
00467 dglSPReport_s *pSPReport;
00468 dglInt32_t nDistance;
00469
00470 G_debug(3, "Vect_net_shortest_path(): from = %d, to = %d", from, to);
00471
00472
00473
00474
00475 if (List != NULL)
00476 Vect_reset_list(List);
00477
00478
00479 if (from == to) {
00480 if (cost != NULL)
00481 *cost = 0;
00482 return 0;
00483 }
00484
00485 From_node = from;
00486
00487 pclip = NULL;
00488 if (List != NULL) {
00489 nRet =
00490 dglShortestPath(&(Map->graph), &pSPReport, (dglInt32_t) from,
00491 (dglInt32_t) to, clipper, pclip, &(Map->spCache));
00492
00493
00494
00495
00496 }
00497 else {
00498 nRet =
00499 dglShortestDistance(&(Map->graph), &nDistance, (dglInt32_t) from,
00500 (dglInt32_t) to, clipper, pclip, &(Map->spCache));
00501
00502
00503
00504
00505 }
00506
00507 if (nRet == 0) {
00508
00509 if (cost != NULL)
00510 *cost = PORT_DOUBLE_MAX;
00511 return -1;
00512 }
00513 else if (nRet < 0) {
00514 G_warning(_("dglShortestPath error: %s"), dglStrerror(&(Map->graph)));
00515 return -1;
00516 }
00517
00518 if (List != NULL) {
00519 for (i = 0; i < pSPReport->cArc; i++) {
00520 line = dglEdgeGet_Id(&(Map->graph), pSPReport->pArc[i].pnEdge);
00521 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,
00522 line, pSPReport->pArc[i].nDistance);
00523 Vect_list_append(List, line);
00524 }
00525 }
00526
00527 if (cost != NULL) {
00528 if (List != NULL)
00529 *cost = (double)pSPReport->nDistance / Map->cost_multip;
00530 else
00531 *cost = (double)nDistance / Map->cost_multip;
00532 }
00533
00534 if (List != NULL) {
00535 cArc = pSPReport->cArc;
00536 dglFreeSPReport(&(Map->graph), pSPReport);
00537 }
00538 else
00539 cArc = 0;
00540
00541 return (cArc);
00542 }
00543
00556 int
00557 Vect_net_get_line_cost(struct Map_info *Map, int line, int direction,
00558 double *cost)
00559 {
00560
00561
00562 G_debug(5, "Vect_net_get_line_cost(): line = %d, dir = %d", line,
00563 direction);
00564
00565 if (direction == GV_FORWARD) {
00566
00567
00568
00569
00570
00571
00572 if (Map->edge_fcosts[line] == -1) {
00573 *cost = -1;
00574 return 0;
00575 }
00576 else
00577 *cost = Map->edge_fcosts[line];
00578 }
00579 else if (direction == GV_BACKWARD) {
00580
00581
00582
00583
00584
00585 if (Map->edge_bcosts[line] == -1) {
00586 *cost = -1;
00587 return 0;
00588 }
00589 else
00590 *cost = Map->edge_bcosts[line];
00591 G_debug(5, "Vect_net_get_line_cost(): edge_bcosts = %f",
00592 Map->edge_bcosts[line]);
00593 }
00594 else {
00595 G_fatal_error(_("Wrong line direction in Vect_net_get_line_cost()"));
00596 }
00597
00598 return 1;
00599 }
00600
00610 int Vect_net_get_node_cost(struct Map_info *Map, int node, double *cost)
00611 {
00612 G_debug(3, "Vect_net_get_node_cost(): node = %d", node);
00613
00614 *cost = Map->node_costs[node];
00615
00616 G_debug(3, " -> cost = %f", *cost);
00617
00618 return 1;
00619 }
00620
00639 int Vect_net_nearest_nodes(struct Map_info *Map,
00640 double x, double y, double z,
00641 int direction, double maxdist,
00642 int *node1, int *node2, int *ln, double *costs1,
00643 double *costs2, struct line_pnts *Points1,
00644 struct line_pnts *Points2, double *distance)
00645 {
00646 int line, n1, n2, nnodes;
00647 int npoints;
00648 int segment;
00649 static struct line_pnts *Points = NULL;
00650 double cx, cy, cz, c1, c2;
00651 double along;
00652 double length;
00653
00654 G_debug(3, "Vect_net_nearest_nodes() x = %f y = %f", x, y);
00655
00656
00657 if (node1)
00658 *node1 = 0;
00659 if (node2)
00660 *node2 = 0;
00661 if (ln)
00662 *ln = 0;
00663 if (costs1)
00664 *costs1 = PORT_DOUBLE_MAX;
00665 if (costs2)
00666 *costs2 = PORT_DOUBLE_MAX;
00667 if (Points1)
00668 Vect_reset_line(Points1);
00669 if (Points2)
00670 Vect_reset_line(Points2);
00671 if (distance)
00672 *distance = PORT_DOUBLE_MAX;
00673
00674 if (!Points)
00675 Points = Vect_new_line_struct();
00676
00677
00678 line = Vect_find_line(Map, x, y, z, Map->graph_line_type, maxdist, 0, 0);
00679
00680 if (line < 1)
00681 return 0;
00682
00683 Vect_read_line(Map, Points, NULL, line);
00684 npoints = Points->n_points;
00685 Vect_get_line_nodes(Map, line, &n1, &n2);
00686
00687 segment =
00688 Vect_line_distance(Points, x, y, z, 0, &cx, &cy, &cz, distance, NULL,
00689 &along);
00690
00691 G_debug(4, "line = %d n1 = %d n2 = %d segment = %d", line, n1, n2,
00692 segment);
00693
00694
00695 G_debug(4, "cx = %f cy = %f first = %f %f last = %f %f", cx, cy,
00696 Points->x[0], Points->y[0], Points->x[npoints - 1],
00697 Points->y[npoints - 1]);
00698
00699 if (Points->x[0] == cx && Points->y[0] == cy) {
00700 if (node1)
00701 *node1 = n1;
00702 if (ln)
00703 *ln = line;
00704 if (costs1)
00705 *costs1 = 0;
00706 if (Points1) {
00707 Vect_append_point(Points1, x, y, z);
00708 Vect_append_point(Points1, cx, cy, cz);
00709 }
00710 G_debug(3, "first node nearest");
00711 return 1;
00712 }
00713 if (Points->x[npoints - 1] == cx && Points->y[npoints - 1] == cy) {
00714 if (node1)
00715 *node1 = n2;
00716 if (ln)
00717 *ln = line;
00718 if (costs1)
00719 *costs1 = 0;
00720 if (Points1) {
00721 Vect_append_point(Points1, x, y, z);
00722 Vect_append_point(Points1, cx, cy, cz);
00723 }
00724 G_debug(3, "last node nearest");
00725 return 1;
00726 }
00727
00728 nnodes = 2;
00729
00730
00731
00732 if (direction == GV_FORWARD) {
00733 Vect_net_get_line_cost(Map, line, GV_BACKWARD, &c1);
00734 Vect_net_get_line_cost(Map, line, GV_FORWARD, &c2);
00735 }
00736 else {
00737 Vect_net_get_line_cost(Map, line, GV_FORWARD, &c1);
00738 Vect_net_get_line_cost(Map, line, GV_BACKWARD, &c2);
00739 }
00740
00741 if (c1 < 0)
00742 nnodes--;
00743 if (c2 < 0)
00744 nnodes--;
00745 if (nnodes == 0)
00746 return 0;
00747
00748 length = Vect_line_length(Points);
00749
00750 if (ln)
00751 *ln = line;
00752
00753 if (nnodes == 1 && c1 < 0) {
00754 if (node1)
00755 *node1 = n2;
00756
00757 if (costs1) {
00758 *costs1 = c2 * (length - along) / length;
00759 }
00760
00761 if (Points1) {
00762 int i;
00763
00764 if (direction == GV_FORWARD) {
00765 Vect_append_point(Points1, x, y, z);
00766 Vect_append_point(Points1, cx, cy, cz);
00767 for (i = segment; i < npoints; i++)
00768 Vect_append_point(Points1, Points->x[i], Points->y[i],
00769 Points->z[i]);
00770 }
00771 else {
00772 for (i = npoints - 1; i >= segment; i--)
00773 Vect_append_point(Points1, Points->x[i], Points->y[i],
00774 Points->z[i]);
00775
00776 Vect_append_point(Points1, cx, cy, cz);
00777 Vect_append_point(Points1, x, y, z);
00778 }
00779 }
00780 }
00781 else {
00782 if (node1)
00783 *node1 = n1;
00784 if (node2)
00785 *node2 = n2;
00786
00787 if (costs1) {
00788 *costs1 = c1 * along / length;
00789 }
00790
00791 if (costs2) {
00792 *costs2 = c2 * (length - along) / length;
00793 }
00794
00795 if (Points1) {
00796 int i;
00797
00798 if (direction == GV_FORWARD) {
00799 Vect_append_point(Points1, x, y, z);
00800 Vect_append_point(Points1, cx, cy, cz);
00801 for (i = segment - 1; i >= 0; i--)
00802 Vect_append_point(Points1, Points->x[i], Points->y[i],
00803 Points->z[i]);
00804 }
00805 else {
00806 for (i = 0; i < segment; i++)
00807 Vect_append_point(Points1, Points->x[i], Points->y[i],
00808 Points->z[i]);
00809
00810 Vect_append_point(Points1, cx, cy, cz);
00811 Vect_append_point(Points1, x, y, z);
00812 }
00813 }
00814
00815 if (Points2) {
00816 int i;
00817
00818 if (direction == GV_FORWARD) {
00819 Vect_append_point(Points2, x, y, z);
00820 Vect_append_point(Points2, cx, cy, cz);
00821 for (i = segment; i < npoints; i++)
00822 Vect_append_point(Points2, Points->x[i], Points->y[i],
00823 Points->z[i]);
00824 }
00825 else {
00826 for (i = npoints - 1; i >= segment; i--)
00827 Vect_append_point(Points2, Points->x[i], Points->y[i],
00828 Points->z[i]);
00829
00830 Vect_append_point(Points2, cx, cy, cz);
00831 Vect_append_point(Points2, x, y, z);
00832 }
00833 }
00834 }
00835
00836 return nnodes;
00837 }
00838
00857 int
00858 Vect_net_shortest_path_coor(struct Map_info *Map,
00859 double fx, double fy, double fz, double tx,
00860 double ty, double tz, double fmax, double tmax,
00861 double *costs, struct line_pnts *Points,
00862 struct ilist *List, struct line_pnts *FPoints,
00863 struct line_pnts *TPoints, double *fdist,
00864 double *tdist)
00865 {
00866 int fnode[2], tnode[2];
00867 double fcosts[2], tcosts[2], cur_cst;
00868 int nfnodes, ntnodes, fline, tline;
00869 static struct line_pnts *APoints, *SPoints, *fPoints[2], *tPoints[2];
00870 static struct ilist *LList;
00871 static int first = 1;
00872 int reachable, shortcut;
00873 int i, j, fn = 0, tn = 0;
00874
00875 G_debug(3, "Vect_net_shortest_path_coor()");
00876
00877 if (first) {
00878 APoints = Vect_new_line_struct();
00879 SPoints = Vect_new_line_struct();
00880 fPoints[0] = Vect_new_line_struct();
00881 fPoints[1] = Vect_new_line_struct();
00882 tPoints[0] = Vect_new_line_struct();
00883 tPoints[1] = Vect_new_line_struct();
00884 LList = Vect_new_list();
00885 first = 0;
00886 }
00887
00888
00889 if (costs)
00890 *costs = PORT_DOUBLE_MAX;
00891 if (Points)
00892 Vect_reset_line(Points);
00893 if (fdist)
00894 *fdist = 0;
00895 if (tdist)
00896 *tdist = 0;
00897 if (List)
00898 List->n_values = 0;
00899 if (FPoints)
00900 Vect_reset_line(FPoints);
00901 if (TPoints)
00902 Vect_reset_line(TPoints);
00903
00904
00905 fnode[0] = fnode[1] = tnode[0] = tnode[1] = 0;
00906
00907 nfnodes =
00908 Vect_net_nearest_nodes(Map, fx, fy, fz, GV_FORWARD, fmax, &(fnode[0]),
00909 &(fnode[1]), &fline, &(fcosts[0]),
00910 &(fcosts[1]), fPoints[0], fPoints[1], fdist);
00911 if (nfnodes == 0)
00912 return 0;
00913
00914 ntnodes =
00915 Vect_net_nearest_nodes(Map, tx, ty, tz, GV_BACKWARD, tmax,
00916 &(tnode[0]), &(tnode[1]), &tline, &(tcosts[0]),
00917 &(tcosts[1]), tPoints[0], tPoints[1], tdist);
00918 if (ntnodes == 0)
00919 return 0;
00920
00921 G_debug(3, "fline = %d tline = %d", fline, tline);
00922
00923 reachable = shortcut = 0;
00924 cur_cst = PORT_DOUBLE_MAX;
00925
00926
00927 if (fline == tline && (nfnodes > 1 || ntnodes > 1)) {
00928 double len, flen, tlen, c, fseg, tseg;
00929 double fcx, fcy, fcz, tcx, tcy, tcz;
00930
00931 Vect_read_line(Map, APoints, NULL, fline);
00932 len = Vect_line_length(APoints);
00933
00934
00935 fseg =
00936 Vect_line_distance(APoints, fx, fy, fz, 0, &fcx, &fcy, &fcz, NULL,
00937 NULL, &flen);
00938 tseg =
00939 Vect_line_distance(APoints, tx, ty, tz, 0, &tcx, &tcy, &tcz, NULL,
00940 NULL, &tlen);
00941
00942 Vect_reset_line(SPoints);
00943 if (flen == tlen) {
00944 cur_cst = 0;
00945 reachable = shortcut = 1;
00946 }
00947 else if (flen < tlen) {
00948 Vect_net_get_line_cost(Map, fline, GV_FORWARD, &c);
00949 if (c >= 0) {
00950 cur_cst = c * (tlen - flen) / len;
00951
00952 Vect_append_point(SPoints, fx, fy, fz);
00953 Vect_append_point(SPoints, fcx, fcy, fcz);
00954 for (i = fseg; i < tseg; i++)
00955 Vect_append_point(SPoints, APoints->x[i], APoints->y[i],
00956 APoints->z[i]);
00957
00958 Vect_append_point(SPoints, tcx, tcy, tcz);
00959 Vect_append_point(SPoints, tx, ty, tz);
00960
00961 reachable = shortcut = 1;
00962 }
00963 }
00964 else {
00965 Vect_net_get_line_cost(Map, fline, GV_BACKWARD, &c);
00966 if (c >= 0) {
00967 cur_cst = c * (flen - tlen) / len;
00968
00969 Vect_append_point(SPoints, fx, fy, fz);
00970 Vect_append_point(SPoints, fcx, fcy, fcz);
00971 for (i = fseg - 1; i >= tseg; i--)
00972 Vect_append_point(SPoints, APoints->x[i], APoints->y[i],
00973 APoints->z[i]);
00974
00975 Vect_append_point(SPoints, tcx, tcy, tcz);
00976 Vect_append_point(SPoints, tx, ty, tz);
00977
00978 reachable = shortcut = 1;
00979 }
00980 }
00981 }
00982
00983
00984 for (i = 0; i < nfnodes; i++) {
00985 for (j = 0; j < ntnodes; j++) {
00986 double ncst, cst;
00987 int ret;
00988
00989 G_debug(3, "i = %d fnode = %d j = %d tnode = %d", i, fnode[i], j,
00990 tnode[j]);
00991
00992 ret =
00993 Vect_net_shortest_path(Map, fnode[i], tnode[j], NULL, &ncst);
00994 if (ret == -1)
00995 continue;
00996
00997 cst = fcosts[i] + ncst + tcosts[j];
00998 if (reachable == 0 || cst < cur_cst) {
00999 cur_cst = cst;
01000 fn = i;
01001 tn = j;
01002 shortcut = 0;
01003 }
01004 reachable = 1;
01005 }
01006 }
01007
01008 G_debug(3, "reachable = %d shortcut = %d cur_cst = %f", reachable,
01009 shortcut, cur_cst);
01010 if (reachable) {
01011 int ret;
01012
01013 if (shortcut) {
01014 if (Points)
01015 Vect_append_points(Points, SPoints, GV_FORWARD);
01016 }
01017 else {
01018 ret =
01019 Vect_net_shortest_path(Map, fnode[fn], tnode[tn], LList,
01020 NULL);
01021 G_debug(3, "Number of lines %d", LList->n_values);
01022
01023 if (Points)
01024 Vect_append_points(Points, fPoints[fn], GV_FORWARD);
01025
01026 if (FPoints)
01027 Vect_append_points(FPoints, fPoints[fn], GV_FORWARD);
01028
01029 for (i = 0; i < LList->n_values; i++) {
01030 int line;
01031
01032 line = LList->value[i];
01033 G_debug(3, "i = %d line = %d", i, line);
01034
01035 if (Points) {
01036 Vect_read_line(Map, APoints, NULL, abs(line));
01037
01038 if (line > 0)
01039 Vect_append_points(Points, APoints, GV_FORWARD);
01040 else
01041 Vect_append_points(Points, APoints, GV_BACKWARD);
01042 }
01043
01044 if (List)
01045 Vect_list_append(List, line);
01046 }
01047
01048 if (Points)
01049 Vect_append_points(Points, tPoints[tn], GV_FORWARD);
01050
01051 if (TPoints)
01052 Vect_append_points(TPoints, tPoints[tn], GV_FORWARD);
01053 }
01054
01055 if (costs)
01056 *costs = cur_cst;
01057 }
01058
01059 return reachable;
01060 }