1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
|
/*!
\file vector/neta/spanningtree.c
\brief Network Analysis library - spanning tree
Computes minimum spanning tree in the network.
(C) 2009-2010 by Daniel Bundala, and the GRASS Development Team
This program is free software under the GNU General Public License
(>=v2). Read the file COPYING that comes with GRASS for details.
\author Daniel Bundala (Google Summer of Code 2009)
*/
#include <stdio.h>
#include <stdlib.h>
#include <grass/gis.h>
#include <grass/vector.h>
#include <grass/glocale.h>
#include <grass/dgl/graph.h>
struct union_find {
int *parent;
};
static int uf_initialize(struct union_find *uf, int size)
{
int i;
uf->parent = (int *)G_calloc(size, sizeof(int));
if (!uf->parent)
return 0;
for (i = 0; i < size; i++)
uf->parent[i] = i;
return 1;
}
static void uf_release(struct union_find *uf)
{
G_free(uf->parent);
}
static int uf_find(struct union_find *uf, int v)
{
int cur = v, tmp;
while (uf->parent[cur] != cur)
cur = uf->parent[cur];
while (uf->parent[v] != v) {
tmp = uf->parent[v];
uf->parent[v] = cur;
v = tmp;
}
return cur;
}
/*TODO: union by rank */
static void uf_union(struct union_find *uf, int u, int v)
{
int parent_u = uf_find(uf, u);
int parent_v = uf_find(uf, v);
if (parent_u != parent_v)
uf->parent[parent_u] = parent_v;
}
typedef struct {
dglInt32_t cost;
dglInt32_t *edge;
} edge_cost_pair;
static int cmp_edge(const void *pa, const void *pb)
{
if (((edge_cost_pair *)pa)->cost < ((edge_cost_pair *)pb)->cost)
return -1;
return (((edge_cost_pair *)pa)->cost > ((edge_cost_pair *)pb)->cost);
}
/*!
\brief Get number of edges in the spanning forest
\param graph input graph
\param[out] list of edges
\return number of edges
\return -1 on failure
*/
int NetA_spanning_tree(dglGraph_s *graph, struct ilist *tree_list)
{
int nnodes, edges, nedges, i, index;
edge_cost_pair *perm; /*permutation of edges in ascending order */
struct union_find uf;
dglEdgesetTraverser_s et;
/* TODO: consider closed nodes / node costs */
nnodes = dglGet_NodeCount(graph);
nedges = dglGet_EdgeCount(graph);
perm = (edge_cost_pair *)G_calloc(nedges, sizeof(edge_cost_pair));
if (!perm || !uf_initialize(&uf, nnodes + 1)) {
G_fatal_error(_("Out of memory"));
return -1;
}
/* dglGetEdge is only supported with graphs version > 1. Therefore this
* complicated enumeration of the edges... */
index = 0;
G_message(_("Computing minimum spanning tree..."));
G_percent_reset();
for (i = 1; i <= nnodes; i++) {
G_percent(i, nnodes + nedges, 1);
dglInt32_t *edge;
dglEdgeset_T_Initialize(
&et, graph,
dglNodeGet_OutEdgeset(graph, dglGetNode(graph, (dglInt32_t)i)));
for (edge = dglEdgeset_T_First(&et); edge;
edge = dglEdgeset_T_Next(&et))
if (dglEdgeGet_Id(graph, edge) > 0) {
perm[index].edge = edge;
perm[index].cost = dglEdgeGet_Cost(graph, edge);
index++;
}
dglEdgeset_T_Release(&et);
}
edges = 0;
qsort((void *)perm, index, sizeof(edge_cost_pair), cmp_edge);
for (i = 0; i < index; i++) {
G_percent(i + nnodes, nnodes + nedges, 1);
dglInt32_t head =
dglNodeGet_Id(graph, dglEdgeGet_Head(graph, perm[i].edge));
dglInt32_t tail =
dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, perm[i].edge));
if (uf_find(&uf, head) != uf_find(&uf, tail)) {
uf_union(&uf, head, tail);
edges++;
Vect_list_append(tree_list, dglEdgeGet_Id(graph, perm[i].edge));
}
}
G_percent(index, index, 1);
G_free(perm);
uf_release(&uf);
return edges;
}
|
