diglib/poly.c

Go to the documentation of this file.

/*****************************************************************************
 *
 * MODULE:       Vector library
 *
 * AUTHOR(S):    Original author CERL, probably Dave Gerdes.
 *               Update to GRASS 5.7 Radim Blazek.
 *               Update to GRASS 7.0 Markus Metz
 *
 * PURPOSE:      Lower level functions for reading/writing/manipulating vectors.
 *
 * SPDX-FileCopyrightText: 2009 GRASS Development Team
 * SPDX-License-Identifier: GPL-2.0-or-later
 *
 *****************************************************************************/
 
#include <math.h>
#include <grass/vector.h>
 
#ifndef HUGE_VAL
#define HUGE_VAL 9999999999999.0
#endif
 
/*
 * fills BPoints (must be inited previously) by points from input
 * array LPoints
 *
 * each input LPoints[i] must have at least 2 points
 *
 * returns number of points or -1 on error
 */
int dig_get_poly_points(int n_lines, struct line_pnts **LPoints,
                        int *direction, /* line direction: > 0 or < 0 */
                        struct line_pnts *BPoints)
{
    register int i, j, point, start, end, inc;
    struct line_pnts *Points;
    int n_points;
 
    BPoints->n_points = 0;
 
    if (n_lines < 1) {
        return 0;
    }
 
    /* Calc required space */
    n_points = 0;
    for (i = 0; i < n_lines; i++) {
        Points = LPoints[i];
        n_points += Points->n_points - 1; /* each line from first to last - 1 */
    }
    n_points++; /* last point */
 
    if (0 > dig_alloc_points(BPoints, n_points))
        return (-1);
 
    point = 0;
    j = 0;
    for (i = 0; i < n_lines; i++) {
        Points = LPoints[i];
        if (direction[i] > 0) {
            start = 0;
            end = Points->n_points - 1;
            inc = 1;
        }
        else {
            start = Points->n_points - 1;
            end = 0;
            inc = -1;
        }
 
        for (j = start; j != end; j += inc) {
            BPoints->x[point] = Points->x[j];
            BPoints->y[point] = Points->y[j];
            point++;
        }
    }
    /* last point */
    BPoints->x[point] = Points->x[j];
    BPoints->y[point] = Points->y[j];
 
    BPoints->n_points = n_points;
 
    return (BPoints->n_points);
}
 
/*
 * calculate signed area size for polygon
 *
 * points must be closed polygon with first point = last point
 *
 * returns signed area, positive for clockwise, negative for
 * counterclockwise, 0 for degenerate
 */
int dig_find_area_poly(struct line_pnts *Points, double *totalarea)
{
    int i, mid;
    double *x, *y, xshift, yshift;
    double tot_area;
 
    x = Points->x;
    y = Points->y;
 
    /* Get a reference point close to the polygon as new origin.
     * The first point would be good enough, particularly for small
     * polygons. The average of the first and the mid-point is a fast
     * approximation of a reference point with reduced distance to the
     * ring's vertices, also for larger rings.
     *
     * Shift coordinates towards this reference point to make
     * calculation of the signed area more robust by increasing the
     * accuracy for given fp precision limits.
     *
     * Considering the basic formula
     *    area += (x2 - x1) * (y2 + y1)
     * the shift is in theory only needed for addition, not subtraction,
     * but does no harm and is kept for symmetry treating x and y coords.
     *
     * Keep in sync with G_planimetric_polygon_area() in lib/gis/area_poly2.c */
 
    mid = Points->n_points / 2;
    xshift = (x[0] + x[mid]) / 2.;
    yshift = (y[0] + y[mid]) / 2.;
 
    /* line integral: *Points do not need to be pruned */
    /* surveyor's formula is more common, but more prone to
     * fp precision limit errors, and *Points would need to be pruned */
    tot_area = 0.0;
    for (i = 1; i < Points->n_points; i++) {
        tot_area += ((x[i] - xshift) - (x[i - 1] - xshift)) *
                    ((y[i] - yshift) + (y[i - 1] - yshift));
    }
    *totalarea = 0.5 * tot_area;
 
    return (0);
}
 
