v.net.alloc: add option to use costs from centers or costs to centers

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@68032 15284696-431f-4ddb-bdfa-cd5b030d7da7

vector/v.net.alloc/Makefile

@@ -3,7 +3,7 @@ MODULE_TOPDIR = ../..
 
 PGM = v.net.alloc
 
-LIBES = $(VECTORLIB) $(GISLIB)
+LIBES = $(VECTORLIB) $(GISLIB) $(GRAPHLIB)
 DEPENDENCIES = $(VECTORDEP) $(GISDEP)
 EXTRA_INC = $(VECT_INC)
 EXTRA_CFLAGS = $(VECT_CFLAGS)

vector/v.net.alloc/alloc.c

@@ -0,0 +1,266 @@
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/dbmi.h>
+#include <grass/dgl/graph.h>
+#include <grass/glocale.h>
+#include "alloc.h"
+
+int alloc_from_centers_loop_tt(struct Map_info *Map, NODE *Nodes,
+                               CENTER *Centers, int ncenters,
+                               int tucfield)
+{
+    int center, line, nlines, i;
+    int node1;
+    int cat;
+    struct line_cats *Cats;
+    struct line_pnts *Points;
+    double cost, n1cost, n2cost;
+    int ret;
+
+    Cats = Vect_new_cats_struct();
+    Points = Vect_new_line_struct();
+    
+    nlines = Vect_get_num_lines(Map);
+
+    for (i = 2; i <= (nlines * 2 + 2); i++) {
+	Nodes[i].center = -1;/* NOTE: first two items of Nodes are not used */
+	Nodes[i].cost = -1;
+	Nodes[i].edge = 0;
+    }
+
+    for (center = 0; center < ncenters; center++) {
+	G_percent(center, ncenters, 1);
+	node1 = Centers[center].node;
+	Vect_net_get_node_cost(Map, node1, &n1cost);
+	G_debug(2, "center = %d node = %d cat = %d", center, node1,
+		Centers[center].cat);
+
+	for (line = 1; line <= nlines; line++) {
+	    G_debug(5, "  node1 = %d line = %d", node1, line);
+	    Vect_net_get_node_cost(Map, line, &n2cost);
+	    /* closed, left it as not attached */
+
+	    if (Vect_read_line(Map, Points, Cats, line) < 0)
+		continue;
+	    if (Vect_get_line_type(Map, line) != GV_LINE)
+		continue;
+	    if (!Vect_cat_get(Cats, tucfield, &cat))
+		continue;
+
+	    for (i = 0; i < 2; i++) {
+		if (i == 1)
+		    cat *= -1;
+
+		ret =
+		    Vect_net_ttb_shortest_path(Map, node1, 0, cat, 1,
+					       tucfield, NULL,
+					       &cost);
+		if (ret == -1) {
+		    continue;
+		}		/* node unreachable */
+
+		/* We must add center node costs (not calculated by Vect_net_shortest_path() ), but
+		 *  only if center and node are not identical, because at the end node cost is add later */
+		if (ret != 1)
+		    cost += n1cost;
+
+		G_debug(5,
+			"Arc nodes: %d %d cost: %f (x old cent: %d old cost %f",
+			node1, line, cost, Nodes[line * 2 + i].center,
+			Nodes[line * 2 + i].cost);
+		if (Nodes[line * 2 + i].center == -1 ||
+		    Nodes[line * 2 + i].cost > cost) {
+		    Nodes[line * 2 + i].cost = cost;
+		    Nodes[line * 2 + i].center = center;
+		}
+	    }
+	}
+    }
+    G_percent(1, 1, 1);
+
+    Vect_destroy_cats_struct(Cats);
+    Vect_destroy_line_struct(Points);
+
+    return 0;
+}
+
+int alloc_from_centers(dglGraph_s *graph, NODE *Nodes, CENTER *Centers, int ncenters)
+{
+    int i, nnodes;
+    dglHeap_s heap;
+    dglEdgesetTraverser_s et;
+    int have_node_costs;
+    dglInt32_t ncost;
+
+    nnodes = dglGet_NodeCount(graph);
+
+    /* initialize nodes */
+    for (i = 1; i <= nnodes; i++) {
+	Nodes[i].cost = -1;
+	Nodes[i].center = -1;
+	Nodes[i].edge = 0;
+    }
+
+    ncost = 0;
+    have_node_costs = dglGet_NodeAttrSize(graph);
+
+    dglHeapInit(&heap);
+
+    for (i = 0; i < ncenters; i++) {
+	int v = Centers[i].node;
+
+	if (Nodes[v].cost == 0)
+	    continue;		/* ignore duplicates */
+	Nodes[v].cost = 0;		/* make sure all centers are processed first */
+	Nodes[v].center = i;
+	dglHeapData_u heap_data;
+
+	heap_data.ul = v;
+	dglHeapInsertMin(&heap, 0, ' ', heap_data);
+    }
+    while (1) {
+	dglInt32_t v, dist;
+	dglHeapNode_s heap_node;
+	dglHeapData_u heap_data;
+	dglInt32_t *edge;
+	dglInt32_t *node;
+
+	if (!dglHeapExtractMin(&heap, &heap_node))
+	    break;
+	v = heap_node.value.ul;
+	dist = heap_node.key;
+	if (Nodes[v].cost < dist)
+	    continue;
+
+	node = dglGetNode(graph, v);
+
+	if (have_node_costs && Nodes[v].edge) {
+	    memcpy(&ncost, dglNodeGet_Attr(graph, node),
+		   sizeof(ncost));
+	    if (ncost > 0)
+		dist += ncost;
+	    /* do not go through closed nodes */
+	    if (ncost < 0)
+		continue;
+	}
+
+	dglEdgeset_T_Initialize(&et, graph,
+				dglNodeGet_OutEdgeset(graph, node));
+
+	for (edge = dglEdgeset_T_First(&et); edge;
+	     edge = dglEdgeset_T_Next(&et)) {
+	    dglInt32_t *to = dglEdgeGet_Tail(graph, edge);
+	    dglInt32_t to_id = dglNodeGet_Id(graph, to);
+	    dglInt32_t d = dglEdgeGet_Cost(graph, edge);
+
+	    if (Nodes[to_id].cost < 0 || Nodes[to_id].cost > dist + d) {
+		Nodes[to_id].cost = dist + d;
+		Nodes[to_id].edge = dglEdgeGet_Id(graph, edge);
+		Nodes[to_id].center = Nodes[v].center;
+		heap_data.ul = to_id;
+		dglHeapInsertMin(&heap, dist + d, ' ', heap_data);
+	    }
+	}
+
+	dglEdgeset_T_Release(&et);
+    }
+
+    dglHeapFree(&heap, NULL);
+
+    return 0;
+}
+
+int alloc_to_centers(dglGraph_s *graph, NODE *Nodes, CENTER *Centers, int ncenters)
+{
+    int i, nnodes;
+    dglHeap_s heap;
+    dglEdgesetTraverser_s et;
+    int have_node_costs;
+    dglInt32_t ncost;
+
+    if (graph->Version < 2) {
+	G_warning("Directed graph must be version 2 or 3 for distances to centers");
+	return -1;
+    }
+
+    nnodes = dglGet_NodeCount(graph);
+
+    /* initialize nodes */
+    for (i = 1; i <= nnodes; i++) {
+	Nodes[i].cost = -1;
+	Nodes[i].center = -1;
+	Nodes[i].edge = 0;
+    }
+
+    ncost = 0;
+    have_node_costs = dglGet_NodeAttrSize(graph);
+
+    dglHeapInit(&heap);
+
+    for (i = 0; i < ncenters; i++) {
+	int v = Centers[i].node;
+
+	if (Nodes[v].cost == 0)
+	    continue;		/* ignore duplicates */
+	Nodes[v].cost = 0;		/* make sure all centers are processed first */
+	Nodes[v].center = i;
+	dglHeapData_u heap_data;
+
+	heap_data.ul = v;
+	dglHeapInsertMin(&heap, 0, ' ', heap_data);
+    }
+    while (1) {
+	dglInt32_t v, dist;
+	dglHeapNode_s heap_node;
+	dglHeapData_u heap_data;
+	dglInt32_t *edge;
+	dglInt32_t *node;
+
+	if (!dglHeapExtractMin(&heap, &heap_node))
+	    break;
+	v = heap_node.value.ul;
+	dist = heap_node.key;
+	if (Nodes[v].cost < dist)
+	    continue;
+
+	node = dglGetNode(graph, v);
+
+	if (have_node_costs && Nodes[v].edge) {
+	    memcpy(&ncost, dglNodeGet_Attr(graph, node),
+		   sizeof(ncost));
+	    if (ncost > 0)
+		dist += ncost;
+	    /* do not go through closed nodes */
+	    if (ncost < 0)
+		continue;
+	}
+
+	dglEdgeset_T_Initialize(&et, graph,
+				dglNodeGet_InEdgeset(graph, node));
+
+	for (edge = dglEdgeset_T_First(&et); edge;
+	     edge = dglEdgeset_T_Next(&et)) {
+	    dglInt32_t *from = dglEdgeGet_Head(graph, edge);
+	    dglInt32_t from_id = dglNodeGet_Id(graph, from);
+	    dglInt32_t d = dglEdgeGet_Cost(graph, edge);
+
+	    if (Nodes[from_id].cost < 0 || Nodes[from_id].cost > dist + d) {
+		Nodes[from_id].cost = dist + d;
+		Nodes[from_id].edge = dglEdgeGet_Id(graph, edge);
+		Nodes[from_id].center = Nodes[v].center;
+		heap_data.ul = from_id;
+		dglHeapInsertMin(&heap, dist + d, ' ', heap_data);
+	    }
+	}
+
+	dglEdgeset_T_Release(&et);
+    }
+
+    dglHeapFree(&heap, NULL);
+
+    return 0;
+}
+

