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 00032 int NetA_articulation_points(dglGraph_s * graph, 00033 struct ilist *articulation_list) 00034 { 00035 int nnodes; 00036 int points = 0; 00037 00038 dglEdgesetTraverser_s *current; /*edge to be processed when the node is visited */ 00039 int *tin, *min_tin; /*time in, and smallest tin over all successors. 0 if not yet visited */ 00040 dglInt32_t **parent; /*parents of the nodes */ 00041 dglInt32_t **stack; /*stack of nodes */ 00042 dglInt32_t **current_edge; /*current edge for each node */ 00043 int *mark; /*marked articulation points */ 00044 dglNodeTraverser_s nt; 00045 dglInt32_t *current_node; 00046 int stack_size; 00047 int i, time; 00048 00049 nnodes = dglGet_NodeCount(graph); 00050 current = 00051 (dglEdgesetTraverser_s *) G_calloc(nnodes + 1, 00052 sizeof(dglEdgesetTraverser_s)); 00053 tin = (int *)G_calloc(nnodes + 1, sizeof(int)); 00054 min_tin = (int *)G_calloc(nnodes + 1, sizeof(int)); 00055 parent = (dglInt32_t **) G_calloc(nnodes + 1, sizeof(dglInt32_t *)); 00056 stack = (dglInt32_t **) G_calloc(nnodes + 1, sizeof(dglInt32_t *)); 00057 current_edge = (dglInt32_t **) G_calloc(nnodes + 1, sizeof(dglInt32_t *)); 00058 mark = (int *)G_calloc(nnodes + 1, sizeof(int)); 00059 if (!tin || !min_tin || !parent || !stack || !current || !mark) { 00060 G_fatal_error(_("Out of memory")); 00061 return -1; 00062 } 00063 00064 for (i = 1; i <= nnodes; i++) { 00065 dglEdgeset_T_Initialize(¤t[i], graph, 00066 dglNodeGet_OutEdgeset(graph, 00067 dglGetNode(graph, i))); 00068 current_edge[i] = dglEdgeset_T_First(¤t[i]); 00069 tin[i] = mark[i] = 0; 00070 } 00071 00072 dglNode_T_Initialize(&nt, graph); 00073 00074 time = 0; 00075 for (current_node = dglNode_T_First(&nt); current_node; 00076 current_node = dglNode_T_Next(&nt)) { 00077 dglInt32_t current_id = dglNodeGet_Id(graph, current_node); 00078 00079 if (tin[current_id] == 0) { 00080 int children = 0; /*number of subtrees rooted at the root/current_node */ 00081 00082 stack[0] = current_node; 00083 stack_size = 1; 00084 parent[current_id] = NULL; 00085 while (stack_size) { 00086 dglInt32_t *node = stack[stack_size - 1]; 00087 dglInt32_t node_id = dglNodeGet_Id(graph, node); 00088 00089 if (tin[node_id] == 0) /*vertex visited for the first time */ 00090 min_tin[node_id] = tin[node_id] = ++time; 00091 else { /*return from the recursion */ 00092 dglInt32_t to = dglNodeGet_Id(graph, 00093 dglEdgeGet_Tail(graph, 00094 current_edge 00095 [node_id])); 00096 if (min_tin[to] >= tin[node_id]) /*no path from the subtree above the current node */ 00097 mark[node_id] = 1; /*so the current node must be an articulation point */ 00098 00099 if (min_tin[to] < min_tin[node_id]) 00100 min_tin[node_id] = min_tin[to]; 00101 current_edge[node_id] = dglEdgeset_T_Next(¤t[node_id]); /*proceed to the next edge */ 00102 } 00103 for (; current_edge[node_id]; current_edge[node_id] = dglEdgeset_T_Next(¤t[node_id])) { //try next edges 00104 dglInt32_t *to = 00105 dglEdgeGet_Tail(graph, current_edge[node_id]); 00106 if (to == parent[node_id]) 00107 continue; /*skip parrent */ 00108 int to_id = dglNodeGet_Id(graph, to); 00109 00110 if (tin[to_id]) { /*back edge, cannot be a bridge/articualtion point */ 00111 if (tin[to_id] < min_tin[node_id]) 00112 min_tin[node_id] = tin[to_id]; 00113 } 00114 else { /*forward edge */ 00115 if (node_id == current_id) 00116 children++; /*if root, increase number of children */ 00117 parent[to_id] = node; 00118 stack[stack_size++] = to; 00119 break; 00120 } 00121 } 00122 if (!current_edge[node_id]) 00123 stack_size--; /*current node completely processed */ 00124 } 00125 if (children > 1) 00126 mark[current_id] = 1; /*if the root has more than 1 subtrees rooted at it, then it is an 00127 * articulation point */ 00128 } 00129 } 00130 00131 for (i = 1; i <= nnodes; i++) 00132 if (mark[i]) { 00133 points++; 00134 Vect_list_append(articulation_list, i); 00135 } 00136 00137 dglNode_T_Release(&nt); 00138 for (i = 1; i <= nnodes; i++) 00139 dglEdgeset_T_Release(¤t[i]); 00140 00141 G_free(current); 00142 G_free(tin); 00143 G_free(min_tin); 00144 G_free(parent); 00145 G_free(stack); 00146 G_free(current_edge); 00147 return points; 00148 }