Geografic-Information-System/lib/vector/Vlib/clean_nodes.c

/*!
   \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);
}