vector/v.net.alloc/alloc.h

@@ -0,0 +1,20 @@
+
+typedef struct
+{
+    int cat;			/* category number */
+    int node;			/* node number */
+} CENTER;
+
+typedef struct
+{
+    int center;			/* neares center, initially -1 */
+    double cost;		/* costs from this center, initially not defined */
+    int edge;			/* edge to follow from this node */
+} NODE;
+
+int alloc_from_centers_loop_tt(struct Map_info *Map, NODE *Nodes,
+                               CENTER *Centers, int ncenters,
+                               int tucfield);
+
+int alloc_from_centers(dglGraph_s *graph, NODE *Nodes, CENTER *Centers, int ncenters);
+int alloc_to_centers(dglGraph_s *graph, NODE *Nodes, CENTER *Centers, int ncenters);

vector/v.net.alloc/main.c

@@ -5,10 +5,11 @@
  *
  * AUTHOR(S):    Radim Blazek
  *               Stepan Turek <stepan.turek seznam.cz> (turns support)
+ *               Markus Metz (costs from/to centers)
  *
  * PURPOSE:      Allocate subnets for nearest centers
  *               
- * COPYRIGHT:    (C) 2001, 2014 by the GRASS Development Team
+ * COPYRIGHT:    (C) 2001, 2016 by the GRASS Development Team
  *
  *               This program is free software under the 
  *               GNU General Public License (>=v2). 
