plus_area.c

Go to the documentation of this file.

/**
 * \file plus_area.c
 *
 * \brief Vector library - update topo for areas (lower level functions)
 *
 * Lower level functions for reading/writing/manipulating vectors.
 *
 * This program is free software under the GNU General Public License
 * (>=v2). Read the file COPYING that comes with GRASS for details.
 *
 * \author CERL (probably Dave Gerdes), Radim Blazek
 *
 * \date 2001-2006
 */

#include <stdlib.h>
#include <grass/vector.h>
#include <grass/glocale.h>

static int debug_level = -1;

/*!
 * \brief Build topo for area from lines
 *
 * Area is built in clockwise order.
 * Take a given line and start off to the RIGHT/LEFT and try to complete
 * an area. 
 * 
 * Possible Scenarios:
 *  - I.    path runs into first line.                              : AREA!
 *  - II.   path runs into a dead end (no other area lines at node) : no area
 *  - III.  path runs into a previous line that is not 1st line or to 1st line but not to start node : no area
 *
 * After we find an area then we call point_in_area() to see if the
 * specified point is w/in the area
 *
 * Old returns  -1:  error   0:  no area    (1:  point in area)
 *              -2: island  !!
 *
 * \param[in] plus pointer to Plus_head structure
 * \param[in] first_line line id of first line
 * \param[in] side side of line to build area on (GV_LEFT | GV_RIGHT)
 * \param[in] lines pointer to array of lines
 *
 * \return  -1 on error   
 * \return   0 no area
 * \return   number of lines
 */
int
dig_build_area_with_line(struct Plus_head *plus, plus_t first_line, int side,
                         plus_t ** lines)
{
    register int i;
    int prev_line, next_line;
    static plus_t *array;
    char *p;
    static int array_size;      /* 0 on startup */
    int n_lines;
    struct P_line *Line;
    struct P_topo_b *topo;
    int node;

    if (debug_level == -1) {
        const char *dstr = G_getenv_nofatal("DEBUG");

        if (dstr != NULL)
            debug_level = atoi(dstr);
        else
            debug_level = 0;
    }

    G_debug(3, "dig_build_area_with_line(): first_line = %d, side = %d",
            first_line, side);

    /* First check if line is not degenerated (degenerated lines have angle -9) 
     *  Following degenerated lines are skip by dig_angle_next_line() */
    Line = plus->Line[first_line];
    if (Line->type != GV_BOUNDARY)
        return -1;

    topo = (struct P_topo_b *)Line->topo;
    node = topo->N1;            /* to check one is enough, because if degenerated N1 == N2 */
    if (dig_node_line_angle(plus, node, first_line) == -9.) {
        G_debug(3, "First line degenerated");
        return (0);
    }

    if (array_size == 0) {      /* first time */
        array_size = 1000;
        array = (plus_t *) dig__falloc(array_size, sizeof(plus_t));
        if (array == NULL)
            return (dig_out_of_memory());
    }

    if (side == GV_LEFT) {
        first_line = -first_line;       /* start at node1, reverse direction */
    }
    array[0] = first_line;
    prev_line = -first_line;    /* start at node2 for direct and node1 for
                                   reverse direction */
    /* angle of first line */
    n_lines = 1;
    while (1) {
        next_line =
            dig_angle_next_line(plus, prev_line, GV_RIGHT, GV_BOUNDARY, NULL);
        G_debug(3, "next_line = %d", next_line);

        if (next_line == 0) {
            G_debug(3, "Cannot build area, no next line for line %d",
                    prev_line);
            return (-1);        /* Not found */
        }

        /* Check if adjacent lines do not have the same angle */
        if (!dig_node_angle_check(plus, next_line, GV_BOUNDARY)) {
            G_debug(3,
                    "Cannot build area, a neighbour of the line %d has the same angle at the node",
                    next_line);
            return 0;
        }

        /*  I. Area closed. This also handles the problem w/ 1 single area line */
        if (first_line == next_line) {
            /* GOT ONE!  fill area struct  and return */
            G_debug(3, "Got one! :");

            /* avoid loop when not debugging */
            if (debug_level > 2) {
                for (i = 0; i < n_lines; i++) {
                    G_debug(3, " area line (%d) = %d", i, array[i]);
                }
            }

            *lines = array;
            return (n_lines);
        }

        /* II. Note this is a dead end */
        /* ( if prev_line != -first_line so it goes after the previous test) ? */
        if (prev_line == next_line) {
            G_debug(3, "Dead_end:");
            return (0);         /* dead end */
        }

        /* III. Unclosed ?, I would say started from free end */
        for (i = 0; i < n_lines; i++)
            if (abs(next_line) == abs(array[i])) {
                G_debug(3, "Unclosed area:");
                return (0);     /* ran into a different area */
            }

        /* otherwise keep going */
        if (n_lines >= array_size) {
            p = dig__frealloc(array, array_size + 100, sizeof(plus_t),
                              array_size);
            if (p == NULL)
                return (dig_out_of_memory());
            array = (plus_t *) p;
            array_size += 100;
        }
        array[n_lines++] = next_line;
        prev_line = -next_line;
    }

    return 0;
}

