GRASS Programmer's Manual
6.4.1(2011)
|
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; /* from node set in SP and used by clipper for first arc */ 00029 00030 static int clipper(dglGraph_s * pgraph, 00031 dglSPClipInput_s * pargIn, 00032 dglSPClipOutput_s * pargOut, void *pvarg) 00033 { /* caller's pointer */ 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) { /* do not clip first */ 00047 if (dglGet_NodeAttrSize(pgraph) > 0) { 00048 memcpy(&cost, dglNodeGet_Attr(pgraph, pargIn->pnNodeFrom), 00049 sizeof(cost)); 00050 if (cost == -1) { /* closed, cannot go from this node except it is 'from' node */ 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 /* TODO int costs -> double (waiting for dglib) */ 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; /* LL */ 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 /* Allocate space for costs, later replace by functions reading costs from graph */ 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 /* Set to -1 initially */ 00147 for (i = 1; i <= nlines; i++) { 00148 Map->edge_fcosts[i] = -1; /* forward */ 00149 Map->edge_bcosts[i] = -1; /* backward */ 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 /* --- Add arcs --- */ 00172 /* Open db connection */ 00173 if (afcol != NULL) { 00174 /* Get field info */ 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 /* Open database */ 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 /* Load costs to array */ 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); /* must be before any continue */ 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; /* Both directions */ 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 { /* DB_C_TYPE_DOUBLE */ 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 { /* DB_C_TYPE_DOUBLE */ 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 /* Set node attributes */ 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 /* Load costs to array */ 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 /* TODO: what happens if we set attributes of not existing node (skipped lines, 00377 * nodes without lines) */ 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)) { /* point with category of field found */ 00390 /* Set costs */ 00391 if (fctype == DB_C_TYPE_INT) { 00392 ret = 00393 db_CatValArray_get_value_int(&fvarr, cat, &cost); 00394 dcost = cost; 00395 } 00396 else { /* DB_C_TYPE_DOUBLE */ 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) { /* closed */ 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 /* init SP cache */ 00435 /* disable to debug dglib cache */ 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 /* Note : if from == to dgl goes to nearest node and returns back (dgl feature) => 00473 * check here for from == to */ 00474 00475 if (List != NULL) 00476 Vect_reset_list(List); 00477 00478 /* Check if from and to are identical, otherwise dglib returns path to neares node and back! */ 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 /* comment out belwo and uncomment above to debug dglib cache */ 00493 /* nRet = 00494 dglShortestPath(&(Map->graph), &pSPReport, (dglInt32_t) from, 00495 (dglInt32_t) to, clipper, pclip, NULL); */ 00496 } 00497 else { 00498 nRet = 00499 dglShortestDistance(&(Map->graph), &nDistance, (dglInt32_t) from, 00500 (dglInt32_t) to, clipper, pclip, &(Map->spCache)); 00501 /* comment out belwo and uncomment above to debug dglib cache */ 00502 /* nRet = 00503 dglShortestDistance(&(Map->graph), &nDistance, (dglInt32_t) from, 00504 (dglInt32_t) to, clipper, pclip, NULL); */ 00505 } 00506 00507 if (nRet == 0) { 00508 /* G_warning( "Destination node %d is unreachable from node %d\n" , to , from ); */ 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, /* this is the cost from clip() */ 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 /* dglInt32_t *pEdge; */ 00561 00562 G_debug(5, "Vect_net_get_line_cost(): line = %d, dir = %d", line, 00563 direction); 00564 00565 if (direction == GV_FORWARD) { 00566 /* V1 has no index by line-id -> array used */ 00567 /* 00568 pEdge = dglGetEdge ( &(Map->graph), line ); 00569 if ( pEdge == NULL ) return 0; 00570 *cost = (double) dglEdgeGet_Cost ( &(Map->graph), pEdge ); 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 pEdge = dglGetEdge ( &(Map->graph), -line ); 00582 if ( pEdge == NULL ) return 0; 00583 *cost = (double) dglEdgeGet_Cost ( &(Map->graph), pEdge ); 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; /* nearest line segment (first is 1) */ 00649 static struct line_pnts *Points = NULL; 00650 double cx, cy, cz, c1, c2; 00651 double along; /* distance along the line to nearest point */ 00652 double length; 00653 00654 G_debug(3, "Vect_net_nearest_nodes() x = %f y = %f", x, y); 00655 00656 /* Reset */ 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 /* Find nearest line */ 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 /* Check first or last point and return one node in that case */ 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 /* c1 - costs to get from/to the first vertex */ 00731 /* c2 - costs to get from/to the last vertex */ 00732 if (direction == GV_FORWARD) { /* from point to net */ 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; /* both directions closed */ 00747 00748 length = Vect_line_length(Points); 00749 00750 if (ln) 00751 *ln = line; 00752 00753 if (nnodes == 1 && c1 < 0) { /* first direction is closed, return node2 as node1 */ 00754 if (node1) 00755 *node1 = n2; 00756 00757 if (costs1) { /* to node 2, i.e. forward */ 00758 *costs1 = c2 * (length - along) / length; 00759 } 00760 00761 if (Points1) { /* to node 2, i.e. forward */ 00762 int i; 00763 00764 if (direction == GV_FORWARD) { /* from point to net */ 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) { /* to node 1, i.e. backward */ 00788 *costs1 = c1 * along / length; 00789 } 00790 00791 if (costs2) { /* to node 2, i.e. forward */ 00792 *costs2 = c2 * (length - along) / length; 00793 } 00794 00795 if (Points1) { /* to node 1, i.e. backward */ 00796 int i; 00797 00798 if (direction == GV_FORWARD) { /* from point to net */ 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) { /* to node 2, i.e. forward */ 00816 int i; 00817 00818 if (direction == GV_FORWARD) { /* from point to net */ 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]; /* nearest nodes, *node[1] is 0 if only one was found */ 00867 double fcosts[2], tcosts[2], cur_cst; /* costs to nearest nodes on the network */ 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 /* Reset */ 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 /* Find nearest nodes */ 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 /* It may happen, that 2 points are at the same line. */ 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 /* distance along the line */ 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 { /* flen > tlen */ 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 /* Find the shortest variant from maximum 4 */ 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; /* not reachable */ 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 }