/*
 * find orientation of polygon (clockwise or counterclockwise)
 * in theory faster than signed area for > 4 vertices, but is not robust
 * against special cases
 * use dig_find_area_poly instead
 *
 * points must be closed polygon with first point = last point
 *
 * this code uses bits and pieces from softSurfer and GEOS
 * (C) 2000 softSurfer (www.softsurfer.com)
 * (C) 2006 Refractions Research Inc.
 *
 * copes with partially collapsed boundaries and 8-shaped isles
 * the code is long and not much faster than dig_find_area_poly
 * it can be written much shorter, but that comes with speed penalty
 *
 * returns orientation, positive for CW, negative for CCW, 0 for degenerate
 */
double dig_find_poly_orientation(struct line_pnts *Points)
{
    unsigned int pnext, pprev, pcur = 0;
    unsigned int lastpoint = Points->n_points - 1;
    double *x, *y, orientation;
 
    x = Points->x;
    y = Points->y;
 
    /* first find leftmost highest vertex of the polygon */
    for (pnext = 1; pnext < lastpoint; pnext++) {
        if (y[pnext] < y[pcur])
            continue;
        else if (y[pnext] == y[pcur]) { /* just as high */
            if (x[pnext] > x[pcur])     /* but to the right */
                continue;
            if (x[pnext] ==
                x[pcur]) { /* duplicate point, self-intersecting polygon ? */
                pprev = (pcur == 0 ? lastpoint - 1 : pcur - 1);
                if (y[pnext - 1] < y[pprev])
                    continue;
            }
        }
        pcur = pnext; /* a new leftmost highest vertex */
    }
 
    /* Points are not pruned, so ... */
    pnext = pcur;
    pprev = pcur;
 
    /* find next distinct point */
    do {
        if (pnext < lastpoint - 1)
            pnext++;
        else
            pnext = 0;
    } while (pnext != pcur && x[pcur] == x[pnext] && y[pcur] == y[pnext]);
 
    /* find previous distinct point */
    do {
        if (pprev > 0)
            pprev--;
        else
            pprev = lastpoint - 1;
    } while (pprev != pcur && x[pcur] == x[pprev] && y[pcur] == y[pprev]);
 
    /* orientation at vertex pcur == signed area for triangle pprev, pcur, pnext
     * rather use robust determinant of Olivier Devillers? */
    orientation = (x[pnext] - x[pprev]) * (y[pcur] - y[pprev]) -
                  (x[pcur] - x[pprev]) * (y[pnext] - y[pprev]);
 
    if (orientation)
        return orientation;
 
    /* orientation is 0, can happen with dirty boundaries, next check */
    /* find rightmost highest vertex of the polygon */
    pcur = 0;
    for (pnext = 1; pnext < lastpoint; pnext++) {
        if (y[pnext] < y[pcur])
            continue;
        else if (y[pnext] == y[pcur]) { /* just as high */
            if (x[pnext] < x[pcur])     /* but to the left */
                continue;
            if (x[pnext] ==
                x[pcur]) { /* duplicate point, self-intersecting polygon ? */
                pprev = (pcur == 0 ? lastpoint - 1 : pcur - 1);
                if (y[pnext - 1] < y[pprev])
                    continue;
            }
        }
        pcur = pnext; /* a new rightmost highest vertex */
    }
 
    /* Points are not pruned, so ... */
    pnext = pcur;
    pprev = pcur;
 
    /* find next distinct point */
    do {
        if (pnext < lastpoint - 1)
            pnext++;
        else
            pnext = 0;
    } while (pnext != pcur && x[pcur] == x[pnext] && y[pcur] == y[pnext]);
 
    /* find previous distinct point */
    do {
        if (pprev > 0)
            pprev--;
        else
            pprev = lastpoint - 1;
    } while (pprev != pcur && x[pcur] == x[pprev] && y[pcur] == y[pprev]);
 
    /* orientation at vertex pcur == signed area for triangle pprev, pcur, pnext
     * rather use robust determinant of Olivier Devillers? */
    orientation = (x[pnext] - x[pprev]) * (y[pcur] - y[pprev]) -
                  (x[pcur] - x[pprev]) * (y[pnext] - y[pprev]);
 
    if (orientation)
        return orientation;
 
    /* orientation is 0, next check */
    /* find leftmost lowest vertex of the polygon */
    pcur = 0;
    for (pnext = 1; pnext < lastpoint; pnext++) {
        if (y[pnext] > y[pcur])
            continue;
        else if (y[pnext] == y[pcur]) { /* just as low */
            if (x[pnext] > x[pcur])     /* but to the right */
                continue;
            if (x[pnext] ==
                x[pcur]) { /* duplicate point, self-intersecting polygon ? */
                pprev = (pcur == 0 ? lastpoint - 1 : pcur - 1);
                if (y[pnext - 1] > y[pprev])
                    continue;
            }
        }
        pcur = pnext; /* a new leftmost lowest vertex */
    }
 
    /* Points are not pruned, so ... */
    pnext = pcur;
    pprev = pcur;
 
    /* find next distinct point */
    do {
        if (pnext < lastpoint - 1)
            pnext++;
        else
            pnext = 0;
    } while (pnext != pcur && x[pcur] == x[pnext] && y[pcur] == y[pnext]);
 
    /* find previous distinct point */
    do {
        if (pprev > 0)
            pprev--;
        else
            pprev = lastpoint - 1;
    } while (pprev != pcur && x[pcur] == x[pprev] && y[pcur] == y[pprev]);
 
    /* orientation at vertex pcur == signed area for triangle pprev, pcur, pnext
     * rather use robust determinant of Olivier Devillers? */
    orientation = (x[pnext] - x[pprev]) * (y[pcur] - y[pprev]) -
                  (x[pcur] - x[pprev]) * (y[pnext] - y[pprev]);
 
    if (orientation)
        return orientation;
 
    /* orientation is 0, last check */
    /* find rightmost lowest vertex of the polygon */
    pcur = 0;
    for (pnext = 1; pnext < lastpoint; pnext++) {
        if (y[pnext] > y[pcur])
            continue;
        else if (y[pnext] == y[pcur]) { /* just as low */
            if (x[pnext] < x[pcur])     /* but to the left */
                continue;
            if (x[pnext] ==
                x[pcur]) { /* duplicate point, self-intersecting polygon ? */
                pprev = (pcur == 0 ? lastpoint - 1 : pcur - 1);
                if (y[pnext - 1] > y[pprev])
                    continue;
            }
        }
        pcur = pnext; /* a new rightmost lowest vertex */
    }
 
    /* Points are not pruned, so ... */
    pnext = pcur;
    pprev = pcur;
 
    /* find next distinct point */
    do {
        if (pnext < lastpoint - 1)
            pnext++;
        else
            pnext = 0;
    } while (pnext != pcur && x[pcur] == x[pnext] && y[pcur] == y[pnext]);
 
    /* find previous distinct point */
    do {
        if (pprev > 0)
            pprev--;
        else
            pprev = lastpoint - 1;
    } while (pprev != pcur && x[pcur] == x[pprev] && y[pcur] == y[pprev]);
 
    /* orientation at vertex pcur == signed area for triangle pprev, pcur, pnext
     * rather use robust determinant of Olivier Devillers? */
    orientation = (x[pnext] - x[pprev]) * (y[pcur] - y[pprev]) -
                  (x[pcur] - x[pprev]) * (y[pnext] - y[pprev]);
 
    return orientation; /* 0 for degenerate */
}