/*!
 * \brief Allocate space for new area and create boundary info from array.
 *
 * Then for each line in area, update line (right,left) info.
 *
 * Neither islands nor centroids are filled.
 *
 * \param[in] plus pointer to Plus_head structure
 * \param[in] n_lines number of lines
 * \param[in] lines array of lines, negative for reverse direction
 * \param[in] box bounding box
 *
 * \return number of new area
 * \return -1 on error
 */
int dig_add_area(struct Plus_head *plus, int n_lines, plus_t * lines,
                 struct bound_box *box)
{
    register int i;
    register int area, line;
    struct P_area *Area;
    struct P_line *Line;
    struct P_topo_b *topo;

    G_debug(3, "dig_add_area():");
    /* First look if we have space in array of pointers to areas
     *  and reallocate if necessary */
    if (plus->n_areas >= plus->alloc_areas) {   /* array is full */
        if (dig_alloc_areas(plus, 1000) == -1)
            return -1;
    }

    /* allocate area structure */
    area = plus->n_areas + 1;
    G_debug(3, "    new area = %d", area);
    Area = dig_alloc_area();
    if (Area == NULL)
        return -1;

    if (dig_area_alloc_line(Area, n_lines) == -1)
        return -1;

    for (i = 0; i < n_lines; i++) {
        line = lines[i];
        Area->lines[i] = line;
        Line = plus->Line[abs(line)];
        topo = (struct P_topo_b *)Line->topo;
        if (line < 0) {         /* reverse direction -> area on left */
            if (topo->left != 0) {
                G_warning(_("Line %d already has area/isle %d to left"), line,
                          topo->left);
                return -1;
            }

            G_debug(3, "  Line %d left set to %d.", line, area);
            topo->left = area;
        }
        else {
            if (topo->right != 0) {
                G_warning(_("Line %d already has area/isle %d to right"),
                          line, topo->right);
                return -1;
            }

            G_debug(3, "  Line %d right set to %d.", line, area);
            topo->right = area;
        }
    }
    Area->n_lines = n_lines;
    Area->centroid = 0;

    plus->Area[area] = Area;

    dig_spidx_add_area(plus, area, box);

    plus->n_areas++;

    return (area);
}

/*!
 * \brief Add isle to area if does not exist yet.
 *
 * \param[in] plus pointer to Plus_head structure
 * \param[in] area area id
 * \param[in] isle isle id
 *
 * \return 0
 */
int dig_area_add_isle(struct Plus_head *plus, int area, int isle)
{
    int i;
    struct P_area *Area;

    G_debug(3, "dig_area_add_isle(): area = %d isle = %d", area, isle);

    if (debug_level == -1) {
        const char *dstr = G_getenv_nofatal("DEBUG");

        if (dstr != NULL)
            debug_level = atoi(dstr);
        else
            debug_level = 0;
    }

    Area = plus->Area[area];
    if (Area == NULL)
        G_fatal_error("Attempt to add isle to dead area");

    if (debug_level > 0) {
        for (i = 0; i < Area->n_isles; i++) {
            if (Area->isles[i] == isle) {
                /* Already exists: bug in vector libs */
                G_warning(_("Isle already registered in area"));
                return 0;
            }
        }
    }

    if (Area->alloc_isles <= Area->n_isles)     /* array is full */
        dig_area_alloc_isle(Area, 1);

    Area->isles[Area->n_isles] = isle;
    Area->n_isles++;
    G_debug(3, "  -> n_isles = %d", Area->n_isles);

    return 0;
}

/*!
 * \brief Delete isle from area.
 *
 * \param[in] plus pointer to Plus_head structure
 * \param[in] area area id
 * \param[in] isle isle id
 *
 * \return 0
 */