@@ -23,27 +24,17 @@
 #include <grass/vector.h>
 #include <grass/dbmi.h>
 #include <grass/glocale.h>
+#include "alloc.h"
 
-typedef struct
-{
-    int cat;			/* category number */
-    int node;			/* node number */
-} CENTER;
-
-typedef struct
-{
-    int center;			/* neares center, initially -1 */
-    double cost;		/* costs from this center, initially not undefined */
-} NODE;
 
 int main(int argc, char **argv)
 {
-    int i, j, ret, center, line, center1, center2;
+    int i, ret, line, center1, center2;
     int nlines, nnodes, type, ltype, afield, nfield, geo, cat, tfield,
 	tucfield;
     int node1, node2;
-    double cost, e1cost, e2cost, n1cost, n2cost, s1cost, s2cost, l, l1, l2;
-    struct Option *map, *output;
+    double e1cost, e2cost, n1cost, n2cost, s1cost, s2cost, l, l1, l2;
+    struct Option *map, *output, *method_opt;
     struct Option *afield_opt, *nfield_opt, *afcol, *abcol, *ncol, *type_opt,
 	*term_opt, *tfield_opt, *tucfield_opt;
     struct Flag *geo_f, *turntable_f;
@@ -55,6 +46,8 @@ int main(int argc, char **argv)
     NODE *Nodes;
     struct line_cats *Cats;
     struct line_pnts *Points, *SPoints;
+    int graph_version;
+    int from_centers;
 
     /* initialize GIS environment */
     G_gisinit(argv[0]);		/* reads grass env, stores program name to G_program_name() */
@@ -74,6 +67,15 @@ int main(int argc, char **argv)
     map = G_define_standard_option(G_OPT_V_INPUT);
     output = G_define_standard_option(G_OPT_V_OUTPUT);
 
+    method_opt = G_define_option();
+    method_opt->key = "method";
+    method_opt->type = TYPE_STRING;
+    method_opt->required = NO;
+    method_opt->options = "from,to";
+    method_opt->answer = "from";
+    method_opt->description = _("Use costs from centers or costs to centers");
+    method_opt->guisection = _("Cost");
+
     term_opt = G_define_standard_option(G_OPT_V_CATS);
     term_opt->key = "center_cats";
     term_opt->required = YES;
@@ -177,13 +179,19 @@ int main(int argc, char **argv)
     tucfield = Vect_get_field_number(&Map, tucfield_opt->answer);
 
     /* Build graph */
+    graph_version = 1;
+    from_centers = 1;
+    if (method_opt->answer[0] == 't' && !turntable_f->answer) {
+	graph_version = 2;
+	from_centers = 0;
+    }
     if (turntable_f->answer)
 	Vect_net_ttb_build_graph(&Map, type, afield, nfield, tfield, tucfield,
 				 afcol->answer, abcol->answer, ncol->answer,
 				 geo, 0);
     else
 	Vect_net_build_graph(&Map, type, afield, nfield, afcol->answer,
-			     abcol->answer, ncol->answer, geo, 0);
+			     abcol->answer, ncol->answer, geo, graph_version);
 
     nnodes = Vect_get_num_nodes(&Map);
     nlines = Vect_get_num_lines(&Map);
