GRASS Programmer's Manual  6.4.1(2011)
net.c
Go to the documentation of this file.
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 }
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Defines