int dig_area_del_isle(struct Plus_head *plus, int area, int isle)
{
    int i;
    struct P_area *Area;

    G_debug(3, "dig_area_del_isle(): area = %d isle = %d", area, isle);

    Area = plus->Area[area];
    if (Area == NULL)
        G_fatal_error(_("Attempt to delete isle from dead area"));

    /* find index of the isle */
    i = 0;
    while (i < Area->n_isles && Area->isles[i] != isle)
        i++;

    if (i == Area->n_isles) {
        G_fatal_error(_("Attempt to delete not registered isle %d from area %d"),
                      isle, area);
    }

    i++;
    while (i < Area->n_isles) {
        Area->isles[i - 1] = Area->isles[i];
        i++;
    }

    Area->n_isles--;

    return 0;
}

/*!
 * \brief Delete area from Plus_head structure
 *
 *  This function deletes area from the topo structure and resets references
 *  to this area in boundaries, isles (within), and the centroid (if any) to 0.
 *  Possible new area is not created by this function, so that
 *  old boundaries participating in this area are left without area information
 *  even if form new area.
 *  Not enabled now: If area is inside other area, area info for islands within 
 *                   deleted area is reset to that area outside.
 *  (currently area info of isles is set to 0)
 *
 * \param[in] plus pointer to Plus_head structure
 * \param[in] area area id
 *
 * \return 0 on error
 * \return 1 on success
 */
int dig_del_area(struct Plus_head *plus, int area)
{
    int i, line;

    struct P_area *Area;
    struct P_line *Line;
    struct P_isle *Isle;
    struct P_topo_b *btopo;
    struct P_topo_c *ctopo;
    
    G_debug(3, "dig_del_area() area =  %d", area);
    Area = plus->Area[area];

    if (Area == NULL) {
        G_warning(_("Attempt to delete dead area"));
        return 0;
    }

    dig_spidx_del_area(plus, area);

    /* Set area for all lines to 0 */
    /* isle = 0; */
    for (i = 0; i < Area->n_lines; i++) {
        line = Area->lines[i];  /* >0 = clockwise -> right, <0 = counterclockwise ->left */
        Line = plus->Line[abs(line)];
        btopo = (struct P_topo_b *)Line->topo;
        if (line > 0) {
            G_debug(3, "  Set line %d right side to 0", line);
            btopo->right = 0;
        }
        else {
            G_debug(3, "  Set line %d left side to 0", line);
            btopo->left = 0;
        }

        /* Find the isle this area is part of (used late below) */
        /*
           if ( line > 0 ) {
           if ( Line->left < 0 ) isle = Line->left; 
           } else {
           if ( Line->right < 0 ) isle = Line->right; 
           }
         */
    }

    /* Unset area information of centroid */
    /* TODO: duplicate centroids have also area information ->
     *        1) do not save such info
     *        2) find all by box and reset info */
    line = Area->centroid;
    if (line > 0) {
        Line = plus->Line[line];
        if (!Line) {
            G_warning(_("Dead centroid %d registered for area (bug in the vector library)"),
                      line);
        }
        else {
            ctopo = (struct P_topo_c *)Line->topo;
            ctopo->area = 0;
        }
    }

    /* Find the area this area is within */
    /*
       area_out = 0;
       if ( isle > 0 ) { 
       Isle =  plus->Isle[abs(isle)];
       area_out = Isle->area;
       }
     */

    /* Reset information about area outside for isles within this area */
    G_debug(3, "  n_isles = %d", Area->n_isles);
    for (i = 0; i < Area->n_isles; i++) {
        Isle = plus->Isle[Area->isles[i]];
        if (Isle == NULL) {
            G_fatal_error(_("Attempt to delete area %d info from dead isle %d"),
                          area, Area->isles[i]);
        }
        else {
            /* Isle->area = area_out; */
            Isle->area = 0;
        }
    }

    /* free structures */
    dig_free_area(Area);
    plus->Area[area] = NULL;
    return 1;
}


/*!
 * \brief Find line number of next angle to follow a line
 *
 * Assume that lines are sorted in increasing angle order and angles
 * of points and degenerated lines are set to -9 (ignored).
 *
 * \param[in] plus pointer to Plus_head structure
 * \param[in] current_line current line id, negative if request for end node
 * \param[in] side side GV_RIGHT or GV_LEFT
 * \param[in] type line type (GV_LINE, GV_BOUNDARY or both)
 * \param[in] angle
 *
 * \return line number of next angle to follow a line (negative if connected by end node)
 *               (number of current line may be returned if dangle - this is used in build)
 * \return 0 on error or not found
 */