@@ -235,101 +243,30 @@ int main(int argc, char **argv)
 
     /* alloc and reset space for all lines */
     if (turntable_f->answer) {
-
 	/* if turntable is used we are looking for lines as destinations, not the intersections (nodes) */
 	Nodes = (NODE *) G_calloc((nlines * 2 + 2), sizeof(NODE));
-	for (i = 2; i <= (nlines * 2 + 2); i++) {
-	    Nodes[i].center = -1;/* NOTE: first two items of Nodes are not used */
-	}
-
     }
     else {
 	Nodes = (NODE *) G_calloc((nnodes + 1), sizeof(NODE));
-	for (i = 1; i <= nnodes; i++) {
-	    Nodes[i].center = -1;
-	}
     }
 
     /* Fill Nodes by nearest center and costs from that center */
-    G_message(_("Calculating costs from centers ..."));
-
-    for (center = 0; center < ncenters; center++) {
-	G_percent(center, ncenters, 1);
-	node1 = Centers[center].node;
-	Vect_net_get_node_cost(&Map, node1, &n1cost);
-	G_debug(2, "center = %d node = %d cat = %d", center, node1,
-		Centers[center].cat);
-
-	if (turntable_f->answer)
-	    for (line = 1; line <= nlines; line++) {
-		G_debug(5, "  node1 = %d line = %d", node1, line);
-		Vect_net_get_node_cost(&Map, line, &n2cost);
-		/* closed, left it as not attached */
-
-		if (Vect_read_line(&Map, Points, Cats, line) < 0)
-		    continue;
-		if (Vect_get_line_type(&Map, line) != GV_LINE)
-		    continue;
-		if (!Vect_cat_get(Cats, tucfield, &cat))
-		    continue;
-
-		for (j = 0; j < 2; j++) {
-		    if (j == 1)
-			cat *= -1;
-
-		    ret =
-			Vect_net_ttb_shortest_path(&Map, node1, 0, cat, 1,
-						   tucfield, NULL,
-						   &cost);
-		    if (ret == -1) {
-			continue;
-		    }		/* node unreachable */
-
-		    /* We must add center node costs (not calculated by Vect_net_shortest_path() ), but
-		     *  only if center and node are not identical, because at the end node cost is add later */
-		    if (ret != 1)
-			cost += n1cost;
-
-		    G_debug(5,
-			    "Arc nodes: %d %d cost: %f (x old cent: %d old cost %f",
-			    node1, line, cost, Nodes[line * 2 + j].center,
-			    Nodes[line * 2 + j].cost);
-		    if (Nodes[line * 2 + j].center == -1 ||
-			cost < Nodes[line * 2 + j].cost) {
-			Nodes[line * 2 + j].cost = cost;
-			Nodes[line * 2 + j].center = center;
-		    }
-		}
-	    }
-	else
-	    for (node2 = 1; node2 <= nnodes; node2++) {
-		G_debug(5, "  node1 = %d node2 = %d", node1, node2);
-		Vect_net_get_node_cost(&Map, node2, &n2cost);
-		if (n2cost == -1) {
-		    continue;
-		}		/* closed, left it as not attached */
 
-		ret = Vect_net_shortest_path(&Map, node1, node2, NULL, &cost);
-		if (ret == -1) {
-		    continue;
-		}		/* node unreachable */
-
-		/* We must add center node costs (not calculated by Vect_net_shortest_path() ), but
-		 *  only if center and node are not identical, because at the end node cost is add later */
-		if (node1 != node2)
-		    cost += n1cost;
-
-		G_debug(5,
-			"Arc nodes: %d %d cost: %f (x old cent: %d old cost %f",
-			node1, node2, cost, Nodes[node2].center,
-			Nodes[node2].cost);
-		if (Nodes[node2].center == -1 || cost < Nodes[node2].cost) {
-		    Nodes[node2].cost = cost;
-		    Nodes[node2].center = center;
-		}
-	    }
+    if (turntable_f->answer) {
+	G_message(_("Calculating costs from centers ..."));
+	alloc_from_centers_loop_tt(&Map, Nodes, Centers, ncenters,
+				   tucfield);
+    }
+    else {
+	if (from_centers) {
+	    G_message(_("Calculating costs from centers ..."));
+	    alloc_from_centers(Vect_net_get_graph(&Map), Nodes, Centers, ncenters);
+	}
+	else {
+	    G_message(_("Calculating costs to centers ..."));
+	    alloc_to_centers(Vect_net_get_graph(&Map), Nodes, Centers, ncenters);
+	}
     }
-    G_percent(1, 1, 1);
 
     /* Write arcs to new map */
     if (Vect_open_new(&Out, output->answer, Vect_is_3d(&Map)) < 0)
@@ -362,8 +299,15 @@ int main(int argc, char **argv)
 	    Vect_net_get_node_cost(&Map, node2, &n2cost);
 	}
 
