[GRASS-SVN] r36401 - in grass/trunk/lib/vector: Vlib diglib
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Mar 17 09:32:50 EDT 2009
Author: mmetz
Date: 2009-03-17 09:32:50 -0400 (Tue, 17 Mar 2009)
New Revision: 36401
Modified:
grass/trunk/lib/vector/Vlib/build_nat.c
grass/trunk/lib/vector/Vlib/map.c
grass/trunk/lib/vector/Vlib/open.c
grass/trunk/lib/vector/diglib/poly.c
Log:
improved topology building; bugfix for Vect_copy
Modified: grass/trunk/lib/vector/Vlib/build_nat.c
===================================================================
--- grass/trunk/lib/vector/Vlib/build_nat.c 2009-03-17 10:35:04 UTC (rev 36400)
+++ grass/trunk/lib/vector/Vlib/build_nat.c 2009-03-17 13:32:50 UTC (rev 36401)
@@ -86,7 +86,6 @@
/* dig_find_area_poly(APoints, &area_size); */
- Vect_line_prune(APoints);
area_size = dig_find_poly_orientation(APoints);
/* area_size is not real area size, we are only interested in the sign */
@@ -210,7 +209,7 @@
poly = Vect_point_in_area_outer_ring(Node->x, Node->y, Map, area);
G_debug(3, " poly = %d", poly);
- if (poly == 1) { /* pint in area, but node is not part of area inside isle (would be poly == 2) */
+ if (poly == 1) { /* point in area, but node is not part of area inside isle (would be poly == 2) */
/* In rare case island is inside more areas in that case we have to calculate area
* of outer ring and take the smaller */
if (sel_area == 0) { /* first */
@@ -220,23 +219,27 @@
if (cur_size < 0) { /* second area */
/* This is slow, but should not be called often */
Vect_get_area_points(Map, sel_area, APoints);
- G_begin_polygon_area_calculations();
+ /* G_begin_polygon_area_calculations();
cur_size =
G_area_of_polygon(APoints->x, APoints->y,
- APoints->n_points);
+ APoints->n_points); */
+ /* this is faster, but there may be latlon problems: the poles */
+ dig_find_area_poly(APoints, &cur_size);
G_debug(3, " first area size = %f (n points = %d)",
cur_size, APoints->n_points);
}
Vect_get_area_points(Map, area, APoints);
- size =
+ /* size =
G_area_of_polygon(APoints->x, APoints->y,
- APoints->n_points);
- G_debug(3, " area size = %f (n points = %d)", cur_size,
+ APoints->n_points); */
+ /* this is faster, but there may be latlon problems: the poles */
+ dig_find_area_poly(APoints, &size);
+ G_debug(3, " area size = %f (n points = %d)", size,
APoints->n_points);
- if (size < cur_size) {
+ if (size > 0 && size < cur_size) {
sel_area = area;
cur_size = size;
}
Modified: grass/trunk/lib/vector/Vlib/map.c
===================================================================
--- grass/trunk/lib/vector/Vlib/map.c 2009-03-17 10:35:04 UTC (rev 36400)
+++ grass/trunk/lib/vector/Vlib/map.c 2009-03-17 13:32:50 UTC (rev 36401)
@@ -174,7 +174,8 @@
if (Vect_legal_filename(out) < 0)
G_fatal_error(_("Vector map name is not SQL compliant"));
- xmapset = G_find_vector2(in, mapset);
+ /* remove mapset from fully qualified name with G_find_vector(), confuses G__file_name() */
+ xmapset = G_find_vector(in, mapset);
if (!xmapset) {
G_warning(_("Unable to find vector map <%s> in <%s>"), in, mapset);
return -1;
Modified: grass/trunk/lib/vector/Vlib/open.c
===================================================================
--- grass/trunk/lib/vector/Vlib/open.c 2009-03-17 10:35:04 UTC (rev 36400)
+++ grass/trunk/lib/vector/Vlib/open.c 2009-03-17 13:32:50 UTC (rev 36401)
@@ -215,7 +215,7 @@
_("Unable to open vector map <%s> on level %d. "
"Try to rebuild vector topology by v.build."),
Vect_get_full_name(Map), level_request);
- G_warning(_("Unable to read head file"));
+ G_warning(_("Unable to read head file of vector <%s>"), Vect_get_full_name(Map));
}
G_debug(1, "Level request = %d", level_request);
Modified: grass/trunk/lib/vector/diglib/poly.c
===================================================================
--- grass/trunk/lib/vector/diglib/poly.c 2009-03-17 10:35:04 UTC (rev 36400)
+++ grass/trunk/lib/vector/diglib/poly.c 2009-03-17 13:32:50 UTC (rev 36401)
@@ -8,11 +8,11 @@
*
* PURPOSE: Lower level functions for reading/writing/manipulating vectors.
*
- * COPYRIGHT: (C) 2001 by the GRASS Development Team
+ * COPYRIGHT: (C) 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.
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
*
*****************************************************************************/
#include <math.h>
@@ -89,6 +89,9 @@
** 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.
*/
int dig_find_area_poly(struct line_pnts *Points, double *totalarea)
{
@@ -96,15 +99,16 @@
double *x, *y;
double tot_area;
+ /* prune first with Vect_line_prune(Points) for speed? */
+
x = Points->x;
y = Points->y;
- tot_area = 0.0;
+ /* first point 0 == point n */
+ tot_area = y[0] * (x[1] - x[n - 1]);
for (i = 1; i < n; i++) {
tot_area += y[i] * (x[i + 1] - x[i - 1]);
}
- /* add last point with i == Points->n_points - 1 */
- tot_area += y[i] * (x[1] - x[i - 1]);
*totalarea = 0.5 * tot_area;
@@ -118,12 +122,15 @@
* return value is positive for CW, negative for CCW, 0 for degenerate
*
* Points must be closed polygon
- * only works with pruned lines (no consecutive duplicate vertices)
+ *
+ * this code uses bits and pieces from softSurfer and GEOS
+ * (C) 2000 softSurfer (www.softsurfer.com)
+ * (C) 2006 Refractions Research Inc.
*/
double dig_find_poly_orientation(struct line_pnts *Points)
{
- unsigned int i, pcur = 0;
- unsigned int n = Points->n_points -1; /* skip last point == first point */
+ unsigned int pnext, pprev, pcur = 0;
+ unsigned int npoints = Points->n_points - 1; /* skip last point == first point */
double *x, *y;
/* first find leftmost highest vertex of the polygon */
@@ -132,26 +139,43 @@
x = Points->x;
y = Points->y;
- for (i = 1; i < n; i++) {
-
- if (y[i] < y[pcur])
+ for (pnext = 1; pnext < npoints; pnext++) {
+ if (y[pnext] < y[pcur])
continue;
- else if (y[i] == y[pcur]) { /* just as high */
- if (x[i] < x[pcur]) /* but to the right */
+ else if (y[pnext] == y[pcur]) { /* just as high */
+ if (x[pnext] < x[pcur]) /* but to the right */
continue;
}
- pcur = i; /* a new leftmost highest vertex */
+ pcur = pnext; /* a new leftmost highest vertex */
}
- /* orientation at vertex pcur == signed area for triangle
- * rather use robust determinant of Olivier Dilliers? */
- if (pcur > 0) {
- return (x[pcur + 1] - x[pcur - 1]) * (y[pcur] - y[pcur - 1])
- - (x[pcur] - x[pcur - 1]) * (y[pcur + 1] - y[pcur - 1]);
+ /* Points are not pruned, so ... */
+
+ /* find next distinct point */
+ pnext = pcur + 1;
+ while (pnext != pcur) {
+ if (pnext == npoints)
+ pnext = 0;
+ if (x[pcur] != x[pnext] || y[pcur] != y[pnext])
+ break;
+ pnext++;
}
- else {
- n -= 1;
- return (x[1] - x[n]) * (y[0] - y[n])
- - (x[0] - x[n]) * (y[1] - y[n]);
+
+ /* find previous distinct point */
+ if (pcur == 0)
+ pprev = npoints - 1;
+ else
+ pprev = pcur - 1;
+ while (pprev != pcur) {
+ if (pprev == 0)
+ pprev = npoints;
+ if (x[pcur] != x[pprev] || y[pcur] != y[pprev])
+ break;
+ pprev--;
}
+
+ /* 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])
+ - (x[pcur] - x[pprev]) * (y[pnext] - y[pprev]);
}
More information about the grass-commit
mailing list