buffer.c

Go to the documentation of this file.

/*!
   \file lib/vector/Vlib/buffer.c
 
   \brief Vector library - nearest, adjust, parallel lines
 
   Higher level functions for reading/writing/manipulating vectors.
 
   See buffer2.c for replacement.
 
   (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 <math.h>
#include <grass/vector.h>
 
#define LENGTH(DX, DY) (sqrt((DX * DX) + (DY * DY)))
#define PI             M_PI
 
/* vector() calculates normalized vector form two points */
static void vect(double x1, double y1, double x2, double y2, double *x,
                 double *y)
{
    double dx, dy, l;
 
    dx = x2 - x1;
    dy = y2 - y1;
    l = LENGTH(dx, dy);
    if (l == 0) {
        /* assume that dx == dy == 0, which should give (NaN,NaN) */
        /* without this, very small dx or dy could result in Infinity */
        dx = dy = 0;
    }
    *x = dx / l;
    *y = dy / l;
}
 
/* find_cross find first crossing between segments from s1 to s2 and from s3 to
 *s4
 ** s5 is set to first segment and s6 to second
 ** neighbours are taken as crossing each other only if overlap
 ** returns: 1 found
 **         -1 found overlap
 **          0 not found
 */
static int find_cross(struct line_pnts *Points, int s1, int s2, int s3, int s4,
                      int *s5, int *s6)
{
    int i, j, ret;
    double *x, *y;
 
    G_debug(5, "find_cross(): npoints = %d, s1 = %d, s2 = %d, s3 = %d, s4 = %d",
            Points->n_points, s1, s2, s3, s4);
 
    x = Points->x;
    y = Points->y;
 
    for (i = s1; i <= s2; i++) {
        for (j = s3; j <= s4; j++) {
            if (j == i) {
                continue;
            }
            ret = dig_test_for_intersection(x[i], y[i], x[i + 1], y[i + 1],
                                            x[j], y[j], x[j + 1], y[j + 1]);
            if (ret == 1 && ((i - j) > 1 || (i - j) < -1)) {
                *s5 = i;
                *s6 = j;
                G_debug(5, "  intersection: s5 = %d, s6 = %d", *s5, *s6);
                return 1;
            }
            if (ret == -1) {
                *s5 = i;
                *s6 = j;
                G_debug(5, "  overlap: s5 = %d, s6 = %d", *s5, *s6);
                return -1;
            }
        }
    }
    G_debug(5, "  no intersection");
    return 0;
}
 
/* point_in_buf - test if point px,py is in d buffer of Points
 ** returns:  1 in buffer
 **           0 not  in buffer
 */
static int point_in_buf(struct line_pnts *Points, double px, double py,
                        double d)
{
    int i, np;
    double sd;
 
    np = Points->n_points;
    d *= d;
    for (i = 0; i < np - 1; i++) {
        sd = dig_distance2_point_to_line(px, py, 0, Points->x[i], Points->y[i],
                                         0, Points->x[i + 1], Points->y[i + 1],
                                         0, 0, NULL, NULL, NULL, NULL, NULL);
        if (sd <= d) {
            return 1;
        }
    }
    return 0;
}
 
/* clean_parallel - clean parallel line created by parallel_line:
 ** - looking for loops and if loop doesn't contain any other loop
 **   and centroid of loop is in buffer removes this loop (repeated)
 ** - optionally removes all end points in buffer
 *    parameters:
 *      Points - parallel line
 *      origPoints - original line
 *      d - offset
 *      rm_end - remove end points in buffer
 ** note1: on some lines (multiply selfcrossing; lines with end points
 **        in buffer of line other; some shapes of ends ) may create nosense
 ** note2: this function is stupid and slow, somebody more clever
 **        than I am should write paralle_line + clean_parallel
 **        better;    RB March 2000
 */
static void clean_parallel(struct line_pnts *Points,
                           struct line_pnts *origPoints, double d, int rm_end)
{
    int i, j, np, npn, sa, sb;
    int sa_max = 0;
    int first = 0, current, last, lcount;
    double *x, *y, px, py, ix, iy;
    static struct line_pnts *sPoints = NULL;
 
    G_debug(4, "clean_parallel(): npoints = %d, d = %f, rm_end = %d",
            Points->n_points, d, rm_end);
 
    x = Points->x;
    y = Points->y;
    np = Points->n_points;
 
    if (sPoints == NULL)
        sPoints = Vect_new_line_struct();
 
    Vect_reset_line(sPoints);
 
    npn = 1;
 
