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

Generated on Thu Jul 16 13:21:17 2009 for GRASS Programmer's Manual by  doxygen 1.5.6