clean_nodes.c

Go to the documentation of this file.

/*!
   \file lib/vector/Vlib/clean_nodes.c

   \brief Vector library - Clean boundaries at nodes

   Higher level functions for reading/writing/manipulating vectors.

   (C) 2001-2009 by 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 Radim Blazek
 */

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

/*!
   \brief Clean small angles at nodes.

   It may happen that even if the angle between 2 boundaries at node
   is very small, the calculated angle is 0 because of representation
   error. The map must be built at least on level GV_BUILD_BASE

   \param Map input map
   \param otype feature type
   \param[out] Err vector map where error line segments are written

   \return number of line modifications
 */
int
Vect_clean_small_angles_at_nodes(struct Map_info *Map, int otype,
                                 struct Map_info *Err)
{
    int node, nnodes;
    int nmodif = 0;
    struct line_pnts *Points;
    struct line_cats *SCats, *LCats, *OCats;

    Points = Vect_new_line_struct();
    SCats = Vect_new_cats_struct();
    LCats = Vect_new_cats_struct();
    OCats = Vect_new_cats_struct();

    nnodes = Vect_get_num_nodes(Map);
    for (node = 1; node <= Vect_get_num_nodes(Map); node++) {
        int i, nlines;

        if (node <= nnodes)
            G_percent(node, nnodes, 1);
        G_debug(3, "node = %d", node);
        if (!Vect_node_alive(Map, node))
            continue;

        while (1) {
            float angle1 = -100;
            int line1 = -999;   /* value not important, just for debug */
            int clean = 1;

            nlines = Vect_get_node_n_lines(Map, node);
            G_debug(3, "nlines = %d", nlines);

            for (i = 0; i < nlines; i++) {
                struct P_line *Line;
                int line2;
                float angle2;

                line2 = Vect_get_node_line(Map, node, i);
                Line = Map->plus.Line[abs(line2)];
                if (!Line)
                    continue;
                G_debug(4, "  type = %d", Line->type);
                if (!(Line->type & (otype & GV_LINES)))
                    continue;

                angle2 = Vect_get_node_line_angle(Map, node, i);
                if (angle2 == -9.0)
                    continue;   /* Degenerated line */

                G_debug(4, "  line1 = %d angle1 = %e line2 = %d angle2 = %e",
                        line1, angle1, line2, angle2);

                if (angle2 == angle1) {
                    int j;
                    double length1, length2;
                    int short_line;     /* line with shorter end segment */
                    int long_line;      /* line with longer end segment */
                    int new_short_line = 0;     /* line number of short line after rewrite */
                    int short_type, long_type, type;
                    double x, y, z, nx, ny, nz;

                    G_debug(4, "  identical angles -> clean");

                    /* Length of end segments for both lines */
                    Vect_read_line(Map, Points, NULL, abs(line1));
                    if (line1 > 0) {
                        length1 =
                            Vect_points_distance(Points->x[0], Points->y[0],
                                                 0.0, Points->x[1],
                                                 Points->y[1], 0.0, 0);
                    }
                    else {
                        int np;

                        np = Points->n_points;
                        length1 =
                            Vect_points_distance(Points->x[np - 1],
                                                 Points->y[np - 1], 0.0,
                                                 Points->x[np - 2],
                                                 Points->y[np - 2], 0.0, 0);
                    }

                    Vect_read_line(Map, Points, NULL, abs(line2));
                    if (line2 > 0) {
                        length2 =
                            Vect_points_distance(Points->x[0], Points->y[0],
                                                 0.0, Points->x[1],
                                                 Points->y[1], 0.0, 0);
                    }
                    else {
                        int np;

                        np = Points->n_points;
                        length2 =
                            Vect_points_distance(Points->x[np - 1],
                                                 Points->y[np - 1], 0.0,
                                                 Points->x[np - 2],
                                                 Points->y[np - 2], 0.0, 0);
                    }

                    G_debug(4, "  length1 = %f length2 = %f", length1,
                            length2);

                    if (length1 < length2) {
                        short_line = line1;
                        long_line = line2;
                    }
                    else {
                        short_line = line2;
                        long_line = line1;
                    }

                    /* Remove end segment from short_line */
                    short_type =
                        Vect_read_line(Map, Points, SCats, abs(short_line));

                    if (short_line > 0) {
                        x = Points->x[1];
                        y = Points->y[1];
                        z = Points->z[1];
                        Vect_line_delete_point(Points, 0);      /* first */
                    }
                    else {
                        x = Points->x[Points->n_points - 2];
                        y = Points->y[Points->n_points - 2];
                        z = Points->z[Points->n_points - 2];
                        Vect_line_delete_point(Points, Points->n_points - 1);   /* last */
                    }

                    /* It may happen that it is one line: node could be deleted,
                     * in that case we have to read the node coords first */
                    Vect_get_node_coor(Map, node, &nx, &ny, &nz);

                    if (Points->n_points > 1) {
                        new_short_line =
                            Vect_rewrite_line(Map, abs(short_line),
                                              short_type, Points, SCats);
                    }
                    else {
                        Vect_delete_line(Map, abs(short_line));
                    }

                    /* It may happen that it is one line, in that case we have to take the new
                     * short line as long line, orientation is not changed */
                    if (abs(line1) == abs(line2)) {
                        if (long_line > 0)
                            long_line = new_short_line;
                        else
                            long_line = -new_short_line;
                    }

                    /* Add new line (must be before rewrite of long_line otherwise node could be deleted) */
                    long_type =
                        Vect_read_line(Map, NULL, LCats, abs(long_line));

                    Vect_reset_cats(OCats);
                    for (j = 0; j < SCats->n_cats; j++) {
                        Vect_cat_set(OCats, SCats->field[j], SCats->cat[j]);
                    }
                    for (j = 0; j < LCats->n_cats; j++) {
                        Vect_cat_set(OCats, LCats->field[j], LCats->cat[j]);
                    }

                    if (long_type == GV_BOUNDARY || short_type == GV_BOUNDARY) {
                        type = GV_BOUNDARY;
                    }
                    else {
                        type = GV_LINE;
                    }

                    Vect_reset_line(Points);
                    Vect_append_point(Points, nx, ny, nz);
                    Vect_append_point(Points, x, y, z);
                    Vect_write_line(Map, type, Points, OCats);

                    if (Err) {
                        Vect_write_line(Err, type, Points, OCats);
                    }

                    /* Snap long_line to the new short_line end */
                    long_type =
                        Vect_read_line(Map, Points, LCats, abs(long_line));
                    if (long_line > 0) {
                        Points->x[0] = x;
                        Points->y[0] = y;
                        Points->z[0] = z;
                    }
                    else {
                        Points->x[Points->n_points - 1] = x;
                        Points->y[Points->n_points - 1] = y;
                        Points->z[Points->n_points - 1] = z;
                    }
                    Vect_line_prune(Points);
                    if (Points->n_points > 1) {
                        Vect_rewrite_line(Map, abs(long_line), long_type,
                                          Points, LCats);
                    }
                    else {
                        Vect_delete_line(Map, abs(long_line));
                    }

                    nmodif += 3;
                    clean = 0;

                    break;
                }

                line1 = line2;
                angle1 = angle2;
            }

            if (clean || !Vect_node_alive(Map, node))
                break;
        }
    }
    G_verbose_message(_("Modifications: %d"), nmodif);

    return (nmodif);
}