[GRASS-SVN] r36533 - grass/trunk/lib/vector/diglib
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Mar 30 12:16:08 EDT 2009
Author: mmetz
Date: 2009-03-30 12:16:06 -0400 (Mon, 30 Mar 2009)
New Revision: 36533
Modified:
grass/trunk/lib/vector/diglib/poly.c
Log:
find poly orientation now deals with partially collapsed boundaries
Modified: grass/trunk/lib/vector/diglib/poly.c
===================================================================
--- grass/trunk/lib/vector/diglib/poly.c 2009-03-30 16:14:50 UTC (rev 36532)
+++ grass/trunk/lib/vector/diglib/poly.c 2009-03-30 16:16:06 UTC (rev 36533)
@@ -5,6 +5,7 @@
*
* 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.
*
@@ -24,15 +25,18 @@
#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
+ * 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)
+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;
@@ -48,7 +52,7 @@
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 += Points->n_points - 1; /* each line from first to last - 1 */
}
n_points++; /* last point */
@@ -86,12 +90,16 @@
}
/*
- ** Calculate signed area size for polygon.
- **
- ** Total area is positive for clockwise and negative for counterclockwise
- ** Formula modified from
- ** Sunday, Daniel. 2002. Fast Polygon Area and Newell Normal Computation.
- ** Journal of Graphics Tools; 7(2):9-13.
+ * calculate signed area size for polygon
+ *
+ * points must be closed polygon with first point = last point
+ *
+ * formula modified from:
+ * Sunday, Daniel. 2002. Fast Polygon Area and Newell Normal Computation.
+ * Journal of Graphics Tools; 7(2):9-13.
+ *
+ * returns signed area, positive for clockwise, negative for
+ * counterclockwise, 0 for degenerate
*/
int dig_find_area_poly(struct line_pnts *Points, double *totalarea)
{
@@ -117,41 +125,44 @@
}
/*
- * find orientation of polygon
+ * find orientation of polygon (clockwise or counterclockwise)
* faster than signed area for > 4 vertices
- *
- * return value is positive for CW, negative for CCW, 0 for degenerate
*
- * Points must be closed polygon
+ * 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.
+ *
+ * and now copes with partially collapsed boundaries
+ * the code is long but fast
+ * 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;
+ double *x, *y, orientation;
- /* first find leftmost highest vertex of the polygon */
- /* could also be leftmost lowest, rightmost highest or rightmost lowest */
-
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 */
+ if (x[pnext] > x[pcur]) /* but to the right */
continue;
}
pcur = pnext; /* a new leftmost highest vertex */
}
/* Points are not pruned, so ... */
- pprev = pnext = pcur;
+ pnext = pcur;
+ pprev = pcur;
/* find next distinct point */
do {
@@ -171,6 +182,131 @@
/* orientation at vertex pcur == signed area for triangle pprev, pcur, pnext
* rather use robust determinant of Olivier Devillers? */
- return (x[pnext] - x[pprev]) * (y[pcur] - y[pprev])
+ 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;
+ }
+ 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;
+ }
+ 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;
+ }
+ 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 */
}
More information about the grass-commit
mailing list