    /* remove loops */
    while (first < np - 2) {
        /* find first loop which doesn't contain any other loop */
        current = first;
        last = Points->n_points - 2;
        lcount = 0;
        while (find_cross(Points, current, last - 1, current + 1, last, &sa,
                          &sb) != 0) {
            if (lcount == 0) {
                first = sa;
            } /* move first forward */
 
            current = sa + 1;
            last = sb;
            lcount++;
            G_debug(5, "  current = %d, last = %d, lcount = %d", current, last,
                    lcount);
        }
        if (lcount == 0) {
            break;
        } /* loop not found */
 
        /* ensure sa is monotonically increasing, so npn doesn't reset low */
        if (sa > sa_max)
            sa_max = sa;
        if (sa < sa_max)
            break;
 
        /* remove loop if in buffer */
        if ((sb - sa) == 1) { /* neighbouring lines overlap */
            j = sb + 1;
            npn = sa + 1;
        }
        else {
            Vect_reset_line(sPoints);
            dig_find_intersection(x[sa], y[sa], x[sa + 1], y[sa + 1], x[sb],
                                  y[sb], x[sb + 1], y[sb + 1], &ix, &iy);
            Vect_append_point(sPoints, ix, iy, 0);
            for (i = sa + 1; i < sb + 1; i++) { /* create loop polygon */
                Vect_append_point(sPoints, x[i], y[i], 0);
            }
            Vect_find_poly_centroid(sPoints, &px, &py);
            if (point_in_buf(origPoints, px, py, d)) { /* is loop in buffer ? */
                npn = sa + 1;
                x[npn] = ix;
                y[npn] = iy;
                j = sb + 1;
                npn++;
                if (lcount == 0) {
                    first = sb;
                }
            }
            else { /* loop is not in buffer */
                first = sb;
                continue;
            }
        }
 
        for (i = j; i < Points->n_points; i++) { /* move points down */
            x[npn] = x[i];
            y[npn] = y[i];
            npn++;
        }
        Points->n_points = npn;
    }
 
    if (rm_end) {
        /* remove points from start in buffer */
        j = 0;
        for (i = 0; i < Points->n_points - 1; i++) {
            px = (x[i] + x[i + 1]) / 2;
            py = (y[i] + y[i + 1]) / 2;
            if (point_in_buf(origPoints, x[i], y[i], d * 0.9999) &&
                point_in_buf(origPoints, px, py, d * 0.9999)) {
                j++;
            }
            else {
                break;
            }
        }
        if (j > 0) {
            npn = 0;
            for (i = j; i < Points->n_points; i++) {
                x[npn] = x[i];
                y[npn] = y[i];
                npn++;
            }
            Points->n_points = npn;
        }
        /* remove points from end in buffer */
        j = 0;
        for (i = Points->n_points - 1; i >= 1; i--) {
            px = (x[i] + x[i - 1]) / 2;
            py = (y[i] + y[i - 1]) / 2;
            if (point_in_buf(origPoints, x[i], y[i], d * 0.9999) &&
                point_in_buf(origPoints, px, py, d * 0.9999)) {
                j++;
            }
            else {
                break;
            }
        }
        if (j > 0) {
            Points->n_points -= j;
        }
    }
}
 
/* parallel_line - remove duplicate points from input line and
 *  creates new parallel line in 'd' offset distance;
 *  'tol' is tolerance between arc and polyline;
 *  this function doesn't care about created loops;
 *
 *  New line is written to existing nPoints structure.
 */
static void parallel_line(struct line_pnts *Points, double d, double tol,
                          struct line_pnts *nPoints)
{
    int i, j, np, na, side;
    double *x, *y, nx, ny, tx, ty, vx, vy, ux, uy, wx, wy;
    double atol, atol2, a, av, aw;
 
    G_debug(4, "parallel_line()");
 
    Vect_reset_line(nPoints);
 
    Vect_line_prune(Points);
    np = Points->n_points;
    x = Points->x;
    y = Points->y;
 
    if (np == 0)
        return;
 
    if (np == 1) {
        Vect_append_point(nPoints, x[0], y[0],
                          0); /* ? OK, should make circle for points ? */
        return;
    }
 
    if (d == 0) {
        Vect_copy_xyz_to_pnts(nPoints, x, y, NULL, np);
        return;
    }
 
    side = (int)(d / fabs(d));
    atol = 2 * acos(1 - tol / fabs(d));
 
    for (i = 0; i < np - 1; i++) {
        vect(x[i], y[i], x[i + 1], y[i + 1], &tx, &ty);
        vx = ty * d;
        vy = -tx * d;
 
        nx = x[i] + vx;
        ny = y[i] + vy;
        Vect_append_point(nPoints, nx, ny, 0);
 
        nx = x[i + 1] + vx;
        ny = y[i + 1] + vy;
        Vect_append_point(nPoints, nx, ny, 0);
 
        if (i <
            np - 2) { /* use polyline instead of arc between line segments */
            vect(x[i + 1], y[i + 1], x[i + 2], y[i + 2], &ux, &uy);
            wx = uy * d;
            wy = -ux * d;
            av = atan2(vy, vx);
            aw = atan2(wy, wx);
            a = (aw - av) * side;
            if (a < 0)
                a += 2 * PI;
 
            /* TODO: a <= PI can probably fail because of representation error
             */
            if (a <= PI && a > atol) {
                na = (int)(a / atol);
                atol2 = a / (na + 1) * side;
                for (j = 0; j < na; j++) {
                    av += atol2;
                    nx = x[i + 1] + fabs(d) * cos(av);
                    ny = y[i + 1] + fabs(d) * sin(av);
                    Vect_append_point(nPoints, nx, ny, 0);
                }
            }
        }
    }
    Vect_line_prune(nPoints);
}
 
