#1898: Nathan Wagner's patch that adds a flag 2 to allow ST_DelaunayTriangles to...

git-svn-id: http://svn.osgeo.org/postgis/trunk@11361 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/doc/reference_processing.xml b/doc/reference_processing.xml
index 92acc3a955c8756527c87617117b36a3c38cacff..af17d0db4ba94d15babcce6250e2a43380c0fd58 100644 (file)
--- a/doc/reference_processing.xml
+++ b/doc/reference_processing.xml
@@ -898,12 +898,13 @@ Return a <ulink
 url="http://en.wikipedia.org/wiki/Delaunay_triangulation">Delaunay
 triangulation</ulink> around the vertices of the input geometry.
 Output is a COLLECTION of polygons (for flags=0) or a MULTILINESTRING
-(for flags=1).  The tolerance, if any, is used to snap input vertices
+(for flags=1) or TIN (for flags=2).  The tolerance, if any, is used to snap input vertices
 togheter.
                        </para>
 
                        <para>Availability: 2.1.0 - requires GEOS &gt;= 3.4.0.</para>
                        <para>&Z_support;</para>
+                       <para>&T_support;</para>
 
                </refsection>
                          <refsection>

diff --git a/liblwgeom/lwgeom_geos.c b/liblwgeom/lwgeom_geos.c
index 64906cfde041ef6b8aded8f4ec94e18dfc8acf8f..50bdba9c5499f584e957184fefe74135e2fc9e15 100644 (file)
--- a/liblwgeom/lwgeom_geos.c
+++ b/liblwgeom/lwgeom_geos.c
@@ -17,6 +17,8 @@
 
 #include <stdlib.h>
 
+LWTIN *lwtin_from_geos(const GEOSGeometry *geom, int want3d);
+
 #undef LWGEOM_PROFILE_BUILDAREA
 
 #define LWGEOM_GEOS_ERRMSG_MAXSIZE 256
@@ -1298,9 +1300,10 @@ lwgeom_offsetcurve(const LWLINE *lwline, double size, int quadsegs, int joinStyl
 #endif /* POSTGIS_GEOS_VERSION < 32 */
 }
 
-LWGEOM*
-lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int edgeOnly)
-{
+/*
+ * output = 1 for edges, 2 for TIN, 0 for polygons
+ */
+LWGEOM* lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int output) {
 #if POSTGIS_GEOS_VERSION < 34
        lwerror("lwgeom_delaunay_triangulation: GEOS 3.4 or higher required");
        return NULL;
@@ -1308,6 +1311,11 @@ lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int edg
        GEOSGeometry *g1, *g3;
        LWGEOM *lwgeom_result;
 
+       if (output < 0 || output > 2) {
+               lwerror("lwgeom_delaunay_triangulation: invalid output type specified %d", output);
+               return NULL;
+       }
+
        initGEOS(lwnotice, lwgeom_geos_error);
 
        g1 = (GEOSGeometry *)LWGEOM2GEOS(lwgeom_in);
@@ -1317,7 +1325,8 @@ lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int edg
                return NULL;
        }
 
-       g3 = GEOSDelaunayTriangulation(g1, tolerance, edgeOnly);
+       /* if output != 1 we want polys */
+       g3 = GEOSDelaunayTriangulation(g1, tolerance, output == 1);
 
        /* Don't need input geometry anymore */
        GEOSGeom_destroy(g1);
@@ -1331,12 +1340,21 @@ lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, int edg
        /* LWDEBUGF(3, "result: %s", GEOSGeomToWKT(g3)); */
 
        GEOSSetSRID(g3, lwgeom_get_srid(lwgeom_in));
-       lwgeom_result = GEOS2LWGEOM(g3, lwgeom_has_z(lwgeom_in));
+
+       if (output == 2) {
+               lwgeom_result = (LWGEOM *)lwtin_from_geos(g3, lwgeom_has_z(lwgeom_in));
+       } else {
+               lwgeom_result = GEOS2LWGEOM(g3, lwgeom_has_z(lwgeom_in));
+       }
+
        GEOSGeom_destroy(g3);
 
-       if (lwgeom_result == NULL)
-       {
-               lwerror("lwgeom_delaunay_triangulation: GEOS2LWGEOM returned null");
+       if (lwgeom_result == NULL) {
+               if (output != 2) {
+                       lwerror("lwgeom_delaunay_triangulation: GEOS2LWGEOM returned null");
+               } else {
+                       lwerror("lwgeom_delaunay_triangulation: lwtin_from_geos returned null");
+               }
                return NULL;
        }
 

diff --git a/liblwgeom/lwtin.c b/liblwgeom/lwtin.c
index d7fef8524db127fc04f5abb28981bfcb476a49de..daf001204a27bd6c76c835de85f5d680c4b577ed 100644 (file)
--- a/liblwgeom/lwtin.c
+++ b/liblwgeom/lwtin.c
@@ -13,9 +13,12 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+
+#include "lwgeom_geos.h"
 #include "liblwgeom_internal.h"
 #include "lwgeom_log.h"
 
+LWTIN *lwtin_from_geos(const GEOSGeometry *geom, int want3d);
 
 
 LWTIN* lwtin_add_lwtriangle(LWTIN *mobj, const LWTRIANGLE *obj)
@@ -23,6 +26,68 @@ LWTIN* lwtin_add_lwtriangle(LWTIN *mobj, const LWTRIANGLE *obj)
        return (LWTIN*)lwcollection_add_lwgeom((LWCOLLECTION*)mobj, (LWGEOM*)obj);
 }
 