int
dig_angle_next_line(struct Plus_head *plus, plus_t current_line, int side,
                    int type, float *angle)
{
    int next;
    int line;
    plus_t node;
    struct P_node *Node;
    struct P_line *Line;

    if (debug_level == -1) {
        const char *dstr = G_getenv_nofatal("DEBUG");

        if (dstr != NULL)
            debug_level = atoi(dstr);
        else
            debug_level = 0;
    }

    G_debug(3, "dig__angle_next_line: line = %d, side = %d, type = %d",
            current_line, side, type);

    Line = plus->Line[abs(current_line)];
    
    if (!(Line->type & GV_LINES)) {
        if (angle)
            *angle = -9.;
        return 0;
    }

    node = 0;
    if (current_line > 0) {
        if (Line->type == GV_LINE) {
            struct P_topo_l *topo = (struct P_topo_l *)Line->topo;
            node = topo->N1;
        }
        else if (Line->type == GV_BOUNDARY) {
            struct P_topo_b *topo = (struct P_topo_b *)Line->topo;
            node = topo->N1;
        }
    }
    else {
        if (Line->type == GV_LINE) {
            struct P_topo_l *topo = (struct P_topo_l *)Line->topo;
            node = topo->N2;
        }
        else if (Line->type == GV_BOUNDARY) {
            struct P_topo_b *topo = (struct P_topo_b *)Line->topo;
            node = topo->N2;
        }
    }

    G_debug(3, " node = %d", node);

    Node = plus->Node[node];
    G_debug(3, "  n_lines = %d", Node->n_lines);
    /* avoid loop when not debugging */
    if (debug_level > 2) {
        int i;

        for (i = 0; i < Node->n_lines; i++) {
            G_debug(3, "  i = %d line = %d angle = %f", i, Node->lines[i],
                    Node->angles[i]);
        }
    }

    /* first find index for that line */
    next = Node->n_lines - 1;
    while (next >= 0 && Node->lines[next] != current_line) {
        next--;
    }

    if (next == -1) {
        /* internal error, should not happen */
        G_fatal_error("dig_angle_next_line(): line %d not found at its own node %d",
                      current_line, node);
        if (angle)
            *angle = -9.;
        return 0;               /* not found */
    }

    G_debug(3, "  current position = %d", next);
    while (1) {
        if (side == GV_RIGHT) { /* go up (greater angle) */
            if (next == Node->n_lines - 1)
                next = 0;
            else
                next++;
        }
        else {                  /* go down (smaller angle) */
            if (next == 0)
                next = Node->n_lines - 1;
            else
                next--;
        }
        G_debug(3, "  next = %d line = %d angle = %f", next,
                Node->lines[next], Node->angles[next]);

        if (Node->angles[next] == -9.) {        /* skip points and degenerated */
            G_debug(3, "  point/degenerated -> skip");
            if (Node->lines[next] == current_line)
                break;          /* Yes, that may happen if input line is degenerated and isolated and this breaks loop */
            else
                continue;
        }

        line = Node->lines[next];
        Line = plus->Line[abs(line)];
        
        if (Line->type & type) {        /* line found */
            G_debug(3, "  this one");
            if (angle)
                *angle = Node->angles[next];
            return line;
        }

        /* input line reached, this must be last, because current_line may be correct return value (dangle) */
        if (line == current_line)
            break;
    }
    G_debug(3, "  No next line for line %d at node %d",
            current_line, (int)node);
    if (angle)
        *angle = -9.;

    return 0;
}

/*!
 * \brief Check if angles of adjacent lines differ.
 *
 * Negative line number for end point. Assume that lines are sorted
 * in increasing angle order and angles of points and degenerated
 * lines are set to 9 (ignored).
 *
 * \param[in] plus pointer to Plus_head structure
 * \param[in] line current line id, negative if request for node 2
 * \param[in] type line type (GV_LINE, GV_BOUNDARY or both)
 *
 * \return 1 angles differ
 * \return 0 angle of a line up or down is identical
 */
