extend.c

Go to the documentation of this file.

/*!
  \file lib/vector/vedit/extend.c

  \brief Vedit library - extend lines (adopted from break.c)

  (C) 2017 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 Huidae Cho <grass4u gmail.com>
*/

#include <math.h>
#include <grass/vedit.h>

#define TOL 1e-9

static int extend_lines(struct Map_info *, int, int, int, int, double,
                        struct ilist *);
static int find_extended_intersection(double, double, double, double, double,
                                      double, double *, double *);
static int check_extended_direction(double, double, double, int, double,
                                    double);

/*!
  \brief Extend lines in given threshold

  \code
  1. Extend first line only
          \                      \
     id1   \             ->       \
                                   \
     id2 ----------            -----+----


  2. Extend both lines
          \                      \
     id1   \             ->       \
                                   \
     id2        ---                 +----


  3. Extend first line when both are on the same line
     id1 ---    --- id2  ->    -----+----


  4. Connect two parallel lines (parallel=1)
     id1 ------                -------
                         ->         /
     id2     ------                +-----


  5. Don't connect two parallel lines (parallel=0)
     id1 ------                ------
                         ->
     id2     ------                ------
  \endcode

  \param Map pointer to Map_info
  \param List list of selected lines
  \param nodes 1 for start node, 2 for end node, other for both
  \param parallel connect parallel lines
  \param thresh threshold value

  \return number of modified lines
 */
int Vedit_extend_lines(struct Map_info *Map, struct ilist *List, int nodes,
                       int parallel, double thresh)
{
    int nlines_modified;
    int i, j, first_node, n_nodes;

    struct ilist *List_exclude, *List_found;

    nlines_modified = 0;

    List_exclude = Vect_new_list();
    List_found = Vect_new_list();

    first_node = 0;
    n_nodes = 2;

    switch (nodes) {
        case 1:
            n_nodes = 1;
            break;
        case 2:
            first_node = 1;
            break;
    }

    /* collect lines to be modified */
    for (i = 0; i < List->n_values; i++) {
        int line, extended, node[2];

        line = List->value[i];

        if (!Vect_line_alive(Map, line))
            continue;

        if (Vect_get_line_type(Map, line) & GV_POINTS)
            continue;

        node[0] = node[1] = -1;
        Vect_get_line_nodes(Map, line, &(node[0]), &(node[1]));
        if (node[0] < 0 || node[1] < 0)
            continue;

        extended = 0;
        Vect_reset_list(List_exclude);
        Vect_list_append(List_exclude, line);
        for (j = first_node; j < n_nodes && !extended; j++) {
            double x, y, z;

            /* for each line node find lines in threshold */
            Vect_get_node_coor(Map, node[j], &x, &y, &z);

            do {
                int found;

                /* find first nearest line */
                found = Vect_find_line_list(Map, x, y, z,
                                            GV_LINES, thresh, WITHOUT_Z,
                                            List_exclude, List_found);

                if (found > 0 && Vect_line_alive(Map, found)) {
                    /* try to extend lines (given node) */
                    G_debug(3, "Vedit_extend_lines(): lines=%d,%d", line, found);
                    if (extend_lines(Map, !j, line, found, parallel, thresh,
                                     List)) {
                        G_debug(3, "Vedit_extend_lines(): lines=%d,%d -> extended",
                                line, found);
                        nlines_modified += 2;
                        extended = 1;
                    }
                }

                Vect_list_append(List_exclude, found);
            } while(List_found->n_values > 0 && !extended);
        }
    }

    Vect_destroy_list(List_exclude);
    Vect_destroy_list(List_found);

    return nlines_modified;
}

