improved topology building; bugfix for Vect_copy

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@36401 15284696-431f-4ddb-bdfa-cd5b030d7da7

lib/vector/Vlib/build_nat.c

@@ -86,7 +86,6 @@ int Vect_build_line_area(struct Map_info *Map, int iline, int side)
 
     /* 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 @@ int Vect_isle_find_area(struct Map_info *Map, int isle)
 	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 @@ int Vect_isle_find_area(struct Map_info *Map, int isle)
 		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;
 		}

lib/vector/Vlib/map.c

@@ -174,7 +174,8 @@ Vect_copy(const char *in, const char *mapset, const char *out)
     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;

lib/vector/Vlib/open.c

@@ -215,7 +215,7 @@ Vect__open_old(struct Map_info *Map, const char *name, const char *mapset,
 		_("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);

lib/vector/diglib/poly.c

@@ -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 @@ int dig_get_poly_points(int n_lines, struct line_pnts **LPoints, int *direction,
  **  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 @@ int dig_find_area_poly(struct line_pnts *Points, double *totalarea)
     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 @@ int dig_find_area_poly(struct line_pnts *Points, double *totalarea)
  * 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 @@ double dig_find_poly_orientation(struct line_pnts *Points)
     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]);
 }