/*!
   \brief Create parallel line
 
   This function is replaced by Vect_line_parallel2().
 
   \param InPoints input line
   \param distance create parallel line in distance
   \param tolerance maximum distance between theoretical arc and polygon
   segments \param rm_end remove end points falling into distance \param[out]
   OutPoints output line
 
   \return
 */
void Vect_line_parallel(struct line_pnts *InPoints, double distance,
                        double tolerance, int rm_end,
                        struct line_pnts *OutPoints)
{
    G_debug(4,
            "Vect_line_parallel(): npoints = %d, distance = %f, tolerance = %f",
            InPoints->n_points, distance, tolerance);
 
    parallel_line(InPoints, distance, tolerance, OutPoints);
 
    clean_parallel(OutPoints, InPoints, distance, rm_end);
 
    return;
}
 
/*!
   \brief Create buffer around the line line.
 
   This function is replaced by Vect_line_buffer().
 
   Buffer is closed counter clockwise polygon. Warning: output line
   may contain loops!
 
   \param InPoints input line
   \param distance create buffer in distance
   \param tolerance maximum distance between theoretical arc and polygon
   segments \param[out] OutPoints output line
 */
void Vect_line_buffer(const struct line_pnts *InPoints, double distance,
                      double tolerance, struct line_pnts *OutPoints)
{
    double dangle;
    int side, npoints;
    static struct line_pnts *Points = NULL;
    static struct line_pnts *PPoints = NULL;
 
    distance = fabs(distance);
 
    dangle = 2 * acos(1 - tolerance / fabs(distance)); /* angle step */
 
    if (Points == NULL)
        Points = Vect_new_line_struct();
 
    if (PPoints == NULL)
        PPoints = Vect_new_line_struct();
 
    /* Copy and prune input */
    Vect_reset_line(Points);
    Vect_append_points(Points, InPoints, GV_FORWARD);
    Vect_line_prune(Points);
 
    Vect_reset_line(OutPoints);
 
    npoints = Points->n_points;
    if (npoints <= 0) {
        return;
    }
    else if (npoints == 1) { /* make a circle */
        double angle, x, y;
 
        for (angle = 0; angle < 2 * PI; angle += dangle) {
            x = Points->x[0] + distance * cos(angle);
            y = Points->y[0] + distance * sin(angle);
            Vect_append_point(OutPoints, x, y, 0);
        }
        /* Close polygon */
        Vect_append_point(OutPoints, OutPoints->x[0], OutPoints->y[0], 0);
    }
    else { /* 2 and more points */
        for (side = 0; side < 2; side++) {
            double angle, sangle;
            double lx1, ly1, lx2, ly2;
            double x, y, nx, ny, sx, sy, ex, ey;
 
            /* Parallel on one side */
            if (side == 0) {
                Vect_line_parallel(Points, distance, tolerance, 0, PPoints);
                Vect_append_points(OutPoints, PPoints, GV_FORWARD);
            }
            else {
                Vect_line_parallel(Points, -distance, tolerance, 0, PPoints);
                Vect_append_points(OutPoints, PPoints, GV_BACKWARD);
            }
 
            /* Arc at the end */
            /* 2 points at theend of original line */
            if (side == 0) {
                lx1 = Points->x[npoints - 2];
                ly1 = Points->y[npoints - 2];
                lx2 = Points->x[npoints - 1];
                ly2 = Points->y[npoints - 1];
            }
            else {
                lx1 = Points->x[1];
                ly1 = Points->y[1];
                lx2 = Points->x[0];
                ly2 = Points->y[0];
            }
 
            /* normalized vector */
            vect(lx1, ly1, lx2, ly2, &nx, &ny);
 
            /* starting point */
            sangle = atan2(-nx, ny); /* starting angle */
            sx = lx2 + ny * distance;
            sy = ly2 - nx * distance;
 
            /* end point */
            ex = lx2 - ny * distance;
            ey = ly2 + nx * distance;
 
            Vect_append_point(OutPoints, sx, sy, 0);
 
            /* arc */
            for (angle = dangle; angle < PI; angle += dangle) {
                x = lx2 + distance * cos(sangle + angle);
                y = ly2 + distance * sin(sangle + angle);
                Vect_append_point(OutPoints, x, y, 0);
            }
 
            Vect_append_point(OutPoints, ex, ey, 0);
        }
 
        /* Close polygon */
        Vect_append_point(OutPoints, OutPoints->x[0], OutPoints->y[0], 0);
    }
    Vect_line_prune(OutPoints);
 
    return;
}