net.c

Go to the documentation of this file.
00001 /*
00002 ****************************************************************************
00003 *
00004 * MODULE:       Vector library 
00005 *               
00006 * AUTHOR(S):    Radim Blazek
00007 *
00008 * PURPOSE:      Higher level functions for reading/writing/manipulating vectors.
00009 *
00010 * COPYRIGHT:    (C) 2001 by the GRASS Development Team
00011 *
00012 *               This program is free software under the GNU General Public
00013 *               License (>=v2). Read the file COPYING that comes with GRASS
00014 *               for details.
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;   /* from node set in SP and used by clipper for first arc */  
00024 
00025 static int clipper ( dglGraph_s    *pgraph ,
00026                      dglSPClipInput_s  * pargIn ,
00027                      dglSPClipOutput_s * pargOut ,
00028                      void *          pvarg )         /* caller's pointer */
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 ) { /* do not clip first */
00043         if ( dglGet_NodeAttrSize(pgraph) > 0 ) {
00044             memcpy( &cost, dglNodeGet_Attr(pgraph, pargIn->pnNodeFrom), sizeof(cost) );
00045             if ( cost == -1 ) { /* closed, cannot go from this node except it is 'from' node */
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,   /* line type for arcs */
00091                         int afield,  /* arc costs field (if 0, use length) */
00092                         int nfield,  /* node costs field (if 0, do not use node costs) */
00093                         char *afcol, /* column with forward costs for arc */
00094                         char *abcol, /* column with backward costs for arc (if NULL, back = forward) */
00095                         char *ncol,  /* column with costs for nodes */
00096                         int geo,     /* use geodesic calculation for length (LL) */
00097                         int algorithm ) /* not used, in future code for 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     /* TODO int costs -> double (waiting for dglib) */
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; /* LL */
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     /* Allocate space for costs, later replace by functions reading costs from graph */
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     /* Set to -1 initialy */
00142     for ( i  = 1; i <= nlines; i++ ) {
00143         Map->edge_fcosts[i] = -1; /* forward */
00144         Map->edge_bcosts[i] = -1; /* backward */
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     /* --- Add arcs --- */
00162     /* Open db connection */
00163     if ( afcol != NULL ) {
00164         /* Get field info */
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         /* Open database */
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         /* Load costs to array */
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 ); /* must be before any continue */
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; /* Both directions */ 
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 { /* DB_C_TYPE_DOUBLE */
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 { /* DB_C_TYPE_DOUBLE */
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     /* Set node attributes */
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         /* Load costs to array */
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             /* TODO: what happens if we set attributes of not existing node (skipped lines,
00315              *       nodes without lines) */
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) ) { /* point with category of field found */
00327                     /* Set costs */
00328                     if ( fctype == DB_C_TYPE_INT ) { 
00329                         ret = db_CatValArray_get_value_int ( &fvarr, cat, &cost );
00330                         dcost = cost;
00331                     } else { /* DB_C_TYPE_DOUBLE */
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 ) { /* closed */
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     /* init SP cache */
00366     /* Disabled because of BUG1 in dglib. Without cache it is terribly slow, but with cache there
00367     *  are too many errors. */
00368     /* dglInitializeSPCache( gr, &(Map->spCache) ); */
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     /* Note : if from == to dgl goes to nearest node and returns back (dgl feature) => 
00396     *         check here for from == to */
00397     
00398     if ( List != NULL ) Vect_reset_list ( List);
00399 
00400     /* Check if from and to are identical, otherwise dglib returns path to neares node and back! */
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         /*nRet = dglShortestPath ( &(Map->graph), &pSPReport, from, to, clipper, pclip, &(Map->spCache));*/
00411         nRet = dglShortestPath ( &(Map->graph), &pSPReport, from, to, clipper, pclip, NULL);
00412     } else {
00413         /*nRet = dglShortestDistance ( &(Map->graph), &nDistance, from, to, clipper, pclip, &(Map->spCache));*/
00414         nRet = dglShortestDistance ( &(Map->graph), &nDistance, from, to, clipper, pclip, NULL);
00415     }
00416 
00417     if ( nRet == 0 ) {
00418         /* G_warning( "Destination node %d is unreachable from node %d\n" , to , from ); */         
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, /* this is the cost from clip() */
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 *  Returns in cost for given direction (GV_FORWARD,GV_BACKWARD) in *cost.
00460 *  cost is set to -1 if closed.
00461 *  Return : 1 : OK
00462 *           0 : does not exist (was not inserted)
00463 */
00464 int
00465 Vect_net_get_line_cost ( struct Map_info *Map, int line, int direction, double *cost )
00466 {
00467     /* dglInt32_t *pEdge; */
00468     
00469     G_debug (5, "Vect_net_get_line_cost(): line = %d, dir = %d", line, direction ); 
00470     
00471     if ( direction == GV_FORWARD ) {
00472         /* V1 has no index by line-id -> array used */
00473         /*
00474         pEdge = dglGetEdge ( &(Map->graph), line );
00475         if ( pEdge == NULL ) return 0;
00476         *cost = (double) dglEdgeGet_Cost ( &(Map->graph), pEdge );
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         pEdge = dglGetEdge ( &(Map->graph), -line );
00485         if ( pEdge == NULL ) return 0;
00486         *cost = (double) dglEdgeGet_Cost ( &(Map->graph), pEdge );
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 *  Returns in cost for given node in *cost.
00502 *  Return : 1 : OK
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; /* nearest line segment (first is 1) */
00549     static struct line_pnts *Points = NULL;
00550     double cx, cy, cz, c1, c2;
00551     double along; /* distance along the line to nearest point */
00552     double length;
00553 
00554     G_debug (3, "Vect_net_nearest_nodes() x = %f y = %f", x, y );
00555     
00556     /* Reset */
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     /* Find nearest line */
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     /* Check first or last point and return one node in that case */
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     /* c1 - costs to get from/to the first vertex */
00611     /* c2 - costs to get from/to the last vertex */
00612     if ( direction == GV_FORWARD ) { /* from point to net */
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; /* both directions closed */
00623     
00624     length = Vect_line_length ( Points );
00625     
00626     if ( ln ) *ln = line;
00627 
00628     if ( nnodes == 1 && c1 < 0 ) { /* first direction is closed, return node2 as node1 */
00629         if ( node1 ) *node1 = n2;
00630         
00631         if ( costs1 ) { /* to node 2, i.e. forward */
00632             *costs1 = c2 * (length - along) / length;
00633         }
00634 
00635         if ( Points1 ) { /* to node 2, i.e. forward */
00636             int i;
00637 
00638             if ( direction == GV_FORWARD ) { /* from point to net */
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 ) { /* to node 1, i.e. backward */
00656             *costs1 = c1 * along / length;
00657         }
00658 
00659         if ( costs2 ) { /* to node 2, i.e. forward */
00660             *costs2 = c2 * (length - along) / length;
00661         }
00662 
00663         if ( Points1 ) { /* to node 1, i.e. backward */
00664             int i;
00665 
00666             if ( direction == GV_FORWARD ) { /* from point to net */
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 ) { /* to node 2, i.e. forward */
00681             int i;
00682             
00683             if ( direction == GV_FORWARD ) { /* from point to net */
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]; /* nearest nodes, *node[1] is 0 if only one was found */
00737     double fcosts[2], tcosts[2], cur_cst; /* costs to nearest nodes on the network */
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     /* Reset */
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     /* Find nearest nodes */
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     /* It may happen, that 2 points are at the same line. */
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         /* distance along the line */
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 {  /* flen > tlen */
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     /* Find the shortest variant from maximum 4 */
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; /* not reachable */
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 

Generated on Sat Jul 22 22:05:57 2006 for GRASS by  doxygen 1.4.7