+LWTIN *lwtin_from_geos(const GEOSGeometry *geom, int want3d) {
+       int type = GEOSGeomTypeId(geom);
+       int hasZ;
+       int SRID = GEOSGetSRID(geom);
+
+       /* GEOS's 0 is equivalent to our unknown as for SRID values */
+       if ( SRID == 0 ) SRID = SRID_UNKNOWN;
+
+       if ( want3d ) {
+               hasZ = GEOSHasZ(geom);
+               if ( ! hasZ ) {
+                       LWDEBUG(3, "Geometry has no Z, won't provide one");
+                       want3d = 0;
+               }
+       }
+
+       switch (type) {
+               LWTRIANGLE **geoms;
+               uint32_t i, ngeoms;
+       case GEOS_GEOMETRYCOLLECTION:
+               LWDEBUG(4, "lwgeom_from_geometry: it's a Collection or Multi");
+
+               ngeoms = GEOSGetNumGeometries(geom);
+               geoms = NULL;
+               if ( ngeoms ) {
+                       geoms = lwalloc(ngeoms * sizeof *geoms);
+                       if (!geoms) {
+                               lwerror("lwtin_from_geos: can't allocate geoms");
+                               return NULL;
+                       }
+                       for (i=0; i<ngeoms; i++) {
+                               const GEOSGeometry *poly, *ring;
+                               const GEOSCoordSequence *cs;
+                               POINTARRAY *pa;
+
+                               poly = GEOSGetGeometryN(geom, i);
+                               ring = GEOSGetExteriorRing(poly);
+                               cs = GEOSGeom_getCoordSeq(ring);
+                               pa = ptarray_from_GEOSCoordSeq(cs, want3d);
+
+                               geoms[i] = lwtriangle_construct(SRID, NULL, pa);
+                       }
+               }
+               return (LWTIN *)lwcollection_construct(TINTYPE, SRID, NULL, ngeoms, (LWGEOM **)geoms);
+       case GEOS_POLYGON:
+       case GEOS_MULTIPOINT:
+       case GEOS_MULTILINESTRING:
+       case GEOS_MULTIPOLYGON:
+       case GEOS_LINESTRING:
+       case GEOS_LINEARRING:
+       case GEOS_POINT:
+               lwerror("lwtin_from_geos: invalid geometry type for tin: %d", type);
+               break;
+
+       default:
+               lwerror("GEOS2LWGEOM: unknown geometry type: %d", type);
+               return NULL;
+       }
+
+       /* shouldn't get here */
+       return NULL;
+}
 
 void lwtin_free(LWTIN *tin)
 {

diff --git a/regress/delaunaytriangles.sql b/regress/delaunaytriangles.sql
index f87cdc66c08db9d22da5f6b4a33461a2a805febc..3f98221cd6b5efbb41d9b5823f0209e9cde9186c 100644 (file)
--- a/regress/delaunaytriangles.sql
+++ b/regress/delaunaytriangles.sql
@@ -3,3 +3,6 @@ SELECT 1,  ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9)'::geometry)
 SELECT 2,  ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9, 8 9)'::geometry));
 SELECT 3,  ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9, 8 9)'::geometry, 2));
 SELECT 4,  ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9, 8 9)'::geometry, 2, 1));
+SELECT 5,  ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9)'::geometry, 0.0, 2));
+SELECT 6,  ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9, 8 9)'::geometry, 0.0, 2));
+SELECT 7,  ST_AsText(ST_DelaunayTriangles('MULTIPOINT(5 5, 6 0, 7 9, 8 9)'::geometry, 2.0, 2));

diff --git a/regress/delaunaytriangles_expected b/regress/delaunaytriangles_expected
index 78e5c6d31f9dcf185bb510e25c1674a61f7d2ae3..d69923b5591c55a4e703b759b8bc8e915f2ff18a 100644 (file)
--- a/regress/delaunaytriangles_expected
+++ b/regress/delaunaytriangles_expected
@@ -2,3 +2,6 @@
 2|GEOMETRYCOLLECTION(POLYGON((5 5,6 0,8 9,5 5)),POLYGON((5 5,8 9,7 9,5 5)))
 3|GEOMETRYCOLLECTION(POLYGON((5 5,6 0,7 9,5 5)))
 4|MULTILINESTRING((5 5,7 9),(5 5,6 0),(6 0,7 9))
+5|TIN(((5 5,6 0,7 9,5 5)))
+6|TIN(((5 5,6 0,8 9,5 5)),((5 5,8 9,7 9,5 5)))
+7|TIN(((5 5,6 0,7 9,5 5)))