int extend_lines(struct Map_info *Map, int first, int line_from, int line_to,
                 int parallel, double thresh, struct ilist *List)
{
    /* TODO: If line_from extends to the i'th segment of line_to but the
     * line_from node is closest to the j'th segment of line_to, this function
     * wouldn't work because it only checks intersection of the start/end
     * segment of line_from and the closest segment of line_to (i'th segment).
     */
    int line_new;
    int type_from, type_to;

    struct line_pnts *Points_from, *Points_to, *Points_final;
    struct line_cats *Cats_from, *Cats_to;

    Points_from = Vect_new_line_struct();
    Points_to = Vect_new_line_struct();
    Points_final = Vect_new_line_struct();
    Cats_from = Vect_new_cats_struct();
    Cats_to = Vect_new_cats_struct();

    type_from = Vect_read_line(Map, Points_from, Cats_from, line_from);
    type_to = Vect_read_line(Map, Points_to, Cats_to, line_to);

    line_new = 0;
    if (!(type_from & GV_LINES) || !(type_to & GV_LINES))
        line_new = -1;

    /* avoid too much indentation */
    do {
        int n_points, seg, is, line_to_extended;
        double x, y, px, py, x1, y1;
        double dist, spdist, lpdist, length;
        double angle_t, angle_f;

        if (line_new == -1)
            break;

        n_points = Points_from->n_points - 1;

        if (first) {
            x = Points_from->x[0];
            y = Points_from->y[0];
        }
        else {
            x = Points_from->x[n_points];
            y = Points_from->y[n_points];
        }
        seg = Vect_line_distance(Points_to, x, y, 0.0, WITHOUT_Z,
                                 &px, &py, NULL, &dist, &spdist, &lpdist);

        if (!(seg > 0 && dist > 0.0 && (thresh < 0. || dist <= thresh)))
            break;

        /* lines in threshold */
        length = first ? 0 : Vect_line_length(Points_from);

        /* find angles */
        if (!Vect_point_on_line(Points_from, length, NULL, NULL, NULL, &angle_f,
                               NULL) ||
            !Vect_point_on_line(Points_to, lpdist, NULL, NULL, NULL, &angle_t,
                               NULL))
            break;

        line_to_extended = 0;

        /* extend both lines and find intersection */
        if (!find_extended_intersection(x, y, angle_f, px, py, angle_t,
                                        &x1, &y1)) {
            /* parallel lines */
            if (!parallel)
                break;

            x1 = px;
            y1 = py;
            if (first)
                Vect_line_insert_point(Points_from, 0, x1, y1, 0.0);
            else
                Vect_append_point(Points_from, x1, y1, 0.0);
        } else {
            /* skip if extended into the wrong direction */
            if (!check_extended_direction(x, y, angle_f, first, x1, y1))
                break;

            /* skip if extended too far from line_from */
            if (!Vect_line_distance(Points_from, x1, y1, 0.0, WITHOUT_Z, NULL,
                                    NULL, NULL, &dist, NULL, NULL) ||
                dist > thresh)
                break;

            Vect_line_distance(Points_to, x1, y1, 0.0, WITHOUT_Z, NULL, NULL,
                               NULL, &dist, NULL, NULL);
            /* if intersection point is not on line_to */
            if (dist > TOL) {
                double x2, y2;

                /* skip if not extended from a line_to node */
                if (seg > 1 && seg < Points_to->n_points - 1)
                    break;

                if (seg == 1) {
                    line_to_extended = 1;

                    x2 = Points_to->x[0];
                    y2 = Points_to->y[0];
                } else {
                    line_to_extended = 2;

                    x2 = Points_to->x[Points_to->n_points - 1];
                    y2 = Points_to->y[Points_to->n_points - 1];
                }

                /* skip if extended into the wrong direction */
                if (!check_extended_direction(x2, y2, angle_t, seg == 1, x1,
                                              y1))
                    break;
            }
            /* otherwise, split line_to later */

            /* lines extended -> extend/split line_to */
            /* update line_from */
            if (first) {
                Points_from->x[0] = x1;
                Points_from->y[0] = y1;
            } else {
                Points_from->x[n_points] = x1;
                Points_from->y[n_points] = y1;
            }
        }

        line_new = Vect_rewrite_line(Map, line_from, type_from, Points_from,
                                     Cats_from);
        /* Vect_list_append(List, line_new); */

        Vect_reset_line(Points_final);
        if (line_to_extended == 1) {
            /* extend line_to start node */
            Vect_append_point(Points_final, x1, y1, 0.0);
            for (is = 0; is < Points_to->n_points; is++)
                Vect_append_point(Points_final, Points_to->x[is],
                                  Points_to->y[is], Points_to->z[is]);
            line_new = Vect_rewrite_line(Map, line_to, type_to, Points_final,
                                         Cats_to);
        } else if (line_to_extended == 2) {
            /* extend line_to end node */
            for (is = 0; is < Points_to->n_points; is++)
                Vect_append_point(Points_final, Points_to->x[is],
                                  Points_to->y[is], Points_to->z[is]);
            Vect_append_point(Points_final, x1, y1, 0.0);
            line_new = Vect_rewrite_line(Map, line_to, type_to, Points_final,
                                         Cats_to);
        } else {
            int n_parts = 0;

            /* break line_to */
            /* update line_to  -- first part */
            for (is = 0; is < seg; is++)
                Vect_append_point(Points_final, Points_to->x[is],
                                  Points_to->y[is], Points_to->z[is]);
            Vect_append_point(Points_final, x1, y1, 0.0);

            if (Vect_line_length(Points_final) > 0) {
                n_parts++;
                line_new = Vect_rewrite_line(Map, line_to, type_to,
                                             Points_final, Cats_to);
                /* Vect_list_append(List, line_new); */
            }

            /* write second part */
            Vect_reset_line(Points_final);
            Vect_append_point(Points_final, x1, y1, 0.0);
            for (is = seg; is < Points_to->n_points; is++)
                Vect_append_point(Points_final, Points_to->x[is],
                                  Points_to->y[is], Points_to->z[is]);

            if (Vect_line_length(Points_final) > 0) {
                if (n_parts > 0)
                    line_new = Vect_write_line(Map, type_to, Points_final,
                                               Cats_to);
                else
                    line_new = Vect_rewrite_line(Map, line_to, type_to,
                                                 Points_final, Cats_to);
                /* Vect_list_append(List, line_new); */
            }
        }
    } while(0);

    Vect_destroy_line_struct(Points_from);
    Vect_destroy_line_struct(Points_to);
    Vect_destroy_line_struct(Points_final);
    Vect_destroy_cats_struct(Cats_from);
    Vect_destroy_cats_struct(Cats_to);

    return line_new > 0 ? 1 : 0;
}

