00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include <stdlib.h>
00018 #include <string.h>
00019 #include "gis.h"
00020 #include "dbmi.h"
00021 #include "Vect.h"
00022
00023 static int From_node;
00024
00025 static int clipper ( dglGraph_s *pgraph ,
00026 dglSPClipInput_s * pargIn ,
00027 dglSPClipOutput_s * pargOut ,
00028 void * pvarg )
00029 {
00030 dglInt32_t cost;
00031 dglInt32_t from;
00032
00033 G_debug ( 3, "Net: clipper()" );
00034
00035 from = dglNodeGet_Id(pgraph, pargIn->pnNodeFrom);
00036
00037 G_debug ( 3, " Edge = %d NodeFrom = %d NodeTo = %d edge cost = %d",
00038 dglEdgeGet_Id (pgraph, pargIn->pnEdge),
00039 from, dglNodeGet_Id(pgraph, pargIn->pnNodeTo),
00040 pargOut->nEdgeCost);
00041
00042 if ( from != From_node ) {
00043 if ( dglGet_NodeAttrSize(pgraph) > 0 ) {
00044 memcpy( &cost, dglNodeGet_Attr(pgraph, pargIn->pnNodeFrom), sizeof(cost) );
00045 if ( cost == -1 ) {
00046 G_debug ( 3, " closed node" );
00047 return 1;
00048 } else {
00049 G_debug ( 3, " EdgeCost += %d (node)", cost );
00050 pargOut->nEdgeCost += cost;
00051 }
00052 }
00053 } else {
00054 G_debug ( 3, " don't clip first node" );
00055 }
00056
00057 return 0;
00058 }
00059
00088 int
00089 Vect_net_build_graph ( struct Map_info *Map,
00090 int ltype,
00091 int afield,
00092 int nfield,
00093 char *afcol,
00094 char *abcol,
00095 char *ncol,
00096 int geo,
00097 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 ] = { 360000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
00107 struct field_info *Fi;
00108 dbDriver *driver;
00109 dbHandle handle;
00110 dbString stmt;
00111 dbColumn *Column;
00112 dbCatValArray fvarr, bvarr;
00113 int fctype, bctype, nrec;
00114
00115
00116 G_debug (1, "Vect_build_graph(): ltype = %d, afield = %d, nfield = %d", ltype, afield, nfield);
00117 G_debug (1, " afcol = %s, abcol = %s, ncol = %s", afcol, abcol, ncol);
00118
00119 fprintf ( stderr, "Building graph:\n");
00120
00121 Map->graph_line_type = ltype;
00122
00123 Points = Vect_new_line_struct ();
00124 Cats = Vect_new_cats_struct ();
00125
00126 ll = 0;
00127 if( G_projection() == 3) ll = 1;
00128
00129 if ( afcol == NULL && ll && !geo ) Map->cost_multip = 1000000;
00130 else Map->cost_multip = 1000;
00131
00132 nlines = Vect_get_num_lines ( Map );
00133 nnodes = Vect_get_num_nodes ( Map );
00134
00135 gr = &(Map->graph);
00136
00137
00138 Map->edge_fcosts = (double*) G_malloc ( (nlines+1) * sizeof(double) );
00139 Map->edge_bcosts = (double*) G_malloc ( (nlines+1) * sizeof(double) );
00140 Map->node_costs = (double*) G_malloc ( (nnodes+1) * sizeof(double) );
00141
00142 for ( i = 1; i <= nlines; i++ ) {
00143 Map->edge_fcosts[i] = -1;
00144 Map->edge_bcosts[i] = -1;
00145 }
00146 for ( i = 1; i <= nnodes; i++ ) {
00147 Map->node_costs[i] = 0;
00148 }
00149
00150 if ( ncol != NULL ) dglInitialize ( gr, 1, sizeof(dglInt32_t), 0, opaqueset );
00151 else dglInitialize ( gr, 1, 0, 0, opaqueset );
00152
00153 if ( gr == NULL ) G_fatal_error ("Cannot build network graph");
00154
00155 db_init_handle (&handle);
00156 db_init_string ( &stmt);
00157
00158 if ( abcol != NULL && afcol == NULL )
00159 G_fatal_error ("Forward costs column not specified");
00160
00161
00162
00163 if ( afcol != NULL ) {
00164
00165 if ( afield < 1 ) G_fatal_error ("Arc field < 1");
00166 Fi = Vect_get_field( Map, afield);
00167 if ( Fi == NULL ) G_fatal_error ("Cannot get field info");
00168
00169
00170 driver = db_start_driver_open_database ( Fi->driver, Fi->database );
00171 if ( driver == NULL )
00172 G_fatal_error("Cannot open database %s by driver %s", Fi->database, Fi->driver) ;
00173
00174
00175 if ( db_get_column ( driver, Fi->table, afcol, &Column ) != DB_OK)
00176 G_fatal_error("Cannot get column info");
00177
00178 fctype = db_sqltype_to_Ctype ( db_get_column_sqltype(Column) );
00179
00180 if ( fctype != DB_C_TYPE_INT && fctype != DB_C_TYPE_DOUBLE )
00181 G_fatal_error ( "Column type not supported" );
00182
00183 db_CatValArray_init ( &fvarr );
00184 nrec = db_select_CatValArray ( driver, Fi->table, Fi->key, afcol, NULL, &fvarr );
00185 G_debug (1, "forward costs: nrec = %d", nrec );
00186
00187 if ( abcol != NULL ) {
00188 if ( db_get_column ( driver, Fi->table, abcol, &Column ) != DB_OK)
00189 G_fatal_error("Cannot get column info");
00190
00191 bctype = db_sqltype_to_Ctype ( db_get_column_sqltype(Column) );
00192
00193 if ( bctype != DB_C_TYPE_INT && bctype != DB_C_TYPE_DOUBLE )
00194 G_fatal_error ( "Column type not supported" );
00195
00196 db_CatValArray_init ( &bvarr );
00197 nrec = db_select_CatValArray ( driver, Fi->table, Fi->key, abcol, NULL, &bvarr );
00198 G_debug (1, "backward costs: nrec = %d", nrec );
00199 }
00200 }
00201
00202 skipped = 0;
00203 fprintf ( stderr, "Registering arcs ... ");
00204 for ( i = 1; i <= nlines; i++ ) {
00205 G_percent ( i, nlines, 1 );
00206 dofw = dobw = 1;
00207 Vect_get_line_nodes ( Map, i, &from, &to );
00208 type = Vect_read_line ( Map, Points, Cats, i );
00209 if ( !(type & ltype & (GV_LINE | GV_BOUNDARY) ) ) continue;
00210
00211 if ( afcol != NULL ) {
00212 if ( !(Vect_cat_get(Cats, afield, &cat) ) ) {
00213 G_debug ( 2, "Category of field %d not attached to the line %d -> line skipped", afield, i);
00214 skipped += 2;
00215 continue;
00216 } else {
00217 if ( fctype == DB_C_TYPE_INT ) {
00218 ret = db_CatValArray_get_value_int ( &fvarr, cat, &cost );
00219 dcost = cost;
00220 } else {
00221 ret = db_CatValArray_get_value_double ( &fvarr, cat, &dcost );
00222 }
00223 if ( ret != DB_OK ) {
00224 G_warning ("Database record for line %d (cat = %d, forward/both direction(s)) not found "
00225 "(forward/both direction(s) of line skipped)", i, cat);
00226 dofw = 0;
00227 }
00228
00229 if ( abcol != NULL ) {
00230 if ( bctype == DB_C_TYPE_INT ) {
00231 ret = db_CatValArray_get_value_int ( &bvarr, cat, &bcost );
00232 bdcost = bcost;
00233 } else {
00234 ret = db_CatValArray_get_value_double ( &bvarr, cat, &bdcost );
00235 }
00236 if ( ret != DB_OK ) {
00237 G_warning ( "Database record for line %d (cat = %d, backword direction) not found"
00238 "(direction of line skipped)", i, cat);
00239 dobw = 0;
00240 }
00241 } else {
00242 if (dofw) bdcost = dcost;
00243 else dobw = 0;
00244 }
00245 }
00246 } else {
00247 if ( ll ) {
00248 if ( geo ) dcost = Vect_line_geodesic_length ( Points );
00249 else dcost = Vect_line_length ( Points );
00250 } else
00251 dcost = Vect_line_length ( Points );
00252
00253 bdcost = dcost;
00254 }
00255 if ( dofw && dcost != -1 ) {
00256 cost = (dglInt32_t) Map->cost_multip * dcost;
00257 G_debug (5, "Add arc %d from %d to %d cost = %d", i, from, to, cost);
00258 ret = dglAddEdge ( gr, from, to, cost, i );
00259 Map->edge_fcosts[i] = dcost;
00260 if ( ret < 0 ) G_fatal_error ("Cannot add network arc");
00261 }
00262
00263 G_debug (5, "bdcost = %f edge_bcosts = %f", bdcost, Map->edge_bcosts[i]);
00264 if ( dobw && bdcost != -1 ) {
00265 bcost = (dglInt32_t) Map->cost_multip * bdcost;
00266 G_debug (5, "Add arc %d from %d to %d bcost = %d", -i, to, from, bcost);
00267 ret = dglAddEdge ( gr, to, from, bcost, -i );
00268 Map->edge_bcosts[i] = bdcost;
00269 if ( ret < 0 ) G_fatal_error ("Cannot add network arc");
00270 }
00271 }
00272
00273 if ( afcol != NULL && skipped > 0 )
00274 G_debug ( 2, "%d lines missing category of field %d skipped", skipped, afield);
00275
00276 if ( afcol != NULL ) {
00277 db_close_database_shutdown_driver ( driver );
00278 db_CatValArray_free ( &fvarr);
00279
00280 if ( abcol != NULL ) {
00281 db_CatValArray_free ( &bvarr);
00282 }
00283 }
00284
00285
00286 G_debug ( 2, "Register nodes");
00287 if ( ncol != NULL ) {
00288 G_debug ( 2, "Set nodes' costs");
00289 if ( nfield < 1 ) G_fatal_error ("Node field < 1");
00290
00291 fprintf ( stderr, "Setting node costs ... ");
00292
00293 Fi = Vect_get_field( Map, nfield);
00294 if ( Fi == NULL ) G_fatal_error ("Cannot get field info");
00295
00296 driver = db_start_driver_open_database ( Fi->driver, Fi->database );
00297 if ( driver == NULL )
00298 G_fatal_error("Cannot open database %s by driver %s", Fi->database, Fi->driver) ;
00299
00300
00301 if ( db_get_column ( driver, Fi->table, ncol, &Column ) != DB_OK)
00302 G_fatal_error("Cannot get column info");
00303
00304 fctype = db_sqltype_to_Ctype ( db_get_column_sqltype(Column) );
00305
00306 if ( fctype != DB_C_TYPE_INT && fctype != DB_C_TYPE_DOUBLE )
00307 G_fatal_error ( "Column type not supported" );
00308
00309 db_CatValArray_init ( &fvarr );
00310 nrec = db_select_CatValArray ( driver, Fi->table, Fi->key, ncol, NULL, &fvarr );
00311 G_debug (1, "node costs: nrec = %d", nrec );
00312
00313 for ( i = 1; i <= nnodes; i++ ) {
00314
00315
00316
00317 nlines = Vect_get_node_n_lines ( Map, i );
00318 G_debug ( 2, " node = %d nlines = %d", i, nlines );
00319 cfound = 0;
00320 dcost = 0;
00321 for ( j = 0; j < nlines; j++ ) {
00322 line = Vect_get_node_line ( Map, i, j );
00323 G_debug ( 2, " line (%d) = %d", j, line );
00324 type = Vect_read_line ( Map, NULL, Cats, line);
00325 if ( !(type & GV_POINT) ) continue;
00326 if ( Vect_cat_get(Cats, nfield, &cat) ) {
00327
00328 if ( fctype == DB_C_TYPE_INT ) {
00329 ret = db_CatValArray_get_value_int ( &fvarr, cat, &cost );
00330 dcost = cost;
00331 } else {
00332 ret = db_CatValArray_get_value_double ( &fvarr, cat, &dcost );
00333 }
00334 if ( ret != DB_OK ) {
00335 G_warning ( "Database record for node %d (cat = %d) not found "
00336 "(cost set to 0)", i, cat);
00337 }
00338 cfound = 1;
00339 break;
00340 }
00341 }
00342 if ( !cfound ) {
00343 G_debug ( 2, "Category of field %d not attached to any points in node %d"
00344 "(costs set to 0)", nfield, i);
00345 }
00346 if ( dcost == -1 ) {
00347 cost = -1;
00348 } else {
00349 cost = (dglInt32_t) Map->cost_multip * dcost;
00350 }
00351 G_debug ( 3, "Set node's cost to %d", cost);
00352 dglNodeSet_Attr(gr, dglGetNode(gr,i), (dglInt32_t *) &cost);
00353 Map->node_costs[i] = dcost;
00354 }
00355 db_close_database_shutdown_driver ( driver );
00356 db_CatValArray_free ( &fvarr);
00357 fprintf ( stderr, "done.\n");
00358 }
00359
00360 fprintf ( stderr, "Flattening the graph ... ");
00361 ret = dglFlatten ( gr );
00362 if ( ret < 0 ) G_fatal_error ("GngFlatten error");
00363 fprintf ( stderr, "done.\n");
00364
00365
00366
00367
00368
00369
00370 fprintf ( stderr, "Graph was built.\n");
00371
00372 return 0;
00373 }
00374
00375
00385 #include<fcntl.h>
00386 int
00387 Vect_net_shortest_path ( struct Map_info *Map, int from, int to, struct ilist *List, double *cost )
00388 {
00389 int i, line, *pclip, cArc, nRet;
00390 dglSPReport_s * pSPReport;
00391 dglInt32_t nDistance;
00392
00393 G_debug (3, "Vect_net_shortest_path(): from = %d, to = %d", from, to );
00394
00395
00396
00397
00398 if ( List != NULL ) Vect_reset_list ( List);
00399
00400
00401 if ( from == to ) {
00402 if ( cost != NULL ) *cost = 0;
00403 return 0;
00404 }
00405
00406 From_node = from;
00407
00408 pclip = NULL;
00409 if ( List != NULL ) {
00410
00411 nRet = dglShortestPath ( &(Map->graph), &pSPReport, from, to, clipper, pclip, NULL);
00412 } else {
00413
00414 nRet = dglShortestDistance ( &(Map->graph), &nDistance, from, to, clipper, pclip, NULL);
00415 }
00416
00417 if ( nRet == 0 ) {
00418
00419 if ( cost != NULL )
00420 *cost = PORT_DOUBLE_MAX;
00421 return -1;
00422 }
00423 else if ( nRet < 0 ) {
00424 fprintf( stderr , "dglShortestPath error: %s\n", dglStrerror( &(Map->graph) ) );
00425 return -1;
00426 }
00427
00428 if ( List != NULL ) {
00429 for( i = 0 ; i < pSPReport->cArc ; i ++ ) {
00430 line = dglEdgeGet_Id(&(Map->graph), pSPReport->pArc[i].pnEdge);
00431 G_debug( 2, "From %ld to %ld - cost %ld user %d distance %ld" ,
00432 pSPReport->pArc[i].nFrom,
00433 pSPReport->pArc[i].nTo,
00434 dglEdgeGet_Cost(&(Map->graph),
00435 pSPReport->pArc[i].pnEdge) / Map->cost_multip,
00436 line,
00437 pSPReport->pArc[i].nDistance );
00438 Vect_list_append ( List, line );
00439 }
00440 }
00441
00442 if ( cost != NULL ) {
00443 if ( List != NULL )
00444 *cost = (double) pSPReport->nDistance / Map->cost_multip;
00445 else
00446 *cost = (double) nDistance / Map->cost_multip;
00447 }
00448
00449 if ( List != NULL ) {
00450 cArc = pSPReport->cArc;
00451 dglFreeSPReport( &(Map->graph), pSPReport );
00452 } else
00453 cArc = 0;
00454
00455 return (cArc);
00456 }
00457
00458
00459
00460
00461
00462
00463
00464 int
00465 Vect_net_get_line_cost ( struct Map_info *Map, int line, int direction, double *cost )
00466 {
00467
00468
00469 G_debug (5, "Vect_net_get_line_cost(): line = %d, dir = %d", line, direction );
00470
00471 if ( direction == GV_FORWARD ) {
00472
00473
00474
00475
00476
00477
00478 if ( Map->edge_fcosts[line] == -1 ) {
00479 *cost = -1;
00480 return 0;
00481 } else *cost = Map->edge_fcosts[line];
00482 } else if ( direction == GV_BACKWARD ) {
00483
00484
00485
00486
00487
00488 if ( Map->edge_bcosts[line] == -1 ) {
00489 *cost = -1;
00490 return 0;
00491 } else *cost = Map->edge_bcosts[line];
00492 G_debug (5, "Vect_net_get_line_cost(): edge_bcosts = %f", Map->edge_bcosts[line] );
00493 } else {
00494 G_fatal_error ("Wrong line direction in Vect_net_get_line_cost()" );
00495 }
00496
00497 return 1;
00498 }
00499
00500
00501
00502
00503
00504 int
00505 Vect_net_get_node_cost ( struct Map_info *Map, int node, double *cost )
00506 {
00507 G_debug (3, "Vect_net_get_node_cost(): node = %d", node );
00508
00509 *cost = Map->node_costs[node];
00510
00511 G_debug (3, " -> cost = %f", *cost );
00512
00513 return 1;
00514 }
00515
00539 int Vect_net_nearest_nodes ( struct Map_info *Map,
00540 double x, double y, double z,
00541 int direction, double maxdist,
00542 int *node1, int *node2, int *ln, double *costs1, double *costs2,
00543 struct line_pnts *Points1, struct line_pnts *Points2,
00544 double *distance )
00545 {
00546 int line, n1, n2, nnodes;
00547 int npoints;
00548 int segment;
00549 static struct line_pnts *Points = NULL;
00550 double cx, cy, cz, c1, c2;
00551 double along;
00552 double length;
00553
00554 G_debug (3, "Vect_net_nearest_nodes() x = %f y = %f", x, y );
00555
00556
00557 if ( node1 ) *node1 = 0;
00558 if ( node2 ) *node2 = 0;
00559 if ( ln ) *ln = 0;
00560 if ( costs1 ) *costs1 = PORT_DOUBLE_MAX;
00561 if ( costs2 ) *costs2 = PORT_DOUBLE_MAX;
00562 if ( Points1 ) Vect_reset_line ( Points1 );
00563 if ( Points2 ) Vect_reset_line ( Points2 );
00564 if ( distance ) *distance = PORT_DOUBLE_MAX;
00565
00566 if ( !Points ) Points = Vect_new_line_struct();
00567
00568
00569 line = Vect_find_line ( Map, x, y, z, Map->graph_line_type, maxdist, 0, 0 );
00570
00571 if ( line < 1 ) return 0;
00572
00573 Vect_read_line ( Map, Points, NULL, line );
00574 npoints = Points->n_points;
00575 Vect_get_line_nodes ( Map, line, &n1, &n2);
00576
00577 segment = Vect_line_distance ( Points, x, y, z, 0, &cx, &cy, &cz, distance, NULL, &along);
00578
00579 G_debug (4, "line = %d n1 = %d n2 = %d segment = %d", line, n1, n2, segment );
00580
00581
00582 G_debug (4, "cx = %f cy = %f first = %f %f last = %f %f", cx, cy, Points->x[0], Points->y[0],
00583 Points->x[npoints-1], Points->y[npoints-1] );
00584
00585 if ( Points->x[0] == cx && Points->y[0] == cy ) {
00586 if ( node1 ) *node1 = n1;
00587 if ( ln ) *ln = line;
00588 if ( costs1 ) *costs1 = 0;
00589 if ( Points1 ) {
00590 Vect_append_point ( Points1, x, y, z );
00591 Vect_append_point ( Points1, cx, cy, cz );
00592 }
00593 G_debug (3, "first node nearest");
00594 return 1;
00595 }
00596 if ( Points->x[npoints-1] == cx && Points->y[npoints-1] == cy ) {
00597 if ( node1 ) *node1 = n2;
00598 if ( ln ) *ln = line;
00599 if ( costs1 ) *costs1 = 0;
00600 if ( Points1 ) {
00601 Vect_append_point ( Points1, x, y, z );
00602 Vect_append_point ( Points1, cx, cy, cz );
00603 }
00604 G_debug (3, "last node nearest");
00605 return 1;
00606 }
00607
00608 nnodes = 2;
00609
00610
00611
00612 if ( direction == GV_FORWARD ) {
00613 Vect_net_get_line_cost ( Map, line, GV_BACKWARD, &c1 );
00614 Vect_net_get_line_cost ( Map, line, GV_FORWARD, &c2 );
00615 } else {
00616 Vect_net_get_line_cost ( Map, line, GV_FORWARD, &c1 );
00617 Vect_net_get_line_cost ( Map, line, GV_BACKWARD, &c2 );
00618 }
00619
00620 if ( c1 < 0 ) nnodes--;
00621 if ( c2 < 0 ) nnodes--;
00622 if ( nnodes == 0 ) return 0;
00623
00624 length = Vect_line_length ( Points );
00625
00626 if ( ln ) *ln = line;
00627
00628 if ( nnodes == 1 && c1 < 0 ) {
00629 if ( node1 ) *node1 = n2;
00630
00631 if ( costs1 ) {
00632 *costs1 = c2 * (length - along) / length;
00633 }
00634
00635 if ( Points1 ) {
00636 int i;
00637
00638 if ( direction == GV_FORWARD ) {
00639 Vect_append_point ( Points1, x, y, z );
00640 Vect_append_point ( Points1, cx, cy, cz );
00641 for ( i = segment; i < npoints; i++ )
00642 Vect_append_point ( Points1, Points->x[i], Points->y[i], Points->z[i] );
00643 } else {
00644 for ( i = npoints - 1; i >= segment; i-- )
00645 Vect_append_point ( Points1, Points->x[i], Points->y[i], Points->z[i] );
00646
00647 Vect_append_point ( Points1, cx, cy, cz );
00648 Vect_append_point ( Points1, x, y, z );
00649 }
00650 }
00651 } else {
00652 if ( node1 ) *node1 = n1;
00653 if ( node2 ) *node2 = n2;
00654
00655 if ( costs1 ) {
00656 *costs1 = c1 * along / length;
00657 }
00658
00659 if ( costs2 ) {
00660 *costs2 = c2 * (length - along) / length;
00661 }
00662
00663 if ( Points1 ) {
00664 int i;
00665
00666 if ( direction == GV_FORWARD ) {
00667 Vect_append_point ( Points1, x, y, z );
00668 Vect_append_point ( Points1, cx, cy, cz );
00669 for ( i = segment - 1; i >= 0; i-- )
00670 Vect_append_point ( Points1, Points->x[i], Points->y[i], Points->z[i] );
00671 } else {
00672 for ( i = 0; i < segment; i++ )
00673 Vect_append_point ( Points1, Points->x[i], Points->y[i], Points->z[i] );
00674
00675 Vect_append_point ( Points1, cx, cy, cz );
00676 Vect_append_point ( Points1, x, y, z );
00677 }
00678 }
00679
00680 if ( Points2 ) {
00681 int i;
00682
00683 if ( direction == GV_FORWARD ) {
00684 Vect_append_point ( Points2, x, y, z );
00685 Vect_append_point ( Points2, cx, cy, cz );
00686 for ( i = segment; i < npoints; i++ )
00687 Vect_append_point ( Points2, Points->x[i], Points->y[i], Points->z[i] );
00688 } else {
00689 for ( i = npoints - 1; i >= segment; i-- )
00690 Vect_append_point ( Points2, Points->x[i], Points->y[i], Points->z[i] );
00691
00692 Vect_append_point ( Points2, cx, cy, cz );
00693 Vect_append_point ( Points2, x, y, z );
00694 }
00695 }
00696 }
00697
00698 return nnodes;
00699 }
00700
00728 int
00729 Vect_net_shortest_path_coor ( struct Map_info *Map,
00730 double fx, double fy, double fz, double tx, double ty, double tz,
00731 double fmax, double tmax,
00732 double *costs, struct line_pnts *Points,
00733 struct ilist *List, struct line_pnts *FPoints, struct line_pnts *TPoints,
00734 double *fdist, double *tdist )
00735 {
00736 int fnode[2], tnode[2];
00737 double fcosts[2], tcosts[2], cur_cst;
00738 int nfnodes, ntnodes, fline, tline;
00739 static struct line_pnts *APoints, *SPoints, *fPoints[2], *tPoints[2];
00740 static struct ilist *LList;
00741 static int first = 1;
00742 int reachable, shortcut;
00743 int i, j, fn, tn;
00744
00745 G_debug (3, "Vect_net_shortest_path_coor()");
00746
00747 if ( first ) {
00748 APoints = Vect_new_line_struct();
00749 SPoints = Vect_new_line_struct();
00750 fPoints[0] = Vect_new_line_struct();
00751 fPoints[1] = Vect_new_line_struct();
00752 tPoints[0] = Vect_new_line_struct();
00753 tPoints[1] = Vect_new_line_struct();
00754 LList = Vect_new_list ();
00755 first = 0;
00756 }
00757
00758
00759 if ( costs ) *costs = PORT_DOUBLE_MAX;
00760 if ( Points ) Vect_reset_line ( Points );
00761 if ( fdist ) *fdist = 0;
00762 if ( tdist ) *tdist = 0;
00763 if ( List ) List->n_values = 0;
00764 if ( FPoints ) Vect_reset_line ( FPoints );
00765 if ( TPoints ) Vect_reset_line ( TPoints );
00766
00767
00768 fnode[0] = fnode[1] = tnode[0] = tnode[1] = 0;
00769
00770 nfnodes = Vect_net_nearest_nodes ( Map, fx, fy, fz, GV_FORWARD, fmax, &(fnode[0]), &(fnode[1]), &fline,
00771 &(fcosts[0]), &(fcosts[1]), fPoints[0], fPoints[1], fdist );
00772 if ( nfnodes == 0 ) return 0;
00773
00774 ntnodes = Vect_net_nearest_nodes ( Map, tx, ty, tz, GV_BACKWARD, tmax, &(tnode[0]), &(tnode[1]), &tline,
00775 &(tcosts[0]), &(tcosts[1]), tPoints[0], tPoints[1], tdist );
00776 if ( ntnodes == 0 ) return 0;
00777
00778 G_debug (3, "fline = %d tline = %d", fline, tline);
00779
00780 reachable = shortcut = 0;
00781 cur_cst = PORT_DOUBLE_MAX;
00782
00783
00784 if ( fline == tline && (nfnodes > 1 || ntnodes > 1) ) {
00785 double len, flen, tlen, c, fseg, tseg;
00786 double fcx, fcy, fcz, tcx, tcy, tcz;
00787
00788 Vect_read_line ( Map, APoints, NULL, fline );
00789 len = Vect_line_length ( APoints );
00790
00791
00792 fseg = Vect_line_distance ( APoints, fx, fy, fz, 0, &fcx, &fcy, &fcz, NULL, NULL, &flen);
00793 tseg = Vect_line_distance ( APoints, tx, ty, tz, 0, &tcx, &tcy, &tcz, NULL, NULL, &tlen);
00794
00795 Vect_reset_line ( SPoints );
00796 if ( flen == tlen ) {
00797 cur_cst = 0;
00798 reachable = shortcut = 1;
00799 } else if ( flen < tlen ) {
00800 Vect_net_get_line_cost ( Map, fline, GV_FORWARD, &c );
00801 if ( c >= 0 ) {
00802 cur_cst = c * (tlen - flen) / len;
00803
00804 Vect_append_point (SPoints, fx, fy, fz);
00805 Vect_append_point (SPoints, fcx, fcy, fcz);
00806 for ( i = fseg; i < tseg; i++ )
00807 Vect_append_point (SPoints, APoints->x[i], APoints->y[i], APoints->z[i]);
00808
00809 Vect_append_point (SPoints, tcx, tcy, tcz);
00810 Vect_append_point (SPoints, tx, ty, tz);
00811
00812 reachable = shortcut = 1;
00813 }
00814 } else {
00815 Vect_net_get_line_cost ( Map, fline, GV_BACKWARD, &c );
00816 if ( c >= 0 ) {
00817 cur_cst = c * (flen - tlen) / len;
00818
00819 Vect_append_point (SPoints, fx, fy, fz);
00820 Vect_append_point (SPoints, fcx, fcy, fcz);
00821 for ( i = fseg - 1; i >= tseg; i-- )
00822 Vect_append_point (SPoints, APoints->x[i], APoints->y[i], APoints->z[i]);
00823
00824 Vect_append_point (SPoints, tcx, tcy, tcz);
00825 Vect_append_point (SPoints, tx, ty, tz);
00826
00827 reachable = shortcut = 1;
00828 }
00829 }
00830 }
00831
00832
00833 for ( i = 0; i < nfnodes; i++ ) {
00834 for ( j = 0; j < ntnodes; j++ ) {
00835 double ncst, cst;
00836 int ret;
00837
00838 G_debug (3, "i = %d fnode = %d j = %d tnode = %d", i, fnode[i], j, tnode[j]);
00839
00840 ret = Vect_net_shortest_path ( Map, fnode[i], tnode[j], NULL, &ncst);
00841 if ( ret == -1 ) continue;
00842
00843 cst = fcosts[i] + ncst + tcosts[j];
00844 if ( reachable == 0 || cst < cur_cst ) {
00845 cur_cst = cst;
00846 fn = i;
00847 tn = j;
00848 shortcut = 0;
00849 }
00850 reachable = 1;
00851 }
00852 }
00853
00854 G_debug (3, "reachable = %d shortcut = %d cur_cst = %f", reachable, shortcut, cur_cst);
00855 if ( reachable ) {
00856 int ret;
00857
00858 if ( shortcut) {
00859 if ( Points )
00860 Vect_append_points ( Points, SPoints, GV_FORWARD );
00861 } else {
00862 ret = Vect_net_shortest_path ( Map, fnode[fn], tnode[tn], LList, NULL);
00863 G_debug (3, "Number of lines %d", LList->n_values);
00864
00865 if ( Points )
00866 Vect_append_points ( Points, fPoints[fn], GV_FORWARD );
00867
00868 if ( FPoints )
00869 Vect_append_points ( FPoints, fPoints[fn], GV_FORWARD );
00870
00871 for ( i = 0; i < LList->n_values; i++ ) {
00872 int line;
00873
00874 line = LList->value[i];
00875 G_debug (3, "i = %d line = %d", i, line);
00876
00877 if ( Points ) {
00878 Vect_read_line ( Map, APoints, NULL, abs(line) );
00879
00880 if ( line > 0 )
00881 Vect_append_points ( Points, APoints, GV_FORWARD );
00882 else
00883 Vect_append_points ( Points, APoints, GV_BACKWARD );
00884 }
00885
00886 if ( List ) Vect_list_append ( List, line );
00887 }
00888
00889 if ( Points )
00890 Vect_append_points ( Points, tPoints[tn], GV_FORWARD );
00891
00892 if ( TPoints )
00893 Vect_append_points ( TPoints, tPoints[tn], GV_FORWARD );
00894 }
00895
00896 if ( costs ) *costs = cur_cst;
00897 }
00898
00899 return reachable;
00900 }
00901