v.buffer: use GEOS

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

vector/v.buffer/geos.c

@@ -0,0 +1,128 @@
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+
+#ifdef HAVE_GEOS
+
+static int ring2pts(const GEOSGeometry *geom, struct line_pnts *Points)
+{
+    int i, ncoords;
+    double x, y, z;
+    const GEOSCoordSequence *seq = NULL;
+
+    G_debug(3, "ring2pts()");
+
+    z = 0.0;
+    ncoords = GEOSGetNumCoordinates(geom);
+    if (!ncoords)
+	G_warning(_("No coordinates in GEOS geometry!"));
+    seq = GEOSGeom_getCoordSeq(geom);
+    for (i = 0; i < ncoords; i++) {
+	GEOSCoordSeq_getX(seq, i, &x);
+	GEOSCoordSeq_getY(seq, i, &y);
+	Vect_append_point(Points, x, y, z);
+    }
+
+    return 1;
+}
+
+static int geom2ring(GEOSGeometry *geom, struct line_pnts **oPoints,
+		       struct line_pnts ***iPoints, int *inner_count)
+{
+    int i, nrings, ngeoms, count;
+    const GEOSGeometry *geom2;
+    struct line_pnts *Points, **arrPoints = *iPoints;
+
+    if (GEOSGeomTypeId(geom) == GEOS_LINESTRING ||
+        GEOSGeomTypeId(geom) == GEOS_LINEARRING) {
+
+	Points = *oPoints;
+	if (Points->n_points == 0)
+	    ring2pts(geom, Points);
+	else {
+	    count = *inner_count + 1;
+	    arrPoints = G_realloc(arrPoints, count * sizeof(struct line_pnts *));
+	    arrPoints[*inner_count] = Vect_new_line_struct();
+	    ring2pts(geom, arrPoints[*inner_count]);
+	    *inner_count = count;
+	    *iPoints = arrPoints;
+	}
+    }
+    else if (GEOSGeomTypeId(geom) == GEOS_POLYGON) {
+	geom2 = GEOSGetExteriorRing(geom);
+	Points = *oPoints;
+	if (Points->n_points == 0)
+	    ring2pts(geom2, Points);
+	else {
+	    count = *inner_count + 1;
+	    arrPoints = G_realloc(arrPoints, count * sizeof(struct line_pnts *));
+	    arrPoints[*inner_count] = Vect_new_line_struct();
+	    ring2pts(geom2, arrPoints[*inner_count]);
+	    *inner_count = count;
+	    *iPoints = arrPoints;
+	}
+
+	nrings = GEOSGetNumInteriorRings(geom);
+	
+	if (nrings > 0) {
+
+	    count = *inner_count + nrings;
+	    arrPoints = G_realloc(arrPoints, count * sizeof(struct line_pnts *));
+
+	    for (i = 0; i < nrings; i++) {
+		geom2 = GEOSGetInteriorRingN(geom, i);
+		arrPoints[*inner_count + i] = Vect_new_line_struct();
+		ring2pts(geom2, arrPoints[*inner_count + i]);
+	    }
+	    *inner_count = count;
+	    *iPoints = arrPoints;
+	}
+    }
+    else if (GEOSGeomTypeId(geom) == GEOS_MULTILINESTRING ||
+             GEOSGeomTypeId(geom) == GEOS_MULTIPOLYGON ||
+	     GEOSGeomTypeId(geom) == GEOS_GEOMETRYCOLLECTION) {
+
+	ngeoms = GEOSGetNumGeometries(geom);
+	for (i = 0; i < ngeoms; i++) {
+	    geom2 = GEOSGetGeometryN(geom, i);
+	    geom2ring((GEOSGeometry *)geom2, oPoints, iPoints, inner_count);
+	}
+    }
+    else
+	G_fatal_error(_("Unknown GEOS geometry type"));
+
+    return 1;
+}
+
+int geos_buffer(struct Map_info *In, int id, int type, double da,
+                GEOSBufferParams *buffer_params, struct line_pnts **oPoints,
+		struct line_pnts ***iPoints, int *inner_count)
+{
+    GEOSGeometry *IGeom = NULL;
+    GEOSGeometry *OGeom = NULL;
+    
+    *oPoints = Vect_new_line_struct();
+    *iPoints = NULL;
+    *inner_count = 0;
+
+    if (type == GV_AREA)
+	IGeom = Vect_read_area_geos(In, id);
+    else
+	IGeom = Vect_read_line_geos(In, id, &type);
+
+    OGeom = GEOSBufferWithParams(IGeom, buffer_params, da);
+    
+    if (!OGeom)
+	G_warning(_("Buffering failed"));
+    
+    geom2ring(OGeom, oPoints, iPoints, inner_count);
+
+    if (IGeom)
+	GEOSGeom_destroy(IGeom);
+    if (OGeom)
+	GEOSGeom_destroy(OGeom);
+
+    return 1;
+}
+
+#endif /* HAVE_GEOS */

vector/v.buffer/main.c

@@ -34,6 +34,11 @@
 #define MAX(X,Y) ((X>Y)?X:Y)
 #endif
 
