1 /*!
2 \file vector/neta/path.c
3
4 \brief Network Analysis library - shortest path
5
6 Shortest paths from a set of nodes.
7
8 (C) 2009-2010 by Daniel Bundala, and the GRASS Development Team
9
10 This program is free software under the GNU General Public License
11 (>=v2). Read the file COPYING that comes with GRASS for details.
12
13 \author Daniel Bundala (Google Summer of Code 2009)
14 \author Markus Metz
15 */
16
17 #include <stdio.h>
18 #include <stdlib.h>
19 #include <grass/gis.h>
20 #include <grass/vector.h>
21 #include <grass/glocale.h>
22 #include <grass/dgl/graph.h>
23 #include <grass/neta.h>
24
25 /*!
26 \brief Computes shortest paths to every node from nodes in "from".
27
28 Array "dst" contains the cost of the path or -1 if the node is not
29 reachable. Prev contains edges from predecessor along the shortest
30 path.
31
32 \param graph input graph
33 \param from list of 'from' positions
34 \param[out] dst array of costs to reach nodes
35 \param[out] prev array of edges from predecessor along the shortest path
36
37 \return 0 on success
38 \return -1 on failure
39 */
NetA_distance_from_points(dglGraph_s * graph,struct ilist * from,int * dst,dglInt32_t ** prev)40 int NetA_distance_from_points(dglGraph_s *graph, struct ilist *from,
41 int *dst, dglInt32_t **prev)
42 {
43 int i, nnodes;
44 dglHeap_s heap;
45 int have_node_costs;
46 dglInt32_t ncost;
47
48 nnodes = dglGet_NodeCount(graph);
49 dglEdgesetTraverser_s et;
50
51 /* initialize costs and edge list */
52 for (i = 1; i <= nnodes; i++) {
53 dst[i] = -1;
54 prev[i] = NULL;
55 }
56
57 ncost = 0;
58 have_node_costs = dglGet_NodeAttrSize(graph);
59
60 dglHeapInit(&heap);
61
62 for (i = 0; i < from->n_values; i++) {
63 int v = from->value[i];
64
65 if (dst[v] == 0)
66 continue; /* ignore duplicates */
67 dst[v] = 0; /* make sure all from nodes are processed first */
68 dglHeapData_u heap_data;
69
70 heap_data.ul = v;
71 dglHeapInsertMin(&heap, 0, ' ', heap_data);
72 }
73 while (1) {
74 dglInt32_t v, dist;
75 dglHeapNode_s heap_node;
76 dglHeapData_u heap_data;
77 dglInt32_t *edge;
78 dglInt32_t *node;
79
80 if (!dglHeapExtractMin(&heap, &heap_node))
81 break;
82 v = heap_node.value.ul;
83 dist = heap_node.key;
84 if (dst[v] < dist)
85 continue;
86
87 node = dglGetNode(graph, v);
88
89 if (have_node_costs && prev[v]) {
90 memcpy(&ncost, dglNodeGet_Attr(graph, node),
91 sizeof(ncost));
92 if (ncost > 0)
93 dist += ncost;
94 /* do not go through closed nodes */
95 if (ncost < 0)
96 continue;
97 }
98
99 dglEdgeset_T_Initialize(&et, graph,
100 dglNodeGet_OutEdgeset(graph, node));
101
102 for (edge = dglEdgeset_T_First(&et); edge;
103 edge = dglEdgeset_T_Next(&et)) {
104 dglInt32_t *to = dglEdgeGet_Tail(graph, edge);
105 dglInt32_t to_id = dglNodeGet_Id(graph, to);
106 dglInt32_t d = dglEdgeGet_Cost(graph, edge);
107
108 if (dst[to_id] < 0 || dst[to_id] > dist + d) {
109 dst[to_id] = dist + d;
110 prev[to_id] = edge;
111 heap_data.ul = to_id;
112 dglHeapInsertMin(&heap, dist + d, ' ', heap_data);
113 }
114 }
115
116 dglEdgeset_T_Release(&et);
117 }
118
119 dglHeapFree(&heap, NULL);
120
121 return 0;
122 }
123
124 /*!
125 \brief Computes shortest paths from every node to nodes in "to".
126
127 Array "dst" contains the cost of the path or -1 if the node is not
128 reachable. Nxt contains edges from successor along the shortest
129 path. This method does reverse search starting with "to" nodes and
130 going backward.
131
132 \param graph input graph
133 \param to list of 'to' positions
134 \param[out] dst array of costs to reach nodes
135 \param[out] nxt array of edges from successor along the shortest path
136
137 \return 0 on success
138 \return -1 on failure
139 */
NetA_distance_to_points(dglGraph_s * graph,struct ilist * to,int * dst,dglInt32_t ** nxt)140 int NetA_distance_to_points(dglGraph_s *graph, struct ilist *to,
141 int *dst, dglInt32_t **nxt)
142 {
143 int i, nnodes;
144 dglHeap_s heap;
145 dglEdgesetTraverser_s et;
146 int have_node_costs;
147 dglInt32_t ncost;
148
149 nnodes = dglGet_NodeCount(graph);
150
151 /* initialize costs and edge list */
152 for (i = 1; i <= nnodes; i++) {
153 dst[i] = -1;
154 nxt[i] = NULL;
155 }
156
157 if (graph->Version < 2) {
158 G_warning("Directed graph must be version 2 or 3 for NetA_distance_to_points()");
159 return -1;
160 }
161
162 ncost = 0;
163 have_node_costs = dglGet_NodeAttrSize(graph);
164
165 dglHeapInit(&heap);
166
167 for (i = 0; i < to->n_values; i++) {
168 int v = to->value[i];
169
170 if (dst[v] == 0)
171 continue; /* ignore duplicates */
172 dst[v] = 0; /* make sure all to nodes are processed first */
173 dglHeapData_u heap_data;
174
175 heap_data.ul = v;
176 dglHeapInsertMin(&heap, 0, ' ', heap_data);
177 }
178 while (1) {
179 dglInt32_t v, dist;
180 dglHeapNode_s heap_node;
181 dglHeapData_u heap_data;
182 dglInt32_t *edge;
183 dglInt32_t *node;
184
185 if (!dglHeapExtractMin(&heap, &heap_node))
186 break;
187 v = heap_node.value.ul;
188 dist = heap_node.key;
189 if (dst[v] < dist)
190 continue;
191
192 node = dglGetNode(graph, v);
193
194 if (have_node_costs && nxt[v]) {
195 memcpy(&ncost, dglNodeGet_Attr(graph, node),
196 sizeof(ncost));
197 if (ncost > 0)
198 dist += ncost;
199 /* do not go through closed nodes */
200 if (ncost < 0)
201 continue;
202 }
203
204 dglEdgeset_T_Initialize(&et, graph,
205 dglNodeGet_InEdgeset(graph, node));
206
207 for (edge = dglEdgeset_T_First(&et); edge;
208 edge = dglEdgeset_T_Next(&et)) {
209 dglInt32_t *from = dglEdgeGet_Head(graph, edge);
210 dglInt32_t from_id = dglNodeGet_Id(graph, from);
211 dglInt32_t d = dglEdgeGet_Cost(graph, edge);
212
213 if (dst[from_id] < 0 || dst[from_id] > dist + d) {
214 dst[from_id] = dist + d;
215 nxt[from_id] = edge;
216 heap_data.ul = from_id;
217 dglHeapInsertMin(&heap, dist + d, ' ', heap_data);
218 }
219 }
220
221 dglEdgeset_T_Release(&et);
222 }
223
224 dglHeapFree(&heap, NULL);
225
226 return 0;
227 }
228
229 /*!
230 \brief Find a path (minimum number of edges) from 'from' to 'to'
231 using only edges flagged as valid in 'edges'. Edge costs are not
232 considered. Closed nodes are not traversed.
233
234 Precisely, edge with id I is used if edges[abs(i)] == 1. List
235 stores the indices of lines on the path. The method returns the
236 number of edges or -1 if no path exists.
237
238 \param graph input graph
239 \param from 'from' position
240 \param to 'to' position
241 \param edges array of edges indicating whether an edge should be used
242 \param[out] list list of edges
243
244 \return number of edges
245 \return -1 on failure
246 */
NetA_find_path(dglGraph_s * graph,int from,int to,int * edges,struct ilist * list)247 int NetA_find_path(dglGraph_s * graph, int from, int to, int *edges,
248 struct ilist *list)
249 {
250 dglInt32_t **prev, *queue;
251 dglEdgesetTraverser_s et;
252 char *vis;
253 int begin, end, cur, nnodes;
254 int have_node_costs;
255 dglInt32_t ncost;
256
257 nnodes = dglGet_NodeCount(graph);
258 prev = (dglInt32_t **) G_calloc(nnodes + 1, sizeof(dglInt32_t *));
259 queue = (dglInt32_t *) G_calloc(nnodes + 1, sizeof(dglInt32_t));
260 vis = (char *)G_calloc(nnodes + 1, sizeof(char));
261 if (!prev || !queue || !vis) {
262 G_fatal_error(_("Out of memory"));
263 return -1;
264 }
265 Vect_reset_list(list);
266
267 ncost = 0;
268 have_node_costs = dglGet_NodeAttrSize(graph);
269
270 begin = 0;
271 end = 1;
272 vis[from] = 'y';
273 queue[0] = from;
274 prev[from] = NULL;
275 while (begin != end) {
276 dglInt32_t vertex = queue[begin++];
277 dglInt32_t *edge, *node;
278
279 if (vertex == to)
280 break;
281
282 /* do not go through closed nodes */
283 if (have_node_costs && prev[vertex]) {
284 memcpy(&ncost, dglNodeGet_Attr(graph, dglEdgeGet_Tail(graph, edge)),
285 sizeof(ncost));
286 if (ncost < 0)
287 continue;
288 }
289
290 node = dglGetNode(graph, vertex);
291
292 dglEdgeset_T_Initialize(&et, graph,
293 dglNodeGet_OutEdgeset(graph, node));
294 for (edge = dglEdgeset_T_First(&et); edge;
295 edge = dglEdgeset_T_Next(&et)) {
296 dglInt32_t edge_id = abs(dglEdgeGet_Id(graph, edge));
297 dglInt32_t node_id =
298 dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
299 if (edges[edge_id] && !vis[node_id]) {
300 vis[node_id] = 'y';
301 prev[node_id] = edge;
302 queue[end++] = node_id;
303 }
304 }
305 dglEdgeset_T_Release(&et);
306 }
307 G_free(queue);
308 if (!vis[to]) {
309 G_free(prev);
310 G_free(vis);
311 return -1;
312 }
313
314 cur = to;
315 while (prev[cur] != NULL) {
316 Vect_list_append(list, abs(dglEdgeGet_Id(graph, prev[cur])));
317 cur = dglNodeGet_Id(graph, dglEdgeGet_Head(graph, prev[cur]));
318 }
319
320 G_free(prev);
321 G_free(vis);
322 return list->n_values;
323 }
324