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
00437
00438
00439 G_message(_("Graph was built"));
00440
00441 return 0;
00442 }
00443
00444
00463 int
00464 Vect_net_shortest_path(struct Map_info *Map, int from, int to,
00465 struct ilist *List, double *cost)
00466 {
00467 int i, line, *pclip, cArc, nRet;
00468 dglSPReport_s *pSPReport;
00469 dglInt32_t nDistance;
00470
00471 G_debug(3, "Vect_net_shortest_path(): from = %d, to = %d", from, to);
00472
00473
00474
00475
00476 if (List != NULL)
00477 Vect_reset_list(List);
00478
00479
00480 if (from == to) {
00481 if (cost != NULL)
00482 *cost = 0;
00483 return 0;
00484 }
00485
00486 From_node = from;
00487
00488 pclip = NULL;
00489 if (List != NULL) {
00490
00491 nRet =
00492 dglShortestPath(&(Map->graph), &pSPReport, (dglInt32_t) from,
00493 (dglInt32_t) to, clipper, pclip, NULL);
00494 }
00495 else {
00496
00497 nRet =
00498 dglShortestDistance(&(Map->graph), &nDistance, (dglInt32_t) from,
00499 (dglInt32_t) to, clipper, pclip, NULL);
00500 }
00501
00502 if (nRet == 0) {
00503
00504 if (cost != NULL)
00505 *cost = PORT_DOUBLE_MAX;
00506 return -1;
00507 }
00508 else if (nRet < 0) {
00509 G_warning(_("dglShortestPath error: %s"), dglStrerror(&(Map->graph)));
00510 return -1;
00511 }
00512
00513 if (List != NULL) {
00514 for (i = 0; i < pSPReport->cArc; i++) {
00515 line = dglEdgeGet_Id(&(Map->graph), pSPReport->pArc[i].pnEdge);
00516 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,
00517 line, pSPReport->pArc[i].nDistance);
00518 Vect_list_append(List, line);
00519 }
00520 }
00521
00522 if (cost != NULL) {
00523 if (List != NULL)
00524 *cost = (double)pSPReport->nDistance / Map->cost_multip;
00525 else
00526 *cost = (double)nDistance / Map->cost_multip;
00527 }
00528
00529 if (List != NULL) {
00530 cArc = pSPReport->cArc;
00531 dglFreeSPReport(&(Map->graph), pSPReport);
00532 }
00533 else
00534 cArc = 0;
00535
00536 return (cArc);
00537 }
00538
00551 int
00552 Vect_net_get_line_cost(struct Map_info *Map, int line, int direction,
00553 double *cost)
00554 {
00555
00556
00557 G_debug(5, "Vect_net_get_line_cost(): line = %d, dir = %d", line,
00558 direction);
00559
00560 if (direction == GV_FORWARD) {
00561
00562
00563
00564
00565
00566
00567 if (Map->edge_fcosts[line] == -1) {
00568 *cost = -1;
00569 return 0;
00570 }
00571 else
00572 *cost = Map->edge_fcosts[line];
00573 }
00574 else if (direction == GV_BACKWARD) {
00575
00576
00577
00578
00579
00580 if (Map->edge_bcosts[line] == -1) {
00581 *cost = -1;
00582 return 0;
00583 }
00584 else
00585 *cost = Map->edge_bcosts[line];
00586 G_debug(5, "Vect_net_get_line_cost(): edge_bcosts = %f",
00587 Map->edge_bcosts[line]);
00588 }
00589 else {
00590 G_fatal_error(_("Wrong line direction in Vect_net_get_line_cost()"));
00591 }
00592
00593 return 1;
00594 }
00595
00605 int Vect_net_get_node_cost(struct Map_info *Map, int node, double *cost)
00606 {
00607 G_debug(3, "Vect_net_get_node_cost(): node = %d", node);
00608
00609 *cost = Map->node_costs[node];
00610
00611 G_debug(3, " -> cost = %f", *cost);
00612
00613 return 1;
00614 }
00615
00634 int Vect_net_nearest_nodes(struct Map_info *Map,
00635 double x, double y, double z,
00636 int direction, double maxdist,
00637 int *node1, int *node2, int *ln, double *costs1,
00638 double *costs2, struct line_pnts *Points1,
00639 struct line_pnts *Points2, double *distance)
00640 {
00641 int line, n1, n2, nnodes;
00642 int npoints;
00643 int segment;
00644 static struct line_pnts *Points = NULL;
00645 double cx, cy, cz, c1, c2;
00646 double along;
00647 double length;
00648
00649 G_debug(3, "Vect_net_nearest_nodes() x = %f y = %f", x, y);
00650
00651
00652 if (node1)
00653 *node1 = 0;
00654 if (node2)
00655 *node2 = 0;
00656 if (ln)
00657 *ln = 0;
00658 if (costs1)
00659 *costs1 = PORT_DOUBLE_MAX;
00660 if (costs2)
00661 *costs2 = PORT_DOUBLE_MAX;
00662 if (Points1)
00663 Vect_reset_line(Points1);
00664 if (Points2)
00665 Vect_reset_line(Points2);
00666 if (distance)
00667 *distance = PORT_DOUBLE_MAX;
00668
00669 if (!Points)
00670 Points = Vect_new_line_struct();
00671
00672
00673 line = Vect_find_line(Map, x, y, z, Map->graph_line_type, maxdist, 0, 0);
00674
00675 if (line < 1)
00676 return 0;
00677
00678 Vect_read_line(Map, Points, NULL, line);
00679 npoints = Points->n_points;
00680 Vect_get_line_nodes(Map, line, &n1, &n2);
00681
00682 segment =
00683 Vect_line_distance(Points, x, y, z, 0, &cx, &cy, &cz, distance, NULL,
00684 &along);
00685
00686 G_debug(4, "line = %d n1 = %d n2 = %d segment = %d", line, n1, n2,
00687 segment);
00688
00689
00690 G_debug(4, "cx = %f cy = %f first = %f %f last = %f %f", cx, cy,
00691 Points->x[0], Points->y[0], Points->x[npoints - 1],
00692 Points->y[npoints - 1]);
00693
00694 if (Points->x[0] == cx && Points->y[0] == cy) {
00695 if (node1)
00696 *node1 = n1;
00697 if (ln)
00698 *ln = line;
00699 if (costs1)
00700 *costs1 = 0;
00701 if (Points1) {
00702 Vect_append_point(Points1, x, y, z);
00703 Vect_append_point(Points1, cx, cy, cz);
00704 }
00705 G_debug(3, "first node nearest");
00706 return 1;
00707 }
00708 if (Points->x[npoints - 1] == cx && Points->y[npoints - 1] == cy) {
00709 if (node1)
00710 *node1 = n2;
00711 if (ln)
00712 *ln = line;
00713 if (costs1)
00714 *costs1 = 0;
00715 if (Points1) {
00716 Vect_append_point(Points1, x, y, z);
00717 Vect_append_point(Points1, cx, cy, cz);
00718 }
00719 G_debug(3, "last node nearest");
00720 return 1;
00721 }
00722
00723 nnodes = 2;
00724
00725
00726
00727 if (direction == GV_FORWARD) {
00728 Vect_net_get_line_cost(Map, line, GV_BACKWARD, &c1);
00729 Vect_net_get_line_cost(Map, line, GV_FORWARD, &c2);
00730 }
00731 else {
00732 Vect_net_get_line_cost(Map, line, GV_FORWARD, &c1);
00733 Vect_net_get_line_cost(Map, line, GV_BACKWARD, &c2);
00734 }
00735
00736 if (c1 < 0)
00737 nnodes--;
00738 if (c2 < 0)
00739 nnodes--;
00740 if (nnodes == 0)
00741 return 0;
00742
00743 length = Vect_line_length(Points);
00744
00745 if (ln)
00746 *ln = line;
00747
00748 if (nnodes == 1 && c1 < 0) {
00749 if (node1)
00750 *node1 = n2;
00751
00752 if (costs1) {
00753 *costs1 = c2 * (length - along) / length;
00754 }
00755
00756 if (Points1) {
00757 int i;
00758
00759 if (direction == GV_FORWARD) {
00760 Vect_append_point(Points1, x, y, z);
00761 Vect_append_point(Points1, cx, cy, cz);
00762 for (i = segment; i < npoints; i++)
00763 Vect_append_point(Points1, Points->x[i], Points->y[i],
00764 Points->z[i]);
00765 }
00766 else {
00767 for (i = npoints - 1; i >= segment; i--)
00768 Vect_append_point(Points1, Points->x[i], Points->y[i],
00769 Points->z[i]);
00770
00771 Vect_append_point(Points1, cx, cy, cz);
00772 Vect_append_point(Points1, x, y, z);
00773 }
00774 }
00775 }
00776 else {
00777 if (node1)
00778 *node1 = n1;
00779 if (node2)
00780 *node2 = n2;
00781
00782 if (costs1) {
00783 *costs1 = c1 * along / length;
00784 }
00785
00786 if (costs2) {
00787 *costs2 = c2 * (length - along) / length;
00788 }
00789
00790 if (Points1) {
00791 int i;
00792
00793 if (direction == GV_FORWARD) {
00794 Vect_append_point(Points1, x, y, z);
00795 Vect_append_point(Points1, cx, cy, cz);
00796 for (i = segment - 1; i >= 0; i--)
00797 Vect_append_point(Points1, Points->x[i], Points->y[i],
00798 Points->z[i]);
00799 }
00800 else {
00801 for (i = 0; i < segment; i++)
00802 Vect_append_point(Points1, Points->x[i], Points->y[i],
00803 Points->z[i]);
00804
00805 Vect_append_point(Points1, cx, cy, cz);
00806 Vect_append_point(Points1, x, y, z);
00807 }
00808 }
00809
00810 if (Points2) {
00811 int i;
00812
00813 if (direction == GV_FORWARD) {
00814 Vect_append_point(Points2, x, y, z);
00815 Vect_append_point(Points2, cx, cy, cz);
00816 for (i = segment; i < npoints; i++)
00817 Vect_append_point(Points2, Points->x[i], Points->y[i],
00818 Points->z[i]);
00819 }
00820 else {
00821 for (i = npoints - 1; i >= segment; i--)
00822 Vect_append_point(Points2, Points->x[i], Points->y[i],
00823 Points->z[i]);
00824
00825 Vect_append_point(Points2, cx, cy, cz);
00826 Vect_append_point(Points2, x, y, z);
00827 }
00828 }
00829 }
00830
00831 return nnodes;
00832 }
00833
00852 int
00853 Vect_net_shortest_path_coor(struct Map_info *Map,
00854 double fx, double fy, double fz, double tx,
00855 double ty, double tz, double fmax, double tmax,
00856 double *costs, struct line_pnts *Points,
00857 struct ilist *List, struct line_pnts *FPoints,
00858 struct line_pnts *TPoints, double *fdist,
00859 double *tdist)
00860 {
00861 int fnode[2], tnode[2];
00862 double fcosts[2], tcosts[2], cur_cst;
00863 int nfnodes, ntnodes, fline, tline;
00864 static struct line_pnts *APoints, *SPoints, *fPoints[2], *tPoints[2];
00865 static struct ilist *LList;
00866 static int first = 1;
00867 int reachable, shortcut;
00868 int i, j, fn = 0, tn = 0;
00869
00870 G_debug(3, "Vect_net_shortest_path_coor()");
00871
00872 if (first) {
00873 APoints = Vect_new_line_struct();
00874 SPoints = Vect_new_line_struct();
00875 fPoints[0] = Vect_new_line_struct();
00876 fPoints[1] = Vect_new_line_struct();
00877 tPoints[0] = Vect_new_line_struct();
00878 tPoints[1] = Vect_new_line_struct();
00879 LList = Vect_new_list();
00880 first = 0;
00881 }
00882
00883
00884 if (costs)
00885 *costs = PORT_DOUBLE_MAX;
00886 if (Points)
00887 Vect_reset_line(Points);
00888 if (fdist)
00889 *fdist = 0;
00890 if (tdist)
00891 *tdist = 0;
00892 if (List)
00893 List->n_values = 0;
00894 if (FPoints)
00895 Vect_reset_line(FPoints);
00896 if (TPoints)
00897 Vect_reset_line(TPoints);
00898
00899
00900 fnode[0] = fnode[1] = tnode[0] = tnode[1] = 0;
00901
00902 nfnodes =
00903 Vect_net_nearest_nodes(Map, fx, fy, fz, GV_FORWARD, fmax, &(fnode[0]),
00904 &(fnode[1]), &fline, &(fcosts[0]),
00905 &(fcosts[1]), fPoints[0], fPoints[1], fdist);
00906 if (nfnodes == 0)
00907 return 0;
00908
00909 ntnodes =
00910 Vect_net_nearest_nodes(Map, tx, ty, tz, GV_BACKWARD, tmax,
00911 &(tnode[0]), &(tnode[1]), &tline, &(tcosts[0]),
00912 &(tcosts[1]), tPoints[0], tPoints[1], tdist);
00913 if (ntnodes == 0)
00914 return 0;
00915
00916 G_debug(3, "fline = %d tline = %d", fline, tline);
00917
00918 reachable = shortcut = 0;
00919 cur_cst = PORT_DOUBLE_MAX;
00920
00921
00922 if (fline == tline && (nfnodes > 1 || ntnodes > 1)) {
00923 double len, flen, tlen, c, fseg, tseg;
00924 double fcx, fcy, fcz, tcx, tcy, tcz;
00925
00926 Vect_read_line(Map, APoints, NULL, fline);
00927 len = Vect_line_length(APoints);
00928
00929
00930 fseg =
00931 Vect_line_distance(APoints, fx, fy, fz, 0, &fcx, &fcy, &fcz, NULL,
00932 NULL, &flen);
00933 tseg =
00934 Vect_line_distance(APoints, tx, ty, tz, 0, &tcx, &tcy, &tcz, NULL,
00935 NULL, &tlen);
00936
00937 Vect_reset_line(SPoints);
00938 if (flen == tlen) {
00939 cur_cst = 0;
00940 reachable = shortcut = 1;
00941 }
00942 else if (flen < tlen) {
00943 Vect_net_get_line_cost(Map, fline, GV_FORWARD, &c);
00944 if (c >= 0) {
00945 cur_cst = c * (tlen - flen) / len;
00946
00947 Vect_append_point(SPoints, fx, fy, fz);
00948 Vect_append_point(SPoints, fcx, fcy, fcz);
00949 for (i = fseg; i < tseg; i++)
00950 Vect_append_point(SPoints, APoints->x[i], APoints->y[i],
00951 APoints->z[i]);
00952
00953 Vect_append_point(SPoints, tcx, tcy, tcz);
00954 Vect_append_point(SPoints, tx, ty, tz);
00955
00956 reachable = shortcut = 1;
00957 }
00958 }
00959 else {
00960 Vect_net_get_line_cost(Map, fline, GV_BACKWARD, &c);
00961 if (c >= 0) {
00962 cur_cst = c * (flen - tlen) / len;
00963
00964 Vect_append_point(SPoints, fx, fy, fz);
00965 Vect_append_point(SPoints, fcx, fcy, fcz);
00966 for (i = fseg - 1; i >= tseg; i--)
00967 Vect_append_point(SPoints, APoints->x[i], APoints->y[i],
00968 APoints->z[i]);
00969
00970 Vect_append_point(SPoints, tcx, tcy, tcz);
00971 Vect_append_point(SPoints, tx, ty, tz);
00972
00973 reachable = shortcut = 1;
00974 }
00975 }
00976 }
00977
00978
00979 for (i = 0; i < nfnodes; i++) {
00980 for (j = 0; j < ntnodes; j++) {
00981 double ncst, cst;
00982 int ret;
00983
00984 G_debug(3, "i = %d fnode = %d j = %d tnode = %d", i, fnode[i], j,
00985 tnode[j]);
00986
00987 ret =
00988 Vect_net_shortest_path(Map, fnode[i], tnode[j], NULL, &ncst);
00989 if (ret == -1)
00990 continue;
00991
00992 cst = fcosts[i] + ncst + tcosts[j];
00993 if (reachable == 0 || cst < cur_cst) {
00994 cur_cst = cst;
00995 fn = i;
00996 tn = j;
00997 shortcut = 0;
00998 }
00999 reachable = 1;
01000 }
01001 }
01002
01003 G_debug(3, "reachable = %d shortcut = %d cur_cst = %f", reachable,
01004 shortcut, cur_cst);
01005 if (reachable) {
01006 int ret;
01007
01008 if (shortcut) {
01009 if (Points)
01010 Vect_append_points(Points, SPoints, GV_FORWARD);
01011 }
01012 else {
01013 ret =
01014 Vect_net_shortest_path(Map, fnode[fn], tnode[tn], LList,
01015 NULL);
01016 G_debug(3, "Number of lines %d", LList->n_values);
01017
01018 if (Points)
01019 Vect_append_points(Points, fPoints[fn], GV_FORWARD);
01020
01021 if (FPoints)
01022 Vect_append_points(FPoints, fPoints[fn], GV_FORWARD);
01023
01024 for (i = 0; i < LList->n_values; i++) {
01025 int line;
01026
01027 line = LList->value[i];
01028 G_debug(3, "i = %d line = %d", i, line);
01029
01030 if (Points) {
01031 Vect_read_line(Map, APoints, NULL, abs(line));
01032
01033 if (line > 0)
01034 Vect_append_points(Points, APoints, GV_FORWARD);
01035 else
01036 Vect_append_points(Points, APoints, GV_BACKWARD);
01037 }
01038
01039 if (List)
01040 Vect_list_append(List, line);
01041 }
01042
01043 if (Points)
01044 Vect_append_points(Points, tPoints[tn], GV_FORWARD);
01045
01046 if (TPoints)
01047 Vect_append_points(TPoints, tPoints[tn], GV_FORWARD);
01048 }
01049
01050 if (costs)
01051 *costs = cur_cst;
01052 }
01053
01054 return reachable;
01055 }