static int find_extended_intersection(double x1, double y1, double angle1,
                                      double x2, double y2, double angle2,
                                      double *x, double *y)
{
    double c1, s1, c2, s2, d, a;

    if (fabs(sin(angle1 - angle2)) <= TOL) {
        /* two lines are parallel */
        double angle;

        angle = atan2(y2 - y1, x2 - x1);
        if (fabs(sin(angle - angle1)) <= TOL) {
            /* they are on the same line */
            *x = x2;
            *y = y2;

            return 1;
        }

        /* two lines don't intersect */
        return 0;
    }

    c1 = cos(angle1);
    s1 = sin(angle1);
    c2 = cos(angle2);
    s2 = sin(angle2);
    d = -c1 * s2 + c2 * s1;
    if (d == 0.0)
        /* shouldn't happen again */
        return 0;

    a = (-s2 * (x2 - x1) + c2 * (y2 - y1)) / d;
    *x = x1 + a * c1;
    *y = y1 + a * s1;

    return 1;
}

static int check_extended_direction(double x, double y, double angle,
                                    int start_node, double extx, double exty)
{
    double tmp;
    int xdir, ydir, xext, yext;

    if (start_node)
        angle += M_PI;

    /* expected directions */
    tmp = cos(angle);
    xdir = (fabs(tmp) <= TOL ? 0 : (tmp > 0 ? 1 : -1));
    tmp = sin(angle);
    ydir = (fabs(tmp) <= TOL ? 0 : (tmp > 0 ? 1 : -1));

    /* extended directions */
    tmp = extx - x;
    xext = (fabs(tmp) <= TOL ? 0 : (tmp > 0 ? 1 : -1));
    tmp = exty - y;
    yext = (fabs(tmp) <= TOL ? 0 : (tmp > 0 ? 1 : -1));

    if (xext != 0 && yext != 0) {
        /* no if extended into the wrong direction */
        if (xdir / xext <= 0 || ydir / yext <= 0)
            return 0;
        /* otherwise, ok */
    } else if (xext == 0 && yext == 0) {
        /* snapped to the node, ok */
    } else if (xext == 0) {
        /* vertical extension */
        /* no if not expected or extended into the wrong direction */
        if (xdir != 0 || ydir / yext <= 0)
            return 0;
        /* otherwise, ok */
    } else {
        /* horizontal extension */
        /* no if not expected or extended into the wrong direction */
        if (ydir != 0 || xdir / xext <= 0)
            return 0;
        /* otherwise, ok */
    }

    return 1;
}