articulation_point.c

Go to the documentation of this file.

/*!
   \file vector/neta/articulation_point.c
 
   \brief Network Analysis library - connected components
 
   Computes network articulation points.
 
   (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>
 
/*!
   \brief Get number of articulation points in the graph
 
   \param graph input graph
   \param[out] articulation_list list of articulation points
 
   \return number of points
   \return -1 on error
 */
int NetA_articulation_points(dglGraph_s *graph, struct ilist *articulation_list)
{
    int nnodes;
    int points = 0;
 
    dglEdgesetTraverser_s
        *current;        /*edge to be processed when the node is visited */
    int *tin, *min_tin;  /*time in, and smallest tin over all successors. 0 if
                            not yet visited */
    dglInt32_t **parent; /*parents of the nodes */
    dglInt32_t **stack;  /*stack of nodes */
    dglInt32_t **current_edge; /*current edge for each node */
    int *mark;                 /*marked articulation points */
    dglNodeTraverser_s nt;
    dglInt32_t *current_node;
    int stack_size;
    int i, time;
 
    nnodes = dglGet_NodeCount(graph);
    current = (dglEdgesetTraverser_s *)G_calloc(nnodes + 1,
                                                sizeof(dglEdgesetTraverser_s));
    tin = (int *)G_calloc(nnodes + 1, sizeof(int));
    min_tin = (int *)G_calloc(nnodes + 1, sizeof(int));
    parent = (dglInt32_t **)G_calloc(nnodes + 1, sizeof(dglInt32_t *));
    stack = (dglInt32_t **)G_calloc(nnodes + 1, sizeof(dglInt32_t *));
    current_edge = (dglInt32_t **)G_calloc(nnodes + 1, sizeof(dglInt32_t *));
    mark = (int *)G_calloc(nnodes + 1, sizeof(int));
    if (!tin || !min_tin || !parent || !stack || !current || !mark) {
        G_fatal_error(_("Out of memory"));
        return -1;
    }
 
    for (i = 1; i <= nnodes; i++) {
        dglEdgeset_T_Initialize(
            &current[i], graph,
            dglNodeGet_OutEdgeset(graph, dglGetNode(graph, i)));
        current_edge[i] = dglEdgeset_T_First(&current[i]);
        tin[i] = mark[i] = 0;
    }
 
    dglNode_T_Initialize(&nt, graph);
 
    time = 0;
    for (current_node = dglNode_T_First(&nt); current_node;
         current_node = dglNode_T_Next(&nt)) {
        dglInt32_t current_id = dglNodeGet_Id(graph, current_node);
 
        if (tin[current_id] == 0) {
            int children =
                0; /*number of subtrees rooted at the root/current_node */
 
            stack[0] = current_node;
            stack_size = 1;
            parent[current_id] = NULL;
            while (stack_size) {
                dglInt32_t *node = stack[stack_size - 1];
                dglInt32_t node_id = dglNodeGet_Id(graph, node);
 
                if (tin[node_id] == 0) /*vertex visited for the first time */
                    min_tin[node_id] = tin[node_id] = ++time;
                else { /*return from the recursion */
                    dglInt32_t to = dglNodeGet_Id(
                        graph, dglEdgeGet_Tail(graph, current_edge[node_id]));
                    if (min_tin[to] >=
                        tin[node_id])      /*no path from the subtree above the
                                              current node */
                        mark[node_id] = 1; /*so the current node must be an
                                              articulation point */
 
                    if (min_tin[to] < min_tin[node_id])
                        min_tin[node_id] = min_tin[to];
                    current_edge[node_id] = dglEdgeset_T_Next(
                        &current[node_id]); /*proceed to the next edge */
                }
                /*try next edges */
                for (; current_edge[node_id];
                     current_edge[node_id] =
                         dglEdgeset_T_Next(&current[node_id])) {
                    dglInt32_t *to =
                        dglEdgeGet_Tail(graph, current_edge[node_id]);
                    if (to == parent[node_id])
                        continue; /*skip parent */
                    int to_id = dglNodeGet_Id(graph, to);
 
                    if (tin[to_id]) { /*back edge, cannot be a
                                         bridge/articualtion point */
                        if (tin[to_id] < min_tin[node_id])
                            min_tin[node_id] = tin[to_id];
                    }
                    else { /*forward edge */
                        if (node_id == current_id)
                            children++; /*if root, increase number of children
                                         */
                        parent[to_id] = node;
                        stack[stack_size++] = to;
                        break;
                    }
                }
                if (!current_edge[node_id])
                    stack_size--; /*current node completely processed */
            }
            if (children > 1)
                mark[current_id] =
                    1; /*if the root has more than 1 subtrees rooted at it, then
                        * it is an articulation point */
        }
    }
 
    for (i = 1; i <= nnodes; i++)
        if (mark[i]) {
            points++;
            Vect_list_append(articulation_list, i);
        }
 
    dglNode_T_Release(&nt);
    for (i = 1; i <= nnodes; i++)
        dglEdgeset_T_Release(&current[i]);
 
    G_free(current);
    G_free(tin);
    G_free(min_tin);
    G_free(parent);
    G_free(stack);
    G_free(current_edge);
    G_free(mark);
    return points;
}