GRASS Programmer's Manual
6.4.2(2012)
|
00001 00016 #include <stdio.h> 00017 #include <stdlib.h> 00018 #include <grass/gis.h> 00019 #include <grass/Vect.h> 00020 #include <grass/glocale.h> 00021 #include <grass/dgl/graph.h> 00022 00023 struct union_find 00024 { 00025 int *parent; 00026 }; 00027 00028 static int uf_initialize(struct union_find *uf, int size) 00029 { 00030 int i; 00031 uf->parent = (int *)G_calloc(size, sizeof(int)); 00032 if (!uf->parent) 00033 return 0; 00034 for (i = 0; i < size; i++) 00035 uf->parent[i] = i; 00036 return 1; 00037 } 00038 00039 static void uf_release(struct union_find *uf) 00040 { 00041 G_free(uf->parent); 00042 } 00043 00044 static int uf_find(struct union_find *uf, int v) 00045 { 00046 int cur = v, tmp; 00047 00048 while (uf->parent[cur] != cur) 00049 cur = uf->parent[cur]; 00050 while (uf->parent[v] != v) { 00051 tmp = uf->parent[v]; 00052 uf->parent[v] = cur; 00053 v = tmp; 00054 } 00055 return cur; 00056 } 00057 00058 /*TODO: union by rank */ 00059 static void uf_union(struct union_find *uf, int u, int v) 00060 { 00061 int parent_u = uf_find(uf, u); 00062 int parent_v = uf_find(uf, v); 00063 00064 if (parent_u != parent_v) 00065 uf->parent[parent_u] = parent_v; 00066 } 00067 00068 typedef struct 00069 { 00070 dglInt32_t cost; 00071 dglInt32_t *edge; 00072 } edge_cost_pair; 00073 00074 static int cmp_edge(const void *pa, const void *pb) 00075 { 00076 return ((edge_cost_pair *) pa)->cost - ((edge_cost_pair *) pb)->cost; 00077 } 00078 00088 int NetA_spanning_tree(dglGraph_s * graph, struct ilist *tree_list) 00089 { 00090 int nnodes, edges, nedges, i, index; 00091 edge_cost_pair *perm; /*permutaion of edges in ascending order */ 00092 struct union_find uf; 00093 dglEdgesetTraverser_s et; 00094 00095 nnodes = dglGet_NodeCount(graph); 00096 nedges = dglGet_EdgeCount(graph); 00097 perm = (edge_cost_pair *) G_calloc(nedges, sizeof(edge_cost_pair)); 00098 if (!perm || !uf_initialize(&uf, nnodes + 1)) { 00099 G_fatal_error(_("Out of memory")); 00100 return -1; 00101 } 00102 /*for some obscure reasons, dglGetEdge always returns NULL. Therefore this complicated enumeration of the edges... */ 00103 index = 0; 00104 G_message(_("Computing minimum spanning tree...")); 00105 G_percent_reset(); 00106 for (i = 1; i <= nnodes; i++) { 00107 G_percent(i, nnodes + nedges, 1); 00108 dglInt32_t *edge; 00109 00110 dglEdgeset_T_Initialize(&et, graph, 00111 dglNodeGet_OutEdgeset(graph, 00112 dglGetNode(graph, 00113 (dglInt32_t) 00114 i))); 00115 for (edge = dglEdgeset_T_First(&et); edge; 00116 edge = dglEdgeset_T_Next(&et)) 00117 if (dglEdgeGet_Id(graph, edge) > 0) { 00118 perm[index].edge = edge; 00119 perm[index].cost = dglEdgeGet_Cost(graph, edge); 00120 index++; 00121 } 00122 00123 dglEdgeset_T_Release(&et); 00124 } 00125 edges = 0; 00126 qsort((void *)perm, index, sizeof(edge_cost_pair), cmp_edge); 00127 for (i = 0; i < index; i++) { 00128 G_percent(i + nnodes, nnodes + nedges, 1); 00129 dglInt32_t head = 00130 dglNodeGet_Id(graph, dglEdgeGet_Head(graph, perm[i].edge)); 00131 dglInt32_t tail = 00132 dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, perm[i].edge)); 00133 if (uf_find(&uf, head) != uf_find(&uf, tail)) { 00134 uf_union(&uf, head, tail); 00135 edges++; 00136 Vect_list_append(tree_list, dglEdgeGet_Id(graph, perm[i].edge)); 00137 } 00138 } 00139 G_free(perm); 00140 uf_release(&uf); 00141 return edges; 00142 }