+#ifdef HAVE_GEOS
+int geos_buffer(struct Map_info *, int, int, double,
+                GEOSBufferParams *, struct line_pnts **,
+		struct line_pnts ***, int *);
+#endif
 
 /* returns 1 if unit_tolerance is adjusted, 0 otherwise */
 int adjust_tolerance(double *tolerance)
@@ -200,11 +205,15 @@ int main(int argc, char *argv[])
 
     /* Attributes if sizecol is used */
     int nrec, ctype;
-    struct field_info *Fi;
+    struct field_info *Fi = NULL;
     dbDriver *Driver;
     dbCatValArray cvarr;
     double size_val, scale;
 
+#ifdef HAVE_GEOS
+    GEOSBufferParams *buffer_params = NULL;
+#endif
+
 
     module = G_define_module();
     G_add_keyword(_("vector"));
@@ -323,7 +332,7 @@ int main(int argc, char *argv[])
     if (scale <= 0.0)
 	G_fatal_error("Illegal scale value");
 
-    da = db = dalpha = 0;
+    da = db = dalpha = unit_tolerance = 0;
     if (dista_opt->answer) {
 	da = atof(dista_opt->answer);
 
@@ -409,7 +418,7 @@ int main(int argc, char *argv[])
     /* Create buffers' boundaries */
     nlines = nareas = 0;
     if ((type & GV_POINTS) || (type & GV_LINES))
-	nlines += Vect_get_num_primitives(&In, type);
+	nlines = Vect_get_num_primitives(&In, type);
     if (type & GV_AREA)
 	nareas = Vect_get_num_areas(&In);
     
@@ -424,12 +433,20 @@ int main(int argc, char *argv[])
 
     Vect_spatial_index_init(&si, 0);
 
+#ifdef HAVE_GEOS
+    initGEOS(G_message, G_fatal_error);
+    buffer_params = GEOSBufferParams_create();
+    GEOSBufferParams_setEndCapStyle(buffer_params, GEOSBUF_CAP_ROUND);
+    GEOSBufferParams_setJoinStyle(buffer_params, GEOSBUF_JOIN_ROUND);
+#endif
+
     /* Lines (and Points) */
-    if ((type & GV_POINTS) || (type & GV_LINES)) {
+    if (nlines > 0) {
 	int ltype;
 
-	if (nlines > 0)
-	    G_message(_("Buffering lines..."));
+	G_message(_("Buffering lines..."));
+
+	nlines = Vect_get_num_lines(&In);
 	for (line = 1; line <= nlines; line++) {
 	    int cat;
 
@@ -492,12 +509,19 @@ int main(int argc, char *argv[])
 
 	    }
 	    else {
+
+#ifdef HAVE_GEOS
+		GEOSBufferParams_setMitreLimit(buffer_params, unit_tolerance);
+		geos_buffer(&In, line, type, da, buffer_params, &(arr_bc_pts.oPoints),
+			    &(arr_bc_pts.iPoints), &(arr_bc_pts.inner_count));
+#else
 		Vect_line_buffer2(Points, da, db, dalpha,
 				  !(straight_flag->answer),
 				  !(nocaps_flag->answer), unit_tolerance,
 				  &(arr_bc_pts.oPoints),
 				  &(arr_bc_pts.iPoints),
 				  &(arr_bc_pts.inner_count));
+#endif
 
 		Vect_write_line(&Out, GV_BOUNDARY, arr_bc_pts.oPoints, BCats);
 		line_id = Vect_write_line(&Buf, GV_BOUNDARY, arr_bc_pts.oPoints, Cats);
@@ -527,11 +551,10 @@ int main(int argc, char *argv[])
     }
 
     /* Areas */
-    if (type & GV_AREA) {
+    if (nareas > 0) {
 	int centroid;
 
-	if (nareas > 0) 
-	    G_message(_("Buffering areas..."));
+	G_message(_("Buffering areas..."));
 	for (area = 1; area <= nareas; area++) {
 	    int cat;
 
@@ -575,12 +598,19 @@ int main(int argc, char *argv[])
 			unit_tolerance);
 	    }
 
+#ifdef HAVE_GEOS
+	    GEOSBufferParams_setSingleSided(buffer_params, 1);
+	    GEOSBufferParams_setMitreLimit(buffer_params, unit_tolerance);
+	    geos_buffer(&In, area, GV_AREA, da, buffer_params, &(arr_bc_pts.oPoints),
+			&(arr_bc_pts.iPoints), &(arr_bc_pts.inner_count));
+#else
 	    Vect_area_buffer2(&In, area, da, db, dalpha,
 			      !(straight_flag->answer),
 			      !(nocaps_flag->answer), unit_tolerance,
 			      &(arr_bc_pts.oPoints),
 			      &(arr_bc_pts.iPoints),
 			      &(arr_bc_pts.inner_count));
+#endif
 
 	    Vect_write_line(&Out, GV_BOUNDARY, arr_bc_pts.oPoints, BCats);
 	    line_id = Vect_write_line(&Buf, GV_BOUNDARY, arr_bc_pts.oPoints, Cats);
@@ -608,6 +638,11 @@ int main(int argc, char *argv[])
 	}
     }
 
+#ifdef HAVE_GEOS
+    GEOSBufferParams_destroy(buffer_params);
+    finishGEOS();
+#endif
+
     verbose = G_verbose();
 
     G_message(_("Cleaning buffers..."));