-	Vect_net_get_line_cost(&Map, line, GV_FORWARD, &e1cost);
-	Vect_net_get_line_cost(&Map, line, GV_BACKWARD, &e2cost);
+	if (from_centers) {
+	    Vect_net_get_line_cost(&Map, line, GV_FORWARD, &e1cost);
+	    Vect_net_get_line_cost(&Map, line, GV_BACKWARD, &e2cost);
+	}
+	else {
+	    /* from node to center */
+	    Vect_net_get_line_cost(&Map, line, GV_FORWARD, &e2cost);
+	    Vect_net_get_line_cost(&Map, line, GV_BACKWARD, &e1cost);
+	}
 
 	G_debug(3, "Line %d:", line);
 	G_debug(3, "Arc centers: %d %d (nodes: %d %d)", center1, center2,
@@ -391,7 +335,7 @@ int main(int argc, char **argv)
 		Vect_write_line(&Out, ltype, Points, Cats);
 	    }
 	    else {		/* each node in different area */
-		/* Check if direction is reachable */
+		/* Check if line is reachable from center */
 		if (center1 == -1 || n1cost == -1 || e1cost == -1) {	/* closed from first node */
 		    G_debug(3,
 			    "    -> arc is not reachable from 1. node -> alloc to 2. node");
@@ -408,7 +352,6 @@ int main(int argc, char **argv)
 		    Vect_write_line(&Out, ltype, Points, Cats);
 		    continue;
 		}
-
 		/* Now we know that arc is reachable from both sides */
 
 		/* Add costs of node to starting costs */

vector/v.net.alloc/v.net.alloc.html

@@ -1,8 +1,8 @@
 <h2>DESCRIPTION</h2>
 
-<em>v.net.alloc</em> allocates subnets for nearest centers 
-(direction from center). center nodes must be opened (costs &gt;= 0).
-Costs of center nodes are used in the calculation. 
+<em>v.net.alloc</em> allocates subnets for nearest centers. Center 
+nodes must be opened (costs &gt;= 0). Costs of center nodes are used in 
+the calculation. 
 <p>
 Costs may be either line lengths, or attributes saved in a 
 database table. These attribute values are taken as costs of whole 
@@ -19,6 +19,10 @@ For areas, costs will be calculated along boundary lines.
 The input vector needs to be prepared with <em>v.net operation=connect</em> 
 in order to connect points representing center nodes to the network.
 <p>
+The nearest center can be determined using either costs from the 
+nearest center or costs to the nearest center with option 
+<b>method</b>. See example below.
+<p>
 Application of flag <b>-t</b> enables a turntable support. 
 This flag requires additional parameters <b>turn_layer</b> and <b>turn_cat_layer</b> 
 that are otherwise ignored.
@@ -151,6 +155,44 @@ d.vect myroads_net_alloc_time cat=3 col=yellow layer=1
 d.vect myroads_net_time col=red icon=basic/triangle fcol=green size=12 type=point layer=2
 </pre></div>
 
+<p>Example 3: <em>Differences between costs from centers and costs to centers</em><br>
+
+Each lane of the two-lane road is a one-way road.
+<p>1. Subnetwork allocation from centers:
+<p><img src="v_net_alloc_from_centers.png" alt="v.net.alloc example from centers" border="1">
+<br>
+A center reaches any point following the one-way lanes.
+<p>2. Subnetwork allocation to centers:
+<p><img src="v_net_alloc_to_centers.png" alt="v.net.alloc example to centers" border="1">
+<br>
+Any node reaches reaches the nearest center following the one-way lanes.
+<p>
+In case of an accident, the ambulance should come from the nearest 
+'from' hospital and go to the nearest 'to' hospital.
+
+<div class="code"><pre>
+# North Carolina
+
+# center nodes are hospitals:
+# connect hospitals to streets as layer 2
+v.net input=streets_wake points=hospitals output=streets_hospitals operation=connect thresh=400 arc_layer=1 node_layer=2
+v.to.db map=streets_hospitals layer=1 type=line option=cat columns=cat
+
+# close oneway roads
+v.db.update map=streets_hospitals column=TF_COST value=-1 where="ONE_WAY = 'FT'"
+v.db.update map=streets_hospitals column=FT_COST value=-1 where="ONE_WAY = 'TF'"
+
+# add costs to newly created lines
+v.db.update map=streets_hospitals column=TF_COST value=0 where="cat > 49746"
+v.db.update map=streets_hospitals column=FT_COST value=0 where="cat > 49746"
+
+# from centers
+v.net.alloc in=streets_hospitals out=streets_hospitals_alloc_from center_cats=1-10000 arc_column=FT_COST arc_backward_column=TF_COST
+
+# to centers
+v.net.alloc in=streets_hospitals out=streets_hospitals_alloc_to method=to center_cats=1-10000 arc_column=FT_COST arc_backward_column=TF_COST
+</pre></div>
+
 <h2>SEE ALSO</h2>
 
 <em><a href="d.path.html">d.path</a></em>,