extend.c

Go to the documentation of this file.

/*!
   \file lib/vector/vedit/extend.c
 
   \brief Vedit library - extend lines (adopted from break.c)
 
   SPDX-FileCopyrightText: 2017 GRASS Development Team
   SPDX-License-Identifier: GPL-2.0-or-later
 
   \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 G_UNUSED)
{
    /* 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;
}