int dig_node_angle_check(struct Plus_head *plus, plus_t line, int type)
{
    int next, prev;
    float angle1, angle2;
    plus_t node = 0;
    struct P_line *Line;

    G_debug(3, "dig_node_angle_check: line = %d, type = %d", line, type);

    Line = plus->Line[abs(line)];
    if (!(Line->type & GV_LINES))
        return 0;

    if (line > 0) {
        if (Line->type == GV_LINE) {
            struct P_topo_l *topo = (struct P_topo_l *)Line->topo;
            node = topo->N1;
        }
        else if (Line->type == GV_BOUNDARY) {
            struct P_topo_b *topo = (struct P_topo_b *)Line->topo;
            node = topo->N1;
        }
    }
    else {
        if (Line->type == GV_LINE) {
            struct P_topo_l *topo = (struct P_topo_l *)Line->topo;
            node = topo->N2;
        }
        else if (Line->type == GV_BOUNDARY) {
            struct P_topo_b *topo = (struct P_topo_b *)Line->topo;
            node = topo->N2;
        }
    }

    angle1 = dig_node_line_angle(plus, node, line);

    /* Next */
    next = dig_angle_next_line(plus, line, GV_RIGHT, type, &angle2);
    /* angle2 = dig_node_line_angle(plus, node, next); */
    if (angle1 == angle2) {
        G_debug(3,
                "  The line to the right has the same angle: node = %d, line = %d",
                node, next);
        return 0;
    }

    /* Previous */
    prev = dig_angle_next_line(plus, line, GV_LEFT, type, &angle2);
    /* angle2 = dig_node_line_angle(plus, node, prev); */
    if (angle1 == angle2) {
        G_debug(3,
                "  The line to the left has the same angle: node = %d, line = %d",
                node, next);
        return 0;
    }

    return 1;                   /* OK */
}

/*!
 * \brief Allocate space for new island and create boundary info from array.
 *
 * The order of input lines is expected to be counter clockwise.
 * Then for each line in isle, update line (right,left) info.
 *
 * Area number the island is within is not filled.
 *
 * \param[in] plus pointer to Plus_head structure
 * \param[in] n_lines number of lines
 * \param[in] lines array of lines, negative for reverse direction 
 * \param[in] box bounding box
 *
 * \return number of new isle
 * \return -1 on error
 */
int dig_add_isle(struct Plus_head *plus, int n_lines, plus_t * lines,
                 struct bound_box *box)
{
    register int i;
    register int isle, line;
    struct P_isle *Isle;
    struct P_line *Line;
    struct P_topo_b *topo;

    G_debug(3, "dig_add_isle():");
    /* First look if we have space in array of pointers to isles
     *  and reallocate if necessary */
    if (plus->n_isles >= plus->alloc_isles) {   /* array is full */
        if (dig_alloc_isles(plus, 1000) == -1)
            return -1;
    }

    /* allocate isle structure */
    isle = plus->n_isles + 1;
    Isle = dig_alloc_isle();
    if (Isle == NULL)
        return -1;

    if ((dig_isle_alloc_line(Isle, n_lines)) == -1)
        return -1;

    Isle->area = 0;

    for (i = 0; i < n_lines; i++) {
        line = lines[i];
        G_debug(3, " i = %d line = %d", i, line);
        Isle->lines[i] = line;
        Line = plus->Line[abs(line)];
        topo = (struct P_topo_b *)Line->topo;
        if (line < 0) {         /* revers direction -> isle on left */
            if (topo->left != 0) {
                G_warning(_("Line %d already has area/isle %d to left"), line,
                          topo->left);
                return -1;
            }
            topo->left = -isle;
        }
        else {
            if (topo->right != 0) {
                G_warning(_("Line %d already has area/isle %d to right"), line,
                          topo->right);
                return -1;
            }

            topo->right = -isle;
        }
    }

    Isle->n_lines = n_lines;

    plus->Isle[isle] = Isle;

    dig_spidx_add_isle(plus, isle, box);

    plus->n_isles++;

    return (isle);
}


/*!
 * \brief Delete island from Plus_head structure
 *
 * Reset references to it in lines and area outside.
 *
 * \param[in] plus pointer to Plus_head structure
 * \param[in] isle isle id
 *
 * \return 1
 */
int dig_del_isle(struct Plus_head *plus, int isle)
{
    int i, line;
    struct P_line *Line;
    struct P_isle *Isle;
    struct P_topo_b *topo;

    G_debug(3, "dig_del_isle() isle =  %d", isle);
    Isle = plus->Isle[isle];

    dig_spidx_del_isle(plus, isle);

    /* Set area for all lines to 0 */
    for (i = 0; i < Isle->n_lines; i++) {
        line = Isle->lines[i];  /* >0 = clockwise -> right, <0 = counterclockwise ->left */
        Line = plus->Line[abs(line)];
        topo = (struct P_topo_b *)Line->topo;
        if (line > 0)
            topo->right = 0;
        else
            topo->left = 0;
    }

    /* Delete reference from area it is within */
    G_debug(3, "  area outside isle = %d", Isle->area);
    if (Isle->area > 0) {
        if (plus->Area[Isle->area] == NULL) {
            G_fatal_error(_("Attempt to delete isle %d info from dead area %d"),
                          isle, Isle->area);
        }
        else {
            dig_area_del_isle(plus, Isle->area, isle);
        }
    }

    /* free structures */
    dig_free_isle(Isle);
    plus->Isle[isle] = NULL;

    return 1;
}